Elevations.Rmd

---
title: "Elevations"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message=FALSE, warning=FALSE)
```

# 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](https://cran.r-project.org/web/packages/tmap/vignettes/tmap-getstarted.html) a lot. Super easy to layer map features together like I want.

### Credits

My example is adapted from [Visualize Spatial Data](http://www.wvview.org/spatial_analytics/Visualizing_Spatial_Data/_site/Visualize_Spatial_Data.html) in a [Open Source Spatial Analytics](http://www.wvview.org/Open_Source_Spatial_Analytics.html) course by [Dr. Aaron Maxwell](http://www.wvview.org/Prof_Maxwell.html) at West Virginia University.


## Step 0 Load the needed packages

```{r}
library(raster)
library(tmap)
library(grid)
```

## 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.

```{r}
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.

```{r eval=FALSE}
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")
```

```{r}
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.

```{r}
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.
```{r}
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.

```{r}
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.

```{r}
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:

```{r eval=FALSE}
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.

```{r fig.width=10}
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.

```{r eval=FALSE}
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.
```{r}
load("raster_data/bbshp.RData")
```

```{r  fig.width=10}
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.

```{r fig.width=10}
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.

```{r eval=FALSE}
togiak <- drawLine()
crs(togiak) <- hill.crs
togiak$ID <- "Togiak"
```

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

```{r eval=FALSE}
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.
```{r}
load("data/bbrivers.RData")
```

```{r fig.width=10}
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.

```{r fig.width=10}
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.

```{r eval=FALSE}
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.

```{r eval=FALSE}
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.
```{r}
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.
```{r fig.width=10}
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.

```{r results="hide", eval=FALSE}
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.

```{r results="hide", eval=FALSE}
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.


```{r}
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)
# Transform the Bristol Bay box and AK outline to Lambert
BB.box2 <- sp::spTransform(BB.box, lambert)
BB.box2$name <- "Bristol Bay"
akshp2 <- sp::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`.
```{r}
xy <- tmap::st_bbox(akshp2)
asp2 <- (xy$xmax - xy$xmin)/(xy$ymax - xy$ymin)
h <- 0.25
w <- asp2 * h
vp <- grid::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")
```

![](BB_sockeye_rivers_inset.png)