GECO style map

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)
library(ggplot2)
library(rgdal)  # for readOGR

## Loading required package: sp

## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-2, (SVN revision 1183)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.3, released 2022/10/21
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/rgdal/gdal
##  GDAL does not use iconv for recoding strings.
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 9.1.0, September 1st, 2022, [PJ_VERSION: 910]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

library(rnaturalearth) # for coastline data

## function for creating the color key
source("../R/plot_discrete_cbar.R")

Example data

Obtain an example raster data as NetCDF file, here the file cwdx80.nc, from Zenodo:

Stocker, Benjamin D. (2021). Global rooting zone water storage capacity and rooting depth estimates (v1.0) [Data set]. Zenodo. https://doi.org/10.5281/zenodo.5515246

Save it in ./data/.

Treating resolution

The file cwdx80.nc contains a map of the estimated root zone water storage capacity at 0.05˚ resolution. This is excessively high for a plot of a global map in a journal article (printed on less than A4 format). Therefore, we regrid the data to 0.5˚ first. We do this with the shell command cdo remapbil using CDO. The grid information file is provided here.

system("cdo remapbil,../data/gridfile_halfdeg.txt ../data/cwdx80.nc ../data/cwdx80_halfdeg.nc")

To facilitate things, the file data/cwdx80_halfdeg.nc is provided as part of this repository.

Read the lower resolution NetCDF file into an R data frame using the function rbeni::nc_to_df().

df <- rbeni::nc_to_df("../data/cwdx80_halfdeg.nc", varnam = "cwdx80") |> 
  mutate(lon = round(lon, digits = 2), lat = round(lat, digits = 2))

This is the starting point of our map plotting.

Plot map

Global map

The workflow below shows how spatial data can be plotted. Here, spatial data is a data frame containing a column lon and a column lat for the longitude and latitude in degrees, respectively. the column cwdx80 contains the values of the variable that is to be displayed and visualized by colors of pixels.

## Reduce extent of map in latitude
latmin <- -60
latmax <- 80
lonmin <- -180 
lonmax <- 180

## Define bins of color scale
breaks <- c(seq(0, 100, by = 20), 150, 200, 300, 500, 700, 900, 1200, Inf)

## define domain object
domain <- c(lonmin, lonmax, latmin, latmax)

## read 110 m resolution coastline from NaturalEarth data (is a shapefile)
coast <- ne_coastline(scale = 110, returnclass = "sf")

# download ocean outlines
ocean <- ne_download(
  scale = 110,
  type = "ocean",
  category = "physical",
  returnclass = "sf")

## resize data frame to limit data volume (geographic clipping is done separately)
df <- df |> 
  dplyr::filter(lon > domain[1] & lon < domain[2] & lat > domain[3] & lat < domain[4])
    
nbin <- length(breaks) - 1
breaks_with <- breaks

## Indicate that color scale extends to minus or plus infinity by a triangle
toptriangle <- FALSE
bottomtriangle <- FALSE
if (is.infinite(breaks[length(breaks)])){
  toptriangle <- TRUE
  breaks <- breaks[-(length(breaks)-1)]
}
if (is.infinite(breaks[1])){
  bottomtriangle <- TRUE
  breaks <- breaks[-2]
}
nbin <- length(breaks) - 1

## add dummy rows to make sure values in layer span the entire range
df <- df |> 
  bind_rows(
    tibble(
      lon = NA,
      lat = NA,
      cwdx80 = breaks[1:(length(breaks)-1)] + 0.5 * (breaks[2]-breaks[1])
    )
  )

## make the data discrete
df$layercut <- as.factor(base::cut(df$cwdx80, breaks = breaks, labels = FALSE, include.lowest = TRUE))

## Define colors
# colorscale <- viridis::cividis(nbin, direction = -1)

## or alternatively
colorscale <- scico::scico(nbin, palette = "batlowK", direction = -1)

## Define color of the top and bottom triangle if required
if (toptriangle){
  colorscale <- c(colorscale, colorscale[length(colorscale)])
}
if (bottomtriangle){
  colorscale <- c(colorscale[1], colorscale)
}

## Create ggplot object
ggmap <- ggplot() +

  ## main raster layer
  ## Note: geom_raster() is a fast special case of geom_tile() used when all the tiles are the same size.
  geom_raster(
    data = df, 
    aes(x = lon, y = lat, fill = layercut, color = layercut), 
    interpolate = TRUE,
    show.legend = FALSE
    ) +

  scale_fill_manual(values = colorscale) +
  scale_color_manual(values = colorscale) +

  # ## add ocean
  # geom_sf(data = ocean,
  #         color = NA,
  #         fill = "lightblue") +
  
  ## add coastline
  geom_sf(data = coast,
          color = "grey10",
          fill = NA,
          size = 0.1) +

  coord_sf(xlim = c(lonmin, lonmax),
           ylim = c(latmin, latmax),
           expand = FALSE   # to draw map strictly bounded by the specified extent
           ) +

  ## some layout modifications
  xlab('') + 
  ylab('') +
  theme_bw() +
  theme(axis.ticks.y.right = element_line(),
        axis.ticks.x.top = element_line(),
        panel.grid = element_blank()) +
  
  ## add a title to the map
  labs(title = expression(paste(italic("S")[CWDX80])))

## Warning: Ignoring unknown aesthetics: colour

## use the function obtained from https://github.com/adrfantini/plot_discrete_cbar
## here, it's in the directory ./R/
gglegend <- plot_discrete_cbar(
  breaks           = breaks_with, # Vector of breaks. If +-Inf are used, triangles will be added to the sides the color bar
  colors           = colorscale,
  legend_title     = "(mm)",
  legend_direction = "vertical",
  spacing = "constant",
  expand_size_y = 0.4,
  width = 0.03,
  font_size = 3  # of color key labels
  )

## Warning in plot_discrete_cbar(breaks = breaks_with, colors = colorscale, :
## Ignoring RColorBrewer palette 'Greys', since colors were passed manually

## Combine map and legend with cowplot
cowplot::plot_grid(ggmap, gglegend, ncol = 2, rel_widths = c(1, 0.15))

## Warning: Removed 12 rows containing missing values (geom_raster).

Zoomed-in map

When zooming in over Europe, we want to display nice high-resolution coast outlines and use these for clipping the coarse resolution raster. This is done as explained nicely by Koen Hufkens on his blog.

## Reduce extent of map in latitude
lonmin = -10
lonmax = 25
latmin = 35
latmax = 55

## Define bins of color scale
breaks <- c(seq(0, 100, by = 20), 150, 200, 300, 500, 700, 900, 1200, Inf)

## define domain object
domain <- c(lonmin, lonmax, latmin, latmax)

## read 110 m resolution coastline and countries from NaturalEarth data (is a shapefile)
coast <- ne_coastline(scale = 50, returnclass = "sf")
countries <- ne_countries(scale = 50, returnclass = "sf")

# download ocean outlines
ocean <- ne_download(
  scale = 50,
  type = "ocean",
  category = "physical",
  returnclass = "sf")

## resize data frame to limit data volume (geographic clipping is done separately)
df <- df |> 
  dplyr::filter(lon > domain[1] & lon < domain[2] & lat > domain[3] & lat < domain[4])
    
nbin <- length(breaks) - 1
breaks_with <- breaks

## Indicate that color scale extends to minus or plus infinity by a triangle
toptriangle <- FALSE
bottomtriangle <- FALSE
if (is.infinite(breaks[length(breaks)])){
  toptriangle <- TRUE
  breaks <- breaks[-(length(breaks)-1)]
}
if (is.infinite(breaks[1])){
  bottomtriangle <- TRUE
  breaks <- breaks[-2]
}
nbin <- length(breaks) - 1

## add dummy rows to make sure values in layer span the entire range
df <- df |> 
  bind_rows(
    tibble(
      lon = NA,
      lat = NA,
      cwdx80 = breaks[1:(length(breaks)-1)] + 0.5 * (breaks[2]-breaks[1])
    )
  )

## make the data discrete
df$layercut <- as.factor(
  base::cut(df$cwdx80, 
            breaks = breaks, 
            labels = FALSE, 
            include.lowest = TRUE
            ))

## Define colors
# colorscale <- viridis::cividis(nbin, direction = -1)

## or alternatively
colorscale <- scico::scico(nbin, palette = "batlowK", direction = -1)

## Define color of the top and bottom triangle if required
if (toptriangle){
  colorscale <- c(colorscale, colorscale[length(colorscale)])
}
if (bottomtriangle){
  colorscale <- c(colorscale[1], colorscale)
}

## Create ggplot object
ggmap <- ggplot() +

  ## main raster layer
  ## Note: geom_raster() is a fast special case of geom_tile() used when all the tiles are the same size.
  geom_raster(data = df, 
              aes(x = lon, y = lat, fill = layercut, color = layercut), 
              show.legend = FALSE) +

  scale_fill_manual(values = colorscale) +
  scale_color_manual(values = colorscale) +
  
  ## add country outlines (in year 2022)
  geom_sf(data = countries,
          color = "gray25",
          fill = NA,
          size = 0.1) +

  ## add ocean
  geom_sf(data = ocean,
          color = NA,
          fill = "grey25") +
  
  ## add coastline
  geom_sf(data = coast,
          color = "grey10",
          fill = NA,
          size = 0.1) +
  
  ## limit longitude and latitude extent
  coord_sf(xlim = c(lonmin, lonmax),
           ylim = c(latmin, latmax),
           expand = FALSE   # to draw map strictly bounded by the specified extent
           ) +

  ## some layout modifications
  xlab('') + 
  ylab('') +
  theme_bw() +
  theme(axis.ticks.y.right = element_line(),
        axis.ticks.x.top = element_line(),
        panel.grid = element_blank()
        ) +
  
  ## add a title to the map
  labs(title = expression(paste(italic("S")[CWDX80])))

## Warning: Ignoring unknown aesthetics: colour

## use the function obtained from https://github.com/adrfantini/plot_discrete_cbar
## here, it's in the directory ./R/
gglegend <- plot_discrete_cbar(
  breaks           = breaks_with, # Vector of breaks. If +-Inf are used, triangles will be added to the sides the color bar
  colors           = colorscale,
  legend_title     = "(mm)",
  legend_direction = "vertical",
  spacing = "constant",
  expand_size_y = 0.4,
  width = 0.03,
  font_size = 3  # of color key labels
  )

## Warning in plot_discrete_cbar(breaks = breaks_with, colors = colorscale, :
## Ignoring RColorBrewer palette 'Greys', since colors were passed manually

## Combine map and legend with cowplot
cowplot::plot_grid(ggmap, gglegend, ncol = 2, rel_widths = c(1, 0.15))

## Warning: Removed 12 rows containing missing values (geom_raster).

Different projection

This follows the example from Koen Hufken’s blog post. Doesn’t run currently.

library(sf)
library(terra)

## Reduce extent of map in latitude
latmin <- -60
latmax <- 85
lonmin <- -179.999 
lonmax <- 179.999

## Define bins of color scale
breaks <- c(seq(0, 100, by = 20), 150, 200, 300, 500, 700, 900, 1200, Inf)

## Read NetCDF data as a raster object
rasta <- terra::rast("../data/cwdx80_halfdeg.nc")

## Define projection object
target_crs <- "+proj=robin +over"
robinson <- CRS(target_crs)

## define domain object
domain <- c(lonmin, lonmax, latmin, latmax)

## read 110 m resolution coastline from NaturalEarth data (is a shapefile)
coast <- ne_coastline(scale = 110, returnclass = "sf") |> 
  dplyr::select(featurecla)

# create a bounding box for the robinson projection
# we'll use this as "trim" to remove jagged edges at
# end of the map (due to the curved nature of the
# robinson projection)
bb <- st_union(st_make_grid(
  st_bbox(c(xmin = lonmin,
            xmax = lonmax,
            ymax = latmax,
            ymin = latmin), 
          crs = st_crs("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")),
  n = 100)) |> 
  st_union()

bb_robinson <- st_transform(bb, as.character(robinson))

# clip countries to bounding box
# and transform
coast <- coast |>
  st_transform(robinson) |>
  st_intersection(bb_robinson)

# convert gridded raster data dataframe
g_df <- project(rasta, bb_robinson) |>
  terra::crop(bb_robinson) |>
  as.data.frame(xy = TRUE)

nbin <- length(breaks) - 1
breaks_with <- breaks

## Indicate that color scale extends to minus or plus infinity by a triangle
toptriangle <- FALSE
bottomtriangle <- FALSE
if (is.infinite(breaks[length(breaks)])){
  toptriangle <- TRUE
  breaks <- breaks[-(length(breaks)-1)]
}
if (is.infinite(breaks[1])){
  bottomtriangle <- TRUE
  breaks <- breaks[-2]
}
nbin <- length(breaks) - 1

## add dummy rows to make sure values in layer span the entire range
df <- df |> 
  bind_rows(
    tibble(
      lon = NA,
      lat = NA,
      cwdx80 = breaks[1:(length(breaks)-1)] + 0.5 * (breaks[2]-breaks[1])
    )
  )

## make the data discrete
df$layercut <- as.factor(base::cut(df$cwdx80, breaks = breaks, labels = FALSE, include.lowest = TRUE))

## Define colors
# colorscale <- viridis::cividis(nbin, direction = -1)

## or alternatively
colorscale <- scico::scico(nbin, palette = "batlowK", direction = -1)

## Define color of the top and bottom triangle if required
if (toptriangle){
  colorscale <- c(colorscale, colorscale[length(colorscale)])
}
if (bottomtriangle){
  colorscale <- c(colorscale[1], colorscale)
}

# # This uses instructions from https://datascience.blog.wzb.eu/2019/04/30/zooming-in-on-maps-with-sf-and-ggplot2/
# disp_win_wgs84 <- st_sfc(st_point(c(lonmin, latmin)),
#                          st_point(c(lonmax, latmax)),
#                          crs = 4326)
# disp_win_trans <- st_transform(disp_win_wgs84, crs = target_crs)
# disp_win_coord <- st_coordinates(disp_win_trans)

## Create ggplot object
ggmap <- ggplot()+

  geom_raster(
    data = df,
    aes(
      x = x,
      y = y,
      fill = layercut
    ),
    interpolate = TRUE,
    show.legend = FALSE
    ) +

  scale_fill_manual(values = colorscale) +
  scale_color_manual(values = colorscale) +
  
  ## bounding box
  geom_sf(data = bb_robinson,
          colour='black',
          linetype='solid',
          fill = NA,
          size = 0.1) +
  
  ## coast
  geom_sf(data = coast,
          colour='grey25',
          linetype='solid',
          fill= NA,
          size = 0.1) +
  
  #coord_sf(crs = st_crs("+proj=robin")) +
  
  # # This uses instructions from https://datascience.blog.wzb.eu/2019/04/30/zooming-in-on-maps-with-sf-and-ggplot2/
  # coord_sf(xlim = disp_win_coord[,'X'], 
  #        ylim = disp_win_coord[,'Y'],
  #        datum = target_crs, 
  #        expand = FALSE) +
  
  # coord_sf(xlim = c(lonmin, lonmax),
  #          ylim = c(latmin, latmax),
  #          crs = st_crs("+proj=robin"),
  #          expand = FALSE   # to draw map strictly bounded by the specified extent
  #          ) +
  
  labs(
    title = expression(paste(italic("S")[CWDX80])),
    caption = "
Some caption
you like
"
  ) +
  theme_void()


## use the function obtained from https://github.com/adrfantini/plot_discrete_cbar
## here, it's in the directory ./R/
gglegend <- plot_discrete_cbar(
  breaks           = breaks_with, # Vector of breaks. If +-Inf are used, triangles will be added to the sides the color bar
  colors           = colorscale,
  legend_title     = "(mm)",
  legend_direction = "vertical",
  spacing = "constant",
  expand_size_y = 0.4,
  width = 0.03,
  font_size = 3  # of color key labels
  )

## Combine map and legend with cowplot
cowplot::plot_grid(ggmap, gglegend, ncol = 2, rel_widths = c(1, 0.15))

Using rbeni

The function plot_map4() from the {rbeni} package implements the steps of plotting a map described above. It has some added flexibility in terms of the type of object to be plotted as a map.

# first read NetCDF into a data frame
df <- rbeni::nc_to_df("../data/cwdx80_halfdeg.nc", varnam = "cwdx80") |> 
  mutate(lon = round(lon, digits = 2), lat = round(lat, digits = 2))

# pass the data frame for plotting
rbeni::plot_map4(df,
          varnam = "cwdx80", 
          lonmin = -180, lonmax = 180, latmin = -60, latmax = 80, 
          breaks = c(seq(0, 100, by = 20), 150, 200, 300, 500, 700, 900, 1200, Inf), 
          spacing = "constant", 
          colorscale = "batlowK",
          combine = TRUE, legend_title = "(mm)",
          expand_size_y = 0.4,
          hillshade = FALSE, rivers = FALSE, lakes = FALSE, ocean = TRUE,
          scale = 110
          )

## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

## Warning: Ignoring unknown aesthetics: colour

## Warning in plot_discrete_cbar(breaks = breaks_with, colors = colorscale, :
## Ignoring RColorBrewer palette 'Greys', since colors were passed manually