Making pretty maps with sdmTMB outputSource:
If the code in this vignette has not been evaluated, a rendered version is available on the documentation site under ‘Articles’.
This vignette illustrates how to make gridded/raster maps of sdmTMB output (or gridded data from any source) on a base map of class sf.
First, we will download coastline data for our region. We will use a
medium resolution version. For applied situations on finer scales, you
may want to use
scale = "large".
map_data <- rnaturalearth::ne_countries( scale = "medium", returnclass = "sf", country = "canada") # Crop the polygon for plotting and efficiency: # st_bbox(map_data) # find the rough coordinates bc_coast <- suppressWarnings(suppressMessages( st_crop(map_data, c(xmin = -134, ymin = 46, xmax = -120, ymax = 57))))
Transform our map into UTM 9 coordinates, which is the equal-area projection we fit in:
utm_zone9 <- 3156 bc_coast_proj <- sf::st_transform(bc_coast, crs = utm_zone9)
If you don’t know the crs number (3156 above), you can look it up by googling it or looking on https://epsg.io/. For example, https://epsg.io/3156.
Look at our base map:
Now we will fit a model and predict on a grid. The coordinates X and Y are already UTM zone 9 coordinates (in km).
d <- subset(pcod, year == 2017) pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 30) #> as(<dgCMatrix>, "dgTMatrix") is deprecated since Matrix 1.5-0; do as(., "TsparseMatrix") instead m <- sdmTMB(density ~ 0 + depth_scaled + depth_scaled2, data = d, mesh = pcod_spde, family = tweedie(link = "log")) #> Warning in checkMatrixPackageVersion(): Package version inconsistency detected. #> TMB was built with Matrix version 1.5.3 #> Current Matrix version is 1.5.1 #> Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package p <- predict(m, newdata = qcs_grid)
We can use
sf::st_boundary() to help find reasonable
limits on our map (output not shown):
Finally, we will combine our gridded predictions with the base map.
We will multiply the
Y columns by 1000
because we worked in UTM km for model fitting (to avoid computational
issues with the scale of the range paramater):
ggplot(bc_coast_proj) + geom_sf() + geom_raster(data = p, aes(x = X * 1000, y = Y * 1000, fill = exp(est))) + xlim(230957.7, 1157991 - 300000) + ylim(5366427, 6353456 - 200000) + scale_fill_viridis_c() + theme_light() + labs(fill = "Predicted\ndensity") + labs(x = "Longitude", y = "Latitude")