inst/examples/mapping.made.easy.r

In AtlanticR/bio.utilities: Package that interacts with bio providing base-level generic routines to make your life easier.

## Mapping made easy
## basic philosophy: KISS
## use simple constructs for mapping as much as possible
## Proposed solution: use the "maps" library
## Remember: this is to do quick and dirty plots/visualization not to do exact analyses
##   for the latter: use "sp", "rproj", "rgdal", "lattice", etc
# Created Wednesday 11 February 2015


library(maps)
library(mapdata)  # high resolution world coastlines/polygons
require(lubridate)
require(lattice)
bioLibrary("bio.groundfish")
bioLibrary("bio.utilities")
bioLibrary("bio.bathymetry")


m = maps::map()   # low resolution map of the world .. default is rectangular coords
str(m)


# dealing with projections in a simple manner: use "mapproj"
graphics.off()
require(mapproj)
maps::map(projection="polyconic")

graphics.off()
maps::map(projection="albers", par=c(30,40) )

graphics.off()
maps::map(projection="mollweide" )

graphics.off()
maps::map("worldHires", projection="gnomonic", orient=c(0,-100,0)) # example of orient

graphics.off()
m = maps::map("worldHires", plot=FALSE)
maps::map("worldHires", projection='orth',orient=c(41,-74,0))
str(m)
maps::map.grid(m, col="lightgray", label=FALSE, lty="dotted")

# to add layers to a projected map, we need to deal with projections
points(mapproject(list(y=42.5,x=-68)),col="red", pch="x", cex=3 )

# if more complex projection-based computations are required then
# use library(sp) and the function "CRS"
# as well as rgdal::project
# .. see bio.spatialmethods::lonlat2planar


graphics.off()
maps::map( "worldHires", projection="ortho", xlim=c(-165,-53), ylim=c(-0, 90 ),  )      # national boundaries
maps::map( "worldHires", projection="azequalarea", xlim=c(-165,-53), ylim=c(-0, 90 ),  )      # national boundaries
maps::map( "worldHires", projection="eisenlohr", xlim=c(-165,-53), ylim=c(-0, 90 ),  )      # national boundaries


# Canada only
graphics.off()
maps::map( "worldHires", regions="Canada", xlim=c(-141,-53), ylim=c(40,85), col="gray90", fill=TRUE)

# oops forgot US
maps::map( "worldHires", regions=c("USA"), resolution=0, add=TRUE )

# or both together
graphics.off()
maps::map( "worldHires", regions=c("Canada", "USA"), xlim=c(-141,-53), ylim=c(40,85), col="gray90", fill=TRUE)


# aside .. how to extract polygons
maps::map( "worldHires", regions=c("Canada", "USA"), xlim=c(-141,-53), ylim=c(40,85))
mCan = maps::map( "worldHires", regions="Canada", plot=FALSE, resolution=0)
lines( mCan,col="green" )


# zoom in and use all data
graphics.off()
maps::map( "worldHires", regions=c("Canada", "USA"), xlim=c(-72, -56 ), ylim=c(42,49), fill=FALSE, resolution=0 )

# add some more things
maps::map.scale()
title('Nova Scotia, High resolution')
points( -60, 44, col="red", pch=21, cex=10 )
text( -60, 44, "Sable Island" )
box()


# bring in real data
gs = groundfish.db( "gsinf" )
gs = gs[ which( year( gs$timestamp) == 2010 ), ]

# locations of sets
points( lat~lon, data=gs, col="green", pch=19, cex=0.2 )


# look at salinity in space
graphics.off()
z = gs$bottom_salinity
ncats = 100 ## number of colours
drange = range(z, na.rm=TRUE)

cols = color.code( "seis", seq( drange[1], drange[2], length.out=ncats) ) [ cut( z, ncats, label = FALSE) ]  ## see color.code.r to contruct your own
zi <- rev(order( z ))
points( lat~lon, data=gs, col=cols[zi], pch=21, cex=1 )

# or to fill in the points
points( lat~lon, data=gs, bg=cols[zi], col=cols[zi], pch=21, cex=1.5 )



# zoom in and use all data .. look at temperature as bubble plots
graphics.off()
maps::map( "worldHires", regions=c("Canada", "USA"), xlim=c(-72, -56 ), ylim=c(42,49), fill=FALSE, resolution=0 )
maps::map.scale()
title('Nova Scotia, High resolution')


z = gs$bottom_temperature
ncats = 100 ## number of colours
drange = range(z, na.rm=TRUE)

sz = seq( drange[1], drange[2], length.out=ncats)
zcat = sz[ cut( z, ncats, label = FALSE) ]
zi <- rev(order( z ))
points( lat~lon, data=gs, col="green", pch=21, cex=zcat[zi]/10 )

# or both colour and size of circles
cols = color.code( "seis", seq( drange[1], drange[2], length.out=ncats) ) [ cut( z, ncats, label = FALSE) ]
points( lat~lon, data=gs,  bg=cols[zi], col=cols[zi], pch=21, cex=zcat[zi]/10 )



# or to use transparencies for temperture in this case
library(scales)  #for transparency
graphics.off()
maps::map( "worldHires", regions=c("Canada", "USA"), xlim=c(-72, -56 ), ylim=c(42,49), fill=FALSE, resolution=0 )
maps::map.scale()
title('Nova Scotia, High resolution')

z = gs$bottom_temperature
ncats = 100 ## number of colours
drange = range(z, na.rm=TRUE)


zcat = seq( drange[1], drange[2], length.out=ncats)[ cut( z, ncats, label = FALSE) ]
zalpha = seq( 0.1, 0.9, length.out=ncats)[ cut( z, ncats, label = FALSE) ]
zi <- rev(order( z ))
points( lat~lon, data=gs, col="green", pchR=21, cex=zcat[zi]/10 )

# or both size of circles and transparancy
cols = color.code( "seis", seq( drange[1], drange[2], length.out=ncats) ) [ cut( z, ncats, label = FALSE) ]
points( lat~lon, data=gs,  bg=alpha("darkgreen", zalpha[zi] ), col=alpha("darkgreen", zalpha[zi] ), pch=21, cex=zcat[zi]/10 )


# or both size of circles and transparency and colour! :)
cols = color.code( "seis", seq( drange[1], drange[2], length.out=ncats) ) [ cut( z, ncats, label = FALSE) ]
points( lat~lon, data=gs,  bg=alpha(cols[zi], zalpha[zi] ), col=alpha(cols[zi], zalpha[zi] ), pch=21, cex=zcat[zi]/10 )


## to get some isobaths
graphics.off()
i1000 = isobath.db( p=p, DS="isobath", depths=c(1000 )
maps::map( "worldHires", regions=c("Canada", "USA"), xlim=c(-72, -56 ), ylim=c(42,49), fill=FALSE, resolution=0 )
lines(i1000)

## to bring in some ARCinfo shape files:
library(maptools)  #for shapefiles .. 100 m isobath
fn = find.bio.gis( "dm100_region" )
dm100 <- readShapePoly( fn)   #
maps::map( "worldHires", regions=c("Canada", "USA"), xlim=c(-72, -56 ), ylim=c(42,49), fill=FALSE, resolution=0 )
plot( dm100, add=TRUE )  # there seem to be some problems with the polygon data ... :)




### an example using lattice plots and passing the map coastline to it
# non-lattice plots do not work well inside of lattice .. but you can pass data
datatoplot = groundfish.db( DS="gsinf" )
datatoplot$year = lubridate::year( datatoplot$timestamp )
datatoplot = datatoplot[ which(datatoplot$year %in% 2000:2014 ), ]

# trellis.par.set(col.whitebg())
xyplot( lat ~ lon| factor(year), data=datatoplot,
  as.table=TRUE, xlab="Longitude", ylab="Latitude", main="Example",
  panel = function(x, y, subscripts, ...) {
    panel.xyplot( x, y, pch=21, cex=0.4, col="orange" , ...)
  }
)

print(pl)


require(maps)
require(mapdata)
coast =  maps::map( "worldHires", regions=c("Canada", "USA"), xlim=c(-72, -56 ), ylim=c(42,49), fill=TRUE, resolution=0, plot=FALSE)


trellis.par.set(col.whitebg())
xyplot( lat ~ lon| factor(year), data=datatoplot, coast=coast,
  as.table=TRUE, xlab="Longitude", ylab="Latitude", main="Example",
  panel = function(x, y, subscripts, ...) {
    panel.lines( coast$x, coast$y, col = "grey", lwd=0.7 )
    panel.xyplot( x, y, pch=21, cex=0.4, col="orange" , ...)
  }
)



ncats = 100
drange = range( z, na.rm=TRUE)
cols = color.code( "seis", seq( drange[1], drange[2], length.out=ncats) ) [ cut( z, ncats, label = FALSE) ]

trellis.par.set(col.whitebg())
xyplot( lat ~ lon| factor(year), data=datatoplot, coast=coast, cols=cols,
  as.table=TRUE, xlab="Longitude", ylab="Latitude", main="Example",
  panel = function(x, y, subscripts, cols=cols, colorkey=colorkey, ...) {
    panel.lines( coast$x, coast$y, col = "grey", lwd=0.7 )
    panel.xyplot( x, y, pch=21, cex=0.4, col=cols , ...)
  }
)

Cairo( file="zooplankton_data_yearly.pdf", type="pdf", bg="white",  units="in", width=6, height=8 )
  print(pl)
dev.off()