Skip to content

Commit

Permalink
add comparison to ERA5 to manual test of climateNA time series.
Browse files Browse the repository at this point in the history
  • Loading branch information
@cmahony
cmahony committed Dec 20, 2024
1 parent 418d452 commit e81a086
Showing 1 changed file with 163 additions and 48 deletions.
211 changes: 163 additions & 48 deletions tests/manual_tests/check_TimeSeriesVsNormals.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -203,7 +203,7 @@ X <- rast(dem) # use the DEM as a template raster
for(e in 1:3){
par(mfrow=c(3,4), mar=c(0,0,0,0))
for(m in 1:12){
X[climatena_anom[, id1]] <- climatena_anom[,get(paste0(elements[e], monthcodes[m]))]
X[climna.anom[, id1]] <- climna.anom[,get(paste0(elements[e], monthcodes[m]))]
if(e==3) X <- log2(X)
X[X>lim] <- lim
X[X < 0-lim] <- 0-lim
Expand Down Expand Up @@ -422,18 +422,18 @@ ahccd.ref <- ahccd.ref[which(ahccd.ref$id %in% unique(ahccd.recent$id))]
ahccd.recent <- ahccd.recent[which(ahccd.recent$id %in% unique(ahccd.ref$id))]

# subtract ref from recent
ahccd.anom <- ahccd.ref
ahccd.anom <- copy(ahccd.ref)
ahccd.anom[, (3:ncol(ahccd.anom)) := lapply(3:ncol(ahccd.recent),
function(i) ahccd.recent[[i]] - ahccd.ref[[i]])]
ppt_columns <- grep("PPT", names(ahccd.anom), value = TRUE)
ahccd.anom[, (ppt_columns) := ahccd.recent[, .SD, .SDcols = ppt_columns] / ahccd.ref[, .SD, .SDcols = ppt_columns]]

## Plot of pairwise comparison of CRU and ClimateNA (with AHCCD stations)
element.names <- c("mean daily minimum temperature (K)", "mean daily maximum temperature (K)", " precipitation (%)")
element.names <- c("mean daily minimum temperature", "mean daily maximum temperature", " precipitation")

X <- rast(dem) # use the DEM as a template raster
e=1
# for(e in 1:3){
for(e in 1:3){

lim.upper <- if(elements[e]=="Pr") 1 else 3
lim.lower <- if(elements[e]=="Pr") -1 else -3
Expand All @@ -446,7 +446,7 @@ e=1
pct <- if(elements[e]=="Pr") 100 else 1

m=2
# for(m in 1:12){
for(m in 1:12){

png(filename=paste("vignettes/CRUvsClimateNA_BCAB", elements[e], monthcodes[m], "png",sep="."), type="cairo", units="in", width=6.5, height=3.3, pointsize=10, res=300)
# pdf(file=paste("results//CMIP6Eval.Fig3", metric,"pdf",sep="."), width=7.5, height=14, pointsize=12)
Expand All @@ -461,8 +461,8 @@ e=1
ahccd.ptData <- ahccd.ptData[hasData]
ahccd.locations <- stn[which(stn$id %in% ahccd.anom[month==m,id][hasData])]
if(e==3) ahccd.ptData <- log2(ahccd.ptData)
ahccd.ptData[ahccd.ptData>lim] <- lim
ahccd.ptData[ahccd.ptData < 0-lim] <- 0-lim
ahccd.ptData[ahccd.ptData>lim.upper] <- lim.upper
ahccd.ptData[ahccd.ptData < lim.lower] <- lim.lower
z_cut <- cut(ahccd.ptData, breaks = breaks, include.lowest = TRUE, labels = colscheme)

for(source in c("cru.gpcc", "climatena")){
Expand All @@ -487,6 +487,8 @@ e=1
lines(bdy.na, lwd=0.4)
points(ahccd.locations$lon, ahccd.locations$lat, bg = as.character(z_cut), pch = 21, cex=1.6, lwd=1.5)

legend("bottomleft", legend="AHCCD stations", bg = as.character(z_cut), pch = 21, pt.cex=1.6, pt.lwd=1.5, bty="n")

print(source)
}

Expand All @@ -495,16 +497,16 @@ e=1
xl <- 0.2; yb <- 0.25; xr <- 0.8; yt <- 0.45
rect(head(seq(xl,xr,(xr-xl)/length(colscheme)),-1), yb, tail(seq(xl,xr,(xr-xl)/length(colscheme)),-1), yt, border=NA, col=colscheme)
rect(xl, yb, xr, yt)
labels <- if(elements[e]=="Pr") paste(round(2^seq(lim.lower,lim.upper,(lim.upper-lim.lower)/2), 2)*pct-pct, "%", sep="") else round(seq(lim.lower,lim.upper,(lim.upper-lim.lower)/2), 2)*pct
labels <- if(elements[e]=="Pr") paste(round(2^seq(lim.lower,lim.upper,(lim.upper-lim.lower)/4), 2)*pct-pct, "%", sep="") else paste0(round(seq(lim.lower,lim.upper,(lim.upper-lim.lower)/4), 2)*pct, "\u00B0", "C")
text(seq(xl,xr,(xr-xl)/(length(labels)-1)),rep(yb,length(labels)),labels,pos=1,cex=1.5,font=1, offset=0.5)
text(mean(c(xl,xr)), yt+0.01, paste("Change in", month.name[m], element.names[e], "\n1961-1990 to 2001-2020"), pos=3, cex=1.5, font=2)

dev.off()

month.abb[m]
# }
}
print(elements[e])
# }
}

#-----------------------------------------
## Compare to AHCCD stations in a small region
Expand Down Expand Up @@ -558,52 +560,165 @@ clim <- downscale(stn.s,

e=1
# for(e in 1:3){
m=2
# for(m in 1:12){
ts.x.climatena <- as.numeric(unique(clim[PERIOD %in% 1901:2020 & DATASET == "climatena", PERIOD]))
ts.y.climatena <- unname(unlist(clim[clim$PERIOD %in% 1901:2020 & DATASET == "climatena",
lapply(.SD, function(x) mean(x, na.rm = TRUE)), by = PERIOD,
.SDcols = paste(c("Tmin", "Tmax", "PPT")[e], monthcodes[m], sep="_")][, 2]))
ts.x.cru <- as.numeric(unique(clim[PERIOD %in% 1901:2020 & DATASET == "cru.gpcc", PERIOD]))
ts.y.cru <- unname(unlist(clim[clim$PERIOD %in% 1901:2020 & DATASET == "cru.gpcc",
lapply(.SD, function(x) mean(x, na.rm = TRUE)), by = PERIOD,
.SDcols = paste(c("Tmin", "Tmax", "PPT")[e], monthcodes[m], sep="_")][, 2]))
ts.ref <- mean(ts.y.climatena[ts.x.climatena%in%1961:1990])

## plot comparing CRU and climateNA time series with AHCCD station data.
png(filename=paste("vignettes/CRUvsClimateNA_ts", elements[e], monthcodes[m], "png",sep="."), type="cairo", units="in", width=6.5, height=4, pointsize=10, res=300)

par(mar = c(2,3,0.5,0.5), mgp = c(1.75, 0.25, 0), mfrow=c(1,1))
ylim <- range(ts.y.climatena, na.rm = T) + c(-1.5,1.5)
plot(ts.x.cru, rep(NA, length(ts.x.cru)), ylim=ylim, type="o", tck=-0.01, col="white", xlab="", ylab=element.names[e])

recent <- vector()
for(station in stations){
data <- ahccd.s[id==station & month==m]
x <- data$year
y <- data[,get(c("Tmin", "Tmax", "PPT")[e])]
if(sum(is.finite(y[x%in%2001:2020]))>16){
y.ref <- mean(y[x%in%1961:1990])
y <- y + (ts.ref-y.ref) # shift to cru baseline
lines(x,y, col=alpha("gray50", .25), lwd=2)
norm <- mean(y[x%in%2001:2020], na.rm=T)
lines(c(2001,2020), rep(norm,2), lty=2, col=alpha("gray50", 1))
recent <- if(station==stations[1]) norm else c(recent, norm)
}
}

lines(ts.x.cru,ts.y.cru, col="dodgerblue", lwd=3)
lines(ts.x.climatena,ts.y.climatena, col="black", lwd=2)
lines(c(1961,1990), rep(mean(ts.y.cru[ts.x.cru%in%1961:1990]),2), lty=2, col=1, lwd=2)
lines(c(2001,2020), rep(mean(ts.y.cru[ts.x.cru%in%2001:2020]),2), lty=2, col="dodgerblue", lwd=2)
lines(c(2001,2020), rep(mean(ts.y.climatena[ts.x.climatena%in%2001:2020]),2), lty=2, lwd=2)
boxplot(recent, add=T, at=2021, width=5, range=0)
legend("topleft", legend=c("ClimateNA", "CRU", "AHCCD stations"), pch=c(NA, NA, NA),
col = c(1,"dodgerblue", "gray80"), lty=c(1,1,1), lwd=c(3,2,2), bty="n")

dev.off()
month.abb[m]
# }
print(elements[e])
# }




#-----------------------------------------
## Compare ERA5, ClimateNA, and AHCCD station data over BC and Alberta
#-----------------------------------------

#calculate the ahccd anomalies
stn <- fread("C:/Users/CMAHONY/OneDrive - Government of BC/Data/AHCCD/AHCCD_location.csv")
ahccd <- fread("C:/Users/CMAHONY/OneDrive - Government of BC/Data/AHCCD/AHCCD.csv")
ahccd <- ahccd[,- "tmean"]
names(ahccd) <- c("id", "year", "month", "Tmin", "Tmax", "PPT")

ahccd.ref <- ahccd[year %in% 1981:2010,
lapply(.SD, function(x) mean(x, na.rm = TRUE)),
by = .(id, month),
.SDcols = c("Tmin", "Tmax", "PPT")]
ahccd.recent <- ahccd[year %in% 2001:2010,
lapply(.SD, function(x) mean(x, na.rm = TRUE)),
by = .(id, month),
.SDcols = c("Tmin", "Tmax", "PPT")]
ahccd.ref <- ahccd.ref[which(ahccd.ref$id %in% unique(ahccd.recent$id))]
ahccd.recent <- ahccd.recent[which(ahccd.recent$id %in% unique(ahccd.ref$id))]

# subtract ref from recent
ahccd.anom <- copy(ahccd.ref)
ahccd.anom[, (3:ncol(ahccd.anom)) := lapply(3:ncol(ahccd.recent),
function(i) ahccd.recent[[i]] - ahccd.ref[[i]])]
ppt_columns <- grep("PPT", names(ahccd.anom), value = TRUE)
ahccd.anom[, (ppt_columns) := ahccd.recent[, .SD, .SDcols = ppt_columns] / ahccd.ref[, .SD, .SDcols = ppt_columns]]

## Plot of pairwise comparison of CRU and ClimateNA (with AHCCD stations)
element.names <- c("mean daily minimum temperature", "mean daily maximum temperature", " precipitation")

X <- rast(dem) # use the DEM as a template raster
e=1
for(e in 1:3){

lim.upper <- if(elements[e]=="Pr") 1 else 2
lim.lower <- if(elements[e]=="Pr") -1 else -2

inc=(lim.upper-lim.lower)/100
breaks=seq(lim.lower, lim.upper+inc, inc)
colscheme <- colorRampPalette(if(elements[e]=="Pr") rev(hcl.colors(5,"Blue-Red 3")) else hcl.colors(5,"Blue-Red 3"))(length(breaks)-1)

anom <- rast(paste0("//objectstore2.nrs.bcgov/ffec/TransferAnomalies/delta.from.1961_1990.to.2001_2020.", elements[e], ".tif"))
pct <- if(elements[e]=="Pr") 100 else 1

m=2
# for(m in 1:12){
ts.x.climatena <- as.numeric(unique(clim[PERIOD %in% 1901:2020 & DATASET == "climatena", PERIOD]))
ts.y.climatena <- unname(unlist(clim[clim$PERIOD %in% 1901:2020 & DATASET == "climatena",
lapply(.SD, function(x) mean(x, na.rm = TRUE)), by = PERIOD,
.SDcols = paste(c("Tmin", "Tmax", "PPT")[e], monthcodes[m], sep="_")][, 2]))
ts.x.cru <- as.numeric(unique(clim[PERIOD %in% 1901:2020 & DATASET == "cru.gpcc", PERIOD]))
ts.y.cru <- unname(unlist(clim[clim$PERIOD %in% 1901:2020 & DATASET == "cru.gpcc",
lapply(.SD, function(x) mean(x, na.rm = TRUE)), by = PERIOD,
.SDcols = paste(c("Tmin", "Tmax", "PPT")[e], monthcodes[m], sep="_")][, 2]))
ts.ref <- mean(ts.y.climatena[ts.x.climatena%in%1961:1990])
for(m in c(1,2,4,7,10)){

png(filename=paste("vignettes/ERA5vsClimateNA_BCAB", elements[e], monthcodes[m], "png",sep="."), type="cairo", units="in", width=6.5, height=3.3, pointsize=10, res=300)

mat <- matrix(c(5,1,2,5,3,4), 3)
layout(mat, widths=c(1,1), heights=c(0.3, .05 ,1))

## plot comparing CRU and climateNA time series with AHCCD station data.
png(filename=paste("vignettes/CRUvsClimateNA_ts", elements[e], monthcodes[m], "png",sep="."), type="cairo", units="in", width=6.5, height=4, pointsize=10, res=300)
par(mar=c(0.1,0.1,0.1,0.1))

par(mar = c(2,3,0.5,0.5), mgp = c(1.75, 0.25, 0), mfrow=c(1,1))
ylim <- range(ts.y.climatena, na.rm = T) + c(-1.5,1.5)
plot(ts.x.cru, rep(NA, length(ts.x.cru)), ylim=ylim, type="o", tck=-0.01, col="white", xlab="", ylab=element.names[e])
ahccd.ptData <- ahccd.anom[month==m, get(c("Tmin", "Tmax", "PPT")[e])]
hasData <- is.finite(ahccd.ptData)
ahccd.ptData <- ahccd.ptData[hasData]
ahccd.locations <- stn[which(stn$id %in% ahccd.anom[month==m,id][hasData])]
if(e==3) ahccd.ptData <- log2(ahccd.ptData)
ahccd.ptData[ahccd.ptData>lim.upper] <- lim.upper
ahccd.ptData[ahccd.ptData < lim.lower] <- lim.lower
z_cut <- cut(ahccd.ptData, breaks = breaks, include.lowest = TRUE, labels = colscheme)

recent <- vector()
for(station in stations){
data <- ahccd.s[id==station & month==m]
x <- data$year
y <- data[,get(c("Tmin", "Tmax", "PPT")[e])]
if(sum(is.finite(y[x%in%2001:2020]))>16){
y.ref <- mean(y[x%in%1961:1990])
y <- y + (ts.ref-y.ref) # shift to cru baseline
lines(x,y, col=alpha("gray50", .25), lwd=2)
norm <- mean(y[x%in%2001:2020], na.rm=T)
lines(c(2001,2020), rep(norm,2), lty=2, col=alpha("gray50", 1))
recent <- if(station==stations[1]) norm else c(recent, norm)
}
for(source in c("era5", "climatena")){

source.name <- if(source=="climatena") "ClimateNA" else "ERA5-Land"
plot(1, type="n", axes=F, xlab="", ylab="")
text(1,1, source.name, font=2,cex=1.35)

par(mar=c(0.1,0.1,0.1,0.1))
if(source=="era5"){
anom <- mean(rast(paste0("C:/Users/CMAHONY/OneDrive - Government of BC/Data/Western North America_",c("tmin", "tmax", "prcp")[e], "_anomaly_",month.abb[m] ,"_2001_2020_data.tif")))
X <- anom
# X <- mask(X, bdy.na)
if(e==3) X <- log2(X/100+1)
} else {
X <- rast(dem) # use the DEM as a template raster
X[climna.anom[, id1]] <- climna.anom[,get(paste0(c("Tmin", "Tmax", "PPT")[e], monthcodes[m]))]
if(e==3) X <- log2(X)
}
X[X>lim.upper] <- lim.upper
X[X < lim.lower] <- lim.lower

image(X, xlim=c(-135, -110), ylim=c(48.4, 59.8), col=colscheme, breaks=breaks, xaxt="n", yaxt="n")
lines(bdy.na, lwd=0.4)
points(ahccd.locations$lon, ahccd.locations$lat, bg = as.character(z_cut), pch = 21, cex=1.6, lwd=1.5)

legend("bottomleft", legend="AHCCD stations", bg = as.character(z_cut), pch = 21, pt.cex=1.6, pt.lwd=1.5, bty="n")

print(source)
}

lines(ts.x.cru,ts.y.cru, col="dodgerblue", lwd=3)
lines(ts.x.climatena,ts.y.climatena, col="black", lwd=2)
lines(c(1961,1990), rep(mean(ts.y.cru[ts.x.cru%in%1961:1990]),2), lty=2, col=1, lwd=2)
lines(c(2001,2020), rep(mean(ts.y.cru[ts.x.cru%in%2001:2020]),2), lty=2, col="dodgerblue", lwd=2)
lines(c(2001,2020), rep(mean(ts.y.climatena[ts.x.climatena%in%2001:2020]),2), lty=2, lwd=2)
boxplot(recent, add=T, at=2021, width=5, range=0)
legend("topleft", legend=c("ClimateNA", "CRU", "AHCCD stations"), pch=c(NA, NA, NA),
col = c(1,"dodgerblue", "gray80"), lty=c(1,1,1), lwd=c(3,2,2), bty="n")
## legend
plot(1, type="n", axes=F, xlab="", ylab="", xlim=c(0,1), ylim=c(0,1))
xl <- 0.2; yb <- 0.25; xr <- 0.8; yt <- 0.45
rect(head(seq(xl,xr,(xr-xl)/length(colscheme)),-1), yb, tail(seq(xl,xr,(xr-xl)/length(colscheme)),-1), yt, border=NA, col=colscheme)
rect(xl, yb, xr, yt)
# legend.labels <- if(e==3) paste(round(2^seq(0-lim,lim,lim/2), 2)*pct-pct, "%", sep="") else paste0(round(seq(0-lim,lim,lim/2), 2)*pct, "\u00B0", "C")
labels <- if(elements[e]=="Pr") paste(round(2^seq(lim.lower,lim.upper,(lim.upper-lim.lower)/4), 2)*pct-pct, "%", sep="") else paste0(round(seq(lim.lower,lim.upper,(lim.upper-lim.lower)/4), 2)*pct, "\u00B0", "C")
text(seq(xl,xr,(xr-xl)/(length(labels)-1)),rep(yb,length(labels)),labels,pos=1,cex=1.5,font=1, offset=0.5)
text(mean(c(xl,xr)), yt+0.01, paste("Change in", month.name[m], element.names[e], "\n1981-2010 to 2001-2020"), pos=3, cex=1.5, font=2)

dev.off()

month.abb[m]
# }
}
print(elements[e])
# }
}

0 comments on commit e81a086

Please sign in to comment.