How to combine SST contour lines and polygons of reef distribution on a world map showing oceans bathymetry?
I would like to create a world map showing the global distribution of coral reefs in combination with the two 20°C isotherms. I have managed to plot them individually, but not in a single graph. The polygons on the coral reef distribution are stored in cd_polygons, but I couldn't add it to the oce plot using map_polygon(). Conversely, the SST data is stored in levitus from the package ocedata, but I was unable to add it to the ggplot using geom_contour(). I prefer the second approach, as I am more familiar with ggplot2 than with base R. Any ideas?
# Load packages
if (!require("pacman")) install.packages("pacman")
#> Lade nötiges Paket: pacman
p_unload("all")
#> The following packages have been unloaded:
#> pacman
pacman::p_load(
marmap, oce, ocedata, sf, terra, tidyterra, ncdf4, ggplot2, dplyr, tidyr, tibble, rnaturalearth, here, pipeR
)
# Download sea surface data from: https://data.unep-wcmc.org/datasets/1; file size: ca. 200 MB
# working, but bloats the reprex with console messages on the download status:
#usethis::use_zip("https://datadownload-production.s3.us-east-1.amazonaws.com/WCMC008_CoralReefs2021_v4_1.zip", destdir = here())
download.file(
"https://datadownload-production.s3.us-east-1.amazonaws.com/WCMC008_CoralReefs2021_v4_1.zip",
destfile = "WCMC008_CoralReefs2021_v4_1.zip",
mode = "wb",
cacheOK = F
)
unzip("WCMC008_CoralReefs2021_v4_1.zip")
unlink("WCMC008_CoralReefs2021_v4_1.zip") # deletes .zip file
# Load data
cd_polygons <- (st_read(here("14_001_WCMC008_CoralReefs2021_v4_1/01_Data/WCMC008_CoralReef2021_Py_v4_1.shp"))) %>>%
select(geometry) %>>%
st_as_sf()
#> Reading layer `WCMC008_CoralReef2021_Py_v4_1' from data source
#> `C:\Users\Marvin Rds\AppData\Local\Temp\RtmpOKlkJI\reprex-3dec40b172a7-brave-mara\14_001_WCMC008_CoralReefs2021_v4_1\01_Data\WCMC008_CoralReef2021_Py_v4_1.shp'
#> using driver `ESRI Shapefile'
#> Simple feature collection with 17504 features and 18 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: -179.9999 ymin: -34.29823 xmax: 179.9999 ymax: 32.51482
#> Geodetic CRS: WGS 84
data(coastlineWorldMedium)
data(levitus, package = "ocedata")
getNOAA.bathy(-180, 180, -90, 90, keep = T) %>>%
as.topo() %>>%
(~bathy)
#> Querying NOAA database ...
#> This may take seconds to minutes, depending on grid size
#> Building bathy matrix ...
#> topo object has data as follows.
#> longitude[1:5400]: -180.00, -179.93, ..., 179.93, 180.00
#> latitude[1:2700]: -90.000, -89.933, ..., 89.933, 90.000
#> z, a 5400x2700 array with value 49.1167679 at [1,1] position
# Plotting attempt with oce
mapPlot(coastlineWorldMedium, col = "lightgray", projection = "+proj=robin", drawBox = F)
mapImage(bathy, col = oceColorsGebco, breaks = seq(-4000, 0, 500)) # bathymetric contours, takes very long to plot
mapGrid()
mapContour(levitus[["longitude"]], levitus[["latitude"]], levitus[["SST"]], levels = 20, lwd = 2)
mapPolygon(coastlineWorldMedium, col = "lightgrey")
## Attempt with ggplot, sf, and tidyterra
# get bathymetric data
getNOAA.bathy(-180, 180, -90, 90, keep = T) %>>%
as.xyz() %>>%
rename(lon = V1, lat = V2, depth = V3) %>>%
filter(depth < .01) %>>%
as_spatraster(xycols = c(1:2), crs = 4326) %>>% # 4326 == EPSG of NOAA bathy
project("+proj=robin") %>>%
(~bat_wld)
#> File already exists ; loading 'marmap_coord_-180;-90;180;90_res_4.csv'
#> class : SpatRaster
#> dimensions : 2700, 5400, 1 (nrow, ncol, nlyr)
#> resolution : 6298.438, 6298.438 (x, y)
#> extent : -17005830, 17005737, -8502985, 8502798 (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
#> source(s) : memory
#> name : depth
#> min value : -1.011325e+04
#> max value : 5.086064e-03
# Import country data
world_map <- ne_countries(scale = "medium", returnclass = "sf")
# Plot using ggplot and sf
world_map %>>%
st_transform(crs = "+proj=robin") %>>%
ggplot() +
geom_spatraster(data = bat_wld, aes(fill = depth)) +
scale_fill_gradientn(
colours = c("#5e24d6", "#22496d", "#042f66", "#054780", "#1074a6", "#218eb7", "#48b5d2", "#72d0e1", "#9ddee7", "#c6edec"),
breaks = c(0, -2500, -5000, -7500),
guide = "none",
na.value = NA
) +
geom_sf() +
geom_sf(data = cd_polygons, aes(col = "red")) +
coord_sf(crs = "+proj=robin", expand = F) +
theme(
panel.background = element_blank(),
panel.grid.major = element_line(
linetype = "solid",
colour = "grey65",
linewidth = .25
),
panel.grid.minor = element_line( # does not show up, why?
linetype = "solid",
colour = "grey65",
linewidth = .25
),
panel.ontop = T,
legend.position = "none"
)
#> <SpatRaster> resampled to 5e+05 cells.
Created on 2024-08-28 with reprex v2.1.1
Solution
This took a minute to figure out! This method combines the data from levitus to the ggplot you created in your question.
This requires the use of a few additional libraries: reshape2 and ggnewscale.
The first step is to manipulate the levitus data. I have a UDF that I know I found from someone else on SO, but I don't remember where (to give them proper credit).
Once the levitus data is processed with the UDF, then it's converted to a SpatRaster (like bat_wld in your question).
sst <- function(dt) { # convert SST into dataframe for raster plotting
dimnames(dt$SST) <- list(long = dt$longitude, lat = dt$latitude)
ret <- melt(dt$SST, value.name = "SST")
cbind(ret, degF = ret$SST * 9/5 + 32)
}
# prep levitus SST for plotting
lsst <- sst(levitus) # manipulate data into data frame
# convert to SpatRaster
lsstr <- rast(lsst[, c(1:2, 4)], type = "xyz", crs = "+init=epsg:4326")
I call the plot exactly as present in your question and appended the call for the levitus data at the end. When I call the levitus data, I use tidyterra's geom_spatraster_contour.
world_map %>>%
st_transform(crs = "+proj=robin") %>>%
ggplot() +
geom_spatraster(data = bat_wld, aes(fill = depth)) +
scale_fill_gradientn(
colours = c("#5e24d6", "#22496d", "#042f66", "#054780", "#1074a6", "#218eb7", "#48b5d2", "#72d0e1", "#9ddee7", "#c6edec"),
breaks = c(0, -2500, -5000, -7500),
guide = "none",
na.value = NA
) +
geom_sf() +
geom_sf(data = cd_polygons, aes(col = "red")) +
coord_sf(crs = "+proj=robin", expand = F) +
theme(
panel.background = element_blank(),
panel.grid.major = element_line(
linetype = "solid",
colour = "grey65",
linewidth = .25
),
panel.grid.minor = element_line( # does not show up, why?
linetype = "solid",
colour = "grey65",
linewidth = .25
),
panel.ontop = T,
legend.position = "none" ### additional layer ###
) + new_scale_color() + # add new color scale & levitus data
# breaks set to match plot created with `oce` plot
tidyterra::geom_spatraster_contour(data = lsstr, breaks = seq(70, 200, 20),
alpha = .7, linewidth = 1, color = "black") +
coord_sf(crs = "+proj=robin", expand = F) # maintain Robinson projection
Here's everything to make this plot altogether
(easier copy + paste)
pacman::p_load(
marmap, oce, ocedata, sf, stringr, terra, tidyterra, ncdf4, ggplot2,
dplyr, tidyr, tibble, rnaturalearth, here, pipeR, usethis, reshape2, ggnewscale
)
sst <- function(dt) { # convert SST into dataframe for raster plotting
dimnames(dt$SST) <- list(long = dt$longitude, lat = dt$latitude)
ret <- melt(dt$SST, value.name = "SST")
cbind(ret, degF = ret$SST * 9/5 + 32)
}
#----------- only need to do one time ----------
## Attempt with oce
# Download sea surface data from: https://data.unep-wcmc.org/datasets/1; file size: ca. 200 MB
# working, but bloats the reprex with console messages on the download status:
#use_zip("https://datadownload-production.s3.us-east-1.amazonaws.com/WCMC008_CoralReefs2021_v4_1.zip", destdir = here())
# download.file(
# "https://datadownload-production.s3.us-east-1.amazonaws.com/WCMC008_CoralReefs2021_v4_1.zip",
# destfile = "WCMC008_CoralReefs2021_v4_1.zip",
# mode = "wb",
# cacheOK = F
# )
#
# unzip("WCMC008_CoralReefs2021_v4_1.zip")
#
# unlink("WCMC008_CoralReefs2021_v4_1.zip") # deletes .zip file
#------------- unchanged from question --------------
# Load data
cd_polygons <- invisible((st_read(here("14_001_WCMC008_CoralReefs2021_v4_1/01_Data/WCMC008_CoralReef2021_Py_v4_1.shp")))) %>>%
select(geometry) %>>%
st_as_sf()
#> Reading layer `WCMC008_CoralReef2021_Py_v4_1' from data source
#> `C:\Users\Marvin Rds\AppData\Local\Temp\RtmpOKlkJI\reprex-3dec730f34d7-milky-cub\14_001_WCMC008_CoralReefs2021_v4_1\01_Data\WCMC008_CoralReef2021_Py_v4_1.shp'
#> using driver `ESRI Shapefile'
#> Simple feature collection with 17504 features and 18 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: -179.9999 ymin: -34.29823 xmax: 179.9999 ymax: 32.51482
#> Geodetic CRS: WGS 84
data(levitus, package = "ocedata")
## Attempt with ggplot, sf, and tidyterra
# get bathymetric data
getNOAA.bathy(-180, 180, -90, 90, keep = T) %>>%
as.xyz() %>>%
rename(lon = V1, lat = V2, depth = V3) %>>%
filter(depth < .01) %>>%
as_spatraster(xycols = c(1:2), crs = 4326) %>>% # 4326 == EPSG of NOAA bathy
project("+proj=robin") %>>%
(~bat_wld)
# Import country data
world_map <- ne_countries(scale = "medium", returnclass = "sf") # package: rnaturalearth
#------------------- new content --------------------
# prep levitus SST for plotting
lsst <- sst(levitus) # manipulate data into data frame
# convert to SpatRaster
lsstr <- rast(lsst[, c(1:2, 4)], type = "xyz", crs = "+init=epsg:4326")
#------------- unchanged from question --------------
world_map %>>%
st_transform(crs = "+proj=robin") %>>%
ggplot() +
geom_spatraster(data = bat_wld, aes(fill = depth)) +
scale_fill_gradientn(
colours = c("#5e24d6", "#22496d", "#042f66", "#054780", "#1074a6", "#218eb7", "#48b5d2", "#72d0e1", "#9ddee7", "#c6edec"),
breaks = c(0, -2500, -5000, -7500),
guide = "none",
na.value = NA
) +
geom_sf() +
geom_sf(data = cd_polygons, aes(col = "red")) +
coord_sf(crs = "+proj=robin", expand = F) +
theme(
panel.background = element_blank(),
panel.grid.major = element_line(
linetype = "solid",
colour = "grey65",
linewidth = .25
),
panel.grid.minor = element_line( # does not show up, why?
linetype = "solid",
colour = "grey65",
linewidth = .25
),
panel.ontop = T,
legend.position = "none" #-------- new layer --------
) + new_scale_color() + # add new color scale & levitus data
# breaks set to match plot created with `oce`
tidyterra::geom_spatraster_contour(data = lsstr, breaks = seq(70, 200, 20),
alpha = .7, linewidth = 1, color = "black") +
coord_sf(crs = "+proj=robin", expand = F) # maintain Robinson projection
- Retrieving expected data.frame for testthat expectation
- How to use an anonymous function with pipe operator?
- R passing RSelenium driver environment as function argument
- Keep file names merging list with lapply?
- replace strings that contains one another
- Stop parsing out zeros after decimals in ggplot2's annotate
- Plot data.tree coloring and labelling by level
- Extracting the date and time out of a date and time format
- How to reshape a dataframe into wide format with specified column pairs
- Show partial segments removed after limit scales for `geom_pointrange`
- Control ggplot2 legend look without affecting the plot
- Concatenate strings over all rows in a single column of a data frame
- Is there some way to keep variable names from.SD+.SDcols together with non .SD variable names in data.table?
- How to create mutually dependent checkbox columns in a DT table?
- A problem with an stacked, grouped and faceted graph y gg plot
- Installing R on Linux: configure: error: libcurl >= 7.28.0 library and headers are required with support for https
- Create a column conditional on multiple other character and numerical columns
- data.frame to array (2 columns)
- Changing options() in a function environment without changing options() in global environment in R?
- pivot(), group_by() and summarise() with nested data
- Why legend is not showing up in plot layout using patchwork in R?
- in dplyr::mutate, dplyr::starts_with works for .before but not .after?
- Finding optimal cut-off points for dividing a variable into equal groups considering ties in which quantiles do not work
- How to add a legend to hline?
- R Shiny Suppress slider handle for sliderInput until click
- brmultinom brglm2 vifs and dredge dont work
- R, how to get current time miliseconds, without extra information ? (Sys.time(), format("%X"))
- Installing multiple versions of R
- Command to see 'R' path that RStudio is using
- Fast pairwise longest common substring from start