A Youtube Video

First, watch this video. It’s not me!

https://www.youtube.com/watch?v=yYWdxExabHo

Here’s the link to my video walkthrough.

https://www.youtube.com/watch?v=RD8ItZBLIvA

Setup

library(sf)
library(raster)
library(tmap)
library(tmaptools)
library(tidyverse)
load("counties_us.Rdata")

Read Rasters

r46_123 = raster("n46_w123_1arc_v3.tif")
r46_124 = raster("n46_w124_1arc_v3.tif")
r47_123 = raster("n47_w123_1arc_v3.tif")
r47_124 = raster("n47_w124_1arc_v3.tif")

Look at Rasters

Combine first.

whole_dem = raster::merge(r46_123,r47_123,r47_124,r46_124)
tm_shape(whole_dem) + tm_raster(midpoint = NA)

## stars object downsampled to 1000 by 1000 cells. See tm_shape manual (argument raster.downsample)

Get Thurston County

Thurston = counties_us %>% filter(STATEFP == 53,
                                  NAME == "Thurston")
thecrs = crs(r46_123)
Thurston = st_transform(Thurston, thecrs)

Compare

Look at the dem data and Thurston County.

tm_shape(whole_dem)  + tm_raster(midpoint = NA) +
tm_shape(Thurston) + 
  tm_borders(lwd = 2, col = "blue")

## stars object downsampled to 1000 by 1000 cells. See tm_shape manual (argument raster.downsample)

Crop it

Trim the dem data down to fit Thurston County.

dem_crop = crop(whole_dem,Thurston)

tm_shape(dem_crop)  + tm_raster(midpoint = NA) +
tm_shape(Thurston) + 
  tm_borders(lwd = 2, col = "blue")

## stars object downsampled to 1530 by 654 cells. See tm_shape manual (argument raster.downsample)

Mask

Eliminate the non-Thurston data.

dem_mask = mask(dem_crop,Thurston) 

tm_shape(dem_mask)  + tm_raster(midpoint = NA) +
tm_shape(Thurston) + 
  tm_borders(lwd = 2, col = "blue")

## stars object downsampled to 1530 by 654 cells. See tm_shape manual (argument raster.downsample)

Continuous Raster

Use style = “cont” in call to tm_raster.

dem_mask = mask(dem_crop,Thurston) 

tm_shape(dem_mask)  + tm_raster(style = "cont",midpoint = NA) +
tm_shape(Thurston) + 
  tm_borders(lwd = 2, col = "blue")

## stars object downsampled to 1530 by 654 cells. See tm_shape manual (argument raster.downsample)

dem_mask = mask(dem_crop,Thurston) 

tm_shape(dem_mask)  + tm_raster(style = "cont",midpoint = NA,
                                palette = "-Blues") +
tm_shape(Thurston) + 
  tm_borders(lwd = 2, col = "black")

## stars object downsampled to 1530 by 654 cells. See tm_shape manual (argument raster.downsample)

Getting DEM Data