coastline_layout = function( p, xy.scalebar=c(-2e5, 1.5e5), depths=c( 100, 200, 300, 400, 500, 600, 700 ), plotmap=FALSE, redo=FALSE ) {
fn = file.path( project.datadirectory("aegis", "polygons", "coastline" ), paste( "coastline_layout", p$spatial.domain, "rdata", sep="." ) )
out = NULL
if (!redo) {
if (file.exists(fn)) load(fn)
if( !is.null(out) ) return(out)
}
bounding_domain = Polygon( matrix( c(
min(p$boundingbox[["xlim"]]), min(p$boundingbox[["ylim"]]),
min(p$boundingbox[["xlim"]]), max(p$boundingbox[["ylim"]]),
max(p$boundingbox[["xlim"]]), max(p$boundingbox[["ylim"]]),
max(p$boundingbox[["xlim"]]), min(p$boundingbox[["ylim"]]),
min(p$boundingbox[["xlim"]]), min(p$boundingbox[["ylim"]])
), ncol = 2, byrow = TRUE ) )
bounding_domain = SpatialPolygons(
list(Polygons(list(bounding_domain), ID = "bb")),
proj4string=sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") # expect lon/lat
)
bounding_domain = spTransform( bounding_domain, sp::CRS(p$internal.crs) )
# coastline for mapping via spplot
p = aegis::aegis_parameters( p=p, DS="coastline" )
pp = aegis::spatial_parameters( spatial.domain="canada.east" ) # largest domain
# coast = aegis::coastline.db( p=p, DS=" gshhg coastline highres", no.clip=TRUE )
coast = coastline.db( p=pp, DS="eastcoast_gadm" )
coast = spTransform( coast, sp::CRS(p$internal.crs) )
coast = gSimplify(coast, tol = 0.01) # simplify the polgons a bit (km)
coast = rgeos::gIntersection( bounding_domain, coast, drop_lower_td=TRUE, byid=TRUE, checkValidity=TRUE ) # crop
# sum(gIsValid(coast, byid=TRUE)==FALSE) # check if any bad polys?
# coast = gBuffer(coast, byid=TRUE, width=0)
polyid = gsub( "^bb[[:space:]]", "", names(coast) )
oo = which(duplicated(polyid) )
if (length(oo)>0) {
for ( i in 1:length(oo) ) {
j = which( polyid == polyid[oo[i]] )
for ( k in 2:length(j)) polyid[j[k]] = paste( polyid[j[k]], k, sep="_")
}
}
coast = sp::spChFIDs( coast, polyid ) #fix id's
#
# depth contours
isobs = aegis::isobath.db( p=pp, depths=depths )
isobs = spTransform( isobs, sp::CRS(p$internal.crs) )
isobs = rgeos::gIntersection( bounding_domain, isobs, drop_lower_td=TRUE, byid=TRUE ) # crop
# sum(gIsValid(isobs, byid=TRUE)==FALSE) # check if any bad polys?
# isobs = gBuffer(isobs, byid=TRUE, width=0)
# plot(isobs)
isoid = gsub( "^bb[[:space:]]", "", names(isobs) )
oo = which(duplicated(isoid) )
if (length(oo)>0) {
for ( i in 1:length(oo) ) {
j = which( isoid == isoid[oo[i]] )
for ( k in 2:length(j)) isoid[j[k]] = paste( isoid[j[k]], k, sep="_")
}
}
isobs = sp::spChFIDs( isobs, isoid ) #fix id's
if (plotmap) {
plot(coast, col="grey95")
plot(isobs, add=TRUE, col="lightblue")
sppoly = aegis::maritimes_groundfish_strata( timeperiod="pre2014", returntype="polygons" )
sppoly = spTransform( sppoly, sp::CRS(p$internal.crs) )
plot(sppoly, col="transparent", add=TRUE)
}
out = list(
coastLayout = list(
list("sp.polygons", coast, fill=FALSE, col="gray", first=FALSE ), # outline of NS for plotting with spplot
list("sp.lines", isobs, col="lightgray" ), # outline of NS for plotting with spplot
# list("SpatialPolygonsRescale", layout.north.arrow(), offset=xy.arrow, scale = 100000),
list("SpatialPolygonsRescale", layout.scale.bar(), offset =xy.scalebar, scale = 100000, fill=c("transparent","gray") ),
list("sp.text", xy.scalebar+c(0,-10000), "0"),
list("sp.text", xy.scalebar+c(100000 ,-10000), "100 km")
),
bounding_domain=bounding_domain
)
save(out, file=fn, compress=TRUE)
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.