Hillshade, colors and marginal plots with tidyterra (I) | One world

[giscus[bot]]

Hillshade, colors and marginal plots with tidyterra (I) | One world

How to overlay SpatRasters - Using shadow effects on relief mappings is a very common technique, that allows to produce informative yet beautiful maps.

https://dieghernan.github.io/202210_tidyterra-hillshade/

[francisvolh]

Hola Diego! Thanks a lot for the tutorial! I am trying to incorporate the use of tidyterra into my workflow, and I am trying to add some elevation hillshade into my maps. I can do it with base R and terra (using the example functions from the terra documentation) but I have already a set of layers in my maps using ggplot2. So I wanted to add the hillshade as a ggplot layer. I tried your method here, but my computer seems to run out of memory when resampling the elevation raster for the colour scale here:

# Use a vector of colors
index <- hill %>%
  mutate(index_col = rescale(shades, to = c(1, length(pal_greys)))) %>%
  mutate(index_col = round(index_col)) %>%
  pull(index_col)

# Error: std::bad_alloc

I am trying a few other options: like using the rescale::scale within the scale_fill_gradient() (I have done this before, but it seems the elevation raster is too large for this, as I also run out of memory), or like not rescaling, and just producing a vector of colours of length(terra:values(raster)) as you did, and using that in the scale_fill_gradient(). But it does take some time to produce (still running, and afraid it may crash). Any tips?

[dieghernan]

Hi @francisvolh and thanks for commenting. I have been testing and I don't find a better solution. It seems that your rasters are too big for this approach, and tidyterra has limitations on handling big rasters.

This is a bit simple solution, but if this is just for mapping I would try to resample the raster as:

# ... code
hill <- shade(slope, aspect, 30, 270)

hill <- terra::spatSample(hill, 500000, method = "regular", as.raster = TRUE)
# normalize names
names(hill) <- "shades"

# ... more code

This should reduce the raster size notably but still have good quality for mapping..

Why 500000 (you may ask)? Well terra::plot() has a parameter maxcell set to 500000, meaning that each raster plot created with terra is resampled before plotting. tidyterra has the same mechanism but we inform about that with a message.

You may want to play with this parameter also when creating your final map, so all the raster may have the same number of cells:

# ....

maxcell <- 1000000
hill <- terra::spatSample(hill, maxcell, method = "regular", as.raster = TRUE)

# ...

ggplot() +
  geom_spatraster(
    data = hill, fill = vector_cols, maxcell = Inf,
    alpha = 1
  ) +
  geom_spatraster(data = r, maxcell = maxcell) +
  scale_fill_hypso_tint_c()

ping @rhijmans so maybe there are other options with terra

[francisvolh]

@dieghernan gracias! So I tried this approach:
https://dominicroye.github.io/en/2022/hillshade-effects/
They suggest to make a dataframe of the hillshade layer (and also the raster, but I just added it with tidyterra as a terra object), and plot it with ggplot2::geom_raster, which works without memory issues.

So it runs nicely, but honestly, I don't see the difference "visually", or it is very subtle. Unless maybe I increase the transparency a lot (reduce the alpha) of my layer of interest, which may "water down" my actual purpose (to show my raster).

Anyway, I will keep playing with the function, and may bug rhijmans for some other issue with the shade function soon! Cheers! and thanks again!

[dieghernan]

Good! @dominicroye is always a great reference for anything, specially related with mapping 😊

[dieghernan]

Hi @francisvolh,

I think I found an alternative approach for big rasters, that consists on plotting the hillshade layer as an RGB layer. This is done mainly in two steps:

1. Reescale the hill object to [0, 255] with terra::stretch().
2. Plotting the result with geom_spatraster_rgb(r = 1, g = 1, b= 1). Basically this is like to create a grey shade for each value of hill depending on the intensity.

I tried with a SpatRaster of China (29.8Mb, 33,741,360 cells), see:

library(terra)
library(tidyterra)
library(ggplot2)
library(dplyr)
library(scales)

# Get the data
library(geodata)

mydir <- "."

# Try with China, 29.8 MB
r_init <- elevation_30s("China", path = mydir)


r_init
#> class       : SpatRaster 
#> dimensions  : 4572, 7380, 1  (nrow, ncol, nlyr)
#> resolution  : 0.008333333, 0.008333333  (x, y)
#> extent      : 73.4, 134.9, 15.6, 53.7  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : CHN_elv_msk.tif 
#> name        : CHN_elv_msk 
#> min value   :        -263 
#> max value   :        8233

label_comma()(ncell(r_init))
#> [1] "33,741,360"

# For better handling we set here the names
names(r_init) <- "alt"

# We don't want values lower than 0
r <- clamp(r_init, lower = 0)

## Create hillshade effect <---  SEE CHANGES!
slope <- terrain(r, "slope", unit = "radians")
aspect <- terrain(r, "aspect", unit = "radians")
hill <- shade(slope, aspect, 30, 270)

# normalize names
names(hill) <- "shades"
hill_stretch <- stretch(hill, minv = 0, maxv = 255) # <------- HERE

mmax <- minmax(r) %>% as.vector()

# Custom breaks
mybrks <- c(0, 250, 500, 750, seq(1000, 8000, by = 1000))

# The plot

# Max cell, 2,000,000 for plotting, this is high-res
mcell <- 2000000
base_text_size <- 35

ggplot() +
  # Hillshade with rgb instead of vector of colors <------ CHANGE HERE
  geom_spatraster_rgb(
    data = hill_stretch, r = 1, g = 1, b = 1, alpha = 1,
    maxcell = mcell
  ) + # <------- UP TO HERE
  # Avoid resampling with maxcell
  geom_spatraster(data = r, maxcell = mcell) +
  scale_fill_hypso_tint_c(
    limits = mmax,
    palette = "dem_poster",
    alpha = 0.4,
    labels = label_comma(),
    # For the legend I use custom breaks
    breaks = mybrks
  ) +
  # Change guide
  guides(fill = guide_legend(
    title = "   m.",
    direction = "horizontal",
    nrow = 1,
    override.aes = list(alpha = 1)
  )) +
  labs(
    title = "Elevation of China",
    subtitle = "Hillshade and hypsometric tint blend",
    caption = paste0(
      "@dhernangomez using tidyterra, ggplot2, geodata R packages.",
      " Data: Shuttle Radar Topography Mission (SRTM)"
    )
  ) +
  theme_minimal(base_family = "serif") +
  theme(
    plot.background = element_rect("grey97", colour = NA),
    plot.margin = margin(20, 20, 20, 20),
    plot.caption = element_text(size = base_text_size * 0.5),
    plot.title = element_text(face = "bold", size = base_text_size * 1.4),
    plot.subtitle = element_text(
      margin = margin(b = 10),
      size = base_text_size
    ),
    axis.text = element_text(size = base_text_size * 0.7),
    legend.position = "bottom",
    legend.title = element_text(size = base_text_size * 0.8),
    legend.text = element_text(size = base_text_size * 0.7),
    legend.key = element_rect("grey50"),
    legend.spacing.x = unit(0, "pt"),
    legend.text.position = "bottom",
    legend.title.position = "top",
    legend.key.width = unit(7, "line"),
    legend.key.height = unit(0.5, "line")
  )
#> SpatRaster resampled to ncells = 2001858
#> SpatRaster resampled to ncells = 2001858

^{Created on 2024-03-25 with reprex v2.1.0}

[francisvolh]

Hola @dieghernan !

Nice! I can give it a go eventually. The method of using a dataframe has worked for me so far. But I can still try this and see how they differ! Thanks for updating me on your results!

[gamamo]

Hi! Would you know how to plot "no data values" in a color different than those in the colour scheme?

[dieghernan]

I guess you refer to NAs. Usually in the ggplot2 system use the na.value arg in the scales, you have several examples in https://dieghernan.github.io/tidyterra/articles/faqs.html addressing this