edit figure files and added percentages

Add_percent.R

@@ -0,0 +1,94 @@
+
+
+
+# what is the ratio between percent and irradiance
+
+dfi2 <- read.csv("F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/march21_all/dfa_2.csv")
+dfi8 <- read.csv("F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/march21_all/dfa_8.csv")
+dfp2 <- read.csv("F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/perc_march/pm_dfa_2.csv")
+dfp8 <- read.csv("F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/perc_march/pm_dfa_8.csv")
+
+plot(dfi2$Combined~dfp2$Combined)
+
+
+### PERCENTAGES FOR FIG6 Anomalies
+
+par(mfrow=c(3,5))
+par(oma = c(2,3,0.8,2))
+par(mar=c(1.5,1,0,1.5))
+
+#plot 5 times
+y_lim = c(-30,90)
+for (i in 1:5){
+  plot(smooth.spline(dfi2$Combined~dfi2$Elevation1), type = 'l', col='slategrey',
+       ylim=y_lim, ylab=expression(Mean~change~'in'~irradiance~(Wm^2)), lwd = 5,
+       lty = 1,xlab='Elevation (m)',cex.lab=1.3, cex.axis=1.8, xaxt='n',yaxt='n')
+  axis(2,cex.axis=1.8, labels=NA)
+  axis(1,cex.axis=1.8, labels=NA)
+  abline(h=0,lty=3)
+  # abline(h=50,lty=3)
+  par(new = T)
+  # y_lim = c(-10,26.6)
+  y_lim2 = y_lim/3
+  plot(smooth.spline(dfp2$Combined~dfp2$Elevation1), type = 'l', col='red',
+       ylim=y_lim2, ylab=NA, lwd = 2,
+       lty = 1,xlab=NA,cex.lab=1.3, cex.axis=1.8, axes=F)
+  axis(4, cex.axis=1.8,cex.lab=1.3)
+}
+# south (plot 5 times)
+for (i in 1:5){
+  plot(smooth.spline(dfi8$Combined~dfi8$Elevation1), type = 'l', col='slategrey',
+       ylim=y_lim, ylab=expression(Mean~change~'in'~irradiance~(Wm^2)), lwd = 5,
+       lty = 1,xlab='Elevation (m)',cex.lab=1.3, cex.axis=1.8, xaxt='n',yaxt='n')
+  axis(2,cex.axis=1.8, labels=NA)
+  axis(1,cex.axis=1.8, labels=NA)
+  abline(h=0,lty=3)
+  # abline(h=50,lty=3)
+  par(new = T)
+  # y_lim = c(-10,26.6)
+  y_lim2 = y_lim/3
+  plot(smooth.spline(dfp8$Combined~dfp8$Elevation1), type = 'l', col='red',
+       ylim=y_lim2, ylab=NA, lwd = 2,
+       lty = 1,xlab=NA,cex.lab=1.3, cex.axis=1.8, axes=F)
+  axis(4, cex.axis=1.8,cex.lab=1.3)
+}
+
+
+
+### FIG 5 PERCENTAGE
+
+dfi2 <- read.csv("F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/march21_all/dfa_2.csv")
+dfp2 <- read.csv("F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/perc_march/pm_dfa_2.csv")
+
+par(mfrow=c(1,2))
+y_lim = c(-150,110) #y_lim = c(-30,100)
+plot(smooth.spline(dfi2$Combined~dfi2$Elevation1), type = 'l', col='slategrey',
+     ylim=y_lim, ylab=expression(Mean~change~'in'~irradiance~(Wm^2)), lwd = 5,
+     lty = 1,xlab='Elevation (m)',cex.lab=1.3, cex.axis=1.8)
+abline(h=0,lty=3)
+# abline(h=50,lty=3)
+par(new = T)
+# y_lim = c(-10,26.6)
+y_lim2 = y_lim/3
+plot(smooth.spline(dfp2$Combined~dfp2$Elevation1), type = 'l', col='red',
+     ylim=y_lim2, ylab=NA, lwd = 2,
+     lty = 1,xlab=NA,cex.lab=1.3, cex.axis=1.8, axes=F)
+axis(4, cex.axis=1.8,cex.lab=1.3)
+
+
+
+### FIG 7 PERCENTAGE (^change to -150,110)
+y_lim = c(-150,110) #y_lim = c(-30,100)
+plot(smooth.spline(dfi2$Combined~dfi2$Elevation1), type = 'l', col='slategrey',
+     ylim=y_lim, ylab=expression(Mean~change~'in'~irradiance~(Wm^2)), lwd = 5,
+     lty = 1,xlab='Elevation (m)',cex.lab=1.3, cex.axis=1.8)
+abline(h=0,lty=3)
+# abline(h=50,lty=3)
+par(new = T)
+# y_lim = c(-10,26.6)
+y_lim2 = y_lim/3
+plot(smooth.spline(dfp2$Combined~dfp2$Elevation1), type = 'l', col='red',
+     ylim=y_lim2, ylab=NA, lwd = 2,
+     lty = 1,xlab=NA,cex.lab=1.3, cex.axis=1.8, axes=F)
+axis(4, cex.axis=1.8,cex.lab=1.3, at=c(-40,-30,-20,-10,0,10,20,30))
+axis(4, cex.axis=1.8,cex.lab=1.3, at=c(-40,-20))

Fig4.R

@@ -177,7 +177,7 @@ p.sw.dem.anom <- function(df, mod = 7, y_lim = c(-100,100), dem_ls = c("ASTER","
 }
 
 
-
+# Fig 7
 # north glacier anomaly
 dfadem2 <- read.csv("F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/dem_comp/dfa_2.csv")
 p.sw.dem.anom(dfadem2, mod = 7, y_lim = c(-150,110), dem_ls = c("ASTER","SRTM", "ALOS", "WV30"), placeLegend = "bottom")

Fig5.R

@@ -109,8 +109,8 @@ par(oma = c(2,3,0.8,0.5))
 par(mar=c(1.5,1,0,0))
 # anomalies for each season
 # march
-p.anom.fig5(dfa2m,y_lim = c(-30,90), drop_res = c(1,4,7), ax_bottom = FALSE)
-p.anom.fig5(dfa8m,y_lim = c(-30,90), drop_res = c(1,4,7), placeLegend = 0)
+p.anom.fig5(dfa2m,y_lim1 = c(-30,90), drop_res = c(1,4,7), ax_bottom = FALSE)
+p.anom.fig5(dfa8m,y_lim2 = c(-30,90), drop_res = c(1,4,7), placeLegend = 0)
 # june
 p.anom.fig5(dfa2j,y_lim = c(-30,90), drop_res = c(1,4,7), ax_bottom = FALSE)
 p.anom.fig5(dfa8j,y_lim = c(-30,90), drop_res = c(1,4,7), placeLegend = 0)
@@ -123,7 +123,7 @@ p.anom.fig5(dfa2d,y_lim = c(-140,200), drop_res = c(1,4,7), placeLegend = 0)
 # change line types!
 
 
-#
+# FIG 6 (now 5) - combined glaciers
 #### PLOT ANOMALY FOR ALL GLACIERS AT ONE RESOLUTION
 dfga_j <- read.csv("F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/june21_all/dfga.csv")
 dfga_m <- read.csv("F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/march21_all/dfga.csv")

toposol2_demcomp.R

@@ -24,10 +24,10 @@ dem_ls = list(wv_dem, aster_dem, srtm_dem, alos_dem, wv_dem)
 # df_alos <- sw.res(alos_dem, shape = glaciers[2,])
 
 #################################################
-sw.dem.comp <- function(dem_ls, shape, gn = NA, date = ISOdate(2017, 3, 21, 0) , savepath = "F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/dem_comp/", isave = "F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/dem_comp/images_tmp/"){
+sw.dem.comp <- function(dem_ls, shape, gn = NA, date = ISOdate(2017, 3, 21, 0), perc = FALSE, savepath = "F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/dem_comp/", isave = "F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/dem_comp/images_tmp/"){
   # creates a dataframe of daily values at each resolution for a given glacier
   # creates anomaly dataframe
-  dem_names = c("WV", "ASTER","SRTM", "ALOS", "WV30")
+  dem_names = c("HiMAT8m", "ASTER","SRTM", "ALOS", "HiMAT30m")
   strt = Sys.time()
 
   # loop through resolutions for single glacier
@@ -40,7 +40,7 @@ sw.dem.comp <- function(dem_ls, shape, gn = NA, date = ISOdate(2017, 3, 21, 0) ,
       location.variables(dem_ls[[dd]], shape, resampleFactor = 1)
     }
 
-    tfstk <- sw.daily(date)
+    tfstk <- sw.daily(date, perc = perc)
     # save output of each resolution to one dataframe
     dft <- sw.res.dataframe(tfstk)
     dft$dem <- dem_names[dd]
@@ -82,15 +82,15 @@ sw.dem.comp <- function(dem_ls, shape, gn = NA, date = ISOdate(2017, 3, 21, 0) ,
 }
 
 ###############################################
-sw.glacier.dems <- function(date = ISOdate(2017, 6, 21, 0), savepath = "F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/dem_comp/"){
+sw.glacier.dems <- function(date = ISOdate(2017, 6, 21, 0), perc=FALSE, savepath = "F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/dem_comp/"){
   # iterates through all glaciers indices given and calls sw.res
   # returns anomaly dataframe for all glaciers
   strt = Sys.time()
   # iterate through all glaciers
   g_num     <- length(glaciers@polygons)
   for (g in 1:g_num){
     # run for single glacier
-    dfg <- sw.dem.comp(dem_ls = dem_ls, shape = glaciers[g,], gn = g, date = date) # anomaly for g
+    dfg <- sw.dem.comp(dem_ls = dem_ls, shape = glaciers[g,], gn = g, date = date, perc = perc) # anomaly for g
 
     # compile dataframe
     if (g == 1){dfga = dfg} else{dfga = rbind(dfga,dfg)}
@@ -106,5 +106,6 @@ sw.glacier.dems <- function(date = ISOdate(2017, 6, 21, 0), savepath = "F:/HiMAT
 }
 
 # run code for all DEMs and glaciers
-dfga <- sw.glacier.dems(date = ISOdate(2017, 3, 21, 0))
+svp = svp = "F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/dem_comp_perc/demp_"
+dfga <- sw.glacier.dems(date = ISOdate(2017, 3, 21, 0), perc=TRUE, savepath = svp)
 

toposol2_main.R

@@ -27,8 +27,8 @@ load.wv.dat(dpath, gpath)
 ## RUN RUN RUN
 #### RUN RUN RUN RUN
 # line 182 of toposol2_functions (Percentage = TRUE/FALSE)
-svp = "F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/NEW_variables/march_irr/mi_"
-dfga <- sw.glacier.res(date = ISOdate(2017, 3, 21, 0), savepath = svp, perc=FALSE)
+svp = "F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/NEW_variables/march_perc/mp_"
+dfga <- sw.glacier.res(date = ISOdate(2017, 3, 21, 0), savepath = svp, perc=TRUE)
 #### RUN RUN RUN RUN
 ## RUN RUN RUN
 
@@ -56,7 +56,7 @@ p.sw.anom(dfa2,y_lim = c(-10,30), plot4 = FALSE, topovar = 7)
 # try to calculate percentage for 90m north glacier
 dftmp <- sw.res(demL, shape = glaciers[2,], gn= 2, resampleFactor = 4, savepath = NULL)
 dftmp2 = read.csv("F:/HiMAT/MATTO/PROJECTS/WV_RESOLUTION/variables/perc_march/pm_dftf_2.csv")
-p.sw.elv(dftmp2,res=8, y_lim = c(-90,30))
+p.sw.elv(dftmp2,res=8, y_lim = c(-110,30))
 
 
 p.sw.anom(dfa,y_lim = c(-50,100))