This is a simple hillshade map in ggplot made from freely available data. However, one can replace this data with their own, change the colors and customise it how they want. It produces a figure that can be published in a scientific journal. We will be using Pitkin County, CO, USA because its topography is complex, so it is a good place to visualise. This method can be applied with any map such as predictive maps, kriging etc. and over larger regions. It adds to the information contained in a map by adding another useful attribute instead of a planar surface. It involves overlaying a desired map with terrain and a background giving more context a visual appeal.

Its recommended that users should have a background in using spatial data in R for this little tutorial. It will allow for more uses and be more meaningful.

library(geodata) #download raster images
library(tigris) #download site as a vector
library(ggmap) #nice background map
library(ggnewscale) #needed for colors
library(ggspatial) #add north and scale bars
library(sf) #manipulate vector data
library(terra) #raster data
library(tidyterra) #more raster and mapping  tools
library(tidyverse) #data tidying and manipulation

First we need to download the data as in region and DEM through tigris and geodata. A lot more can be downloaded through these packages. We need to first get Pitkin County through tigris and then find the centroid of the country to download a 3 arc seconds (~90 m) DEM through geodata

#download Pitkin county, Colorado
pit = counties(state = "Colorado", cb = T)%>% #counties for CO
  filter(NAME == "Pitkin")%>% #filter Pitkin County
  select(NAME)%>% #only select properties needed
  st_transform(4326) #Need to transform coordinates to align with other data

#get centroid for DEM download
mid = st_centroid(pit)%>% #Gets centroid
  st_coordinates()%>% #extract coordinates
  as.vector() #make into a vector (lon, lat)

#download dem for the county
dem = elevation_3s(lon = mid[1], lat = mid[2], path = tempdir())%>% 
  crop(., pit)%>%
  mask(., pit) #Crop and mask to area

#See how we did
plot(dem)

Lets see where we are going to be making the map. Its always good to actually know where and what you are mapping. So, with a simple line of code we can see exactly where Pitkin County is in the lower 48 of the USA, which will give us a better understanding with elevation.

usa = gadm("USA", path = tempdir())%>%
  st_as_sf()%>%
  filter(NAME_1 != "Alaska")%>%
  filter(NAME_1 != "Hawaii")%>%
  select(COUNTRY)

elv = elevation_30s("USA", path = tempdir())%>%
  rename("Elevation" = "USA_elv_msk")

ggplot()+
  geom_spatraster(data = elv)+
  scale_fill_whitebox_c(na.value = NA)+
  geom_sf(data = usa, fill = NA, color = "black", show.legend = F, inherit.aes = F)+
  geom_sf(data = pit, aes(color = NAME), fill = "black")+
  scale_color_manual(values = "black")+
  labs(title = "Lower 48", fill = "Elevation (m)", color = "")+
  theme_bw()

Although this is not a part of the tutorial. Its simple to see we are in the interior at a high elevation of the rocky mountains. In fact, we are close to the second highest peak in the lower 48, Mt Elbert at ~4,400 m and are surrounded by the highest concentration of mountains over 4000 m in the lower US.

Now download background map cropped to area (not masked). Need to register with stadiamaps and run register_stadiamaps(“your key”) first and choose zoom level and maptype (?get_stadiamap). In this case stamen terrain background maps work well, especially for hillshades but there are many options.

#Download map using the region of the county
background = get_stadiamap(as.vector(st_bbox(pit)), zoom = 12,
                   maptype = "stamen_terrain_background")

#Check to see
ggmap(background)+
  labs(x ="", y = "", title = "Background map")

Need to make a hillshade from the DEM. This has a function for a multi angle hillshade, which works well when using color palettes for bathymetric tints. It is not totally necessary as we will not use such a color palette but its nice to have. We need to calculate slope and aspect (in radians with terrain function) and then calculate from what angles do we want shadows for the hillshade (shade function). This is all done with the terra package but can be done with stars as well, just not as straight forward.

#Function for the multi angle hillshade
hillmulti = function(x, angles = c(270, 15, 60, 330), units = "radians"){
  
  top = terrain(x, v = c('slope',"aspect"), unit = units) #slope and aspect
  
  hst = sapply(angles, function(dir){ #hillshade 
    shade(top[[1]], top[[2]],
          angle = 45,
          direction = dir,
          normalize = TRUE)
  })
  sum(rast(hst)) #sum all angles
}

#Run function on dem
hillSlope = hillmulti(dem)

#See how it looks
plot(hillSlope, col = grey(1:100/100), legend = F, main = "Hillshade")
plot(dem, col = rev(rainbow(25, alpha = 0.4)), add = T)

The hillshade looks good so now we make a figure thats publishable using ggmap and ggplot. However, the code starts with ggmap instead of ggplot and this is where almost everything can be customised but I just used the same colours as before. We add a scale bar and north arrow and overlay the hillshade with the background map. Note, the DEM is added for color aesthetics and for the legend as we dont want a hillshade legend. This step requires knowledge of ggplot.

#plot for publication
ggmap(background)+
  geom_spatraster(data = hillSlope, show.legend = F, inherit.aes = F)+
  scale_fill_distiller(palette = "Greys", na.value = NA, guide = "none")+
  new_scale_fill()+
  geom_spatraster(data = dem, alpha = 0.4, show.legend = T, inherit.aes = F)+
  scale_fill_gradientn(colours = rev(rainbow(25)), na.value = NA)+
  geom_sf(data = pit, fill = NA, color = "black", lwd = 1, inherit.aes = F)+
  annotation_north_arrow(location = "bl", pad_y = unit(0.75, "cm"),
                         height = unit(1, "cm"), width = unit(1, "cm"))+
  annotation_scale(location = "bl")+
  coord_sf(crs = 4326)+
  labs(x = "", y = "", fill = "m", title = "Elevation of Pitkin County")+
  theme_bw()

Remember this is can all be custom, so if you dont like the colors, look or size then play around with it and make it your own. For more advanced maps, check out the blog “Adding knowledge to maps” (for bivariate and multivariate mapping, coming soon!)