R Notebook

Rippin’ Dope GIS shit

load packagess

library(dplyr)
library(purrr)
library(raster)
library(lidR)

You can’t pass a list of rasters to raster::merge because it expects you to list out the names of all the rasters. Enter do.call(), which, given a function name, f, and a list of arguments, L will generate the R code you want, then run it. See below tweaked example taken from the examples section of the raster::merge docs.

# assume r1 and r2 are rasters
x <- list(r1, r2)
# this will create the string:
  "merge(r1, r2)", then execute it.
m <- do.call(merge, x)

I also wrapped this stuff into a loop so we don’t have too many variables dirtying up our workspace

# pass the integers of the files we want
allrasters <-c(3,8,9,10) %>% 
  # make a list of the files
  paste0("LTEP_file", ., ".TIF") %>%
  # make a list where each item is one of the files I want
  map(., raster) %>%
  # merge them all into one object
  # see note above about `do.call`
  do.call(merge, .)

# Reprojecting the one large raster
# NOTE: ANDERS isn't sure what is supposed to be happening here...
print(paste("PROJECTION BEFORE:", projection(allrasters)))

[1] "PROJECTION BEFORE: +proj=lcc +lat_0=45.3333333333333 +lon_0=-120.5 +lat_1=47.3333333333333 +lat_2=45.8333333333333 +x_0=500000.0001016 +y_0=0 +ellps=GRS80 +units=us-ft +no_defs"

projection(allrasters)<- 2927
print(paste("PROJECTION AFTER:", projection(allrasters)))

[1] "PROJECTION AFTER: +proj=lcc +lat_0=45.3333333333333 +lon_0=-120.5 +lat_1=47.3333333333333 +lat_2=45.8333333333333 +x_0=500000.0001016 +y_0=0 +ellps=GRS80 +units=us-ft +no_defs"

st_crs(allrasters)$input

[1] "NAD83(HARN) / Washington South (ftUS)"

st_crs(allrasters)

Error in st_crs(allrasters) : could not find function "st_crs"

# read shpfile but readOGR doesn't work for anders??
# shpfile<- readOGR(dsn="~/Downloads/court/", layer= "EXU_MappingPolygons")
shpfile<- readRDS("./shpfile.RDS")

# give an error if `shpfile` isn't the right type of Class
if (class(shpfile) != "SpatialPolygonsDataFrame") {
  print("shpfile isn't a SpatialPolygonsDataFramem, but it should be!")
}
plot(shpfile)

head(shpfile)

crs(allrasters, asText=TRUE)

[1] "+proj=lcc +lat_0=45.3333333333333 +lon_0=-120.5 +lat_1=47.3333333333333 +lat_2=45.8333333333333 +x_0=500000.0001016 +y_0=0 +ellps=GRS80 +units=us-ft +no_defs"

st_crs(allrasters)$input

[1] "NAD83(HARN) / Washington South (ftUS)"

crs(shpfile, asText=TRUE)

[1] "+proj=lcc +lat_0=45.3333333333333 +lon_0=-120.5 +lat_1=47.3333333333333 +lat_2=45.8333333333333 +x_0=500000.0001016 +y_0=0 +ellps=GRS80 +units=us-ft +no_defs"

st_crs(shpfile)$input

[1] "NAD83(HARN) / Washington South (ftUS)"

crs(allrasters, asText=TRUE)

[1] "+proj=lcc +lat_0=45.3333333333333 +lon_0=-120.5 +lat_1=47.3333333333333 +lat_2=45.8333333333333 +x_0=500000.0001016 +y_0=0 +ellps=GRS80 +units=us-ft +no_defs"

crs(shpfile, asText=TRUE)

[1] "+proj=lcc +lat_0=45.3333333333333 +lon_0=-120.5 +lat_1=47.3333333333333 +lat_2=45.8333333333333 +x_0=500000.0001016 +y_0=0 +ellps=GRS80 +units=us-ft +no_defs"

save out objects for later

saveRDS(shpfile, "shpfile.RDS")
saveRDS(allrasters, "allrasters.RDS")

Part II

load old objects

shpfile <- readRDS("shpfile.RDS")

map city, bitch

filter_shp <- function(shpfile, val_list) {
  shpfile[shpfile@data$SppComp %in% val_list,]
}
# filter for only two desired conditions
# make a named list so we can access the indv. SpationPolygondsDataFrames as needed
types <- c("Douglas Fir", "Pioneer")
names(types) <- types
types

  Douglas Fir       Pioneer 
"Douglas Fir"     "Pioneer" 

shp <- map(types, filter_shp, shpfile=shpfile)

plot(out$`Douglas Fir`, add=T, col= "blue")

Error in polypath(x = mcrds[, 1], y = mcrds[, 2], border = border, col = col,  : 
  plot.new has not been called yet

OMG! ggspatial! here’s a ggpstial tutorial

ggplot() +
    layer_spatial(data = dougfir)

shp_pnr <- filter_shp(shpfile, shpfile@data$SppComp == "Pioneer")