refined visuals and change table of temporal adjustment effect

main.R

@@ -23,7 +23,7 @@ outDir <- "D:/BiomassCCI_2019/results"
 dataDir <- "D:/BiomassCCI_2019/data"
 plotsFile <- 'SamplePlots.csv'
 plotsFile1 <- 'SamplePoly.csv'
-agbTilesDir <- "E:/Glob2010/" #*
+agbTilesDir <- "E:/GlobBiomass2017/" #*
 treeCoverDir <- 'E:/treecover2010_v3' #*
 flDir <- 'E:/GFCFolder' 
 forestTHs <- 10 
@@ -73,11 +73,11 @@ plots2 <- BiomePair(plots1)
 # apply growth data to whole plot data by identifying AGB map year
 gez <- sort(as.vector((unique(plots2$GEZ)))) #get unique gez and without NA (sorting removes it also)
 plots.tf <- ldply(lapply (1:length(gez), function(x) 
-  TempApply(plots2, gez[[x]], 2003)), data.frame) #change the year!
+  TempApply(plots2, gez[[x]], 2000)), data.frame) #change the year!
 
 #tree growth data uncertainty estimate
 plots.var <- ldply(lapply (1:length(gez), function(x) 
-  TempVar(plots2, gez[[x]], 2003)), data.frame) 
+  TempVar(plots2, gez[[x]], 2000)), data.frame) 
 
 #get absolute uncertainty of temporally adjusted plots 
 plots.tf$sdGrowth <- abs(plots.tf$AGB_T_HA - plots.var$SD)
@@ -87,8 +87,8 @@ plots3 <- plots2[with(plots2, order(GEZ)), ]
 plots.tf$AGB_T_HA_ORIG <- plots3$AGB_T_HA
 
 #histogram of temporal fix effect
-HistoTemp(plots.tf, 2003)
-HistoShift(plots.tf, 2003)
+HistoTemp(plots.tf, 2000)
+HistoShift(plots.tf, 2000)
 rm(plots1, plots2, plots3, plots.var) 
 
   # export new AGB data according to date generated (optional)
@@ -141,7 +141,7 @@ for(biome in biomes){
 for(continent in continents){
   cat("Processing: ",continent,"\n")
   
-  AGBdata <- invDasymetry("ZONE", continent, wghts = TRUE, is_poly = TRUE)
+  AGBdata <- invDasymetry("ZONE", continent, wghts = TRUE, is_poly = F)
   
   save(AGBdata, file = file.path(outDir,
                                  paste0("InvDasyPlot_", continent, ".Rdata")))

results/TF_pre_post_change_2003.csv

-"agb_Mgha_bins","n_pre","n_post","agb_Mgha_pre","agb_Mgha_post"
-"(0,20]",2,2,4.621,4.621
-"(60,80]",1,2,64.3169,50.60075
-"(80,100]",1,1,97.6846,77.2021277
-"(100,120]",2,1,102.5531915,97.1042553
-"(120,140]",1,1,124.2553191,117.9553191
-"(160,180]",1,1,175.091,169.241
-"(200,300]",2,2,281.6745,273.8745

scripts/InvDasymetry.R

@@ -56,7 +56,7 @@ invDasymetry <- function(clmn = "ZONE", value = "Europe", aggr = NULL,
                      .packages='raster', .export=c('MakeBlockPolygon', 'SRS',
                                                    'sampleTreeCover', 'TCtileNames',
                                                    'AGBtileNames', 'sampleTreeCover',
-                                                   'sampleAGBmap',  'plots',
+                                                   'sampleAGBmap',  'plots.tf',
                                                    'agbTilesDir', 'treeCoverDir',
                                                    'forestTHs')) %dopar% {
                                                      

scripts/TempVis.R

@@ -13,25 +13,25 @@ HistoTemp <- function(df, year){
   # create a bar graph with fixed agb bins
   h1 <- hist(df$AGB_T_HA_ORIG, plot=F, breaks=25)
   h2 <- hist(df$AGB_T_HA, plot=F, breaks=25)
-  png (filename=paste0(outDir,paste0('/histogram_tempfixed_',year,'.png')))
+  png (filename=paste0(outDir,paste0('/histogram_tempfixed_',year,'.png')),  width = 800, height = 600)
   y.ax <- nrow(df) / 4
  # y.ax <-25000
   
-  plot(h1, xaxt="n", col=rgb(0,0,1,1/4), main=main, xlab='AGB_T_HA',ylab='Frequency',
-       xlim = c(0,600), ylim=c(0,y.ax))
-  axis(1,at=0:6*100, labels=c(0:6*100))
+  plot(h1, xaxt="n", col=rgb(0,0,1,1/4), main=main, xlab='AGB(Mg/ha)',ylab='n',
+       xlim = c(0,600), ylim=c(0,y.ax),cex.lab=2, cex.axis=1.5, cex.main=2, cex.sub=2)
+  axis(1,at=0:6*100, labels=c(0:6*100),cex.lab=2, cex.axis=1.5, cex.main=2, cex.sub=2)
   plot(h2,col=rgb(1,0,0,1/4),add=T)
-  legend("topright", c("Before", "After", "Overlap"), 
-         col=c(rgb(0,0,1,1/4), rgb(1,0,0,1/4), rgb(0.5,0,0.5,1/4)), lwd=10)
+  legend("topright", c("Before", "After", "Overlap"),
+         col=c(rgb(0,0,1,1/4), rgb(1,0,0,1/4), rgb(0.5,0,0.5,1/4)), lwd=10, cex=2, bty='n')
   
   dev.off()
   
   
-  plot(h1, xaxt="n", col=rgb(0,0,1,1/4), main=main, xlab='AGB_T_HA',ylab='Frequency',
+  plot(h1, xaxt="n", col=rgb(0,0,1,1/4), main=main, xlab='AGB(Mg/ha)',ylab='n',
        xlim = c(0,600), ylim=c(0,y.ax))
   axis(1,at=0:6*100, labels=c(0:6*100))
   plot(h2,col=rgb(1,0,0,1/4),add=T)
-  legend("topright", c("Before", "After", "Overlap"), 
+  legend("topright", c("Before", "After", "Overlap"),
          col=c(rgb(0,0,1,1/4), rgb(1,0,0,1/4), rgb(0.5,0,0.5,1/4)), lwd=10)
 }
 
@@ -56,11 +56,21 @@ HistoShift <- function(df, year){
   new2 <- as.data.frame(new2)
 
   #calculate change in AGB
-  agg.old <- aggregate(old1["AGB_T_HA_ORIG"], by=old1["group"], mean)
-  agg.new <- aggregate(new1["AGB_T_HA"], by=new1["group"], mean)
+  old3 <- aggregate(old1["AGB_T_HA_ORIG"], by=old1["group"], mean)
+  new3 <- aggregate(new1["AGB_T_HA"], by=new1["group"], mean)
+  
+  if (nrow(old2) != nrow(new2)){
+    fj1 <- full_join(old2, new2, by=c('group'='group'))
+    fj2 <-  full_join(old3, new3, by=c('group'='group'))
+    outs <- cbind(fj1,fj2)
+    outs <- outs[,c(1,2,3,5,6)]
+    
+  }
+  else{
+    outs <- do.call(cbind, list(old2,new2,old3,new3))
+    outs <- outs[,c(1,2,4,6,8)]
+  }
   
-  outs <- do.call(cbind, list(old2,new2,agg.old,agg.new))
-  outs <- outs[,c(1,2,4,6,8)]
   names(outs) <- c('agb_Mgha_bins', 'n_pre', 'n_post', 'agb_Mgha_pre', 'agb_Mgha_post')
   setwd(outDir)
   write.csv(outs, paste0('TF_pre_post_change_',year,'.csv'), row.names = F)