Setup

library(raster)

## Loading required package: sp

library(tmap)

## Registered S3 methods overwritten by 'stars':
##   method             from
##   st_bbox.SpatRaster sf  
##   st_crs.SpatRaster  sf

library(sf)

## Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.5     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.0.2     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x tidyr::extract() masks raster::extract()
## x dplyr::filter()  masks stats::filter()
## x dplyr::lag()     masks stats::lag()
## x dplyr::select()  masks raster::select()

load("counties_us.Rdata")

Use the Viewer

Go to https://www.mrlc.gov/viewer/

Pick out something like Thurston County and download. It gives you a tiff and an xml file.

tclc = raster("NLCD2019.tiff")
is.na(tclc) = 0
tclc[tclc > 21 & tclc < 31] = 1
tclc[tclc != 1] = 2
tm_shape(tclc) + tm_raster(style = "cat")

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

unique(tclc)

## [1] 1 2

NLCD Database

Go to https://www.mrlc.gov/data?f%5B0%5D=category%3Aland%20cover and download the conus 2019 land cover data. It’s about 25 GB.

There are several files. Point to the .img.

nl = raster("nlcd_2019_land_cover_l48_20210604.img")
tm_shape(nl) + tm_raster()

## stars_proxy object shown at 1243 by 805 cells.

## Warning: Duplicated levels found. They have been omitted

Crop and Mask to Thurston County

Note the difference in the two maps.

thurston = counties_us %>% filter(STATEFP ==53, NAME=="Thurston")

# Fix CRS

thecrs = crs(nl)
thurston = st_transform(thurston,thecrs)

thlcc = crop(nl,thurston)

## Warning in matrix(ratvalues, nrow = length(ratvalues)/length(ratnames)): data
## length [1791] is not a sub-multiple or multiple of the number of rows [255]

thlcm = mask(thlcc,thurston)

tm_shape(thlcc) + tm_raster(style = "cat")

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

tm_shape(thlcm) + tm_raster(style = "cat")

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