bcgov / climr-mosaic-arctic

Development of 1981-2010 gridded temperature and precipitation normals for the Yukon and NWT.
Apache License 2.0
0 stars 0 forks source link

12-month plot layout #8

Open cmahony opened 4 months ago

cmahony commented 4 months ago

@sbeale007 here is a plot layout that you can adapt if you feel inspired

WRFCCRN prClimatology

here is the code for the plot. its from a few years ago so you might want to refactor to terra.

## ==============================
## Precipitation plots
## ==============================

library(raster)
library(ncdf4)
library(RColorBrewer)
library(scales)
library(rworldmap)
library(rworldxtra)
data(coastsCoarse)
data(countriesHigh)

monthdays <- c(31, 28.25, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)
monthcodes <- c("01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12")

png(filename=paste("results\\WRFCCRN.prClimatology","png",sep="."), type="cairo", units="in", width=6.5, height=7, pointsize=12, res=300)

mat <- matrix(c(2,3,4,5,13,1,1,6,12,1,1,7,11,10,9,8),4, byrow=T)   #define the plotting order
layout(mat, widths=c(1,1,1,1), heights=c(1,1,1,1))   #set up the multipanel plot

X <- raster(paste("inputs/WRFCCRN_pr_2000-2015.tif", sep=""))
studyarea <- extent(X)
studyarea[c(2,4)] <- 250000
X <- crop(X, studyarea)
values(X) <- log2(values(X))
values(X)[!is.finite(values(X))] <- NA
top <- 0.99
lim1 <- quantile(X, top)
values(X)[which(values(X)>lim1)] <- lim1
bottom <- 0.01
lim0 <- quantile(X, bottom)
values(X)[which(values(X)<lim0)] <- lim0

inc=0.01
breaks=seq(lim0-inc, lim1+inc, inc)
ColScheme=colorRampPalette(brewer.pal(9,"YlGnBu"))(length(breaks)-1)

## ANNUAL map
par(mar=c(0.2, 0.2, 0.2, 0.2))
image(X, col=ColScheme, breaks=breaks, maxpixels=ncell(X), yaxt="n", xaxt="n")
mtext("Annual", line=-1.5, adj=0.95, side=3, cex=1)
plot(bdy.bc.LCC, add=T, lwd=0.5, border=alpha("white", 0.25))
box()

## ANNUAL legend
xl <- studyarea[1]+150000; yb <- studyarea[3]+70000; xr <- studyarea[1]+250000; yt <- studyarea[3]+500000
rect(xl-120000,  yb-50000,  xr+120000,  yt+50000, col=alpha("white", 1), border="black")
rect(xl,  head(seq(yb,yt,(yt-yb)/length(ColScheme)),-1),  xr,  tail(seq(yb,yt,(yt-yb)/length(ColScheme)),-1),  border=NA, col=ColScheme)
rect(xl,  yb,  xr,  yt)
labels <- round(2^(quantile(breaks, c(bottom, 0.25,0.5,0.75, top), na.rm=T)),0) 
text(rep(xr,length(labels)),seq(yb,yt,(yt-yb)/(length(labels)-1)),labels,pos=4,cex=0.9,font=1, offset=0.3)
text(xl-25000, mean(c(yb,yt))-25000, "Mean precipitation\n(mm/day)", srt=90, pos=3, cex=0.85, font=1)

## MONTHLY maps
for(i in 1:12){

  X <- raster(paste("inputs/WRFCCRN_pr_2000-2015_", monthcodes[i], ".tif", sep=""))
  X <- crop(X, studyarea)
  values(X) <- log2(values(X))
  values(X)[!is.finite(values(X))] <- NA
  values(X)[which(values(X)>lim1)] <- lim1
  values(X)[which(values(X)<lim0)] <- lim0

  image(X, col=ColScheme, breaks=breaks, maxpixels=ncell(X), yaxt="n", xaxt="n")
  mtext(month.name[i], line=-1.5, adj=0.95, side=3, cex=0.8)
  plot(bdy.bc.LCC, add=T, lwd=0.5, border=alpha("white", 0.25))
  box()

}

dev.off()