Elevations

Create a map using elevations as the base layer

It used to be that you could make cool maps of study areas with ggmap. But now you need a Google API with a credit card attached to it to do that. So while you could do that personally if you wanted, it’s useless for teaching purposes. Also although there is a lot of free/cheap use allowed with the API, at somepoint you can hit real costs and I don’t want to risk making a mistake with my Google Maps API getting out.

So this example uses elevation data to make a study area map. I use the raster and tmap packages. I like tmap a lot. Super easy to layer map features together like I want.

Credits

My example is adapted from Visualize Spatial Data in a Open Source Spatial Analytics course by Dr. Aaron Maxwell at West Virginia University.

Step 0 Load the needed packages

library(raster)
library(tmap)

Step 1 define the area

In this case, I want Bristol Bay. I just looked at a map. The numbers are longitude and latitude. Note, I know that the elevation data is in lat-lon so this will work. Sometimes your coordinate data are not in lat-lon so you’ll have to convert coordinate systems.

BB <- raster::extent(-162, -154, 56, 60)

Step 2 download elevation data

Download using getData function in raster. There is no documentation on the altitude object that you get. It’s a list and [[1]] is the mainland and [[2]] is Alaska. My guess is [[3]] is Hawaii.

This is multiple files. It’ll download once, put into the folder in path and then load from there next time you run the code. Unfortunately Rmarkdown seems to always download the files again. So I put eval=FALSE on this chunk and load from a file.

USAelevation <- raster::getData("alt", country = "USA", mask=FALSE, download = TRUE, path="raster_data")
BB.elevation <- raster::crop(USAelevation[[2]], BB)
save(BB.elevation, file="raster_data/BB_elevation.rds")

load("raster_data/BB_elevation.rds")

Step 3 Make a hillshade layer

This makes your elevation plots look a lot better by adding hill shading. The hill shading is what will make your map pop so you will need to futz with the numbers in the hillShade function to get the effect that you want. The values you see in ?hillShade make for a drab map in my opinion.

BB.slope <- raster::terrain(BB.elevation, opt = "slope")
BB.aspect <- raster::terrain(BB.elevation, opt = "aspect")
hill <- raster::hillShade(BB.slope, BB.aspect, 1, 55)

Save the coordinate system for the elevation data as I will need this later.

hill.crs <- crs(hill)

Now you can make a plot of that using default raster plotting. It is ugly but makes sure we have the right data and location. At this point, you won’t have an idea what the final plot will look like.

plot(hill, col = grey(0:100/100), legend = FALSE, 
     main = "Bristol Bay")
plot(BB.elevation, col = terrain.colors(7, alpha = 0.35), add = TRUE)

## Step 4. Make pretty with tmap package

The tmap package is the easiest one I found. This map pops the way I want. I futzed with hillShade, the pallettes, and style=... in tm_raster() until I got the effect that I wanted.

pal <- terrain.colors(7)
p <- tm_shape(hill)+
  tm_raster(palette="-Greys", style="cont", legend.show=FALSE)+
  tm_shape(BB.elevation)+
  tm_raster(alpha=.5, style="quantile", n=7, palette=pal)
p

Some other pallettes you might try:

pal <- get_brewer_pal("-Greys", n =12, plot=FALSE)
pal <- get_brewer_pal("cividis", n =12, plot=FALSE)
pal <- topo.colors(7)

Step 5. Add on map components

I want a scale, compass, lat-lon axes, and I don’t want the elevation legend. I set fig.width=10 in this chunk so that the x-axis tick labels (longitude). If plot is small, tmap will not show them.

pplus <- p +
  tm_compass(type="rose", position=c(0.7, 0.1), size=6) +
  tm_scale_bar(position = c(0.6, .005), text.size=.8) +
  tm_layout(legend.show=FALSE) +
  tm_graticules(lines=FALSE, projection=as.character(hill.crs), labels.size=1)
pplus

Step 5. Add a coastline

Here is how to add a coastline if you wanted. Note the cropping takes a long time.

usashp <- raster::getData("GADM", country = "USA", level = 1, path = "data")
akshp <- subset(usashp, NAME_1 %in% c("Alaska"))
akshp <- raster::crop(akshp, raster::extent(-180, -140, 50, 80))
bbshp <- raster::crop(akshp, BB)
save(akshp, bbshp, file="raster_data/bbshp.RData")

So I don’t run the code above over and over. I save the shapefile.

load("raster_data/bbshp.RData")

pplus + tm_shape(bbshp) + tm_borders("black", lwd=3)

Step 6 Make my river lines

This part I did manually because I couldn’t figure a good way to do it otherwise. I make a raster plot where I have futzed with the elevation colors so I can see the rivers.

plot(hill, col = grey(0:100/100), legend = FALSE)
plot(BB.elevation, col = topo.colors(25, alpha = 0.35), add = TRUE, legend=FALSE, breaks=c(0,2,3,4,5,7,10,13,15,20,25,30,35,40,45,50,60,800,1000,2000,3000))

Then I use drawLine() to manually add the lines onto the map. That is going to save a shapefile with all the lines. This is tedious and doesn’t look great but is good enough for my purposes. I have to do each river individually.

Here is the code for Togiak River. I go through this one by one for each river.

togiak <- drawLine()
crs(togiak) <- hill.crs
togiak$ID <- "Togiak"

Once I have made all the lines I bind them together in one SpatialLines object.

bbrivers <- raster::bind(togiak, igushik,
                         wood, nushagak,
                         kvichak, naknek,
                         egegik, ugashik)

Step 7 Add rivers to my plot

I put eval=FALSE on the chunks above and load from a saved data file. Save command is lower after I make labels.

load("data/bbrivers.RData")

ppluswrivers <- pplus +
  tm_shape(bbrivers) +
  tm_lines(lwd=2)
ppluswrivers

Step 8 Add labels to rivers

I could just do this but it is hard to get the labels where I want.

ppluswrivers +
  tm_text("ID", size=0.75, auto.placement=TRUE,just="left",xmod=.2)

Instead I will create a SpatialPoints object with the ends of the rivers. I was careful to always create my rivers starting at the ocean.

fun <- function(x){
  SP<-SpatialPoints(coords = cbind(rev(x[,1])[1],rev(x[,2])[1]))
  crs(SP) <- hill.crs
  SP
  }
aaa=lapply(coordinates(bbrivers),function(x){fun(x[[1]])})
for(i in 1:length(aaa)) aaa[[i]]$ID <- bbrivers[[1]][i]

Now bind together and save.

bbrivers.lab <- aaa[[1]]
for(i in 2:length(aaa))
bbrivers.lab <- raster::bind(bbrivers.lab, aaa[[i]])
save(bbrivers, bbrivers.lab, file="data/bbrivers.RData")

I put eval=FALSE on the chunks above and load from a saved data file.

load("data/bbrivers.RData")

Now I can add the labels and they will be at the ends of the rivers. I tweek xmod and ymod (the jitter) until the labels don’t overlap my rivers.

finalp <- ppluswrivers +
  tm_shape(bbrivers.lab) +
  #tm_text("ID", size=0.75, just="left",xmod=0.2, ymod=0.2, bg.color="white", bg.alpha=0.5)
  tm_text("ID", scale=1, root=4, size.lowerbound = .6, 
        bg.color="white", bg.alpha = .4, 
        just="left",xmod=0.2, ymod=0.2)
finalp

Step 9 Make a high dpi version of the figure

I had to add outer.margins so that it didn’t cut off the y-axis tick labels. You need to pass in fig.width and dpi to the get the effect you want with your image.

Big tiff image (dpi=300) for print publication.

tiff(filename="BB_sockeye_rivers.tiff", units="in", width=6, height=5, res=300)
finalp + tm_layout(outer.margins=c(.02,.05,.01,.03))
dev.off()

Small png image for rmarkdown docs.

png(filename="BB_sockeye_rivers.png", units="in", width=6, height=5, res=96)
finalp + tm_layout(outer.margins=c(.02,.05,.01,.03))
dev.off()

Adding an inset of the location

Create the Alaska inset map on a non-deforming coordinate system.

lambert <- "+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
BB.box <- as(BB, "SpatialPolygons")
crs(BB.box) <- crs(akshp)
BB.box2 <- spTransform(BB.box, lambert)
BB.box2$name <- "Bristol Bay"
akshp2 <- spTransform(akshp, lambert)
insetmap <- tm_shape(akshp2) + 
  tm_polygons() +
  tm_text("NAME_1")+
   tm_shape(BB.box2) + tm_borders(lwd=2, col="blue") +
  tm_text("name", size=0.7, xmod=-2, ymod=2, col="blue")

Add that as an inset. We use grid::viewport.

library(grid)
xy <- st_bbox(akshp2)
asp2 <- (xy$xmax - xy$xmin)/(xy$ymax - xy$ymin)
h <- 0.25
w <- asp2 * h
vp <- viewport(x=.084, y=0.086, width = w, height=h, 
               just=c("left", "bottom"))
tmap_save(
  finalp + tm_layout(outer.margins=c(.02,.05,.01,.03)),
  dpi=100,
          insets_tm = insetmap,
          insets_vp=vp,
  filename="BB_sockeye_rivers_inset.png")