## R-script that reads area-based climate indicators from Met Norway's thredds server
## Rasmus Benestad, 2023-06-06
## Ref: Benestad, R.E., Lussana C., Lutz J., Dobler A. Landgren O.A., Haugen, J.E. Mezghani A., Casati B. Parding K.M. (2022)
## "Global hydro-climatological indicators and changes in the global hydrological cycle and rainfall patterns",
## PLOS Climate, PCLM-D-21-00079R1, DOI: 10.1371/journal.pclm.0000029

## esd - freely available fro https://github.com/metno/esd
library(esd)

## Precipitation fractional surface area
a1mm <- AreaClimateIndicator(param='precip',is='50S-50N',threshold=1)
y <- AreaClimateIndicator(plot=FALSE)
lines(y,lwd=2,col='red')
grid()
legend(1940,0.41,c("global","50S-50N"),col=c('red','black'),lty=1,bty='n')

## Area with heavy precipitation
a10mm <- AreaClimateIndicator(param='precip',threshold=10,FUN='max',plot=FALSE)
a20mm <- AreaClimateIndicator(param='precip',threshold=20,FUN='max',plot=FALSE)
a30mm <- AreaClimateIndicator(param='precip',threshold=30,FUN='max',plot=FALSE)
a50mm <- AreaClimateIndicator(param='precip',threshold=50,FUN='max',plot=FALSE)
par(bty='n')
plot(100*merge(a1mm/mean(subset(a1mm,it=c(1961,1990))),a10mm/mean(subset(a10mm,it=c(1961,1990))),
               a20mm/mean(subset(a20mm,it=c(1961,1990))),a30mm/mean(subset(a30mm,it=c(1961,1990))),
               a50mm/mean(subset(a50mm,it=c(1961,1990)))),
     plot.type='single',main='ERA5: % change in global area with precipitation exceeding',
     col=c("black","darkgreen","lightblue","blue","darkblue"),lwd=c(2,2,2,2,2),xlab='',ylab='%')
grid()
polygon(c(1961,1990,1990,1961,1961),c(80,80,180,180,80),border=NA,col=rgb(0,0,0,0.05))
legend(1940,140,c('1 mm/day','10 mm/day','20 mm/day','30 mm/day','50 mm/day'),cex=0.7,
       lty=1,col=c("black","darkgreen","lightblue","blue","darkblue"),lwd=2,bty='n')

## Artic fractional area with temperatures above freezing
tas <- AreaClimateIndicator(param='t2m',is='Arctic',threshold=0,FUN='mean',plot=FALSE)
tax <- AreaClimateIndicator(param='tmax',is='Arctic',threshold=0,FUN='mean',plot=FALSE)
tan <- AreaClimateIndicator(param='tmin',is='Arctic',threshold=0,FUN='mean',plot=FALSE)
par(bty='n')
plot(100*merge(tas,tax,tan),plot.type='single',main='ERA5: % change in Arctic area with above freezing temperature',
     col=c("black","red","blue"),lty=c(1,2,2),lwd=2,xlab='',ylab='%')
grid()
legend(1940,48,c('TAS','TAX','TAN'),cex=0.7,
       col=c("black","red","blue"),lty=c(1,2,2),lwd=2,bty='n')


## Global surface area with extremely high temperatures
q95 <- function(x) return(quantile(probs=0.95))
tgt40 <- AreaClimateIndicator(param='t2m',FUN='q95',threshold=40)
tgt30 <- AreaClimateIndicator(param='t2m',FUN='q95',threshold=30)