## R-based analysis script for generating the datasets used in this study. This script processes the large postgreSQL query ## file, which joins and flattens the FIA database "PLOT","COND", and "TREE" tables. To avoid more technical and error-prone #SQL querying,basic steps of data filtering are performed here rather than during the query phase. This script also prepares calculated ## data fields such as competition indices and relative density. The final data objects from this script are directly used ## in generalized linear mixed modeling. #Please note that the FIA database maintains fields in U.S. units rather than metric. For the purposes #of this study, units were generally inconsequential due to centering and scaling of predictors #Where necessary, volume increment response units were manually converted in scripts used to generate graphics ##import flattened SQL query of the FIADB: library(readr) september2019_fia_base <- read_csv("E:/Dropbox/ASCC/tree vigor and mortality/september2019_fia_base_4.csv", col_types = cols(carbon_ag = col_double(), dia = col_double(), growcfal = col_number(), ht = col_number(), prevdia = col_double(), totage = col_integer(), tpa_unadj = col_number(), volcfnet = col_number(), dstrbcd1=col_fact(),locale = locale(encoding = "WINDOWS-1252"))) View(september2019_query_3) library(data.table) ##Main plot filtering criteria, excluding PLOT.remper bounds. Dead trees were dropped from the analysis. Microplot trees are excluded ##so that queried plots are likely to contain trees with growth fields, as well as making the calculation of point competition indices ##simpler## aspen_fia=subset(september2019_fia_base, trtcd1==0 #no anthropogenic disturbance & plot_status_cd==1 # & cond_status_cd==1 # forested condition ("condition" syn. with stand) & condprop_unadj==1 #single condition class plots & statuscd!=0 #exclude dropped trees & statuscd!=3) #exclude harvest trees aspen_fia=subset(aspen_fia,qa_status==1 #regular production plots |qa_status==7) #FIA quality assurance/ quality control check plots aspen_fia=subset(aspen_fia,dstrbcd1==0 #undisturbed |dstrbcd1==10 #insect |dstrbcd1==12 #insect |dstrbcd1==20 #disease |dstrbcd1==22 #disease |dstrbcd1==54 #drought |dstrbcd1==60) #competition aspen_fia=as.data.table(aspen_fia) aspen_fia[,plot_id:=as.factor(make.names(paste(statecd,plot)))] ##extract PLOT.remper field (remeasurement period), #and filter via table join out plots where remper is >8 years and < 12 years aspen_fia[,":="(dia=as.numeric(dia),ht=as.numeric(ht),cclcd=as.integer(cclcd),cr=as.integer(cr),tpa_unadj=as.numeric(tpa_unadj),growcfal=as.numeric(growcfal),aspect=as.integer(aspect),slope=as.integer(slope),stdage=as.integer(stdage),measyear=as.integer(measyear),remper=as.integer(remper)),] remper_by_plot=aggregate(remper~plot_id,data=aspen_fia,mean) remper_by_plot$remper_by_plot=remper_by_plot$remper remper_by_plot$remper=NULL aspen_fia=as.data.table(merge(remper_by_plot,aspen_fia,by="plot_id")) aspen_fia=as.data.table(subset(aspen_fia,remper_by_plot>8 & remper_by_plot<12)) rm(remper_by_plot) ##identify ingrowth trees based on absence of TREE.prevdia entry and presence of TREE.growcfal entry--not critical to analysis## ##calculate tree-level density statistics for later use in plot and subplot-level summaries. Not all variables have been of use, ## and planned rescaling of continuous variables to Z-scores prior to analysis reduces the need for unit conversion aspen_fia[,ba:=(as.numeric(dia)^2)*0.005454,by=list(plot_id,kindcd,subp,tree)] aspen_fia[,dia_cm:=dia*2.54,by=list(plot_id,kindcd,subp,tree)] aspen_fia[,ba_met:=((dia_cm)^2)*0.00007854,by=list(plot_id,kindcd,subp,tree,statuscd)] aspen_fia[,baac:=ba*tpa_unadj,by=list(plot_id,kindcd,subp,tree,statuscd)] aspen_fia[,sdiprime_met:=(((dia_cm)/25)^1.605),by=list(plot_id,kindcd,subp,tree,statuscd)] aspen_fia[,sdiprime_us:=(dia/10)^1.605,by=list(plot_id,kindcd,subp,tree,statuscd)] aspen_fia[,sdiprime_ac:=sdiprime_us*tpa_unadj,by=list(plot_id,kindcd,subp,tree,statuscd)] aspen_fia[,sdiprimeha:=sdiprime_met*(tpa_unadj*2.47),by=list(plot_id,kindcd,subp,tree,statuscd)] aspen_fia[,baha:=ba_met*(tpa_unadj*2.47),by=list(plot_id,kindcd,subp,tree,statuscd)] aspen_fia[,tph:=(tpa_unadj*2.47),by=list(plot_id,kindcd,subp,tree,statuscd)] aspen_fia[,sdiha_by_subplot:=sum(sdiprimeha),by=list(plot_id,kindcd,subp)] ##generate point (semi-distance independent) competition indices by subplot## ##competition indices are semi-distance independent, based on tree relative sizes #within a given FIA subplot (as opposed to stand-level aggregate data) aspen_fia_competition_indices=melt(aspen_fia,id.vars=c('plot_id','subp','tree',"kindcd","statuscd"),measure.vars=c("ba_met","dia","dia_cm","tpa_unadj","ht")) aspen_fia_competition_indices=dcast(aspen_fia_competition_indices,kindcd+plot_id+subp+tree+statuscd~variable,mean) aspen_fia_competition_indices_live=subset(aspen_fia_competition_indices,statuscd=="1") aspen_fia_competition_indices=merge(aspen_fia_competition_indices,aspen_fia_competition_indices_live,by=c('plot_id','subp','kindcd'),allow.cartesian=TRUE,all.x=TRUE) aspen_fia_competition_indices=subset(aspen_fia_competition_indices,tree.x!=tree.y) aspen_fia_competition_indices=as.data.table(aspen_fia_competition_indices) aspen_fia_competition_indices[,mean_ba:=mean(ba_met.y),by=list(kindcd,plot_id,subp)] aspen_fia_competition_indices$tree=aspen_fia_competition_indices$tree.x #glover_hool is CI-1 in the manuscript #lorimer is CI-2 in the manuscript, using relative tree diameters #ci_assym is CI-3 in the manuscript #variants of these three CIs with an "ht" prefix use relative tree heights instead #of tree diameters aspen_fia_competition_indices=aspen_fia_competition_indices[,list(glover_hool=sum(ba_met.x/mean_ba),ci_assym=sum(((dia_cm.y/dia_cm.x)^2)*tpa_unadj.y),lorimer=sum((dia.y/dia.x)*tpa_unadj.y),ht_lorimer=sum((ht.y/ht.x)*tpa_unadj.y),ht_ci_assym=sum(((ht.y/ht.x)^2)/tpa_unadj.y)),by=list(kindcd,plot_id,subp,tree)] #rejoin computed competition indices to individual tree records aspen_fia=merge(aspen_fia,aspen_fia_competition_indices,by=c('kindcd','plot_id','subp','tree')) ##create plot and plot x species level summaries for generating data matrices and stand structural statistics## aspen_fia_live=subset(aspen_fia,statuscd=="1") density_by_plot_species=aspen_fia_live[order(plot_id,measyear,common_name,aspect,slope,lat,lon,elev),list(sdipct_rmrs=mean(sdipct_rmrs),sdiprime_ha=sum(na.omit(sdiprimeha)),ht=mean(na.omit(ht)),sdiprimeac=sum(na.omit(sdiprime_ac)),baac=sum(na.omit(baac)),baha=sum(na.omit(baha)),tph=sum(na.omit(tpa_unadj)),tpa=sum(tpa_unadj)),by=list(plot_id,kindcd,common_name)] density_by_plot=aspen_fia_live[order(plot_id),list(sdiprime_ha=sum(na.omit(sdiprimeha)),physclcd=mean(physclcd),sdipct=mean(sdipct_rmrs),sdimax=mean(sdimax_rmrs),sdiprimeac=sum(na.omit(sdiprime_ac)),baac=sum(na.omit(baac)),baha=sum(na.omit(baha)),tph=sum(tpa_unadj*2.47),tpa=sum(tpa_unadj)),by=list(plot_id,kindcd,fortypcd,statecd)] density_by_plot=melt(density_by_plot,id.vars=c("plot_id",'kindcd'),measure.vars=c('sdiprime_ha','baac','baha','growcfal_ac',"ht_over_age")) density_by_plot_wide=dcast(density_by_plot,plot_id~variable+kindcd,mean) #plot basal area density_by_plot_species[,plot_baha:=sum(na.omit(baha)),by=list(plot_id,measyear,statuscd)] #proportion of species by basal area density_by_plot_species[,prop_baha:=baha/plot_baha,by=list(plot_id,statuscd,common_name,measyear)] #proportion of species by summation method SDI #plot biomass. not used. density_by_plot_species[order(plot_id,measyear),plot_biomass:=sum(na.omit(carbon_ag_ac)),by=list(plot_id,statuscd,measyear)] #plot summation methond stand density index density_by_plot_species[order(plot_id,measyear),plot_sdiprimeac:=sum(na.omit(sdiprimeac)),by=list(plot_id,statuscd,measyear)] #proportion of species by plot biomass density_by_plot_species[order(plot_id,measyear),prop_biomass:=carbon_ag_ac/plot_biomass,by=list(plot_id,statuscd,measyear,common_name)] #proportion by stand density index density_by_plot_species[order(plot_id,measyear),prop_sdiprimeac:=sdiprimeac/plot_sdiprimeac,by=list(plot_id,statuscd,measyear,common_name)] density_by_plot_species=merge(density_by_plot_species,REF_SPECIES,by="common_name") melt_density_by_plot_species=melt(density_by_plot_species,id.vars=c("PLOT","symbol","KINDCD","STATUSCD","MEASYEAR"),measure.vars=c("prop_sdiprimeac")) ##create matrices of species proportional abundance by plot and plot measurement date. sdiproportion_matrix=as.data.table(dcast(melt_density_by_plot_species,plot_id+kindcd+measyear~symbol,mean)) sdiproportion_matrix[is.na(sdiproportion_matrix)]=0 ##drop plots with negligeable or 100% aspen relative abundance sdiproportion_matrix_qa=subset(sdiproportion_matrix,POTR5>0 & POTR5<1) ##create a remeasurement indicator variable to only retain plots with both initial and remeasurement data sdiproportion_matrix_qa[,remeasured:=sum(kindcd),by="plot_id"] sdiproportion_matrix_qa=subset(sdiproportion_matrix_qa,remeasured=="3") rownames(sdiproportion_matrix_qa)=make.names(paste(sdiproportion_matrix_qa$plot_id,sdiproportion_matrix_qa$kindcd)) sdiproportion_matrix_qa$measyear=NULL sdiproportion_matrix_qa$remeasured=NULL ##additional proportion aspen by plot table for use in example analysis sdiproportion_matrix_qa_avg=as.data.table(sdiproportion_matrix_qa) sdiproportion_qa_avg=sdiproportion_matrix_qa_avg[,list(POTR5=mean(POTR5)),by="plot_id"] ##convert long format data to wide format aspen_fia_analysis_melt=melt(aspen_fia,id.vars=c("plot_id","tree","subp","lat","lon","kindcd","common_name","statecd","countycd"),measure.vars=c("dia","prevdia","cclcd", "statuscd","remper", "dstrbcd1","measyear","invyr","growcfal","cr","ht","tpa_unadj","volcfnet","glover_hool","sdiha_by_subplot","lorimer","ci_assym","ht_lorimer","ht_ci_assym")) aspen_fia_cast=as.data.table(dcast(aspen_fia_analysis_melt,plot_id+lat+lon+subp+tree+common_name+statecd+countycd~variable+kindcd,mean)) ##calculation of average (initial plus plot remeasurement) values of compeition indices and subplot-level SDI aspen_fia_cast[,glover_hool_avg:=colSums(as.data.frame(glover_hool_1, glover_hool_2))/2,by=list(plot_id,subp,tree)] #aspen_fia_cast[,wykoff_avg:=sum(na.omit(wykoff_1),na.omit(wykoff_2))/2,by=list(plot_id,subp,tree)] aspen_fia_cast[,lorimer_avg:=colSums(as.data.frame(lorimer_1,lorimer_2))/2,by=list(plot_id,subp,tree)] aspen_fia_cast[,ht_lorimer_avg:=colSums(as.data.frame(ht_lorimer_1,ht_lorimer_2))/2,by=list(plot_id,subp,tree)] aspen_fia_cast[,ht_ci_assym_avg:=colSums(ht_ci_assym_1,ht_ci_assym_2)/2,by=list(plot_id,subp,tree)] aspen_fia_cast[,sdiha_by_plot_avg:=colSums(as.data.frame(sdiha_by_subplot_1,sdiha_by_subplot_2))/2,by=list(plot_id,subp)] aspen_fia_cast[,ci_assym_avg:=colSums(as.data.frame(ci_assym_1,ci_assym_2))/2,by=list(plot_id,subp,tree)] #aspen_fia_cast[,ht_wykoff_avg:=sum(na.omit(ht_wykoff_1),na.omit(ht_wykoff_2))/2,,by=list(plot_id,subp,tree)] #ielev_fia_cast[,ht_wykoff_avg:=sum(na.omit(ht_wykoff_1),na.omit(ht_wykoff_2))/2,,by=list(plot_id,subp,tree)] aspen_fia_cast[,scale_log_lorimer_avg:=scale(log(lorimer_avg))] aspen_fia_cast[,scale_log_ci_assym_avg:=scale(log(ci_assym_avg))] aspen_fia_cast[,scale_log_glover_hool_avg_avg:=scale(log(glover_hool_avg))] aspen_fia_cast[,scale_cr:=scale(cr_1)] aspen_fia_cast[,scale_log_dia:=scale(log(dia_1))] aspen_fia_cast[,scale_log_ht:=scale(log(ht_1))] aspen_fia_cast[,scale_log_ht_lorimer:=scale(log(ht_lorimer_avg))] aspen_fia_cast[,scale_log_ht_ci_assym:=scale(log(ht_ci_assym_avg))] ##generate summaries of plot physical characteristics, stand height, and site index based on plot remeasurement data library(pracma) aspen_fia_remeasure=as.data.table(subset(aspen_fia,kindcd=="2")) plot_summary_statistics=aspen_fia_remeasure[order(plot_id),list(siteclcd=mean(siteclcd),physclcd=mean(physclcd),ht=mean(na.omit(ht)),stdage=mean(stdage),sitebase=mean(sibase),lat=mean(lat),lon=mean(lon),elev=mean(elev),aspect=mean(aspect),slope=mean(slope)),by=list(plot_id)] plot_summary_stats_plus_density=merge(density_by_plot_wide, plot_summary_statistics,by=c("plot_id")) #Beers et al. (1966) heat load (not used) plot_summary_stats_plus_density[,heat_load:=abs(180-abs(aspect-225))] #convert slope and aspect to radians plot_summary_stats_plus_density[,":="(lat_rad=deg2rad(lat),slope_rad=deg2rad(slope),aspect_rad=deg2rad(aspect)),by=plot_id] #heat load from Mccune and Keon (2002), incorporating aspect,slope, and latitude plot_summary_stats_plus_density[,mccune_heat_load:=(0.339+0.808*(cos(lat_rad)*cos(slope_rad))-0.196*(sin(lat_rad)*sin(slope_rad))-0.482*(cos(aspect_rad)*sin(slope_rad))),by=plot_id] plot_summary_stats_plus_density[,physclcd:=as.factor(physclcd)] plot_summary_stats_plus_density[,scale_stdage:=scale(stdage)] plot_summary_stats_plus_density[,elev_m:=elev/3.28,] plot_summary_stats_plus_density[,scale_siteclcd:=scale(siteclcd)] ##generate second-to-final files (only ommitting climate data) for use in linear mixed modeling## analysis_aspen_trees=merge(sdiproportion_matrix_qa_avg,aspen_fia_cast,by=c("plot_id")) analysis_aspen_trees=merge(plot_summary_stats_plus_density,analysis_aspen_trees,by="plot_id") analysis_aspen_trees=as.data.table(analysis_aspen_trees) analysis_aspen_trees[,scale_siteclcd:=scale(siteclcd)] analysis_aspen_trees[,min_lat:=min(lat)] analysis_aspen_trees[,latitude_correction:=round(lat-min_lat,0)] analysis_aspen_trees[,cclcd_1:=droplevels(as.factor(cclcd_1))] analysis_aspen_trees[,elev_m_lat_correction:=elev_m+129.4*(round(lat-min_lat,0))] analysis_aspen_trees[,scale_elev_m_lat_correction:=scale(elev_m_lat_correction)] analysis_aspen_trees[,statuscode_2:=statuscd_2-1,] analysis_aspen_trees[,log_growcfal_by_plot:=log(sum(growcfal_2,na.rm=TRUE)),by=list(plot_id,common_name)] analysis_aspen_trees[,plot_id:=as.factor(plot_id)] ##contruct plot-level climate variable averages of VPD and PPT #VPD and PPT were constucted based on import multiband raster of PRISM VPD and PPT files generated in GIS #These data are in monthly format library(raster) library(tidyr) library(stringr) library(prism) library(raster) library(magrittr) ##extract multiband raster from raster stack generated via QGIS gdal "merge" command to fia plots setwd("D:/Dropbox/ASCC/ASCC_GIS/Terraclim") merged_prism_ppt_and_vpd=stack("merged_prism_ppt_and_vpd.tif") mycrs <- merged_prism_ppt_and_vpd@crs@projargs aspen_plot_coordinates=analysis_aspen_trees[,list(lat=mean(lat),lon=mean(lon)),by=plot_id] coordinates(aspen_plot_coordinates) <- c('lon', 'lat') proj4string(aspen_plot_coordinates) <- CRS(mycrs) aspen_points_with_ppt_and_vpd <- data.table(coordinates(aspen_plot_coordinates), aspen_plot_coordinates$plot_id, raster::extract(merged_prism_ppt_and_vpd, aspen_plot_coordinates,method="bilinear")) aspen_points_with_ppt_and_vpd_melt=melt(aspen_points_with_ppt_and_vpd,id.vars=c("aspen_plot_coordinates.plot_id","lat","lon")) ##create table intepreting climate variables based on raster band order climate_data_variables=data.table(levels(as.factor(aspen_points_with_ppt_and_vpd_melt$variable))) climate_data_variables$variable=climate_data_variables$V1 climate_data_variables[,variable_sequence:=as.integer(substr(variable,26,28))] climate_data_variables[,climate_variable:=ifelse(variable_sequence<229,"ppt",ifelse(variable_sequence<457 & variable_sequence>228,"vpd_max","vpd_min"))] climate_data_variables[,year:=rep(2000:2018,times=3,each=12)] climate_data_variables[,month:=rep(1:12,times=57,each=1)] climate_data_variables[,year_shift:=ifelse(month>7 & month<11 | month>10,year+1,year)] climate_data_variables[,season:=ifelse(month>5 & month <11,"warm_season",ifelse(month>10 |month<4, "winter",NA))] ##rejoin to points climate matrix aspen_points_with_ppt_and_vpd_melt=merge(climate_data_variables,aspen_points_with_ppt_and_vpd_melt,by="variable") aspen_points_with_ppt_and_vpd_melt=as.data.table(aspen_points_with_ppt_and_vpd_melt) aspen_points_with_ppt_and_vpd_melt[,plot_id:=aspen_plot_coordinates.plot_id] aspen_points_with_ppt_and_vpd_cast=dcast(aspen_points_with_ppt_and_vpd_melt,plot_id+year_shift+month+season~climate_variable,mean,na.rm=TRUE) aspen_points_with_ppt_and_vpd_cast[,vpd_avg:=vpd_max-vpd_min,] ##arrange aspen climate data by plot and measurement period aspen_plot_measyear_and_remper=analysis_aspen_trees[,list(measyear=ifelse(is.na(measyear_2),mean(invyr_2),mean(measyear_2,na.rm=TRUE)),remper=mean(remper_2,na.rm=TRUE)),by=plot_id] aspen_plot_measyear_and_remper[,plot_id:=as.factor(plot_id)] aspen_points_with_ppt_and_vpd_cast=merge(aspen_plot_measyear_and_remper,aspen_points_with_ppt_and_vpd_cast,by="plot_id",allow.cartesian=TRUE) aspen_points_with_ppt_and_vpd_cast[,timediff:=as.numeric(measyear-2*round(remper,0))] aspen_points_with_ppt_and_vpd_cast[,measyear_plus_1:=measyear+1] aspen_points_with_ppt_and_vpd_cast=subset(aspen_points_with_ppt_and_vpd_cast, year_shift>timediff & year_shift7 & month<11 | month>10,year+1,year)] climate_data_cwd[,season:=ifelse(month>5 & month <11,"warm_season",ifelse(month>10 |month<4, "winter",NA))] aspen_points_with_cwd=merge(climate_data_cwd,aspen_points_with_cwd_melt, by="variable") aspen_plot_measyear_and_remper=unique(analysis_aspen_trees[,list(measyear_2=ifelse(is.na(measyear_2),mean(invyr_2,na.rm=TRUE),mean(invyr_2,na.rm=TRUE)),remper=mean(remper_2,na.rm=TRUE)),by="plot_id"]) aspen_plot_measyear_and_remper[,plot_id:=as.factor(plot_id)] aspen_points_with_cwd=as.data.table(merge(aspen_plot_measyear_and_remper,aspen_points_with_cwd,by="plot_id",allow.cartesian=TRUE)) aspen_points_with_cwd[,timediff:=as.numeric(measyear_2-2*round(remper,0))] aspen_points_with_cwd[,timediff:=as.numeric(measyear_2-2*round(remper,0))] aspen_points_with_cwd[,measyear_plus_1:=measyear_2+1] aspen_points_with_cwd=subset(aspen_points_with_cwd, year_shift>timediff & year_shift4.9,1,0)] analysis_aspen_trees[,live_tally_1:=ifelse(statuscd_1=="1" & dia_1>4.9,1,0)] analysis_aspen_trees[,live_tally_2:=ifelse(statuscode_2=="0" & dia_1>4.9,1,0)] analysis_aspen_trees[,scale_remper:=scale(remper_2)] analysis_aspen_trees=merge(analysis_aspen_trees,aspen_points_with_ppt_and_mean_temp_cast2,by="plot_id") analysis_aspen_trees[,scale_mccune_heat_load:=scale(mccune_heat_load)] analysis_aspen_trees[,PDSI:=(PDSI_avg)] analysis_aspen_trees[,scale_PDSI:=scale(PDSI)] ##create final data file for plot-level recruitment analysis_aspen_recruitment=analysis_aspen_trees[,list(ingrowth=sum(ingrowth,na.rm=TRUE),live_tally_1=sum(live_tally_1,na.rm=TRUE),live_tally_2=sum(live_tally_2,na.rm=TRUE),sdi_by_species_1=sum((dia_1/25)^1.605,na.rm=TRUE),sdi_by_species_2=sum((dia_2/25)^1.605,na.rm=TRUE)),by=list(plot_id,scale_mccune_heat_load,remper_2,scale_elev_m_lat_correction,scale_POTR,scale_log_ppt_cold,scale_log_ppt_total,scale_log_ppt_warm,scale_vpd_total,scale_vpd_warm,common_name_sub_fir,POTR5,cwd_avg)] analysis_aspen_recruitment[,sdi_by_species_mean:=(sdi_by_species_1+sdi_by_species_2)/2] analysis_aspen_recruitment[,change:=(live_tally_1-live_tally_2)] analysis_aspen_recruitment=merge(analysis_aspen_recruitment,plot_summary_stats_plus_density,by=c("plot_id")) analysis_aspen_recruitment[,scale_log_sdiprime_ha:=scale((sdiprime_ha_1+sdiprime_ha_2)/2)] analysis_aspen_recruitment[,scale_log_live_tally_1:=scale(log(live_tally_1+min(live_tally_1+0.001))),common_name_sub_fir] analysis_aspen_recruitment[,ingrowth_present:=ifelse(ingrowth>0,1,0),by=list(plot_id,common_name_sub_fir)] analysis_aspen_recruitment[,scale_remper:=scale(remper_2)] analysis_aspen_recruitment[,null:=1]