# Rasmus Benestad & # DrittvĂ¦r: # Function used to tell how many valid data points there are: library(esd) # Prepare the predictand data: Do this only once since it takes long time: if (!file.exists('Tx.spain.rda')) { # Get the rain gauge data: select stations with at least 50 years Y <- station(src='ecad',cntr='Spain',param='tmax',nmin=50) Y <- subset(Y,it=as.Date(c('1950-01-01','2013-12-31'))) # 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 good <- (apply(Y,2,min,na.rm=TRUE)> -20) Y <- subset(Y,is=(1:length(good))[good]) Y <- subset(Y,is=list(lat=c(35,45))) # Get the anomalies: # Estimate annual aggregates: wet-day mean and frequency: Ma <- annual(Ya) # wet-day mean Sa <- annual(Ya,FUN='sd') # wet-day frequency save(file='Tx.spain.rda',Y,Ma,Sa) } else load('Tx.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('SST.rda')) { sst <- retrieve('data/ncep/sst.mnmean.ncep.ext-reconstr.nc', lon=c(-70,30),lat=c(30,70),param='sst') SST <- annual(sst) save(file='SST.rda',SST) } else load('SST.rda') # Use a smaller subregion for SLP: #SLP <- subset(SLP,is=list(lon=c(-40,30),lat=c(45,70))) # Transform to percentiles to give equal weight to stations # Use anomaly to compute anomalies # plot the anomalies with transparent colours to emphasize common # [col=rgb(1,0,0,0.1)] # PCA to emphasis the common behaviour: compute and plot # EOFs to reduce data size of predictors: compute and plot # CCA for diagnostics: compute and plot # Use loc() to find Santander: for (i in ) { # Check to see if data is close to being normal # Downscale the mean and the standard deviation respectively # Plot the results from the downscaling: # Estimate the annual probability for anomalous T > 3deg C (Tc) # use zoo() and pnorm() Tc <- 3 # Some fixes to make the results inherit correct attributes class(pr) <- class(Ma) pr <- attrcp(subset(Ma,is=1),pr) # fix some house-keeping attributes: attr(pr,'variable') <- 'Pr' attr(pr,'unit') <- 'fraction' attr(pr,'size') <- zoo(annual(Ya,FUN='iv')) # Estimate number of warm events and plot (use NE() from examples) # Compare with the observed number of events: find the actual number of events # use annual and count() }