<jupyter_start><jupyter_text>*Technology transfer for rapid family forest assessments and stewardship planning*  Generate forest indicators from remotely sensed and climate data This notebook facilitates the first steps in developing a predictive model to characterize forest conditions using lidar, high-resolution imagery, and field-based inventory data. The follow methods generate a wide variety forest indicators that will be used in predicting forest type maps for a sample area within Clackamas County Oregon.  This is a small subset of our larger study region and is intended for methodological demonstration.<jupyter_code>#load libraries
library(sp)
library(raster)
library(rgdal)
library(rgeos)
library(dismo)
library(maptools)
library(ggmap)
#set plot size
options(repr.plot.width  = 5,
        repr.plot.height = 4)<jupyter_output><empty_output><jupyter_text>Create spatial layer of FIA plots within the area of interest The Forest Inventory and Analysis (FIA) program of the USDA Forest Service collects inventory data for all forested lands in the U.S., regardless of ownership or availability for harvesting. Data are collected on permanent field plots and consist of a wide variety of attributes, spanning ecological, topographic, and structural characteristics.  Due to privacy requirements specified in the amendments of the FOod Security act of 1985, publicily available FIA plot locations are approximate and can be downloaded via __[FIA Data Mart](https://apps.fs.usda.gov/fia/datamart/datamart_access.html)__. The metrics generated for the entire study area will be based on non-fuzzed plot coordinates, made available through a partnership with FIA researchers at the PNW Forest Service Research Station. For demonstration purposes, fuzzed FIA plots are utilized for this sample case. The PLOT table of the Oregon FIA inventory database was exported as a text file is used to plot the locations of forested plots in Clackamas County.<jupyter_code>#read in fia plot data text filea
fia <- read.csv("./Data/FIA/FIA_OregonStatePlots.csv", header=TRUE, sep=",")
#plot points using lat/lon coordinates
pts <- SpatialPointsDataFrame(fia[,c("LON", "LAT")], fia[c(5:9, 20:21)])
pts_m <- data.frame(fia[c(5:9,20:21)])
#print first five rows
head(pts,n=5)
#set coordinate system
crs(pts) <- NA
crs(pts) <- CRS("+proj=longlat +datum=NAD83")
#reproject to NAD83_2011 Oregon Lambert
utm <- "+proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
nad83 <- "+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=399999.9999984 +y_0=0 +ellps=GRS80 +units=ft +no_defs"
pts <- spTransform(pts, nad83)<jupyter_output><empty_output><jupyter_text>Create plot polygons In this step we remove plots that fall outside of our study region, plots that were inventoried after the lidar acquistion, and plots that had some time of disturbance event occur since the inventory date.  We also add a buffer around the plot center to match the plot dimensions and create a hexagonal grid within each buffer for imputation purposes.   Remove plots that fall outside the study area.<jupyter_code>#read in study area polygon
#reproject to NAD83_2011 Oregon Lambert
area_o <- readOGR("./Data/Misc/Study_Area.shp")
area <- spTransform(area_o, nad83)
area_m <-spTransform(area_o, "+proj=longlat +datum=NAD83")
area_df <- fortify(area_m)

#plot spatial layers
plot(area)
plot(pts, add=T)
#clip points to study area and attach lidar acquistion year
pts_clip_lidar <- intersect(pts, area)
names(pts_clip_lidar)[10]<-"lidarYear"
pts_clip_lidar$lidarYear <-as.numeric(as.character(pts_clip_lidar$lidarYear))
head(pts_clip_lidar)
#plot spatial layers
plot(area)
plot(pts_clip_lidar, add=T)<jupyter_output><empty_output><jupyter_text>Delete plots that were inventoried *after* the lidar acquisition date<jupyter_code>#delete plots with inventory year > lidar acquistion year
pts_clip_lidar_filter <- pts_clip_lidar[pts_clip_lidar$lidarYear > pts_clip_lidar$INVYR, ]
paste(nrow(pts_clip_lidar_filter), "plots remain after filter")
head(pts_clip_lidar_filter, n=1)<jupyter_output><empty_output><jupyter_text>Delete duplicate plots - remove earliest inventory instance<jupyter_code>plot_year <- aggregate(INVYR ~ PLOT, pts_clip_lidar_filter, max)
paste(nrow(plot_year), "plots remain after filter")
head(plot_year, n=1)
pts_clip_lidar_filter <- merge(pts_clip_lidar_filter, plot_year, by.x=c("INVYR", "PLOT"), by.y=c("INVYR", "PLOT"), all.x = FALSE)
head(pts_clip_lidar_filter, n=1)
paste(nrow(pts_clip_lidar_filter), "plots remain after filter")
class(pts_clip_lidar_filter)<jupyter_output><empty_output><jupyter_text>Add a buffer of 180 ft to each point center<jupyter_code>#add buffer of 180 feet to points 
pts_clip_lidar_filter_buff <- gBuffer(pts_clip_lidar_filter, width=180, byid=TRUE)
pts_clip_lidar_filter_buff <- SpatialPolygonsDataFrame( pts_clip_lidar_filter_buff, data=pts_clip_lidar_filter_buff@data )
head(pts_clip_lidar_filter_buff, n=1)
print(nrow(pts_clip_lidar_filter_buff))<jupyter_output>[1] 274<jupyter_text>Remove plots that have experienced a disturbance event since inventory Disturbance data is sourced from the __[Global Forest Watch](http://data.globalforestwatch.org/)__ and is a product of the alogrithms that capture, label, and map land cover change.<jupyter_code>#load landtrendr distrubance (year of disturbance) raster
disturb <- raster("./Data/Landtrendr/ForestLoss_year.tif")
#project raster to nad83
beginCluster(4) 
disturb<- projectRaster(disturb, crs=nad83)
endCluster()
#extract max disturbance year for each plot polygon
beginCluster(4) 
pts_clip_lidar_filter_buff_disturb <- extract(disturb, pts_clip_lidar_filter_buff, fun=max, df=TRUE)
endCluster()
#add plot id
pts_clip_lidar_filter_buff_disturb$plot <- pts_clip_lidar_filter_buff$PLOT
#clean up column names
names(pts_clip_lidar_filter_buff_disturb) <- c('ID', 'distYear', 'plot')
#convert NA distYear values to 0
pts_clip_lidar_filter_buff_disturb$distYear[is.na(pts_clip_lidar_filter_buff_disturb$distYear)] <- 0 
print(nrow(pts_clip_lidar_filter_buff_disturb))
head(pts_clip_lidar_filter_buff_disturb, n=5)
#merge back with rest of dataframe
pts_clip_lidar_filter_buff_disturb <- merge(pts_clip_lidar_filter_buff, pts_clip_lidar_filter_buff_disturb, by.x='PLOT', by.y="plot")
#print header
head(pts_clip_lidar_filter_buff_disturb)
#delete plots with disturbance > FIA inv year
pts_clip_lidar_filter_buff_disturb_filter <- pts_clip_lidar_filter_buff_disturb[pts_clip_lidar_filter_buff_disturb$distYear < pts_clip_lidar_filter_buff_disturb$INVYR, ]<jupyter_output><empty_output><jupyter_text>Create hexagon grid within each FIA buffered plot Make_grid function sourced from __[<NAME>](http://strimas.com/spatial/hexagonal-grids/)__<jupyter_code>make_grid <- function(x, cell_diameter, cell_area, clip = FALSE) {
  if (missing(cell_diameter)) {
    if (missing(cell_area)) {
      stop("Must provide cell_diameter or cell_area")
    } else {
      cell_diameter <- sqrt(2 * cell_area / sqrt(3))
    }
  }
  ext <- as(extent(x) + cell_diameter, "SpatialPolygons")
  projection(ext) <- projection(x)
  # generate array of hexagon centers
  g <- spsample(ext, type = "hexagonal", cellsize = cell_diameter, 
                offset = c(0.5, 0.5))
  # convert center points to hexagons
  g <- HexPoints2SpatialPolygons(g, dx = cell_diameter)
  # clip to boundary of study area
  if (clip) {
    g <- gIntersection(g, x, byid = TRUE)
  } else {
    g <- g[x, ]
  }
  # clean up feature IDs
  row.names(g) <- as.character(1:length(g))
  return(g)
}
#make grid - 1/4 acre hexagons
fia_grid <- make_grid(pts_clip_lidar_filter_buff_disturb_filter, cell_area = 10890, clip = FALSE)
#intersect with plot buffer layer to get attributes
fia_grid <- intersect(fia_grid, pts_clip_lidar_filter_buff_disturb_filter)
#add unique id field to fia_grid layer
fia_grid$hex_id <- 1:nrow(fia_grid)<jupyter_output><empty_output><jupyter_text>Export Data<jupyter_code>#export csv of plot ids, lidar years, and disturbance years to be used to create FVS input
write.csv(pts_clip_lidar_filter_buff_disturb_filter, "./FVS/FVS_input.csv",  row.names=FALSE)
#export shapefile of plots (w/buffer)
writeOGR(obj=pts_clip_lidar_filter_buff_disturb_filter, dsn="./Data/FIA", layer="FIA_plots", driver="ESRI Shapefile",
        overwrite_layer=TRUE)
#export shapefile of plots with hexes
#writeOGR(obj=fia_grid, dsn="./Data/FIA", layer="FIA_plots_hex", driver="ESRI Shapefile", overwrite_layer=TRUE)<jupyter_output><empty_output><jupyter_text>Calculate predictor variables from LANDSAT remotely sensed imagery For over 40 years, a series of __[LANDSAT](https://earthexplorer.usgs.gov/)__ satellites have recorded a temporal record of multispectral data of the Earth’s surface at a global level. At a moderate resolution, 30m x 30m, Landsat imagery is a valuable resource for regional forestry applications.  It is expected that a substantial effort will be required to gather suitable cloud-free LANDSAT images for a given year and to follow procedures such as histogram equalization to normalize images across years. However, Analysis Ready Data (ARD) -- imagery processed to significantly reduce the magnitude of data processing by users, will be utilized whenever possible. Remotely sensed imagery will be summarized within areal units defined by a hexagonal tessellation overlaid on the FIA plots. To respect USFS concerns to conceal FIA plot locations, no data layers with a resolution finer than this tessellation will be taken off-site. This hexagonal tessellation will also be reproduced across the entire LiDAR footprints, and these hexagons will be the level at which forest type predictions are made.<jupyter_code>#create a list of all the landsat bands from landsat directory
landsat_names <-  Sys.glob(file.path("./Data/Landsat/sampleArea", "*.tif"))
images <- sapply(landsat_names, raster) 
#create raster brick from .tif files
ls8<-stack(images)
## tasseled cap coefficients
refB = c(0.3029, 0.2786, 0.4733, 0.5599, 0.508, 0.1872) 
refG = c(-0.2941, -0.243, -0.5424, 0.7276, -0.0713,-0.1608)
refW = c(0.1511, 0.1973, 0.3283, 0.3407, -0.7117, -0.4559)
#function to calculate tasseled cap and band ratio metrics 'calc_landsMets'
calc_landMets <- function(r) {
    b <- r[[2]] * refB[1] + r[[3]] * refB[2] + r[[4]] * refB[3] + r[[5]] * refB[4] + r[[6]] * refB[5] + r[[7]] * refB[6]
    g <- r[[2]] * refG[1] + r[[3]] * refG[2] + r[[4]] * refG[3] + r[[5]] * refG[4] + r[[6]] * refG[5] + r[[7]] * refG[6]
    w <- r[[2]] * refW[1] + r[[3]] * refW[2] + r[[4]] * refW[3] + r[[5]] * refW[4] + r[[6]] * refW[5] + r[[7]] * refW[6]
    ndvi <- (r[[5]] - r[[4]])/(r[[5]] + r[[4]])
    gndvi <- (r[[5]] - r[[3]])/(r[[5]] + r[[3]])
    nbr <- (r[[5]] - r[[7]])/(r[[5]] + r[[7]])
    g_ratio <- r[[5]] - r[[3]]
    b1 <- r[[1]]
    b2 <- r[[2]]
    b3 <- r[[3]]
    b4 <- r[[4]]
    b5 <- r[[5]]
    b7 <- r[[7]]
    l_mets <- stack(b, g, w, ndvi, gndvi, nbr, g_ratio, b1, b2, b3, b4, b5, b7)
    names(l_mets) <- c('brightness', 'greenness', 'wetness', 'ndvi', 'gndvi', 'nbr', 'g_ratio', 'b1', 'b2', 'b3', 'b4', 'b5', 'b7')  
return(l_mets)
}
#calculate metrics using 'calc_landsMets' function on landsat stack
l_mets <- calc_landMets(ls8)
beginCluster(4)
l_mets<- projectRaster(l_mets, crs=nad83)
endCluster()
#calculate mean zonal statistics for tasseled cap and indices for each plot hexagon
ext_mets <- extract(l_mets, fia_grid, fun = mean, na.rm = T, df=T)
ext_mets <- merge(fia_grid, ext_mets, by.x='hex_id', by.y="ID")
head(ext_mets)
#write landsat variables to csv
write.csv(ext_mets, "./Classification/land_mets.csv",  row.names=FALSE)<jupyter_output><empty_output><jupyter_text>Create predictor variables from Climate Data To estimate climate attributes for our study area, we will be using __[ClimateWNA](https://sites.ualberta.ca/~ahamann/data/climatewna.html)__ software. This application provides historical and projected climate data for North America, including monthly, seasonal, and annual attributes, as well as derived metrics such as Growing Degree Days, Potential Evapotranspiration, etc. that are useful for vegetation modeling. A ClimateWNA executable file is included within this repository. In the next couple steps, we create the input needed to derive our selected attributes.  Once that input table is created, the ClimateWNA executable is run and will pop up in an external window.  Set the created text file as input and run the application to create an additional text file containing attributes for our study area.<jupyter_code>#create csv for climateWNA application
head(fia, n=1)
fia_vars <- c("PLOT", "LON", "LAT", "ELEV")
fia_coords <- fia[fia_vars]
fia_coords <- merge(fia_coords, pts_clip_lidar_filter_buff_disturb_filter)
fia_vars2 <- c("PLOT", "STATECD", "LAT", "LON", "ELEV")
fia_coords <- fia_coords[fia_vars2]
names(fia_coords) <- c("ID1", "ID2", "lat", "lon", 'el')
write.csv(fia_coords, "./Data/Climate/climate_input.csv",  row.names=FALSE)<jupyter_output><empty_output><jupyter_text>This next cell will open the ClimateWNA application. Select the 'climate_input.csv" created in the last cell as input.  Save the output as 'climate_input_year_2014Y.csv".<jupyter_code>#Open climateWNA application
wd<-getwd()
system(paste(file.path(wd, "Data/Climate/ClimateWNA_v5.51.exe")))

#select climate_input.csv as input
#output file to data/climate
#read in output created in last step from CLIMATEWNA
climate <- read.csv("./Data/Climate/climate_input_year_2014Y.csv", header=TRUE, sep=",")
head(climate)
#keep variables we are intrested in
climate_vars <- c("ID1", "ID2", "MAT", "MAP", "MWMT", "MSP", "MCMT", "PAS", "DD5", "Eref", "NFFD", "CMD")
climate_mets <- climate[climate_vars]
#rename ID1 and ID2 fields
names(climate_mets)[1]<-"plot"
names(climate_mets)[2]<-"state"
#write to csv file
write.csv(climate_mets, "./Classification/climate_mets.csv",  row.names=FALSE)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11250269920103649, 'fraction_numerical': 0.02317713956668826, 'mean_word_length': 5.857847976307996, 'pattern_counts': {'":': 0, '<': 125, '<?xml version=': 0, '>': 52, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '248874', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5145, 'n_tokens_neox': 4773, 'n_words': 1409, 'path': 'notebooks/ForestMappingMetrics_GIS.ipynb', 'repo_name': 'Ecotrust/ForestMapping', 'size': '13906'}