1. Michael Spencer
  2. SSGB_trends

Source

SSGB_trends / SSGB_snowdays_monthly.R

# ----------------------------------------
# ----------------------------------------
# Plots days of snow cover per elevation
# One plot per month (Nov - Apr)
# Annual plot large plot
# ----------------------------------------
# ----------------------------------------

# Set up
library(RSQLite)
library(dplyr)
library(reshape2)
library(raster)
library(rasterVis)
library(lattice)
library(grid)
library(RColorBrewer)
db = dbConnect(SQLite(), "~/Cloud/Michael/Uni_temp/SSGB/SSGB.sqlite")

# ----------------------------------------
# Read SSGB data

# Get list of complete stations per winter
# Looking at Nov-Apr as gives many more stations and years
st = dbGetQuery(db, "
                SELECT HydroYear, Station, COUNT(Station) AS num FROM SSGB_obs
                WHERE strftime('%m', Date) BETWEEN '11' AND '12'
                AND SnowlineElev !='m'
                OR strftime('%m', Date) BETWEEN '01' AND '04'
                AND SnowlineElev !='m'
                GROUP BY HydroYear, Station
                HAVING num > 180
                ")

# Temp table with Nov-Apr data only
dbSendQuery(db, "
            CREATE TEMP TABLE step AS
            SELECT * FROM SSGB_obs
            WHERE strftime('%m', Date) BETWEEN '11' AND '12'
            AND SnowlineElev !='m'
            OR strftime('%m', Date) BETWEEN '01' AND '04'
            AND SnowlineElev !='m'
            ")

df = lapply(1:nrow(st), function(i){
   dbGetQuery(db, paste0("SELECT HydroYear, Station, SnowlineElev, COUNT(SnowlineElev) AS Days FROM step WHERE Station='", st[i, 2], "' AND HydroYear='", st[i, 1], "' GROUP BY HydroYear, Station, SnowlineElev"))
})

df.mon = lapply(c(11, 12, "01", "02", "03", "04"), function(j){
   print(j)
   lapply(1:nrow(st), function(i){
      dbGetQuery(db, paste0("SELECT HydroYear, Station, SnowlineElev, COUNT(SnowlineElev) AS Days FROM step WHERE Station='", st[i, 2], "' AND HydroYear='", st[i, 1], "' AND strftime('%m', Date)='", j, "' GROUP BY HydroYear, Station, SnowlineElev"))
   })
})
   

# Tidy
dbSendQuery(db, "DROP TABLE step")
dbDisconnect(db)
rm(db)


# ----------------------------------------
# Prep data

prep = function(dl){
   # Make cumulative snowline
   x = lapply(dl, function(i){
      # Remove n and 99
      i = i[i$SnowlineElev != "n" & i$SnowlineElev != "99", ]
      # Fix factor order
      i$SnowlineElev = factor(i$SnowlineElev, levels=seq(0, 1200, 150))
      # Order data frame
      i = with(i, i[order(SnowlineElev),])
      # Get days of snow lying from snowline
      #i$Lying[i$SnowlineElev != "n" & i$SnowlineElev != "99"] = NA
      i$Lying = cumsum(i$Days)
      i
   })
   
   x = do.call("rbind.data.frame", x)
   
   # Convert elevation to numeric
   x$SnowlineElev = as.numeric(as.character(x$SnowlineElev))
   
   # Get median per year/elevation
   x = data.frame(summarise(group_by(x, HydroYear, SnowlineElev), Lying=median(Lying)))
   
   # Make time series regular
   y = data.frame(HydroYear=rep(1946:2006, each=9), SnowlineElev=seq(0, 1200, by=150))
   x = merge(x, y, all.y=T)
   
   # Make wide
   dcast(x, HydroYear ~ SnowlineElev)
}

df1 = prep(df)
df.mon1 = lapply(df.mon, prep)

# ----------------------------------------
# Make plot

# Annual
png("~/Cloud/Michael/Uni_temp/SSGB_trends/figures/SSGB_snowdays_annual.png", width=1000)
par(fig=c(0, 1, 0.25, 1), mar=c(5, 5, 4, 2) + 0.1) # c(bottom, left, top, right)
pal = colorRampPalette(c("pink", "slateblue2", "purple4", "midnightblue"))
image(as.matrix(df1[, -1]), col=pal(20), xaxt="n", yaxt="n", main="Winter (Nov - Apr)", xlab="Winter beginning", ylab="Elevation (m)")
axis(1, seq(0, 1, 1/60), 1946:2006)
axis(2, seq(0, 1, 1/8), seq(0, 1200, 150), las=1)
# Legend
# To add a legend see levelplot and bp/sp from NAO paper
par(fig=c(0, 1, 0, 0.25), mar=c(3, 4.5, 0, 2) + 0.1, new=T)
plot.new()
legend("topleft", legend=paste(round(quantile(df1[, -1], c(0, .33, .66, 1), na.rm=T)), "days of snow cover"), pch=15, pt.cex=2, col=c("pink", "slateblue2", "purple4", "midnightblue"))
dev.off()

# Monthly
x = lapply(df.mon1, function(i){
   y = flip(raster(t(as.matrix(i[, -1]))), "y")
   extent(y) = c(0, 10, 0, 1)
   y
})
x = stack(x)
names(x) = month.name[c(11, 12, 1:4)]

png("~/Cloud/Michael/Uni_temp/SSGB_trends/figures/SSGB_snowdays_monthly.png", width=930, height=700)
myTheme=rasterTheme(region=brewer.pal('Blues', n=9))
levelplot(x, par.settings=myTheme, layout=c(1, 6), xlab="Year", ylab="Elevation (m)", scales=list(
   x=list(
      at=seq(0+(10/120), 10-(10/120), by=((10-(10/120))-(0+(10/120)))/6),
      labels=seq(1946, 2006, by=10)),
   y=list(
      at=seq(0+(1/18), 1-(1/18), by=((1-(1/18))-(0+(1/18)))/8),
      labels=seq(0, 1200, by=150))))
trellis.focus("legend", side="right", clipp.off=T, highlight=F)
grid.text("Days of\nsnow cover", 0.25, 0, hjust=0.5, vjust=2.1)
trellis.unfocus()
dev.off()

# ----------------------------------------
# Table for blog

y = lapply(df.mon1, function(i){
   round(colMeans(i[, -1], na.rm=T))
})
y = do.call("rbind.data.frame", y)
names(y) = names(df1[, -1])
y = t(y[, ncol(y):1])

z = apply(y, 1, function(i){
   paste(i, collapse="</td><td>")
})

paste0("<tr><td>", z, "</td></tr>")