exec/test_predict.R

In bmcnellis/RSFIA: FIA Tree Mortality Analysis

library(ncdf4)
library(raster)
library(ncdf4.helpers)
df0 <- FIA_mortality_with_explanatory

for (i in c('MAP', 'MCMT', 'MWMT')) {
  v0 <- i
  mapNC <- ncdf4::nc_open(paste0(DATA_DIRECTORY, '/', v0, '.nc'))
  MAP <- ncdf4::ncvar_get(mapNC, varid = v0)
  nm <- ncdf4.helpers::nc.get.dim.names(mapNC)
  east <- ncdf4::ncvar_get(mapNC, 'easting')
  north <- ncdf4::ncvar_get(mapNC, 'northing')

  r0 <- raster(t(MAP),
               xmn=min(east),
               xmx=max(east),
               ymn=min(north),
               ymx=max(north),
               crs=CRS("+proj=lcc +ellps=WGS84 +datum=WGS84 +units=m +x_0=0 +y_0=0 +lon_0=-95.0 +lat_0=0 +no_defs+ towgs84=0,0,0")
  )
  r0 <- projectRaster(from = r0, crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
  nc_close(mapNC)

  t0 <- raster::extract(r0, data.frame(df0$LON, df0$LAT))
  var <- data.frame(t0, df0$LON, df0$LAT)
  colnames(var) <- c(v0, 'LON', 'LAT')
  df0 <- dplyr::left_join(df0, var, by = c('LON', 'LAT'))

}
save(mort_df_with_CMIP5, file = paste0(DATA_DIRECTORY, '/mort_df_with_CMIP5.rda'))
rm(list = ls())
DATA_DIRECTORY <- '/media/bem/hot_storage/RSFIA_data'
load(paste0(DATA_DIRECTORY, '/mort_df_with_CMIP5.rda'))
load(paste0(DATA_DIRECTORY, '/stratified_mort_rf.rda'))
load(paste0(DATA_DIRECTORY, '/stratified_tree_training.rda'))
colnames(mort_df_with_CMIP5)[which(colnames(mort_df_with_CMIP5) == 'LON')] <- 'lon'
colnames(mort_df_with_CMIP5)[which(colnames(mort_df_with_CMIP5) == 'LAT')] <- 'lat'
mort_df_with_CMIP5 <- mort_df_with_CMIP5[, c('lon', 'lat', 'MAP', 'MCMT', 'MWMT')]

pred_trees <- dplyr::left_join(trees, mort_df_with_CMIP5, by = c('lon', 'lat'))
pred_trees <- pred_trees[which(pred_trees$STATUSCD == 1), ]
pred_trees$MAP <- pred_trees$MAP / 365
cc <- c('PLT_CN', 'SPCD', 'DIA', 'HT',
        'BLDFIE', 'CECSOL', 'CLYPPT', 'CRFVOL', 'PHIHOX', 'SNDPPT', 'BDTICM',
        'ELEV', 'is_public', 'SLOPE', 'ASPECT', 'MAP', 'MCMT', "MWMT")
pred0 <- pred_trees[, cc]
rm(trees, mort_df_with_CMIP5, pred_trees, cc)
colnames(pred0)[which(colnames(pred0) == 'DIA')] <- 'PREVDIA'
colnames(pred0)[which(colnames(pred0) == 'MAP')] <- 'prcp'
colnames(pred0)[which(colnames(pred0) == 'MCMT')] <- 'tmin'
colnames(pred0)[which(colnames(pred0) == 'MWMT')] <- 'tmax'

new_data <- predict(stratified_mort_rf, newdata = pred0, type = 'response')
new_data <- as.logical(new_data)
pred0 <- cbind(pred0, new_data)
load(paste0(DATA_DIRECTORY, '/stratified_tree_training.rda'))