# Rasmus Benestad # 2014-09-12 SPECS hands on workshop, Santander # Demonstration script for the esd package # Clear the memory rm(list=ls()) # Function used to tell how many valid data points there are: iv <- function(x) sum(is.finite(x)) # Function that computes probabilies for precio, assuming exponential distr. Pr <- function(x,x0=10) { # Pr(X > x) mu <- annual(x,FUN='wetmean') fw <- annual(x,FUN='wetfreq') size <- annual(x,FUN='iv') Pr <- zoo(coredata(fw)*exp(-x0/coredata(mu)), order.by=index(mu)) Pr <- attrcp(mu,Pr) attr(Pr,'variable') <- paste('Pr(X>',x0,'mm/day)') attr(Pr,'unit') <- 'fraction' attr(Pr,'x0') <- x0 attributes(size) <- NULL attr(Pr,'size') <- size class(Pr) <- class(mu) Pr } # Function predicting number of events. NE <- function(Pr) { nel <- zoo(qbinom(p=0.05, size=coredata(attr(Pr,'size')),prob=coredata(Pr)), order.by=index(Pr)) nem <- zoo(qbinom(p=0.5, size=coredata(attr(Pr,'size')),prob=coredata(Pr)), order.by=index(Pr)) neh <- zoo(qbinom(p=0.95, size=coredata(attr(Pr,'size')),prob=coredata(Pr)), order.by=index(Pr)) ne <- merge(nel,nem,neh) ne <- attrcp(Pr,ne) attr(ne,'variable') <- paste('N(X>',x0,'mm/day)') attr(ne,'unit') <- 'days' attr(ne,'location') <- rep(loc(Pr),3) attr(ne,'longitude') <- rep(lon(Pr),3) attr(ne,'latitude') <- rep(lat(Pr),3) attr(ne,'altitude') <- rep(alt(Pr),3) class(ne) <- class(Pr) ne } library(esd) # Prepare the predictand data: Do this only once since it takes long time: if (!file.exists('Pr.spain.rda')) { # Get the rain gauge data: select stations with at least 50 years Y <- station(src='ecad',cntr='Spain',param='precip',nmin=50) Y <- subset(Y,it=as.Date(c('1950-01-01','2013-12-31'))) Y <- subset(Y,is=list(lat=c(35,45))) # If we want to only include the warm months: ngood <- apply(Y,2,iv) # Assume more than 50*365.25 = 18262.5 valid data points Y <- subset(Y,is=(1:length(ngood))[ngood > 22000]) # Y <- subset(Y,it='jjas') # Also need to do the same for the predictors # Estimate annual aggregates: wet-day mean and frequency: Mu <- annual(Y,FUN='wetmean') # wet-day mean Fw <- annual(Y,FUN='wetfreq') # wet-day frequency Nd <- annual(Y,FUN='count',threshold=1) # number of data points save(file='Pr.spain.rda',Y,Mu,Fw,Nd) } else load('Pr.spain.rda') if (!file.exists('slp.ncep.rda')) { slp <- retrieve('data/ncep/slp.mon.mean.nc',lon=c(-70,30),lat=c(30,70)) SLP <- annual(slp) save(file='slp.ncep.rda',SLP) } else load('slp.ncep.rda') if (!file.exists('es.rda')) { sst <- retrieve('data/ncep/sst.mnmean.ncep.ext-reconstr.nc', lon=c(-70,30),lat=c(30,70),param='sst') es <- annual(C.C.eq(sst)) save(file='es.rda',es) } else load('es.rda') # regridding: put SLP on the same grid as es: use 'is' for space dimension slp <- regrid(SLP,is=es) # Visualisation: map(slp) # Extract a timeslice: use 'it' for temporal dimension # Error when using 'as'Date' in this function :-( will be fixed :-) slp <- subset(SLP,it=c(1961,1990)) # Plot the area mean plot(slp) # Use a smaller subregion for SLP: SLP <- subset(SLP,is=list(lon=c(-40,20),lat=c(40,50))) # Transform to percentiles to give equal weight to stations MU <- 100*Mu/colMeans(Mu,na.rm=TRUE) plot(MU,col=rgb(0,0,1,0.1)) # PCA to emphasis the common behaviour pca.mu <- PCA(MU) plot(pca.mu) pca.fw <- PCA(Fw) plot(pca.fw) # EOFs to reduce data size: eof.slp <- EOF(SLP) plot(eof.slp) eof.es <- EOF(es) plot(eof.es) # CCA for diagnostics cca <- CCA(pca.mu,eof.slp) plot(cca) #ns <- dim(Mu)[2] #r2 <- matrix(rep(NA,ns*4),ns,4) #colnames(r2) <- c('mu','fw','ncwd','ncdd') #xc <- r2 # Fix a glitch... attr(Fw,'unit') <- rep('fraction',length(loc(Fw))) for (i in 44) { # downscale the annual wet-day mean: z.mu <- DS(subset(Mu,is=i),eof.es,eofs=1:12) # downscale the annual wet-day frequency: z.fw <- DS(subset(Fw,is=i),eof.slp) # The following commented lines give statistics on the spell length: # nc <- annual(spell(subset(Y,is=i),threshold=1)) # z.ncwd <- DS(subset(nc,is=1),eof.slp,eofs=1:12) #rm z.ncdd <- DS(subset(nc,is=2),eof.slp,eofs=1:12) # The following line collects statistics for this station if one # wants to make a table for a group of stations # r2[i,] <- c(attr(z.mu,'quality'),attr(z.mu,'quality'), # attr(z.ncwd,'quality'),attr(z.ncdd,'quality')) # Plot the results from the downscaling: plot(z.fw) # The wet-day 95-th percentile: # q95 <- aggregate(subset(Y,is=i),year,FUN='exceedance', # fun='quantile',probs=0.95) q95 <- annual(subset(Y,is=i),FUN='exceedance',fun='quantile',probs=0.95) plot(q95) # compare with the 95-percentile for an exponential distribution: lines(-log(1-0.95)*subset(Mu,is=i),lwd=2) # Estimate the annual probability for X > 10mm/day pr <- Pr(subset(Y,is=i),x0=10) plot(pr) # Compare with the downscaled probabilities: pr.ds <- as.station(z.fw)*exp(-10/coredata(z.mu)) # Need to fix a few attributes to get things to work: attr(pr.ds,'variable') <- 'Pr' # variable name index(pr.ds) <- year(pr.ds) # use same time units lines(pr.ds) # Estimate number of heavy rain events per year based on probability: # plotting the 5, 50, and 95 percentiles from a binomial distribution: ne <- NE(pr) # Compare with the observed number of events: no <- annual(subset(Y,is=i),FUN='count',threshold=10) plot(ne) lines(no,lwd=3) }