Author

Carolyn Koehn

Review from session 16

Load libraries:

Code 
library(sf)
library(terra)
library(tidyverse)
library(tmap)
library(tree)
library(randomForest)

Load data:

Code 
download_unzip_read <- function(link){
  tmp <- tempfile()
  download.file(link, tmp)
  tmp2 <- tempfile()
  unzip(zipfile=tmp, exdir=tmp2)
  shapefile.sf <- read_sf(tmp2)
}

### FS Boundaries
fs.url <- "https://data.fs.usda.gov/geodata/edw/edw_resources/shp/S_USA.AdministrativeForest.zip"
fs.bdry <- download_unzip_read(link = fs.url)

### CFLRP Data
cflrp.url <- "https://data.fs.usda.gov/geodata/edw/edw_resources/shp/S_USA.CFLR_HPRP_ProjectBoundary.zip"
cflrp.bdry <- download_unzip_read(link = cflrp.url)

wildfire_haz <- rast("/opt/data/data/assignment01/wildfire_hazard_agg.tif")

cejst <- st_read("/opt/data/data/assignment01/cejst_nw.shp", quiet=TRUE) %>%
  filter(!st_is_empty(.))

Check validity:

Code 
all(st_is_valid(fs.bdry))
[1] FALSE
Code 
all(st_is_valid(cflrp.bdry))
[1] FALSE
Code 
fs.bdry <- st_make_valid(fs.bdry)
cflrp.bdry <- st_make_valid(cflrp.bdry)

Check alignment:

Code 
st_crs(wildfire_haz) == st_crs(fs.bdry)
[1] FALSE
Code 
st_crs(wildfire_haz) == st_crs(cflrp.bdry)
[1] FALSE
Code 
st_crs(wildfire_haz) == st_crs(cejst)
[1] FALSE
Code 
fs.bdry_proj <- st_transform(fs.bdry, crs = st_crs(wildfire_haz))
cflrp.bdry_proj <- st_transform(cflrp.bdry, crs = st_crs(wildfire_haz))
cejst_proj <- st_transform(cejst, crs = st_crs(wildfire_haz))

Subset to relevant geographies:

Code 
fs.bdry_sub <- fs.bdry_proj[cejst_proj, ]
cflrp.bdry_sub <- cflrp.bdry_proj[cejst_proj, ]

cejst_sub <- cejst_proj[fs.bdry_sub, ]

Select relevant attributes:

Code 
cejst_sub <- cejst_sub %>%
  select(GEOID10, LMI_PFS, LHE, HBF_PFS)

Extract wildfire risk:

Code 
wf_risk <- terra::extract(wildfire_haz, cejst_sub, fun=mean)

cejst_sub$WHP_ID <- wf_risk$WHP_ID

CFLRP T or F:

Code 
cflrp <- apply(st_intersects(cejst_sub, cflrp.bdry_sub, sparse = FALSE), 1, any)

cejst_sub$CFLRP <- cflrp

Comparing/Predicting CFLRP tracts

Data preparation

Code 
cejst_mod <- cejst_sub %>%
  st_drop_geometry(.) %>%
  na.omit(.)

cejst_mod[, c("LMI_PFS", "LHE", "HBF_PFS", "WHP_ID")] <- scale(cejst_mod[, c("LMI_PFS", "LHE", "HBF_PFS", "WHP_ID")])

Logistic regression

Code 
logistic.global <- glm(CFLRP ~ LMI_PFS + LHE + HBF_PFS + WHP_ID,
                       family = binomial(link = "logit"),
                       data = cejst_mod)
summary(logistic.global)

Call:
glm(formula = CFLRP ~ LMI_PFS + LHE + HBF_PFS + WHP_ID, family = binomial(link = "logit"), 
    data = cejst_mod)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.66462    0.14347  -4.632 3.62e-06 ***
LMI_PFS      0.17284    0.16881   1.024    0.306    
LHE         -0.25551    0.15791  -1.618    0.106    
HBF_PFS     -0.01511    0.16081  -0.094    0.925    
WHP_ID       0.88902    0.16785   5.297 1.18e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 330.74  on 255  degrees of freedom
Residual deviance: 290.03  on 251  degrees of freedom
AIC: 300.03

Number of Fisher Scoring iterations: 4

Classification Tree

Code 
library(tree)
cejst_mod$CFLRP <- as.factor(ifelse(cejst_mod$CFLRP == 1, "Yes", "No"))
tree.model <- tree(CFLRP ~ LMI_PFS + LHE + HBF_PFS + WHP_ID, cejst_mod)
plot(tree.model, type = "uniform")
text(tree.model, pretty=0)

Random Forest

Code 
library(randomForest)
class.model <- CFLRP ~ .
rf2 <- randomForest(formula = class.model, cejst_mod[,-1])
varImpPlot(rf2)