acs<-function(execesses){ # Compute average cluster size for a given three-dimensional array of excesses # # Description: # # Returns average cluster size # # Usage: acs(execesses) # # Input: # # excesses: p (space) x q (space) x n (time) array of excesses # # Output: # # A p x q matrix containing the average cluster size # # Authors: # # Caio Coelho 28 Feb 2006 # Chris Ferro avenumbexceed<-function(y){ x<-matrix(y,ncol=3,nrow=length(y)/3,byrow=T) sum(x > 0,na.rm=T) / sum(apply(x > 0, 1, any,na.rm=T),na.rm=T) } apply(execesses,c(1,2),avenumbexceed) } anomalyfield <- function(data,firstmonth,firstyear,calc="mac") { # Calculate anomalies by subtracting out mean annual cycle or long-term mean of # a three-dimensional p x q x n array of monthly mean values. The structure of this # array must be as follows: - first dimension p is longitude points # - second dimension q is latitude points # - third dimension n is time (months of the year in # ascending order, e.g. Jan 1948, Feb 1948, ..., Sep 2005) # # Description: # # Returns a p x q x n array containing anomalies # # Usage: # # anomalyfield(data,firstmonth,firstyear,calc="mac") # # Input: # # data: p x q x n array (see data structure description above) # firstmonth: number of the first month containing data in the # data array (e.g. 1 (for January)) # firstyear: first year of data in the data array (e.g. 1948) # calc: string of characters that can be "mac" (Default) to calculate # anomalies by subtracting the mean annual cycle, or "ltm" to # calculate anomalies by subtracting the long-term mean. # # Output: # # $anomaly: p x q x n array containing anomalies # # Authors: # # Caio Coelho 28 October 2005 # Chris Ferro # # Examples: # # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 100) # data <- array(rnorm(prod(dim)), dim) # out<-anomalyfield(data,1,1948) # out<-anomalyfield(data,5,1970) # out<-anomalyfield(data,5,1970,calc="ltm") if(dim(data)[3] < 12) {stop("Not possible to compute either mean annual cycle or long-term mean. Less than a year of data available")} if(calc=="mac"){ mean.month <- function(data) { apply(matrix(data, nrow = 12), 1, mean) } anomaly <- data if (firstmonth == 1) { initialmonth <- firstmonth finalmonth <- floor(dim(data)[3]/12)*12 } if (firstmonth != 1) { initialmonth <- 12-firstmonth+2 finalmonth <- (floor((dim(data)[3]-(12-firstmonth+1))/12)*12)+initialmonth-1 } m <- aperm(apply(data[,,initialmonth:finalmonth],c(1,2),mean.month),c(2,3,1)) anom <- data[,,initialmonth:finalmonth] for(i in 1:((finalmonth-initialmonth+1)/12)) { if(firstmonth == 1) anom[,,(((i-1)*12)+1):(12*i)] <- data[,,(((i-1)*12)+1):(12*i)] - m else anom[,,(((i-1)*12)+1):(12*i)] <- data[,,(((i-1)*12)+initialmonth):((12*i)+initialmonth-1)] - m } if (firstmonth == 1 & finalmonth == dim(data)[3]) { anomaly <- anom } if (firstmonth == 1 & finalmonth < dim(data)[3]) { anomaly[,,firstmonth:finalmonth] <- anom j <- 1 for(i in (finalmonth+1):dim(data)[3]){ anomaly[,,i] <- data[,,i] - m[,,j] j <- j+1 } } if (firstmonth != 1 & finalmonth == dim(data)[3]) { j <- 1 for(i in firstmonth:12){ anomaly[,,j] <- data[,,j] - m[,,i] j <- j+1 } anomaly[,,initialmonth:finalmonth] <- anom } if (firstmonth != 1 & finalmonth < dim(data)[3]){ j <- 1 for(i in firstmonth:12){ anomaly[,,j] <- data[,,j] - m[,,i] j <- j+1 } anomaly[,,initialmonth:finalmonth] <- anom j <- 1 for(i in (finalmonth+1):dim(data)[3]){ anomaly[,,i] <- data[,,i] - m[,,j] j <- j+1 } } if(firstmonth == 1) y1 <- firstyear if(firstmonth != 1) y1 <- firstyear+1 y2<- y1+((dim(anom)[3])/12)-1 cat(paste("Mean annual cycle was computed with",as.character(dim(anom)[3]),"months of data","\n",sep=" ")) cat(paste("from Jan",as.character(y1),"to Dec",as.character(y2),"(i.e.",as.character(((dim(anom)[3])/12)),"years)","\n",sep=" ")) } if(calc=="ltm"){ anomaly <- data m<-apply(data,c(1,2),mean) for(i in 1:dim(data)[3]) { anomaly[,,i] <- data[,,i] - m } cat(paste("Long-term mean was computed with",as.character(dim(data)[3]),"months of data","\n",sep=" ")) } invisible(list(anomaly=anomaly)) } applyfield <- function(array3d,fun = mean, ...) { # Compute basic statistics (e.g. mean, variance, min, max, skewness, quantile) # for a given three-dimensional array with first two space dimensions (e.g. # longitude and latitude) and third time dimension # # Description: # # Returns a matrix with the statistics of interest, specified by the user # (see below the list of statistics that can be computed by this funtion), # which is computed over the time dimension of a three-dimensional array # for all space points of the array. # # Usage: # # applyfield(array3d,fun) # # Input: # # array3d: a three-dimensional array with p longitude points and q latitude # points as the first two dimensions and n as the third time # dimension # # fun: Name of the statistics to be computed # that must be one of the following options: mean, var, # min, max, momentskew, ykskew or quantil, # where momentskew is a moment measure of skewness and # ykskew is the Yule-Kendall skewness statistics # # ...: Additional argument passed to `fun' (e.g. the quantile p to be # computed). The default quantile p to be computed is 0.5 (i.e. # the median). # # Output: # # A matrix of the computed statistics # # Authors: # # Dag Johan Steinskog 9 June 2005 # Caio Coelho # # Examples: # # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 100) # data <- array(rnorm(prod(dim)), dim) # applyfield(data, mean) # applyfield(data, var) # applyfield(data, min) # applyfield(data, max) # applyfield(data, ykskew) # applyfield(data, quantil) # compute the median field # applyfield(data, quantil,p=0.9) # compute the 90th quantile field apply(array3d,c(1,2),fun, ...) } momentskew <- function(z) {mean(((z-mean(z,na.rm = T))/sd(z,na.rm = T))^3,na.rm = T)} ykskew <- function(z) {q<-quantile(z,na.rm = T) (q[2]-2*q[3]+q[4])/(q[4]-q[2]) } quantil <- function(w,p=0.5) {quantile(w,p,na.rm = T)} applyfieldnew <- function(lon,lat,array3d,fun = mean, ...,nonmissing=0.5) { # Compute basic statistics (e.g. mean, variance, min, max, skewness, quantile) # for a given three-dimensional array with first two space dimensions (e.g. # longitude and latitude) and third time dimension. # # Note: This function allows the user to specify through the parameter # nonmissing (defaul is 0.5) the acceptable percentage of nonmissing # values in the time series of each grid point of the three-dimensional # array. This is the main difference between this function and the # function applyfield.r # # Description: # # Returns a matrix with the statistics of interest, specified by the user # (see below the list of statistics that can be computed by this funtion), # which is computed over the time dimension of a three-dimensional array # for all space points of the array. # # Usage: # # applyfieldnew(lon,lat,array3d,fun,nonmissing) # # Input: # # lon: vector with p longitude values # # lat: vector with q latitude values # # array3d: a three-dimensional array with p longitude points and q latitude # points as the first two dimensions and n as the third time # dimension # #nonmissing: Only grid points with fraction given by 'nonmissing' # (between 0 and 1) of non-missing values are used to compute # the statistics below. Default is 0.5. # # fun: Name of the statistics to be computed # that must be one of the following options: mean, var, # min, max, momentskew, ykskew or quantil, # where momentskew is a moment measure of skewness and # ykskew is the Yule-Kendall skewness statistics # # ...: Additional argument passed to `fun' (e.g. the quantile p to be # computed). The default quantile p to be computed is 0.5 (i.e. # the median). # # Output: # # A matrix of the computed statistics # # Authors: # # Dag Johan Steinskog 9 June 2005 # Caio Coelho # # Examples: # # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 100) # data <- array(rnorm(prod(dim)), dim) # applyfieldnew(x,y,data, mean) # applyfieldnew(x,y,data, var) # applyfieldnew(x,y,data, min) # applyfieldnew(x,y,data, max) # applyfieldnew(x,y,data, ykskew) # applyfieldnew(x,y,data, quantil) # compute the median field # applyfieldnew(x,y,data, quantil,p=0.9) # compute the 90th quantile field # reshape three-dimensional array into a data matrix with # time as first dimension and space as sencond dimension data <- reshapefield(lon,lat,array3d)$out # compute percentage of non-missing values at each grid point aux <- apply(data,2,function(x)sum(!is.na(x))/(length(x))) # identify grid points with less than 50% missing values index <- (1:length(aux))[aux >= nonmissing] out<-apply(data[,index],2,fun, ...) out1 <- rep(NA,dim(data)[2]) out1[index]<-out reshapefield(lon,lat,t(as.matrix(out1)))$out[,,1,drop=T] } momentskew <- function(z) {mean(((z-mean(z,na.rm = T))/sd(z,na.rm = T))^3,na.rm = T)} ykskew <- function(z) {q<-quantile(z,na.rm = T) (q[2]-2*q[3]+q[4])/(q[4]-q[2]) } quantil <- function(w,p=0.5) {quantile(w,p,na.rm = T)} bluered <- function(data, scale = c("equal", "linear", "split"), white = median(data), yellow = 0, cyan = 0, invert = FALSE, format = c("hex", "rgb")) { # Generate a Colour Scale for Plotting # # Description: # # Returns a colour scale corresponding to specified data values # or a specified number of colours. # Usage: # # bluered(data, scale = "equal", white = median(data), yellow = 0, # cyan = 0, invert = FALSE, format = "hex") # # Arguments: # # data: A vector of distinct data values corresponding to the colours # or the number of colors (>= 1) to be created. # # scale: Controls the colour gradients (see Details). Must be either # `"equal"' (the default), `"linear"' or `"split"' (or any # unique partial match). # # white: The data value corresponding to the colour white. # # yellow: A number between 0 and 1 indicating the amount of yellow to # insert between white and red. If 0 (the default) then no yellow # is used. Larger values indicate more yellow. # # cyan: A number between 0 and 1 indicating the amount of cyan to # insert between white and blue. If 0 (the default) then no cyan # is used. Larger values indicate more cyan. Not implemented. # # invert: Logical: if `TRUE', negative data values are red and positive # values are blue; if `FALSE' (the default) then vice-versa. # # format: The format in which the colours are returned (see Value). Must # be either `"hex"' (the default) or `"rgb"' (or any unique # partial match). # # Details: # # One colour is created for each value in `data'. Colours range from # blue for data values less than `white' to red for data values # greater than `white'. If `scale' equals `"linear"' then colour # intensity scales linearly with absolute data value: highest # intensity corresponding to the maximum absolute value of `data'. # If `scale' equals `"split"' then linear intensity scales are # constructed separately for blue and red colours so that the # maximum intensity is attained by both blue and red. If `scale' # equals `"equal"' then intensities are equally spaced regardless of # the data values. # # Value: # # If `format' is `"hex"' then returns a character vector of # hexadecimal RGB numbers. If `format' is `"rgb"' then returns a # matrix with three columns containing RGB intensities. # # See Also: # # `col2rgb', `colours', `rgb' # # Author: # # Chris Ferro 4 May 2005 # # Examples: # # image(volcano, col = bluered(seq(100, 200, 10), "linear")) # image(volcano, col = bluered(11)) # image(volcano, col = bluered(11,yellow=1)) if (length(data)==1){data<-seq(0,data-1,by=1)} scale <- match.arg(scale) format <- match.arg(format) data <- unique(data) - white if(invert) data <- -data data <- sort(data) nb <- sum(data < 0) nr <- sum(data > 0) n <- length(data) if(scale == "equal") { if(nb == 0) r <- rep(1, n) else r <- pmin(1:n - 1, nb) / nb if(nr == 0) b <- rep(1, n) else b <- pmin(n - 1:n, nr) / nr } else if(scale == "linear") { dmax <- max(abs(data)) if(nb == 0) r <- rep(1, n) else r <- pmin(data + dmax, dmax) / dmax if(nr == 0) b <- rep(1, n) else b <- pmin(dmax - data, dmax) / dmax } else { if(nb == 0) r <- rep(1, n) else r <- pmax(data[1] - data, data[1]) / data[1] if(nr == 0) b <- rep(1, n) else b <- pmin(data[n] - data, data[n]) / data[n] } g <- pmin(r, b) y <- min(max(yellow, 0), 1) g[b < 1] <- 1 - (1 - y) * (1 - b[b < 1]) - y * (1 - b[b < 1])^2 b <- 1 - (1 + y) * (1 - b) + y * (1 - b)^2 b[b < 1] <- b[b < 1] * (1 - y) # c <- min(max(cyan, 0), 1) # g[r < 1] <- 1 - (1 - c) * (1 - r[r < 1]) - c * (1 - r[r < 1])^2 # r <- 1 - (1 + c) * (1 - r) + c * (1 - r)^2 if(invert) { r <- rev(r) g <- rev(g) b <- rev(b) } if(format == "hex") return(rgb(r, g, b)) cbind(r, g, b) } boundexcesses <- function(paretoparam){ # Compute upper bound of excesses for a given array of Generalized Pareto # distribution parameters # # Description: # # Returns upper bound of excesses # # Usage: boundexcesses(paretoparam) # # Input: # # paretoparam: 2 x p x q array with scale (first level of the array) and # shape (second level of the array) parameters # # Output: # # A matrix containing the upper bound of excesses # # Authors: # # Caio Coelho 28 Feb 2006 paretoparam[2,,][paretoparam[2,,]>=-0.008]<-0 b<--(paretoparam[1,,]/paretoparam[2,,]) b[b==-Inf]<-10000 b } corfield <- function(array3d,ts, ...) { # Compute correlation between each point of a three-dimensional array, with # first two space dimensions and third time dimension, and a given vector of # the same length as the time dimension of the three-dimentional array. # # Description: # # Returns a space matrix with the correlations # # Usage: # # corfield(array3d,ts, ...) # # Input: # # array3d: three-dimensional array with first two space dimensions # (e.g. longitude and latitude) and third time dimension # # ts: vector containing the time series to be correlated with the # three-dimensional array # # ...: Additional argument passed to the correlation function # (e.g. method = "pearson", method = "kendall", method = "spearman", # use = "all.obs", use = "complete.obs", use = "pairwise.complete.obs") # Default computes pearson correlation coefficient (method = "pearson") # using all data available (use = "all.obs") # # Output: # # Matrix of correlations # # Authors: # # Dag Johan Steinskog 9 June 2005 # Caio Coelho # # Examples: # # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 100) # data <- array(rnorm(prod(dim)), dim) # timeseries<-rnorm(100) # corfield(data,timeseries) # corfield(data,timeseries,method = "pearson") # corfield(data,timeseries,method = "kendall") # corfield(data,timeseries,method = "spearman") correl <- function(z) {cor(z,ts,use="complete.obs", ...)} apply(array3d, c(1, 2), correl) } covfield <- function(array3d,ts, ...) { # Compute covariance between each point of a three-dimensional array, with # first two space dimensions and third time dimension, and a given vector of # the same length as the time dimension of the three-dimentional array. # # Description: # # Returns a space matrix with the covariances # # Usage: # # covfield(array3d,ts, ...) # # Input: # # array3d: three-dimensional array with first two space dimensions # (e.g. longitude and latitude) and third time dimension # # ts: a vector containing the time series to be covariate with the # three-dimensional array # # ...: Additional argument passed to the covariance function # (e.g. method = "pearson", method = "kendall", method = "spearman", # use = "all.obs", use = "complete.obs", use = "pairwise.complete.obs") # Default computes pearson covariance (method = "pearson") # using all data available (use = "all.obs") # # Output: # # Matrix of covariances # # Authors: # # Dag Johan Steinskog 25 October 2005 # Caio Coelho # # Examples: # # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 100) # data <- array(rnorm(prod(dim)), dim) # timeseries<-rnorm(100) # covfield(data,timeseries) # covfield(data,timeseries,method = "pearson") # covfield(data,timeseries,method = "kendall") # covfield(data,timeseries,method = "spearman") covar <- function(z) {cov(z,ts,use="complete.obs", ...)} apply(array3d, c(1, 2), covar) } detrend <- function(y,smooth=F,...){ # Description: # # Removes time trend from a time series # # Usage: # # detrend(y,smooth=FALSE,...) # # Arguments: # # y: a vector containing a time series of data # # smooth: Logical. If FALSE (the default) removes linear time trend by # fitting ordinary least squares regression to the time series y. # If TRUE, then time trend is removed by fitting a local weighted # regression (lowess) to the time series y # # ...: Additional arguments passed to lowess # # Value: # # Vector with the detrended data. # # Author: # # Caio Coelho 9 Jan 2006 # # Examples: # # ts <- rnorm(100,25,2) # detrend(ts) # detrend(ts,smooth=TRUE) # detrend(ts,smooth=TRUE,f=1/3) if(!smooth){ x<-seq(1:length(y)) lm(y ~ x)$res } else{ y-lowess(y,...)$y } } detrendfield <- function(array3d,smooth=F,...){ # Description: # # Removes time trend at each grid point of a three-dimensional array of data # with first two space dimensions (e.g. longitude and latitude) and # third time dimension # # Usage: # # detrendfield(array3d,smooth=FALSE,...) # # Arguments: # # array3d: a three-dimensional array, with first two space dimensions (e.g. # longitude and latitude) and third time dimension # # smooth: Logical. If FALSE (the default) removes linear time trend by # fitting ordinary least squares regression at each grid. If TRUE, # then time trend is removed by fitting a local weighted regression # (lowess) at each grid point # # ...: Additional arguments passed to lowess # # Value: # # Three-dimensional array with the detrended data. # # Author: # # Caio Coelho 9 Jan 2006 # # Examples: # # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 100) # data <- array(rnorm(prod(dim),25,2), dim) # detrendfield(data) # detrendfield(data,smooth=TRUE) # detrendfield(data,smooth=TRUE,f=1/3) aperm(apply(array3d, c(1,2), detrend,smooth,...),c(2,3,1)) } eof <- function(lon,lat,array3d, ...) { # Compute EOFs and PCs of a three-dimensional p x q x n array, with first # two space dimensions p (longitude) and q (latitude), and third time # dimension n. # # Description: # # Returns principal component time series, spatial patterns (loadings), # standard deviation of each principal component, percentage of total # variance accounted by each principal component, a vector of longitudes # and a vector of latitudes # # Usage: # # eofs(lon,lat,array3d) # # Input: # # array3d: a three-demensional array with p longitude points and q latitude # points as the first two dimensions and n as the third time # dimension # # lon: vector of longitudes # # lat: vector of latitudes # # ...: Additional arguments to be passed for prcomp function such as, # retx: logical value indicating whether the rotated variables # should be returned. # center: logical value indicating whether the variables should # be zero centered. Alternately, a vector of # length equal the number of columns of x can be supplied. # The value is passed to scale. # scale.: logical value indicating whether the variables should # be scaled to have unit variance before the analysis takes # place. The default is FALSE for consistency with S, but # in general scaling is advisable. Alternatively, a vector # of length equal the number of columns of x can be supplied. # The value is passed to scale. # tol: value indicating the magnitude below which components # should be omitted. (Components are omitted if their # standard deviations are less than or equal to tol times the # standard deviation of the first component). With the default # null setting, no components are omitted. Other settings for # tol could be tol = 0 or tol = sqrt(.Machine$double.eps), # which would omit essentially constant components. # Output: # # $pcs: matrix of principal components (first column contains # first principal component and so on) # $stdev: vector containing the standard deviation of each principal # component, starting from the first. # $eof: three-dimensional array containing the loadings. # $accountedvar: vector containing the amount of variance accounted for # each principal component # $lon: vector of longitudes # $lat: vector of latitudes # # # Authors: # # Dag Johan Steinskog 10 June 2005 # Caio Coelho # # Examples: # # lon <- seq(50,70,20) # lat <- seq(-10,10,20) # t <- 20 # x <- array(rnorm(length(lon)*length(lat)*t,sd=10),dim=c(2,2,20)) # eofs(lon,lat,x) data <- reshapefield(lon,lat,array3d)$out aux <- prcomp(data, ...) eofs<-reshapefield(lon,lat,t(aux$r))$out varsquared <- (aux$sdev)^2 sumvarsquared <- sum(varsquared) accountedvar<-round((varsquared/sumvarsquared)*100,2) invisible(list(pcs=aux$x,stdev=aux$sdev,eofs=eofs,accountedvar=accountedvar,lon=lon,lat=lat)) } extract <- function(x, y, data, xmin, xmax, ymin, ymax, t=1) { # Extract a subset of data from either a two-dimensional map (matrix) with p # longitude and q latitude points, respectively, or a three-dimensional # array (field) with p longitude and q latitude as the first two space # dimensions and n (time) as the third dimension, by providing the latitude and # longitude range of interest and a subset of time slices # # Description: # # Returns a subset of a two-dimensional map (matrix) or a subset # of a three-dimensinal array (field) # # Usage: # # extract(x, y, data, xmin, xmax, ymin, ymax, t=1) # # Arguments: # # x: Vector of longitudes # # y: Vector of latitudes # # data: Two-dimensional map (matrix) or three-dimensional array of data. # Space dimensions of `data' must be `length(x)' by `length(y)' # respectively. # Three-dimensional array must have the following dimension order: # first dimension for longitude, second dimension for latitude, # and third dimension for time # # xmin: Lower bound longitude from where to start extracting data # # xmax: Upper bound longitude up to where to extract data # # ymin: Lower bound latitude from where to start extracting data # # ymax: Upper bound latitude up to where to extract data # # t: Subset of times slices to be extracted. Default is 1 (i.e. extract # only the first time slice # # Output: # # A list with four components is returned invisibly: # # x: Vector of longitudes within the selected range # # y: Vector of latitudes within the selected range # # data: Two-dimensional (matrix) or three-dimensional (array) # containing the selected data # # t: Extracted times slices # # # Author: # # Caio Coelho 28 October 2005 # Chris Ferro # # Examples: # # # Extract data from a three-dimensional array (field) # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 100) # data <- array(rnorm(prod(dim)), dim) # output <- extract(x, y,data,xmin=-12,xmax=12,ymin=32,ymax=58,t=c(1:3)) # output <- extract(x, y,data,xmin=-12,xmax=6,ymin=32,ymax=50,t=c(2,4,6)) # output <- extract(x, y,data,xmin=-18,xmax=11,ymin=36,ymax=59,t=c(10,20,30,40)) # # # Extract data from a two-dimensional map (matrix) # matr <- data[,,1] # output <- extract(x, y,matr,xmin=-12,xmax=12,ymin=32,ymax=58) # output <- extract(x, y,matr,xmin=-12,xmax=6,ymin=32,ymax=50) if (length(dim(data))==2) { xnew <- x ynew <- y xnew <- x[x >= xmin & x <= xmax] ynew <- y[y >= ymin & y <= ymax] datanew <- data[match(xnew, x), match(ynew, y), drop = FALSE] } if(length(dim(data))==3) { xnew <- x ynew <- y tnew <- t xnew <- x[x >= xmin & x <= xmax] ynew <- y[y >= ymin & y <= ymax] datanew <- data[match(xnew, x), match(ynew, y), t, drop = FALSE] } invisible(list(x = xnew, y = ynew, data = datanew, t = t)) } extractmap <- function(x, y, data, imonth, iyear, tmonth, tyear ,xmin=min(x), xmax=max(x), ymin=min(y), ymax=max(y)) { # Extract a matrix (map) on a specific date from a three-dimensional array # containing monthly data with first dimension p (longitude), second dimension # q (latitude) and third dimension n (time) # # Description: # # Returns the data of the specific date # # Usage: # # extractmap(x, y, data, imonth, iyear, tmonth, tyear , # xmin=min(x),xmax=max(x), ymin=min(y), ymax=max(y)) # # Arguments: # # x: Vector of longitudes # # y: Vector of latitudes # # data: Matrix or three-dimensional array of data. Space dimensions of # `data' must be `length(x)' by `length(y)' respectively. # Three-dimensional array must have the following dimension order: # first dimension for longitude, second dimension for latitude, # and third dimention for time # # imonth: Number of the first month of the time series of 'data' # # iyear: Number of the first year of the time series of 'data' # # tmonth: Number of the target month to be extracted from the # time series of 'data' # # tyear: Number of the target year to be extracted from the # time series of 'data' # # xmin: Lower bound longitude from where to start extracting data # # xmax: Upper bound longitude up to where to extract data # # ymin: Lower bound latitude from where to start extracting data # # ymax: Upper bound latitude up to where to extract data # # Output: # # A list with three components is returned invisibly: # # lon: Vector of longitudes within the selected range # # lat: Vector of latitudes within the selected range # # data: Matrix or three-dimentional array containing the selected data # # Authors: # # Caio Coelho 3 November 2005 # Chris Ferro # # Examples: # # # Extract data from a three-dimentional array (field) # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 100) # data <- array(rnorm(prod(dim)), dim) # output <- extractmap(x, y,data,1,1948,9,1951) # output <- extractmap(x, y,data,1,1948,9,1951,xmin=-12,xmax=12,ymin=32,ymax=58) # if (imonth == 1 & tmonth == 1) {t<-((tyear-iyear)*12)+1} if (imonth == 1 & tmonth != 1) {t<-((tyear-iyear)*12)+tmonth} if (imonth != 1) {t<-12-imonth+1+((tyear-iyear)*12)+tmonth-12} xnew <- x ynew <- y tnew <- t xnew <- x[x >= xmin & x <= xmax] ynew <- y[y >= ymin & y <= ymax] datanew <- data[match(xnew, x), match(ynew, y), t, drop = TRUE] cat(paste("Extracted data for month",as.character(tmonth),"of",as.character(tyear),"\n",sep=" ")) invisible(list(lon = xnew, lat = ynew, data = datanew)) } extractwindow <- function(x, y, data, imonth, iyear, month1, year1, month2, year2 ,xmin=min(x), xmax=max(x), ymin=min(y), ymax=max(y)) { # Extract field of monthly data for a given time period from a three-dimensional # array with first dimension p (longitude), second dimension q (latitude) and third # dimension n (time) # # Description: # # Returns field of data for a given time period # # Usage: # # extractwindow(x, y, data, imonth, iyear, month1, year1, # month2, year2 ,xmin=min(x), xmax=max(x), # ymin=min(y), ymax=max(y)) # # Arguments: # # x: Vector of longitudes # # y: Vector of latitudes # # data: Matrix or three-dimensional array of data. Space dimensions of # `data' must be `length(x)' by `length(y)' respectively. # Three-dimensional array must have the following dimension order: # first dimension for longitude, second dimension for latitude, # and third dimention for time # # imonth: Number of the first month of the time series of 'data' # # iyear: Number of the first year of the time series of 'data' # # month1: Number of the month from where to start extracting data # # year1: Number of the year from where to start extracting data # # month2: Number of the month up to where to extract data # # year2: Number of the year up to where to extract data # # xmin: Lower bound longitude from where to start extracting data # Default is the minimum value of x # # xmax: Upper bound longitude up to where to extract data # Default is the maximum value of x # # ymin: Lower bound latitude from where to start extracting data # Default is the minimum value of y # # ymax: Upper bound latitude up to where to extract data # Default is the maximum value of y # # # Output: # # A list with three components is returned invisibly: # # x: Vector of longitudes for the extracted data # # y: Vector of latitudes for the extracted data # # data: Extracted three-dimensional array (field) # # # Authors: # # Caio Coelho 3 November 2005 # Chris Ferro # # Examples: # # # Extract data from a three-dimensional array (field) # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 100) # data <- array(rnorm(prod(dim)), dim) # output <- extractwindow(x, y,data,1,1948,3,1948,5,1949) # output <- extractwindow(x, y,data,1,1948,3,1948,5,1949,xmin=-12,xmax=12,ymin=32,ymax=58) # output <- extractwindow(x, y,data,9,1948,11,1948,5,1949,xmin=-12,xmax=12,ymin=32,ymax=58) # if (imonth == 1 & month1 == 1) { t1<-((year1-iyear)*12)+1 if (month2 ==1){ t2<-((year2-iyear)*12)+1 } else{ t2<-((year2-iyear)*12)+month2 } } if (imonth == 1 & month1 != 1) { t1<-((year1-iyear)*12)+month1 if (month2 ==1){ t2<-((year2-iyear)*12)+1 } else{ t2<-((year2-iyear)*12)+month2 } } if (imonth != 1) { t1<-12-imonth+1+((year1-iyear)*12)+month1-12 t2<-12-imonth+1+((year2-iyear)*12)+month2-12 } xnew <- x ynew <- y tnew <- t xnew <- x[x >= xmin & x <= xmax] ynew <- y[y >= ymin & y <= ymax] datanew <- data[match(xnew, x), match(ynew, y), t1:t2, drop = FALSE] cat(paste("Extracted data from month",as.character(month1),"of",as.character(year1),"to month",as.character(month2),"of",as.character(year2),"\n",sep=" ")) invisible(list(x = xnew, y = ynew, data = datanew)) } fa<-function(x,y,reduction=F,nmodes=3){ # Bayesian forecast assimilation for cross-validated calibration and combination # of forecasts # # Description: This function calls functions faprior, falikelihood, # and faprediction and estimates the mean (ya) and the # square root of the diagonal elements of the # covariance matrix (d) of the posterior distribution # y|x ~ N(ya,d) # # Usage: fa(x,y,reduction,nmodes) # # Arguments: # x: n x p data matrix containing time series of p forecast variables # y: n x q data matrix containing time series of q observable # variables # reduction: Logical. If TRUE performs maximum covariance analysis # on the cross-covariance matrix (y'*x) and uses the # selected number of modes ('nmodes', see below) to model the # likelihood function. Default is FALSE. # nmodes: number of modes used in the likelihood modelling. Default # is 3. Argument only used if 'reduction' is TRUE. # # Outputs: # $ya: q x n matrix containing q predicted observable state time series # $d : q x q matrix containing predicted observable error variance # (i.e. the square root of the diagonal elementes of observable # error covariance matrix d.) # # Authors: Caio Coelho 13 Feb 2006 # David Stephenson # # Reference: Stephenson D. B., C.A.S. Coelho, F. J. Doblas-Reyes and # M. Balmaseda, 2005: Forecast Assimilation: A unified # framework for the combination of multi-model weather and # climate predictions. Tellus A . Vol. 57, 253-264. # # Example 1: Calibration and combination of forecasts for a single location # # xfa<-read.table("xfa.txt") # xfa.txt contains 3 time series of forecasts # yfa<-read.table("yfa.txt") # yfa.txt contains 1 time series of observations # out<-fa(xfa,yfa) # out$ya # combined and calibrated forecast mean # out$d # combined and calibrated forecast standard deviation # # Example 2: Calibration and combination of grided coupled model # forecasts # # xfa<-read.table("xfasst.txt") # xfa.txt contains 168 time series of forecasts # # Columns 1-56 of xfa.txt are forecasts produced by model A # # Columns 57-112 of xfa.txt are forecasts produced by model B # # Columns 113-168 of xfa.txt are forecasts produced by model C # yfa<-read.table("yfasst.txt") # yfa.txt contains 56 time series of observations # out<-fa(xfa,yfa,reduction=T,3) # out$ya # combined and calibrated forecast mean # out$d # combined and calibrated forecast standard deviation y <- as.matrix(y) x <- as.matrix(x) n <- dim(x)[1] ya<-matrix(nrow=n,ncol=ncol(y)) d <-matrix(nrow=n,ncol=ncol(y)) if(!reduction){ # Cross-validation loop nn <- 1 while(nn < n + 1) { xcross <- as.matrix(x[-nn,]) ycross <- as.matrix(y[-nn,]) # Estimate prior parameters yb (1 x q) and c (q x q) prior <- faprior(ycross) # Estimate likelihood parameters g (p x q), yzerogt (n-1 x p) and # s (p x q) lik <- falikelihood(xcross, ycross) # Estimate posterior parameters ya (1 x q) and d (q x q) out <- faprediction(x[nn,], ycross, prior$yb, prior$c, lik$g, lik$yzerogt, lik$s) ya[nn,]<-out$ya d[nn,]<-sqrt(diag(out$d)) nn <- nn + 1 } # end of cross-validation loop } else{ # Cross-validation loop nn <- 1 while(nn < n + 1) { xcross <- as.matrix(x[-nn,]) ycross <- as.matrix(y[-nn,]) mca <- svd(t(ycross) %*% xcross) xtrunc <- xcross %*% mca$v[,1:nmodes] ytrunc <- ycross %*% mca$u[,1:nmodes] # Estimate prior parameters yb (1 x q) and c (q x q) prior <- faprior(ytrunc) # Estimate likelihood parameters g (p x q), yzerogt (n-1 x p) and # s (p x q) lik <- falikelihood(xtrunc, ytrunc) # Estimate posterior parameters ya (1 x q) and d (q x q) out <- faprediction(x[nn,] %*% mca$v[,1:nmodes], ytrunc, prior$yb, prior$c, lik$g, lik$yzerogt, lik$s) ya[nn,] <- out$ya %*% t(mca$u[,1:nmodes]) d[nn,] <- t(sqrt(diag(mca$u[,1:nmodes] %*% out$d %*% t(mca$u[,1:nmodes])))) nn <- nn + 1 } # end of cross-validation loop } list(ya=ya,d=d) } falikelihood <- function(x, y){ # This function estimates the parameters G, yoG' and S # of the likelihood distribution x|y ~ N([y-yo]G',S), # where the symbol ' denotes transpose matrix. # # NOTE: in the code all parameters appear # in lower case letters. # # USAGE: falikelihood(x,y) # # INPUTS: # x -> n x p data matrix containing time series of p forecast variables # y -> n x q data matrix containing time series of q observable var. # # OUTPUTS: # g -> p x q matrix containing the calibration operator # yzerogt -> 1 x q matrix containing the intercept of the # regression between x and y # s -> p x p forecast error covariance matrix # # Authors: Caio Coelho # David Stephenson # # # find out p, q and n dimentions # p <- ncol(x) q <- ncol(y) n <- nrow(y) # # Define g, yzerogt and s # g <- matrix(nrow = p, ncol = q) yzerogt <- matrix(nrow = n, ncol = p) s <- matrix(nrow = p, ncol = p) # # Calculate column means of x and y # # xmean is a 1 x p matrix # xmean <- t(apply(x, 2, mean)) # # ymean is a 1 x q matrix # ymean <- t(apply(y, 2, mean)) # # Estimate covariance matrices # Sxx <- var(x)*((n-1)/n) Syy <- var(y)*((n-1)/n) Sxy <- var(x, y)*((n-1)/n) Syx <- t(Sxy) # # Inverte Syy # Syyi <- solve(Syy) # # Estimate forecast error covariance (s) # s <- Sxx - Sxy %*% Syyi %*% Syx # # Estimate intercept (yzerogt) and slope (g) from # the regression between x and y # # g is a p x q matrix # g <- Sxy %*% Syyi # # yzerogt is a 1 x q matrix. # yzerogt <- - (xmean - ymean %*% t(g)) list(g = g, yzerogt = yzerogt, s = s) } faprediction<-function(x, y, yb, c, g, yzerogt, s){ # This function estimates the mean ya and the covariance d # of the posterior distribution y|x ~ N(ya,d), # where ya = yb + [x - yb g' + yo g'] L' # d = (I - Lg) c # L' = (gcg' +s)^(-1) gc # # The symbol ' indicates transpose matrix # # USAGE: faprediction(x, y, yb, c, g, yzerogt, s) # # INPUTS: # x -> n x p data matrix containing time series of p forecast variables # y -> n x q data matrix containing time series of q observable vars. # yb -> 1 x q matrix containing column means of y # c -> q x q observable state covariance matrix # g -> p x q matrix containing the calibration operator # yzerogt -> 1 x q matrix containing the intercept of the # regression between x and y # s -> p x q forecast error covariance matrix # # OUTPUTS: # # ya -> 1 x q matrix containing q predicted observable state time series # d -> q x q matrix containing predicted observable error covaraiance # # Authors: Caio Coelho # David Stephenson # # Find out p, q and n # p <- ncol(x) q <- ncol(y) n <- nrow(y) # # Define matrices ya and d # ya <- matrix(nrow = n, ncol = q) d <- matrix(nrow = q, ncol = q) # # Transpose g # gt <- t(g) # # Calculate gain/weight matrix transpose (Lt) # Lt <- solve(g %*% c %*% gt + s) %*% g %*% c # # Calculate posterior mean (ya) # ya <- yb + (x - yb %*% gt + yzerogt) %*% Lt # # Calculate the posterior covariance matrix (d) # d <- (diag(q) - t(Lt) %*% g) %*% c # list(ya = ya, d = d) } faprior<-function(y){ # This function estimates the mean (yb) and the background # observable error covariance (c) of the prior distribution # of y given by y~N(yb,c) # # USAGE: faprior(y) # # INPUT: # y -> n x q data matrix containing q observable time series # # OUTPUTS: # yb -> 1 x q matrix containing the column means of y # c -> q x q observable state covariance matrix # # Authors: Caio Coelho # David Stephenson # # Calculate column means of y and store in yb (1 x q) yb <- t(apply(y, 2, mean)) # # # Estimates observable state covariance matrix (c) - biased estimate 1/N # c <- var(y)*(nrow(y)-1)/nrow(y) list(yb = yb, c = c) } ############################################################################## # Load libraries ############################################################################## # load library for reading netCDF data library(RNetCDF) # load library for performing extreme value analysis library(evd) library(ismev) # load libraries for ploting maps library(maps) library(mapdata) library(mapproj) flipmap <- function(lon,lat,data,dimension="lat") { # Flip a two-dimensional map (matrix) or a three-dimensional array up side down, # i.e. first line becomes last, last line becomes first and so on. In other # words, flip either latitude or longitude or both dimensions of a map or a # three-dimensional array (field), with first two space dimensions p for # longitude and q for latitude, and third n time dimension. # # Description: # # Returns a flipped map (matrix) or three-dimensional array # # Usage: # # flipmap(lon,lat,data,dimension) # # Arguments: # # lon: vector of longitudes # # lat: vector of latitudes # # data: A two-dimensional map (matrix) or three-dimensional array of data # values. If data is a matrix, it must # have the structure of longitude for columns and latitude for # rows. If data is a three-dimensional array, it must have the # longitude as the first dimension, latitude as the second # dimension and time as the third dimension # # dimension: A string indicating which dimension(s) will be flipped. # It must be either "lat" (default), "lon" or "both" # # # Output: # # $out: flipped map (matrix) or three-dimensional array # # $lon: longitude vector (flipped appropriately if necessary) # # $lat: latitude vector (flipped appropriately if necessary) # # Details: # # A matrix or array is flipped into a format suitable for plotting. # # Author: # # Dag Johan Steinskog 25 October 2005 # Caio Coelho # Chris Ferro # # Examples: # # #Flip a map (matrix) # x <- array(1:8,dim=c(2,4)) # lon <- seq(50,110,20) # lat <- seq(-10,10,20) # flipmap(lon,lat,x)$out # Defaults flips lat # flipmap(lon,lat,x)$lat # flipmap(lon,lat,x)$lon # flipmap(lon,lat,x,"lon")$out # flipmap(lon,lat,x,"both")$out # # #Flip a three-dimensional array # y <- array(1:8,dim=c(2,2,2)) # lonnew <- seq(50,70,20) # latnew <- seq(-10,10,20) # flipmap(lonnew,latnew,y)$out # Defaults flips lat # flipmap(lonnew,latnew,y)$lat # flipmap(lonnew,latnew,y)$lon # flipmap(lonnew,latnew,y,"lon")$out # flipmap(lonnew,latnew,y,"both")$out index<-NULL if(length(dim(data))==2){ if(dimension == "lon") {flip <- c(F,T)} if(dimension == "lat") {flip <- c(T,F)} if(dimension == "both") {flip <- c(T,T)} for (i in 1:length(dim(data))){ index[[i]] <- 1:dim(data)[i] if (flip[i]) {index[[i]] <- rev(index[[i]])} } out<-data[index[[1]],index[[2]]] lon<-lon[index[[2]]] lat<-lat[index[[1]]] } if(length(dim(data))==3){ if(dimension == "lon") {flip <- c(F,T,F)} if(dimension == "lat") {flip <- c(T,F,F)} if(dimension == "both") {flip <- c(T,T,F)} for (i in 1:length(dim(data))){ index[[i]] <- 1:dim(data)[i] if (flip[i]) {index[[i]] <- rev(index[[i]])} } out<-data[index[[1]],index[[2]],index[[3]]] lon<-lon[index[[2]]] lat<-lat[index[[1]]] } invisible(list(lon=lon,lat=lat,out=out)) } identifyfield <- function(x, y, data, fun=NULL, ..., mark = "*") { # Identify points by clicking on an image plot of the first time slice of # a given three-dimensional data array (i.e. a field) with first two space # dimensions (first dimension longitude, and second dimension latitude) and # third time dimension, or an image plot of a two-dimensional # data matrix (i.e. a map), with first dimension longitude, and second dimension # latitude. Also apply a user specified function to each of # the selected points of the three-dimensional data array. # If a two-dimensional data matrix (i.e. a map) is provided, there is no need # to specify a function to apply at each selected point, and # the function returns the data values for the points selected by clicking on # the image plot. # # Description: # # Applies a function to data associated with cells selected by # clicking on an image plot if data is a three-dimentional array. # If data is a matrix (i.e. a map), the function returns the data # values for the points selected by clicking on the image plot. # Also returns the location of the selected points. # # Usage: # # identifyfield(x, y, data, fun, ..., mark = "*") # # Arguments: # # x: Longitude vector to plot the image. # # y: Latitude vector to plot the image. # # data: A three-dimensional array (i.e. a field) or matrix (i.e. a map) # containing data associated with each cell of the image. The first # two dimensions must be `length(x)' and `length(y)' respectively. # # fun: The function to be applied to the three-dimensional array. Its # first argument must be the data vector associated with a cell. # # ...: Additional arguments passed to `fun'. # # mark: A vector of numbers or strings with which selected cells are # labelled. If shorter than the number of grid cells, values are # recycled; if `NULL', selections are not marked. # # Details: # # Clicking on an image plot created with grid lines at locations `x' # and `y' causes `fun' to be called repeatedly with each of the # vectors of `data' corresponding to selected grid cells. Returned # values from `fun' are stored and returned invisibly. # # Value: # # A list with an element for each grid-cell selection is returned # invisibly. Each element is a list with two or three components: # # loc: The (longitude,latiture)-coordinates of the selected grid cell. # # pos: The positions of the (x,y)-coordinates in `x' and `y'. # # data: The vector of `data' associated with that cell. # # value: The object (if any) returned by `fun' for the three-dimensional # array data ata each selected point. # # See Also: # # `apply', `identify', `image' # # Author: # # Chris Ferro 27 October 2005 # Caio Coelho # # Examples: # # # Three-dimensional array example (field) # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 100) # data <- array(rnorm(prod(dim)), dim) # fun <- function(x, ...) { x11(); plot(x, ...); mean(x) } # plot time # series and compute mean value of this time series # out <- identifyfield(x, y, data, fun, type = "l", mark = 1:prod(dim[1:2])) # fun <- function(x, index, ...) { x11(); plot(index, x, ...) } # out <- identifyfield(x, y, data, fun, index = rexp(100)) # # # Matrix example (map) # out <- identifyfield(x, y, data[,,1]) if (length(dim(data))==3) { fun <- match.fun(fun) image(x,y,data[,,1]) if (min(x)<0){map("world",add=TRUE)} else{map("world2",add=TRUE)} cat("Left-click to select grid cells; right-click to finish.\n") grid <- expand.grid(y = y, x = x)[, 2:1] if(is.null(mark)) lab <- 1:nrow(grid) else lab <- rep(mark, length.out = nrow(grid)) pts <- identify(grid, labels = lab, plot = !is.null(mark), offset = 0) series <- NULL loc <- grid[pts, ] ix <- match(loc[, 1], x) iy <- match(loc[, 2], y) #fun <- match.fun(fun) for(i in 1:length(pts)) { series[[i]] <- list(loc = loc[i, ]) series[[i]]$pos <- c(x = ix[i], y = iy[i]) series[[i]]$data <- data[ix[i], iy[i], ] series[[i]]$value <- fun(series[[i]]$data, ...) } } if (length(dim(data))==2) { image(x,y,data) if (min(x)<0){map("world",add=TRUE)} else{map("world2",add=TRUE)} cat("Left-click to select grid cells; right-click to finish.\n") grid <- expand.grid(y = y, x = x)[, 2:1] if(is.null(mark)) lab <- 1:nrow(grid) else lab <- rep(mark, length.out = nrow(grid)) pts <- identify(grid, labels = lab, plot = !is.null(mark), offset = 0) series <- NULL loc <- grid[pts, ] ix <- match(loc[, 1], x) iy <- match(loc[, 2], y) for(i in 1:length(pts)) { series[[i]] <- list(loc = loc[i, ]) series[[i]]$pos <- c(x = ix[i], y = iy[i]) series[[i]]$data <- data[ix[i], iy[i]] } } invisible(series) } moviefield<-function(lon,lat,array3d,start=1,end=10,delay=0.2, n=11,...){ # Produce a movie of a given three-dimensional array with p longitude points # and q latitude points as the first two dimensions and n as the third time # dimension # # Description: # # Animates a sequence of maps # # Usage: # # moviefield(lon,lat,array3d,start=1,end=10,delay=0.2,n=11,...) # # Input: # # array3d: a three-dimensional array with p longitude points and q latitude # points as the first two dimensions and n as the third time # dimension # # lon: vector of longitudes # # lat: vector of latitudes # # start: first time step of the animation # # end: last time step of the animation # # delay: number of waiting seconds between maps # # n: number of colours/intervals. # # ...: Additional arguments passed to `image' # # Authors: # # David Stephenson 16 June 2005 # Christopher Ferro # Caio Coelho # Dag Johan Steinskog # # Examples: # # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 1000) # data <- array(rnorm(prod(dim)), dim) # moviefield(x,y,data) # moviefield(x,y,data,start=1,end=50,delay=0.5,zlim=c(0.5,0.8)) # rg<-range(array3d, na.rm=T) lowerlim <- rg[1] upperlim <- rg[2] maximum <- max(abs(c(lowerlim,upperlim))) if(lowerlim<0) breaks <- seq(-maximum,maximum,(maximum-(-maximum))/n) else breaks <- seq(lowerlim,upperlim,(upperlim-(lowerlim))/n) wait <- function(time) { now <- start <- proc.time()[3] while(now < (start + time)) now <- proc.time()[3] } # check if lon is from 0 to 360 or -180 to 180 to use appropriate world map if (min(lon)<0){ coast <- map('world',interior = FALSE,resolution = 1) # Low resolution world map (lon -180 to 180) } else{ coast <- map('world2',interior = FALSE,resolution = 1) # Low resolution world map (lon 0 360) } # check if range includes negative values to use appropriate colour scale if (min(array3d,na.rm=T) <0) { colours = bluered(seq(0,n-1,by=1), "linear", yellow =TRUE) } else { #colours = rev(heat.colors(n+1)) colours = rev(heat.colors(n)) } par(mai=c(2, 1, 2, 0.5)) #par(mfrow = c(1, 1), mar = c(1, 1, 1, 1), pty = "m") for(i in start:end) { image(lon, lat, array3d[,,i], col = colours, xlab='', ylab='', axes = FALSE,breaks=breaks,...) map(coast, add = TRUE) title(paste("Frame",i,"of",end,sep=" ")) box() wait(delay) } } mygpd.fit<-function (xdat, threshold, npy = 365, ydat = NULL, sigl = NULL, shl = NULL, siglink = identity, shlink = identity, show = TRUE, method = "Nelder-Mead", maxit = 10000, ...) { # Same as gpd.fit function from ismev package but with standard error # calculation disactivated z <- list() npsc <- length(sigl) + 1 npsh <- length(shl) + 1 n <- length(xdat) z$trans <- FALSE if (is.function(threshold)) stop("`threshold' cannot be a function") u <- rep(threshold, length.out = n) if (length(unique(u)) > 1) z$trans <- TRUE xdatu <- xdat[xdat > u] xind <- (1:n)[xdat > u] u <- u[xind] in2 <- sqrt(6 * var(xdat))/pi in1 <- mean(xdat, na.rm = TRUE) - 0.57722 * in2 if (is.null(sigl)) { sigmat <- as.matrix(rep(1, length(xdatu))) siginit <- in2 } else { z$trans <- TRUE sigmat <- cbind(rep(1, length(xdatu)), ydat[xind, sigl]) siginit <- c(in2, rep(0, length(sigl))) } if (is.null(shl)) { shmat <- as.matrix(rep(1, length(xdatu))) shinit <- 0.1 } else { z$trans <- TRUE shmat <- cbind(rep(1, length(xdatu)), ydat[xind, shl]) shinit <- c(0.1, rep(0, length(shl))) } init <- c(siginit, shinit) z$model <- list(sigl, shl) z$link <- deparse(substitute(c(siglink, shlink))) z$threshold <- threshold z$nexc <- length(xdatu) z$data <- xdatu gpd.lik <- function(a) { sc <- siglink(sigmat %*% (a[seq(1, length = npsc)])) xi <- shlink(shmat %*% (a[seq(npsc + 1, length = npsh)])) y <- (xdatu - u)/sc y <- 1 + xi * y if (min(sc) <= 0) l <- 10^6 else { if (min(y) <= 0) l <- 10^6 else { l <- sum(log(sc)) + sum(log(y) * (1/xi + 1)) } } l } x <- optim(init, gpd.lik, hessian = TRUE, method = method, control = list(maxit = maxit, ...)) sc <- siglink(sigmat %*% (x$par[seq(1, length = npsc)])) xi <- shlink(shmat %*% (x$par[seq(npsc + 1, length = npsh)])) z$conv <- x$convergence z$nllh <- x$value z$vals <- cbind(sc, xi, u) if (z$trans) { z$data <- -log(as.vector((1 + (xi * (xdatu - u))/sc)^(-1/xi))) } z$mle <- x$par z$rate <- length(xdatu)/n #z$cov <- solve(x$hessian) #z$se <- sqrt(diag(z$cov)) z$n <- n z$npy <- npy z$xdata <- xdat if (show) { if (z$trans) print(z[c(2, 3)]) if (length(z[[4]]) == 1) print(z[4]) print(z[c(5, 7)]) if (!z$conv) print(z[c(8, 10, 11, 13)]) } invisible(z) } netcdfinfo <- function(file) { # Extract data structure, names of variables and dimension information # form netcdf file # # Description: # # Returns a summary of the data structure, a list of dimensions, # a list of variables and the lengths (size) of each dimension. # # Usage: # # netcdfinfo(file) # # Arguments: # # file: String containing the full path and netcdf file name. # # Output: # # $dims: dimension names # $vars: Variable names # $lenghts: dimension lenghts (size) # # # Author: # # Chris Ferro 9 June 2005 # Dag Johan Steinskog # Caio Coelho # # Examples: # # netcdfinfo("netcdf_file_name.nc") # netcdfinfo("/home/username/netcdf_file_name.nc") # # Note: The user should use a netcdf file to be able to run these examples # open a netcdf file ncfile <- open.nc(file) # print details print.nc(ncfile) # store dimension names dims <- character(file.inq.nc(ncfile)$ndims) for(i in 1:length(dims)) dims[i] <- dim.inq.nc(ncfile, i - 1)$name # store variables names vars <- character(file.inq.nc(ncfile)$nvars) for(i in 1:length(vars)) vars[i] <- var.inq.nc(ncfile, i - 1)$name # store lenghts of each dimension lenghts <- as.vector(seq(1:file.inq.nc(ncfile)$ndims)) for(i in 1:length(dims)) lenghts[i] <- dim.inq.nc(ncfile, i - 1)$length list(dims=dims,vars=vars,lenghts=lenghts) } netcdfread<-function(file,lonname='longitude',latname='latitude',timename='time',dataname,unpack=F) { # # Extract longitude vector, latitude vector, time vector and # data array from a netcdf file. # # Description: # # Returns a vector of longitudes, a vector of latitudes, a vector of times # and an array of data from a netcdf file. # # Usage: # # netcdfread(file,lonname,latname,timename,dataname,unpack=F) # # Arguments: # # file: String containing the full path and netcdf file name # # lonname: String containing the name of the longitude variable. # Default name is "longitude". # # latname: String containing the name of the latitude variable # Default name is "latitude". # # timename: String containing the name of the time variable. # Default name is "time" # # dataname: String containing the name of the data variable # # unpack: Logical. If TRUE unpacks data contained in variable # dataname by performing the following transformation: # unpacked_data_value=packed_data_value*scale_factor+add_offset, # where, scale_factor and add_offset are constant netcdf attributes # of variable dataname. Default is FALSE. # # Output: # # $lonncfile: Vector of longitudes # $latncfile: Vector of latitudes # $timencfile: Vector of times # $datancfile: Array of data # # Author: # # Chris Ferro 9 June 2005 # Dag Johan Steinskog # Caio Coelho # # Examples: # # # First run the line below to extract the data structure/information # # so that you are able to know the names of variables # netcdfinfo("netcdf_file_name.nc") # # netcdfread("netcdf_file_name.nc","lon_name","lat_name","time_name","data_name") # netcdfread("/home/username/netcdf_file_name.nc","lon_name","lat_name","time_name","data_name") # # Note: The user should use a netcdf file to be able to run these examples # open a netcdf file ncfile <- open.nc(file) # get latitude and longitude lonncfileaux <- var.get.nc(ncfile, lonname) latncfileaux <- var.get.nc(ncfile, latname) # get data and time datancfile <- var.get.nc(ncfile, dataname) timencfile <- var.get.nc(ncfile, timename) if(unpack){ datancfile<-datancfile*att.get.nc(ncfile,dataname,"scale_factor")+att.get.nc(ncfile,dataname,"add_offset") } # sort latitude and longitude in ascending order lonncfile <- sort(lonncfileaux) latncfile <- sort(latncfileaux) # rearange data datancfile<-datancfile[order(lonncfileaux),order(latncfileaux),] invisible(list(lonncfile=lonncfile,latncfile=latncfile,timencfile=timencfile,datancfile=datancfile)) } netcdfwrite <- function(lon,lat,data,filename="data.nc",time=1,mv=-999) { # Writes two-dimensional map or three-dimensional array of data, longitude vector, # latitude vector and time vector into a netcdf file. # # Description: # # Returns a netcdf file that contains a two-dimensional map or a # three-dimensional array of data, a longitude vector, a latitude vector # and a time vector # # Usage: # # netcdfwrite(lon,lat,data,filename,time,mv) # # Arguments: # # data: two-dimensional map with first dimension longitude and # second dimension latitude or three-dimensional array # with first dimension longitude, second dimension latitude and # third dimension time, containing the data to be written # to the netcdf file # lon: a vector containing the longitudes of data # lat: a vector containing the latitudes of data # time: a vector containing the times of data in hours # since 1900-1-1 00:00:0.0 . Default is 1, which is # appropriate for writing a two-dimensional map # filename: string with the name of the netcdf file to be created. # Default is "data.nc" # mv: Value attributed to missing data. Default is -999 # # Output: # # netcdf-file containing the data # # Author: # # Dag Johan Steinskog 5 January 2006 # Caio Coelho # # Example: # # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dims <- c(length(x), length(y), 100) # data <- array(rnorm(prod(dims)), dims) # time <- 1:100 # # write a two-dimensional map of data to the netcdf file # netcdfwrite(x,y,data[,,3],filename="data.nc") # # write a three-dimensional array of data to the netcdf file # netcdfwrite(x,y,data,filename="data.nc",time) if(length(time)==1){data<-array(data,c(dim(data)[1],dim(data)[2],1))} #create netcdf file #ncfile <- create.nc(filename,clobber) ncfile <- create.nc(filename) #give dimensions to variables dim.def.nc(ncfile,"lon",length(lon)) dim.def.nc(ncfile,"lat",length(lat)) dim.def.nc(ncfile,"time",length(time)) # create variables var.def.nc(ncfile,"lon","NC_FLOAT","lon") var.def.nc(ncfile,"lat","NC_FLOAT","lat") var.def.nc(ncfile,"time","NC_DOUBLE","time") var.def.nc(ncfile,"data","NC_FLOAT",c("lon","lat","time")) #var.def.nc(ncfile,"data","NC_DOUBLE",c("lon","lat")) # add data to variables var.put.nc(ncfile,"lat",lat) var.put.nc(ncfile,"lon",lon) var.put.nc(ncfile,"time",time) var.put.nc(ncfile,"data",data) # add description information to variables att.put.nc(ncfile,"data","missing_value","NC_FLOAT",mv) att.put.nc(ncfile,"lon","long_name","NC_FLOAT","Longitude") att.put.nc(ncfile,"lon","units","NC_FLOAT","degrees_east") att.put.nc(ncfile,"lon","axis","NC_FLOAT","X") att.put.nc(ncfile,"lat","long_name","NC_FLOAT","Latitude") att.put.nc(ncfile,"lat","units","NC_FLOAT","degrees_north") att.put.nc(ncfile,"lat","axis","NC_FLOAT","Y") att.put.nc(ncfile,"time","long_name","NC_DOUBLE","Time") att.put.nc(ncfile,"time","units","NC_DOUBLE","hours since 1900-1-1 00:00:0.0") att.put.nc(ncfile,"time","delta_t","NC_DOUBLE","0000-01-00 00:00:00") # syncronize and close the netcdf file sync.nc(ncfile) close.nc(ncfile) } plotmap <- function(lon,lat,data,maintit='',legtit='',equi=TRUE,bw=FALSE, cont=FALSE,reg=FALSE,..., lonlim=c(0,360), latlim=c(50,90),orientation=NULL, mapdat="world",xmaplim=c(-180,180),ymaplim=c(50,90), longrds=NULL,latgrds=NULL, breaks=NULL,n=11,colours=NULL) { # Map two dimensional data matrix, with first dimension p (longitude points) # and second dimension q (latitude points), in either cylindrical equidistant # latitude and longitude projection or stereographic projection # # Description: # # Plots data in either cylindrical equidistant latitude and longitude # projection or stereographic projection with colorbar and contours # if requested # # Usage: # # plotmap(lon,lat,data,maintit='',legtit='',equi=TRUE,bw=FALSE, # cont=FALSE,reg=FALSE,..., lonlim=c(0,360), # latlim=c(50,90),orientation=NULL, # mapdat="world",xmaplim=c(-180,180),ymaplim=c(50,90), # longrds=NULL,latgrds=NULL,breaks=NULL,n=11,colours=NULL) # # Arguments: # # lon: Vector of p longitude coordinates in ascending order. # # lat: Vector of q latitude coordinates in ascending order. # # data: Matrix of p rows (longitudes) and q columns (latitudes) that # contains the data values to be plotted. # #maintit: String with the main title of the plot. If not provided # no title is displayed. # # legtit: String with the title of the colourbar. If not provided # no title is displayed. # # cont: Logical. If TRUE adds contours to the image plot. Default is FALSE. # # breaks: Vector of values for the colour bar. Must always be of # length(n)+1 (see definition of n below). Default is NULL, meaning # that the function produces an automatic vector of breaks based on # n and the range of values of the map to be plotted. # # n: number of colours/intervals for the colour bar. Default is 11. # For the default plot produced by this function n must be an odd number. # #colours: Vector of colours of length n. Defauls is NULL, meaning that an # automatically chosen colour (or gray) scale is produced. # # bw: Logical. If TRUE produces black, gray and white plot. Default is FALSE. # # equi: Logical. If FALSE produces stereographic projection plot. Default is TRUE. # # reg: Logical. To be used only if plotting data in equidistant projection. # If TRUE produces a plot using as much ares as possible of the display # window. Default is FALSE (addequate to plot global data). # # ...: Additional arguments passed to `contour'. # # latlim: Range of latitudes where data will be plotted. e.g. latlim=c(50,90) # or latlim=c(-90,-50) # # lonlim: Range of longitudes where data will be plotted. e.g. lonlim=c(0,360) # or lonlim=c(-180,180) # # orientation: Orientation parameter of mapproject. For Europe # as a whole a reasonable setting is: orientation=c(45,0,7.5) # The latter places the pole of projection to latitude 45 N. # For the South Pole use, for example, orientation=c(-90,0,0) # # mapdat: Name of the dataset from which map data is taken from either maps or mpdata # packages for producing stereographic plots. Default is "world". Other options # could be for example mapdat = "worldHires" plots high resolution map, # mapdat = "none" no map, . # #xmaplim and ymaplim: Range of longitudes/latitudes where the map data will be drawn. # In default, the map is drawn for the area of longitude and latitude provided # in vectors lon and lat. Sometimes it is useful to define e.g. xmaplim=c(-180,180), # ymaplim=c(0,90) i.e. whole hemisphere. But this will be very slow if invoked # with madat="worldHires". # # Authors: # # Caio Coelho 2 November 2005 # Dag Johan Steinskog # # Examples: # # #Cylindrical equidistant global plots with UK on the right hand side # plotmap(seq(0,360,2.5),seq(-90,90,2.5),matrix(rnorm(145*73,sd=10),ncol=73,nrow=145),'Example','Random data') # plotmap(seq(0,360,2.5),seq(-90,90,2.5),matrix(rnorm(145*73,sd=10),ncol=73,nrow=145),'Example','Random data',bw=T) # plotmap(seq(0,360,2.5),seq(-90,90,2.5),abs(matrix(rnorm(145*73,sd=10),ncol=73,nrow=145)),'Example','Random data') # # #Cylindrical equidistant plots with UK in the centre # plotmap(seq(-180,180,2.5),seq(-90,90,2.5),matrix(rnorm(145*73,sd=10),ncol=73,nrow=145),'Example','Random data') # plotmap(seq(-180,180,2.5),seq(-90,90,2.5),matrix(rnorm(145*73,sd=10),ncol=73,nrow=145),'Example','Random data',bw=T) # # #Cylindrical equidistant regional plots # plotmap(seq(-50,50,1),seq(-50,50,1),matrix(rnorm(101*101,sd=10),ncol=101,nrow=101),'Example','Random data',reg=T) # plotmap(seq(-50,50,1),seq(-50,50,1),matrix(rnorm(101*101,sd=10),ncol=101,nrow=101),'Example','Random data',reg=T,bw=T) # plotmap(seq(0,100,1),seq(-50,50,1),matrix(rnorm(101*101,sd=10),ncol=101,nrow=101),'Example','Random data',reg=T) # # #Stereographic plot for the North Pole # plotmap(seq(-180,180,2.5),seq(60,90,2.5),matrix(rnorm(145*13,sd=10),ncol=13,nrow=145),'Example','Random data',equi=F,lonlim=c(-180,180),latlim=c(60,90),xmaplim=c(-180,180),ymaplim=c(60,90)) # plotmap(seq(-180,180,2.5),seq(60,90,2.5),matrix(rnorm(145*13,sd=10),ncol=13,nrow=145),'Example','Random data',equi=F,bw=T,lonlim=c(-180,180),latlim=c(60,90),xmaplim=c(-180,180),ymaplim=c(60,90)) # # #Stereographic plot for the South Pole # plotmap(seq(-180,180,2.5),seq(-90,-60,2.5),matrix(rnorm(145*13,sd=10),ncol=13,nrow=145),'Example','Random data',equi=F,lonlim=c(-180,180),latlim=c(-90,-60),orientation=c(-90,0,0),xmaplim=c(-180,180),ymaplim=c(-90,-60)) # plotmap(seq(-180,180,2.5),seq(-90,-60,2.5),matrix(rnorm(145*13,sd=10),ncol=13,nrow=145),'Example','Random data',equi=F,bw=T,lonlim=c(-180,180),latlim=c(-90,-60),orientation=c(-90,0,0),xmaplim=c(-180,180),ymaplim=c(-90,-60)) rg<-range(data, na.rm=T) lowerlim <- rg[1] upperlim <- rg[2] maximum <- max(abs(c(lowerlim,upperlim))) if(is.null(breaks)){ if(lowerlim<0) breaks <- seq(-maximum,maximum,(maximum-(-maximum))/n) else breaks <- seq(lowerlim,upperlim,(upperlim-(lowerlim))/n) } if(is.null(colours)){ if(!bw){ # check if range includes negative values to use appropriate colour scale if (rg[1] <0) { colours = bluered(seq(0,n-1,by=1), "linear", yellow =TRUE) } else { colours = bluered(seq(0,n-1,by=1), "linear", white=0, invert=T) } } else{ if (rg[1] <0) { colours <- grey(seq(0, 1, length = length(breaks)-1)) colours <- c(colours[1:((n-1)/2)],rev(colours[(((n-1)/2)+1):n])) #if ((n/2)==floor(n/2)) {colours <- c(colours[1:((n-1)/2)],rev(colours[(((n-1)/2)+1):n-1]))} } else { colours <- grey(seq(0, 1, length = length(breaks)-1)) } } } # end if is.null(colours) if (equi){ #par(mfrow = c(2, 1),mai=c(0.1, 0.5, 1, 0.5)) layout(matrix(1:2, ncol = 1, nrow=2), heights = c(5,1)) #par(mar = c(6, 2, 6.5, 0.5)) par(mar = c(6, 2, 6.5, 1)) if (reg) { layout(matrix(1:2, ncol = 1, nrow=2), heights = c(9,1)) #par(mar = c(5, 2, 1.5, 0.5)) par(mar = c(5, 2, 1.5, 1)) } # plot data image(lon,lat,data,axes=F, col = colours, xlab = '', ylab = '',breaks=breaks) if(cont) {contour(lon,lat,data,levels=round(breaks,1),labels=round(breaks,1),add=T,...)} # check if lon is from 0 to 360 or -180 to 180 to use appropriate world map if (min(lon)<0){ map('world',interior = FALSE,add = T, lwd=1) # Low resolution world map (lon -180 to 180) } else{ map('world2',interior = FALSE,add = T, lwd=1) # Low resolution world map (lon 0 360) } box() map.axes() title(maintit) # adding colorbar #par(mai = c(2.8, 0.5, 0.5, 0.5)) #par(mar = c(4.5, 2, 0, 0.5), mgp = c(1.5, 0.3, 0), las = 1) par(mar = c(4.5, 2, 0, 1), mgp = c(1.5, 0.3, 0), las = 1) if (reg) { #par(mar = c(2.5, 2, 0, 0.5), mgp = c(1.5, 0.3, 0), las = 1) par(mar = c(2.5, 2, 0, 1), mgp = c(1.5, 0.3, 0), las = 1) } image(c(1:n),1, t(t(c(1:n))), axes = F, col = colours,xlab = '', ylab = '') box() par(cex = 1.1) ## add the tick marks to the plot without plotting labels axis(1, at = seq(0.5, length(breaks)-0.5), labels = F) #axis(1, at = seq(1, length(breaks)), labels = F) if(maximum>1){ ## add the labels to the plot, without plotting tick marks #axis(1, at = seq(1.5, length(breaks)-1.5), tick = F, axis(1, at = seq(0.5, length(breaks)-0.5), tick = F, # labels = round(breaks[2:(length(breaks)-1)],1)) labels = round(breaks[1:(length(breaks))],1)) } else{ #axis(1, at = seq(1.5, length(breaks)-1.5), tick = F, axis(1, at = seq(0.5, length(breaks)-0.5), tick = F, # labels = round(breaks[2:(length(breaks)-1)],2)) labels = round(breaks[1:(length(breaks))],2)) } ## redifine font size par(cex = 1.2) ## add the title title(xlab = legtit) } #end if equi else{ layout(matrix(1:2, ncol = 1, nrow=2), heights = c(9,1)) #par(mar = c(5, 2, 1.5, 0.5)) par(mar = c(5, 2, 1.5, 1)) #image.map(lon,lat,data,latlim=c(60,90),projection="stereographic",xmaplim=c(-180,180),ymaplim=c(60,90),col=colours,breaks=breaks,longrds=seq(-180,180,by=60),latgrds=seq(60,80,by=10)) ##image.map(lon,lat,data,latlim=latlim,lonlim=lonlim,projection="stereographic",orientation=orientation,mapdat=mapdat,xmaplim=xmaplim,ymaplim=ymaplim,col=colours,breaks=breaks,longr="no",latgr="no") if(cont){projmap(lon,lat,data,latlim=latlim,lonlim=lonlim,projection="stereographic",orientation=orientation,mapdat=mapdat,xmaplim=xmaplim,ymaplim=ymaplim,col=colours,breaks=breaks,longr="no",latgr="no",con=T) } else {projmap(lon,lat,data,latlim=latlim,lonlim=lonlim,projection="stereographic",orientation=orientation,mapdat=mapdat,xmaplim=xmaplim,ymaplim=ymaplim,col=colours,breaks=breaks,longr="no",latgr="no") } title(maintit) # adding colorbar #par(mai = c(2.8, 0.5, 0.5, 0.5)) #par(mar = c(2.5, 2, 0, 0.5), mgp = c(1.5, 0.3, 0), las = 1) par(mar = c(2.5, 2, 0, 1), mgp = c(1.5, 0.3, 0), las = 1) image(c(1:n),1, t(t(c(1:n))), axes = F, col = colours,xlab = '', ylab = '') box() par(cex = 1.1) ## add the tick marks to the plot without plotting labels axis(1, at = seq(0.5, length(breaks)-0.5), labels = F) #axis(1, at = seq(1, length(breaks)), labels = F) if(maximum>1){ ## add the labels to the plot, without plotting tick marks #axis(1, at = seq(1.5, length(breaks)-1.5), tick = F, axis(1, at = seq(0.5, length(breaks)-0.5), tick = F, # labels = round(breaks[2:(length(breaks)-1)],1)) labels = round(breaks[1:(length(breaks))],1)) } else{ #axis(1, at = seq(1.5, length(breaks)-1.5), tick = F, axis(1, at = seq(0.5, length(breaks)-0.5), tick = F, # labels = round(breaks[2:(length(breaks)-1)],2)) labels = round(breaks[1:(length(breaks))],2)) } ## redifine font size par(cex = 1.2) ## add the title title(xlab = legtit) } # end else equi } plotpc <- function(x,i=1,n=11,bw=FALSE,cont=FALSE,...){ # Plot EOFs and Principal Components (PCs, time series) that are output of the # eof function # # Description: # # Plots the principal component time series and the map of loadings # of a given mode # # Usage: # # plotpc(x,i,n=11,bw=FALSE,cont=FALSE,...) # # Input: # # x: the output of the function eofs # # i: the number of the mode to be plotted # # # cont: Logical. If TRUE adds contours to the image plot. Default is FALSE. # # n: An odd number of colours/intervals for the colourbar. Default is 11. # # bw: Logical. If TRUE produces black, gray and white plot. Default is FALSE. # # ...: Additional arguments passed to `contour'. # # Output: # # # Authors: # # Dag Johan Steinskog 10 June 2005 # Caio Coelho # # Example: # # lon <- seq(0,360,5) # lat <- seq(-90,90,5) # t <- 365 # x <- array(rnorm(length(lon)*length(lat)*t,sd=10),dim=c(length(lon),length(lat),t)) # sss<-eof(lon,lat,x) # plotpc(sss,2,cont=T) # plots the second mode # plotpc(sss,10,bw=T,cont=T) # plots the tenth mode rg<-range(x$eof[,,i]*100, na.rm=T) lowerlim <- rg[1] upperlim <- rg[2] maximum <- max(abs(c(lowerlim,upperlim))) if(lowerlim<0) breaks <- seq(-maximum,maximum,(maximum-(-maximum))/n) else breaks <- seq(lowerlim,upperlim,(upperlim-(lowerlim))/n) if(!bw){ # check if range includes negative values to use appropriate colour scale if (rg[1] <0) { colours = bluered(seq(0,n-1,by=1), "linear", yellow =TRUE) } else { colours = bluered(seq(0,n-1,by=1), "linear", white=0, invert=T) } } else{ if (rg[1] <0) { colours <- grey(seq(0, 1, length = length(breaks)-1)) colours <- c(colours[1:((n-1)/2)],rev(colours[(((n-1)/2)+1):n])) } else { colours <- grey(seq(0, 1, length = length(breaks)-1)) } } layout(matrix(1:3, ncol = 1, nrow=3), heights = c(5,8,1)) par(mar = c(6, 4, 4, 0.5)) # Plot PC time series plot(x$pcs[,i],type='l',xlab='Time',ylab=paste('PC',as.character(i))) title(paste('Principal component',as.character(i),sep=" ")) # Plot loaddings par(mar = c(6, 4, 2, 0.5)) image(x$lon,x$lat,x$eof[,,i]*100,col=colours,xlab = '', ylab = '',breaks=breaks) if(cont) {contour(x$lon,x$lat,x$eof[,,i]*100,levels=round(breaks,1),labels=round(breaks,1),add=T,...)} title(paste('EOF',as.character(i),' ',as.character(x$accountedvar[i]),'%')) # check if lon is from 0 to 360 or -180 to 180 to use appropriate world map if (min(x$lon)<0){ map('world',interior = FALSE,add = T, lwd=1) # Low resolution world map (lon -180 to 180) } else{ map('world2',interior = FALSE,add = T, lwd=1) # Low resolution world map (lon 0 360) } box() map.axes() #title(maintit) # adding colorbar par(mar = c(2.5, 2, 0, 0.5), mgp = c(1.5, 0.3, 0), las = 1) image(c(1:n),1, t(t(c(1:n))), axes = F, col = colours,xlab = '', ylab = '') box() par(cex = 1.1) ## add the tick marks to the plot without plotting labels axis(1, at = seq(1.5, length(breaks)-1.5), labels = F) if(maximum>1){ ## add the labels to the plot, without plotting tick marks axis(1, at = seq(1.5, length(breaks)-1.5), tick = F, labels = round(breaks[2:(length(breaks)-1)],1)) } else{ axis(1, at = seq(1.5, length(breaks)-1.5), tick = F, labels = round(breaks[2:(length(breaks)-1)],2)) } ## redifine font size par(cex = 1.2) ## add the title title(xlab = 'loadings') } projmap <- function(x,y,z,rotpol=c(0.,90.), lonlim,latlim, projection="",parameters=NULL,orientation=NULL, mapdat="world",xmaplim,ymaplim,thin.map=0, col=heat.colors(12),breaks,na.col="white", longrds=seq(-180,180,by=20),longr="yes", latgrds=seq(-80,+80,by=20),latgr="yes", ltygrds=3,lwdgrds=1,con=FALSE,...) { # Generate a pixel plot of matrix z on a specified map projection. # This routine deals with sterepgraphic and lambert projections from the mapproj # package. # # Description: # # Returns a pixel plot of matrix z # # Usage: # # projmap(x,y,z,rotpol=c(0.,90.),lonlim,latlim,projection="",parameters=NULL,orientation=NULL, # mapdat="world",xmaplim,ymaplim,thin.map=0,col=heat.colors(12),breaks,na.col="white", # longrds=seq(-180,180,by=20),longr="no",latgrds=seq(-80,+80,by=20),latgr="no", # ltygrds=3,lwdgrds=1,con=T,...) # # Arguments: # # x: A vector of longitudes # # y: A vector of latitudes # # z: A matrix # # rotpol: defines the coordinates of the rotated pole. Default is # geographical north pole (rotpol=c(0.,90.)). # # lonlim: defines the longitude window to be plotted. If # specified the window of the projection is chosen such that the entire # longitude window is comprised in the plot. # # latlim: defines the latitude window to be plotted. If # specified the window of the projection is chosen such that the entire # longitude window is comprised in the plot. # # projection: either lambert or stereographic provided by mapproject(). If # an empty string is used, the projection used in the last application # of mapproject will be used. The last projection settings are stored # in variable .Last.projection # # orientation: is the orientation parameter of mapproject. # # mapdat: is the name of the dataset from which map data is taken. # Default is "world". # # xmaplim: defines a longitude window from which map data shall # be plotted. In default, all map data comprised in the area of x,y is # used for plotting. In some cases this does not comprise all the map # data of the projection window. # # ymaplim: define a latitude window from which map data shall # be plotted. In default, all map data comprised in the area of x,y is # used for plotting. # # thin.map: whether to thin out details of the map (political, continental # outlines). thin.map = 1 means that every second point on the map is # omitted. thin.map = 2 means that this thinning out is applied twice, etc. # By default thin.map = 0, which means no thinning. # # col: breaks the color palette and breaks to be used for the pixels. # length of breaks needs to be one larger than length of colors. If # breaks is not specified equidistant breaks from min to max are used. # # na.col: the color to use for pixels with NA (in image). # Useful values are "white" and "grey80". # # longrds: array of longitude circles for which grid-lines shall be plotted. # # longr: string indicating if the longitude circles should be plotted. Default is "no". # # latgrds: array of latitude circles for which grid-lines shall be plotted. # # longr: string indicating if the latitude circles should be plotted. Default is "no". # # ltygrds: the type of line for longitude and latitude circles. # # lwdgrds: the line thickness for longitude and latitude circles. # # Details: # # Plots a pixel plot with either the lambert or stereographical projection. # This script is used in the function plotmap.r. # # Value: # # A pixel plot of the chosen area and projection. # # See Also: # # `plotmap' # # Author: # # Dag Johan Steinskog 3 November 2005 # Christoph Frei # # Examples: # # lon<-seq(-180,180,2.5) # lat<-seq(60,90,2.5) # data<-matrix(rnorm(length(lon)*length(lat),sd=10),ncol=length(lat),nrow=length(lon)) # projmap(lon,lat,data,projection="stereographic", xmaplim=c(-180,180), ymaplim=c(60,90),longrds=seq(-180,180,by=60),mapdat="world",latgrds=seq(60,80,by=10)) # # lon<-seq(-20,40,2.5) # lat<-seq(50,70,2.5) # data<-matrix(rnorm(length(lon)*length(lat),sd=10),ncol=length(lat),nrow=length(lon)) # projmap(lon,lat,data,projection="lambert",lonlim=c(-20,40),latlim=c(50,70),xmaplim=c(-20,40),ymaplim=c(49,70),longr="no",latgr="no",param=c(-20,60)) # rotated coordinates or not? if (rotpol[2] != 90.) {rot <- TRUE} else {rot <- FALSE} # define generally useful variables nx <- length(x) ny <- length(y) plon <- rotpol[1] plat <- rotpol[2] dx <- diff(x)[1] dy <- diff(y)[1] rgb.miss <- col2rgb(na.col) col.miss <- rgb(rgb.miss[1],rgb.miss[2],rgb.miss[3], maxColorValue = 255) # convert positions/values into single structured arrays xmat <- matrix(rep(x,ny),nrow=nx) ymat <- t(matrix(rep(y,nx),nrow=ny)) xx <- array(xmat) yy <- array(ymat) zz <- array(z) rm(xmat,ymat) # define a lon/lat frame to define the window limits and the # details of the projection # case 1 : if lonlim and latlim are specified if (!missing(lonlim) & !missing(latlim)) { lon <- seq(lonlim[1],lonlim[2],length=nx) lat <- seq(latlim[1],latlim[2],length=ny) fram <- matrix(rep(0,2*2*(nx+ny)),ncol=2) fram[,1] <- c(lon,rep(lon[nx],ny),rev(lon),rep(lon[1],ny)) fram[,2] <- c(rep(lat[1],nx),lat,rep(lat[ny],nx),rev(lat)) } else { # case 2 : frame from field lon <- x; lat <- y fram <- matrix(rep(0,2*2*(nx+ny)),ncol=2) fram[,1] <- c(lon,rep(lon[nx],ny),rev(lon),rep(lon[1],ny)) fram[,2] <- c(rep(lat[1],nx),lat,rep(lat[ny],nx),rev(lat)) if (rot) { fram <- rot2lonlat(fram,pollam=plon,polphi=plat) } } # define the details of the projection and the window limits in # plotting coordinates. all subsequent calls # to mapproject will use the same setting. ll <- mapproject(x=fram[,1],y=fram[,2],projection=projection, parameters=parameters,orientation=orientation) xxlim <- range(ll$x,finite=TRUE); yylim <- range(ll$y,finite=TRUE) rm(ll,fram) # calculate transformed pixel corner positions # a) calculate corner points xxll <- xx-dx/2 xxlr <- xx+dx/2 xxur <- xx+dx/2 xxul <- xx-dx/2 yyll <- yy-dy/2 yylr <- yy-dy/2 yyur <- yy+dy/2 yyul <- yy+dy/2 # b) rotate coordinates if necessary if (rot) { # centers fram <- matrix(rep(0,2*nx*ny),ncol=2) fram[,1] <- xx; fram[,2] <- yy ll <- rot2lonlat(fram,pollam=plon,polphi=plat) xx <- ll[,1]; yy <- ll[,2] # lower left corner fram <- matrix(rep(0,2*nx*ny),ncol=2) fram[,1] <- xxll; fram[,2] <- yyll ll <- rot2lonlat(fram,pollam=plon,polphi=plat) xxll <- ll[,1]; yyll <- ll[,2] # lower right corner fram <- matrix(rep(0,2*nx*ny),ncol=2) fram[,1] <- xxlr; fram[,2] <- yylr ll <- rot2lonlat(fram,pollam=plon,polphi=plat) xxlr <- ll[,1]; yylr <- ll[,2] # upper right corner fram <- matrix(rep(0,2*nx*ny),ncol=2) fram[,1] <- xxur; fram[,2] <- yyur ll <- rot2lonlat(fram,pollam=plon,polphi=plat) xxur <- ll[,1]; yyur <- ll[,2] # upper left corner fram <- matrix(rep(0,2*nx*ny),ncol=2) fram[,1] <- xxul; fram[,2] <- yyul ll <- rot2lonlat(fram,pollam=plon,polphi=plat) xxul <- ll[,1]; yyul <- ll[,2] # clean up rm(fram,ll) } # c) convert to the desired projection # corners xxx <- c(xxll,xxlr,xxur,xxul) yyy <- c(yyll,yylr,yyur,yyul) nnn <- length(xxll) ll <- mapproject(x=xxx,y=yyy) xxll <- ll$x[1:nnn]; yyll <- ll$y[1:nnn] xxlr <- ll$x[(nnn+1):(2*nnn)]; yylr <- ll$y[(nnn+1):(2*nnn)] xxur <- ll$x[(2*nnn+1):(3*nnn)]; yyur <- ll$y[(2*nnn+1):(3*nnn)] xxul <- ll$x[(3*nnn+1):(4*nnn)]; yyul <- ll$y[(3*nnn+1):(4*nnn)] rm(ll,xxx,yyy,nnn) # determine the map window in lon/lat if (!missing(xmaplim) & !missing(ymaplim)) { # case 1 : if xmaplim and ymaplim are specified xlim.map <- xmaplim; ylim.map <- ymaplim } else { # case 2 : general lon <- x; lat <- y fram <- matrix(rep(0,2*2*(nx+ny)),ncol=2) fram[,1] <- c(lon,rep(lon[nx],ny),rev(lon),rep(lon[1],ny)) fram[,2] <- c(rep(lat[1],nx),lat,rep(lat[ny],nx),rev(lat)) if (rot) { fram <- rot2lonlat(fram,pollam=plon,polphi=plat) } xlim.map <- range(fram[,1],finite=TRUE) ylim.map <- range(fram[,2],finite=TRUE) if (!missing(lonlim) & !missing(latlim)) { # if lonlim and latlim are specified xlim.map <- range(c(xlim.map,lonlim),finite=TRUE) ylim.map <- range(c(ylim.map,latlim),finite=TRUE) } } # increments for plotting the grids bylongrd <- diff(xlim.map)/200 bylatgrd <- diff(ylim.map)/200 # get the map coordinates in the appropriate projection if (mapdat != "none") { if (mapdat == "continents") { # continents must be read from external file. Are not available #Êin maps package, data was extracted from package clim.pact load("~/work/R/data/continents.Rdata") mapcors <- mapproject(x=continents) } else { mapcors <- map(mapdat,xlim=xlim.map, ylim=ylim.map,plot=FALSE,interior=FALSE) if(projection == "stereographic" & ylim.map[1]>0) { outside <- mapcors$y < (ylim.map[1]-1) & !is.na(mapcors$y) mapcors$x <- mapcors$x[-which(outside)] mapcors$y <- mapcors$y[-which(outside)] } if(projection == "stereographic" & ylim.map[1]<0) { outside <- mapcors$y > (ylim.map[2]-1) & !is.na(mapcors$y) mapcors$x <- mapcors$x[-which(outside)] mapcors$y <- mapcors$y[-which(outside)] } mapcors <- mapproject(x=mapcors$x,y=mapcors$y) } # thin the map if required if (thin.map > 0) { for (tt in 1:thin.map) { llll <- length(mapcors$x) keep <- rep(c(TRUE,FALSE),length=llll) keep[(is.na(mapcors$x) | is.na(mapcors$y))] <- TRUE keep[llll] <- TRUE keep[c(is.na(mapcors$x[-1]),TRUE)] <- TRUE keep[c(TRUE,is.na(mapcors$x[-llll]))] <- TRUE mapcors <- list(x=mapcors$x[keep],y=mapcors$y[keep]) } rm(keep) } } # determine default breaks zmax <- max(zz); zmin <- min(zz) if (missing(breaks)) { breaks <- seq(zmin,zmax,length=length(col)+1) } else { if (length(breaks) != length(col)+1) { stop("** ERROR ** length of breaks must be length of colors plus 1") } } # determine colors for each pixel val2col <- function(val) { if (is.na(val)) {return(col.miss)} hh <- which(val <= c(breaks,val))[1] if (hh==(length(breaks)+1)) {return(col.miss)} if (hh==1) {return(col.miss)} col[hh-1] } zcols <- sapply(zz,FUN=val2col) if(con){ res<-contourLines(x,y,z) contours_x <- unlist(sapply(res, function(x) c(x$x, NA))) contours_y <- unlist(sapply(res, function(x) c(x$y, NA))) contours <- list(x=contours_x, y=contours_y) contu<-mapproject(contours) } # AND HERE WE GO ---------------------------------------------- # prepare the plot with the frame plot(x=c(1),y=c(1),type="n",axes=FALSE, xlim=xxlim,ylim=yylim,xlab="",ylab="") # plot the colored pixel polygons for (i in seq(1,length(xx))) { polygon(x=c(xxll[i],xxlr[i],xxur[i],xxul[i]), y=c(yyll[i],yylr[i],yyur[i],yyul[i]), col=zcols[i],border=zcols[i],xpd=FALSE) } if (con) { lines(x=contu$x,y=contu$y,xpd=FALSE) } # add the mapoutlines on top if (mapdat != "none") { lines(x=mapcors$x,y=mapcors$y,xpd=FALSE) } # add grid-lines for lon and lat circles if(longr=="yes"){ for (i in longrds) { #ygr <- seq(-90,90,by=bylatgrd) if(ylim.map[1]>0) ygr <- seq(ymaplim[1]-0.5*(diff(y[1:2])),ymaplim[2],by=bylatgrd) if(ylim.map[1]<0) ygr <- seq(ymaplim[2]+0.5*(diff(y[1:2])),ymaplim[1],by=-bylatgrd) grd <- mapproject(x=rep(i,length(ygr)),y=ygr) lines(x=grd$x,y=grd$y,lty=ltygrds,lwd=lwdgrds,xpd=FALSE) }} if(latgr=="yes"){ for (i in latgrds) { xgr <- seq(-180,180,by=bylongrd) grd <- mapproject(x=xgr,y=rep(i,length(xgr))) lines(x=grd$x,y=grd$y,lty=ltygrds,lwd=lwdgrds,xpd=FALSE) }} # add a box # box() } reshapefield <- function(lon,lat,matr) { # Reshape a n x p*q matrix in a three-dimensional n x p x q array or reshapes # a three-dimensional n x p x q array into a n x p*q matrix. # # Description: # # Returns a n x p x q array or a n x (p*q) matrix # # Usage: # # reshapefield(lon,lat,matr) # # Input: # # matr: a n x p*q matrix or a n x p x q array with n time dimension # and p and q space dimensions # lon: vector of longitudes # lat: vector of latitudes # # Output: # # $out: a n x p x q array or a n x p x q matrix # # Authors: # # Dag Johan Steinskog 7 June 2005 # Caio Coelho # # Examples: # # # Reshape a n x p*q matrix in a three-dimensional n x p x q array # x <- array(1:8,dim=c(2,4)) # lon <- seq(50,70,20) # lat <- seq(-10,10,20) # reshapefield(lon,lat,x)$out # # # Reshape a three-dimensional n x p x q array in a n x p*q matrix # y <- array(1:8,dim=c(2,2,2)) # lon <- seq(50,70,20) # lat <- seq(-10,10,20) # reshapefield(lon,lat,y)$out if (length(dim(matr))==2) { dims <- dim(matr) n <- dims[1] spacedim <- dims[2] if(spacedim != length(lat)*length(lon)) stop("Incompatible dimensions") out<-array(t(matr),c(length(lon),length(lat),n)) } if (length(dim(matr))==3) { dims <- dim(matr) n <- dims[3] p <- dims[1] q <- dims[2] if(p*q != length(lat)*length(lon)) stop("Incompatible dimensions") out <- t(array(as.vector(matr),c(p*q,n)))} invisible(list(out=out)) } returnperiod<-function(excess,parameter){ # Compute return period for a given matrix of excesses and a given array of # Generalized Pareto distribution parameters # # Description: # # Returns return period of excesses # # Usage: returnperiod(excess,parameter) # # Input: # # excess: p x q matrix of excesses # # parameter: 2 x p x q array with scale (first level of the array) and # shape (second level of the array) parameters # # Output: # # A matrix containing the return period of excesses # # Authors: # # Caio Coelho 28 Feb 2006 rp<-excess for(i in 1:dim(excess)[1]) for(j in 1:dim(excess)[2]) if (any(is.na(c(excess[i,j],parameter[1,i,j],parameter[2,i,j])))){ rp[i,j]<-NA } else{ rp[i,j]<-1/(1-pgpd(excess[i,j],loc=0,scale=parameter[1,i,j],shape=parameter[2,i,j])) } rp } tvt<-function(y,span,index,prob){ # Compute time-varying threshold for a given monthly time series # # Description: # # Returns smooth long term mean, time-varying threshold and fitted values # given by long term mean plus mean seasonal cycle # # Usage: svt(y,span,index,prob) # # Input: # # y: vector containing a monthly time series of data # # span: fraction of the total number of points of 'y' to be used to compute # the long term mean # # index: index vector composed of logicals (TRUE and FALSE) of the same # length as 'y' defining the months to be considered when computing the # time-varying threshold # # prob: a value between 0 and 1 that informs the percentage of point # that will be left above the time-varying threshold # # Outputs # # $mfit: fitted data (i.e. long term mean + mean annual cycle) # # $ltm: long term mean # #$theshold: time-varying threshold # # Authors: # # Caio Coelho 28 Feb 2006 # Chris Ferro x<-1:length(y) lo<-loess(y~x,span=span) resid<-rep(NA,length(y)) resid[!is.na(y)]<-lo$residuals meanresid<-rep(NA,12) for(i in 1:12){ meanresid[i]<-mean(resid[seq(i,length(y),12)],na.rm=T) } fitted<-rep(NA,length(y)) fitted[!is.na(y)]<-lo$fitted mfit<-fitted+rep(meanresid,length(y)/12) step <- sort(y[index] - mfit[index]) step <- unique(step[step > 0]) th <- mfit[index] p <- mean(y[index] > th, na.rm = TRUE) k <- 0 while(p > prob) { k <- k + 1 th <- mfit[index] + step[k] p <- mean(y[index] > th, na.rm = TRUE) } list(threshold=th,ltm=fitted,mfit=mfit) } xdependence <- function(array3d,ts,u,fun) { # Compute extreme dependence measures between a given three-dimensional array # with first two space dimensions (e.g. longitude and latitude) and third # time dimension, and a given time series with the same length as the time # dimension of the three dimensional array # # Description: # # Returns a map (matrix) of the extreme dependence measures described in # Coles et al. (1999) and in section 8.4 of Coles (2001) # # Usage: # # xdependence(array3d,ts,u,fun) # # Input: # # array3d: three-dimensional array with p longitude points and q latitude # points as the first two dimensions and n as the third time # dimension # # ts: time series of data with the same length as the time dimension # of array3d for which extreme dependence with each grid point of # array3d will be computed # # u: high threshold between 0 and 1 that will define the extreme # level to compute dependence. Must be a high quantile (e.g. 0.95) # # fun: String of characters that defines which extreme dependence # measure is computed (e.g. "chi", "chibar" or "chicount"). # If fun is "chi", computes chi(u) as in Eq.(3.2) of # Coles et al. 1999 and in section 8.4 of Coles (2001). # If fun is "chibar", computes chibar(u) as in section 3.3.2 of # Coles et al. 1999 and in section 8.4 of Coles (2001). # If fun is "chicount", computes the conditional probability of # each grid point variable of array3d to be large (above a large # threshold u) conditional on (given that) the variable ts is # also large (with value above the same large threshold u) # by counting values above threshold and computing relative # frequencies # # Output: # # $out: an map (matrix) of the extreme dependence measure for each grid # point of array3d. # # Authors: # # Caio Coelho 21 Dec 2005 # Christopher Ferro # # References: Coles, S., Heffernan, J., and Tawn, J., 1999: Dependence measures # of extreme value analyses. Extremes 2:4, 339-365. # # Coles, S., 2001: An introduction to statistical modeling of # extreme values. Spring series in statistics. 208pp. # # Examples: # # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 1000) # data <- array(rnorm(prod(dim)), dim) # xdependence(data,data[4,5,],u=0.95,fun="chi")$out # xdependence(data,data[4,5,],u=0.95,fun="chibar")$out # xdependence(data,data[4,5,],u=0.95,fun="chicount")$out out1<-aperm(apply(array3d,c(1,2),rank),c(2,3,1))/dim(array3d)[3] if (fun == "chi"){ ts<-rank(ts)/length(ts) chi <- function(y){ 2-(log((sum(ts<=u & y<=u)+0.5)/(length(ts)+1))/(log((sum(ts<=u)+0.5)/(length(ts)+1)))) } out<-apply(out1,c(1,2),chi) out[out<0]<-0 } if (fun == "chibar"){ ts<-rank(ts)/length(ts) chibar <- function(y){ ((2*log((sum(ts>=u)+0.5)/(length(ts)+1)))/log((sum(ts>=u & y>=u)+0.5)/(length(ts)+1)))-1 } out<-apply(out1,c(1,2),chibar) } if (fun == "chicount"){ chicount <- function(y){ length(y[ts>=quantile(ts,u)][y[ts>=quantile(ts,u)]>=quantile(y,u)])/length(y[ts>=quantile(ts,u)]) } out<-apply(array3d,c(1,2),chicount) } invisible(list(out=out)) } xdependence1 <- function(lon,lat,array3d,ts,u,fun,nonmissing=0.5) { # Compute extreme dependence measures between a given three-dimensional array # with first two space dimensions (e.g. longitude and latitude) and third # time dimension, and a given time series with the same length as the time # dimension of the three dimensional array # # Note: This function allows the user to specify through the parameter # nonmissing (defaul is 0.5) the acceptable percentage of nonmissing # values in the time series of each grid point of the three-dimensional # array. This is the main difference between this function and # xdependence.r # # Description: # # Returns a map (matrix) of the extreme dependence measures described in # Coles et al. (1999) and in section 8.4 of Coles (2001) # # Usage: # # xdependence1(lon,lat,array3d,ts,u,fun,nonmissing) # # Input: # # lon: vector with p longitude values # # lat: vector with q latitude values # # array3d: three-dimensional array with p longitude points and q latitude # points as the first two dimensions and n as the third time # dimension # # ts: time series of data with the same length as the time dimension # of array3d for which extreme dependence with each grid point of # array3d will be computed # # u: high threshold between 0 and 1 that will define the extreme # level to compute dependence. Must be a high quantile (e.g. 0.95) # # fun: String of characters that defines which extreme dependence # measure is computed (e.g. "chi", "chibar" or "chicount"). # If fun is "chi", computes chi(u) as in Eq.(3.2) of # Coles et al. 1999 and in section 8.4 of Coles (2001). # If fun is "chibar", computes chibar(u) as in section 3.3.2 of # Coles et al. 1999 and in section 8.4 of Coles (2001). # If fun is "chicount", computes the conditional probability of # each grid point variable of array3d to be large (above a large # threshold u) conditional on (given that) the variable ts is # also large (with value above the same large threshold u) # by counting values above threshold and computing relative # frequencies # #nonmissing: Only grid points with fraction given by this fraction # (between 0 and 1) of non-missing values are used to compute # the statistics specified in fun. Default is 0.5. # # Output: # # $out: an map (matrix) of the extreme dependence measure for each grid # point of array3d. # # Authors: # # Caio Coelho 21 Dec 2005 # Christopher Ferro # # References: Coles, S., Heffernan, J., and Tawn, J., 1999: Dependence measures # of extreme value analyses. Extremes 2:4, 339-365. # # Coles, S., 2001: An introduction to statistical modeling of # extreme values. Spring series in statistics. 208pp. # # Examples: # # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 1000) # data <- array(rnorm(prod(dim)), dim) # xdependence1(x,y,data,data[4,5,],u=0.95,fun="chi")$out # xdependence1(x,y,data,data[4,5,],u=0.95,fun="chibar")$out # xdependence1(x,y,data,data[4,5,],u=0.95,fun="chicount")$out if (all(is.na(ts))) stop("all missing values for ts") data <- reshapefield(lon,lat,array3d)$out # compute percentage of non-missing values at each grid point aux <- apply(data,2,function(x)sum(!is.na(x))/(length(x))) # identify grid points with more than 50% missing values index <- (1:length(aux))[aux < nonmissing] data[,index]<-NA array3d<-reshapefield(lon,lat,as.matrix(data))$out out1<-aperm(apply(array3d,c(1,2),rank, na.last = "keep"),c(2,3,1)) aux<-apply(array3d,c(1,2),function(x)sum(!is.na(x))) aux1<-array(NA,c(dim(out1))) for (i in 1:(dim(array3d)[3])){ aux1[,,i] <- out1[,,i]/aux } if (fun == "chi"){ ts<-rank(ts, na.last = "keep")/sum(!is.na(ts)) chi <- function(y){ index <- !(is.na(y) | is.na(ts)) if(sum(index)==0) return(NA) 2-(log((sum(ts[index]<=u & y[index]<=u)+0.5)/(length(ts[index])+1))/(log((sum(ts[index]<=u)+0.5)/(length(ts[index])+1)))) } out<-apply(aux1,c(1,2),chi) out[out<0]<-0 } if (fun == "chibar"){ ts<-rank(ts, na.last = "keep")/sum(!is.na(ts)) chibar <- function(y){ index <- !(is.na(y) | is.na(ts)) if(sum(index)==0) return(NA) ((2*log((sum(ts[index]>=u)+0.5)/(length(ts[index])+1)))/log((sum(ts[index]>=u & y[index]>=u)+0.5)/(length(ts[index])+1)))-1 } out<-apply(aux1,c(1,2),chibar) } if (fun == "chicount"){ chicount <- function(y){ index <- !(is.na(y) | is.na(ts)) if(sum(index)==0) return(NA) length(y[index][ts[index]>=quantile(ts[index],u)][y[index][ts[index]>=quantile(ts[index],u)]>=quantile(y[index],u)])/length(y[index][ts[index]>=quantile(ts[index],u)]) } out<-apply(array3d,c(1,2),chicount) } invisible(list(out=out)) } xexcess <- function(array3d,fun = 'meanexcess',p=0.9,threshold = FALSE, upper=TRUE) { # Compute mean excess and variance of excess for a given three-dimensional # array with first two space dimensions and third time dimension. # # Description: # # Returns a matrix with the mean excess or the variance of excess for a # given quantile or user defined threshold. # # Usage: # # xexcess(array3d,fun,p=0.9,threshold = FALSE, upper=TRUE) # # Input: # # array3d: three-dimensional array with p longitude points and # q latitude points as the first 2 dimensions and n as the # third time dimension # # fun: String containing the name of the statistics to be computed # that must be one of the following options "meanexcess", # "varexcess", or "medianexcess". # # p: The quantile to be computed (Default is 0.9) or a particular threshold # chosen by the user. If p is a threshold then the argument threshold # (see below) must be set to TRUE. # # threshold: Logic argument. Default is FALSE. If TRUE p must be a given # threshold value. # # upper: Logic argument. Defalts is TRUE (examines upper tail of the # distribution).If FALSE lower tail is examined. # # Output: # # $out: a matrix of the computed statistics # # Authors: # # Dag Johan Steinskog 14 June 2005 # Caio Coelho # Christopher Ferro # # Examples: # # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 100) # data <- array(rnorm(prod(dim)), dim) # xexcess(data,"meanexcess")$out # xexcess(data,"varexcess")$out # xexcess(data,"meanexcess",p=1.5,threshold =TRUE)$out # xexcess(data,"meanexcess",p=1.5,threshold =TRUE,upper=FALSE)$out # if(fun=="meanexcess"){ mean.excess <- function(y,p) { if(threshold) u<- p else u<-quantile(y,p,na.rm=T) if(upper) mean(y[y>u]-u,na.rm=T) else mean(y[yu]-u,na.rm=T) else var(y[yu]-u,na.rm=T) else median(y[y 14 Dec 2005 # Christopher Ferro # # Examples: # # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 1000) # data <- array(rnorm(prod(dim)), dim) # xgev(data) # xgev(data,n=50) # xgev(data,upper=FALSE) # d<-dim(array3d) block.max <- function(x){ apply(matrix(x, n, d[3]/n), 2, max) } bmax<-aperm(apply(array3d,c(1,2),block.max),c(2,3,1)) estim.par <- function(y) { if(upper) fgev(y, std.err = FALSE)$est else fgev(-y, std.err = FALSE)$est } out<-apply(bmax,c(1,2),estim.par) invisible(list(out=out)) } xindex <- function(y,u,threshold = FALSE) { # Description: # # Evaluates the intervals estimator for the extremal index, an index for # time clusters, for a given time series and threshold. # # Usage: # # xindex(y,u,threshold=FALSE) # # Arguments: # # y: a time series for which the extremal index will be computed # # u: Threshold. Can be either a quantile (e.g. 0.8) or a threshold value # (see definition of theshold below). # # threshold: Logical. If FALSE (the default) u must be a quantile value. # If TRUE, then u must be a threshold value. # # Value: # # Estimate of the extremal index. # # Details: # # Warning: # # The estimator is not constrained to lie in [0,1] and a # default value of 1 is returned if there are fewer than two # extreme values. # # Authors: # # Christopher Ferro 21 Dec 2005 # Caio Coelho # # References: # # Ferro CAT & Segers J (2003) Inference for clusters of extreme # values. Journal of the Royal Statistical Society B 65, 545-556. # # See Also: # # Examples: # # y<-rnorm(1000,25,2) # xindex(y,0.9) # xindex(y,28,threshold=TRUE) # generates a logical vector indicating which positions correspond to # extreme values. if(threshold) z <- (y >= u) else z <- (y >= quantile(y, u,na.rm=T)) if(sum(z,na.rm=T) <= 1) { warning("estimator undefined: too few exceedances") return(1) } else { nz <- length(z) # length of sequence s <- c(1:nz)[z] # exceedance times t <- diff(s) # interexceedance times if(max(t,na.rm=T) <= 2) { t1 <- mean(t,na.rm=T) t2 <- mean(t^2,na.rm=T) } else { t1 <- mean(t-1,na.rm=T) t2 <- mean((t-1)*(t-2),na.rm=T) } } min(1,2*(t1^2)/t2) } xindexfield <- function(y,u,threshold = FALSE) { # Description: # # Computes the intervals estimator for the extremal index, an index for # time clusters, at each grid point of a three-dimensional array of data # with first two space dimensions (e.g. longitude and latitude) and # third time dimension # # Usage: # # xindexfield(y,u,threshold=FALSE) # # Arguments: # # y: a three-dimensional array, with first two space dimensions (e.g. # longitude and latitude) and third time dimension # # u: Threshold. Can be either a quantile (e.g. 0.8) or a threshold value # (see definition of theshold below). # # threshold: Logical. If FALSE (the default) u must be a quantile value. # If TRUE, then u must be a threshold value. # # Value: # # Estimate of the extremal index at each grid point. # # Details: # # Warning: # # The estimator is not constrained to lie in [0,1] and a # default value of 1 is returned if there are fewer than two # extreme values. # # Authors: # # Christopher Ferro 3 Jan 2006 # Caio Coelho # # References: # # Ferro CAT & Segers J (2003) Inference for clusters of extreme # values. Journal of the Royal Statistical Society B 65, 545-556. # # Examples: # # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 100) # data <- array(rnorm(prod(dim),25,2), dim) # xindexfield(data,0.9) # xindexfield(data,28,threshold=TRUE) apply(y,c(1,2),xindex,u,threshold) } xpareto <- function(array3d,p=0.9,upper=TRUE) { # Compute shape and scale parameters of a Generalized Pareto Distribution for # a given three-dimensional array with first two space dimensions and third # time dimension. # # Description: # # Returns an array with the shape and scale parameters of a Generalized # Pareto Distribution for the exceedances above a given quantile that # defines the threshold. # # Usage: # # xpareto(array3d,p=0.9, upper=TRUE) # # Input: # # array3d: three-dimensional array with p longitude points and q latitude # points as the first 2 dimensions and n as the third time # dimension # # p: quantile to be computed (Default is 0.9). # # upper: Logic argument. Default is TRUE (examines upper tail of the # distribution). If FALSE lower tail is examined. # # Output: # # $out: an array of the computed statistics. First dimension contains # the estimated parameters. (e.g. $out[1,,] is the scale parameter # and $out[2,,] is the shape parameter). Second and third dimensions # are the same space dimensions as of array3d. # # # Authors: # # Dag Johan Steinskog 16 June 2005 # Caio Coelho # Christopher Ferro # # Examples: # # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 1000) # data <- array(rnorm(prod(dim)), dim) # xpareto(data) # xpareto(data,p=0.95) # xpareto(data,p=0.1,upper=FALSE) # estim.par <- function(y,p) { if(upper) fpot(y, quantile(y, p,na.rm=T), "gpd",std.err = FALSE)$est else fpot(-y, quantile(-y, 1-p,na.rm=T), "gpd",std.err = FALSE)$est } out<-apply(array3d,c(1,2),estim.par,p) invisible(list(out=out)) } xparetotvt <- function(lon,lat,array3d,span,p=0.9,index,nonmissing=0.5) { # Fit Generalized Pareto distribution with time-varying threshold at each grid # point for a given three-dimensional array of montly data with first two space # dimensions (e.g. longitude and latitude) and third time dimension. # # Description: # # Returns maps of KS test p-values, fitted parameters (constant scale and # constant shape), map of fraction of points above the time-varying # threshold, and the time varying threshold. # # Usage: # # xparetotvt(lon,lat,array3d,span,p,index,nonmissing) # # Input: # # lon: vector with p longitude values # # lat: vector with q latitude values # # array3d: a three-dimensional array of monthly data with p longitude points # and q latitude points as the first two dimensions and n as the # third time dimension # # span: fraction of the total number of points 'n' to be used to compute # the long term mean # # p: a value between 0 and 1 that informs the percentage of point # that will be left below the time-varying threshold # # index: index vector of length n composed of logicals (TRUE and FALSE) # defining the months to be considered when computing the # time-varying threshold # #nonmissing: Only grid points with fraction given by 'nonmissing' # (between 0 and 1) of non-missing values are used to estimate the # Generalized Pareto distribution parameters. Default is 0.5. # # Output: # # $pvaloutput: p x q map of KS test p-values # # $output: 2 x p x q array with scale (first level of the array) and # shape (second level of the array) parameters # # $frac: p x q map with fraction of points above the time-varying threshold # # $th: p x q x n' array containing the time varying threshold, where # n' is the number of months considered when computing the # time-varying threshold # # Authors: # # Caio Coelho 28 Feb 2006 # Chris Ferro # reshape three-dimensional array into a data matrix with # time as first dimension and space as sencond dimension data <- reshapefield(lon,lat,array3d)$out # compute percentage of non-missing values at each grid point aux <- apply(data[index,],2,function(x)sum(!is.na(x))/(length(x))) # identify grid points with less than 50% missing values indexgrid <- (1:length(aux))[aux >= nonmissing] out1 <- rep(NA,dim(data)[2]) out2 <- rep(NA,dim(data)[2]) pval1 <- rep(NA,dim(data)[2]) frac1 <- rep(NA,dim(data)[2]) mthresh <- matrix(nrow=(dim(data)[1])/(12/sum(index[1:12])),ncol=dim(data)[2]) # function to perform KS test ff <- function(x,loc,scale,shape){ if(all(is.na(c(scale,shape)))) {return(NA)} ks.test(x,"pgpd",loc,scale,shape)$p } # Estimate Generalized Pareto distribution parameters (scale and shape) for # exceedances above the threshold p estim.par <- function(y,p,index) { threshold<-tvt(y,span,index,1-p)$threshold yy<-y[index] fraction<-sum(!is.na(yy[yy>threshold]))/sum(!is.na(yy)) params<-gpd.fit(yy[!is.na(yy)], threshold[!is.na(yy)], show=FALSE)$mle pval<-ff((yy-threshold)[yy>threshold],0,params[1],params[2]) list(params=params,pval=pval,fraction=fraction,threshold=threshold) } out<-apply(data[,indexgrid],2,estim.par,p,index) for (i in 1:length(indexgrid)){ #scale parameter out1[indexgrid[i]]<-out[[i]]$params[1] mthresh[,indexgrid[i]]<-out[[i]]$threshold #shape parameter out2[indexgrid[i]]<-out[[i]]$params[2] #p-value pval1[indexgrid[i]]<-out[[i]]$pval frac1[indexgrid[i]]<-out[[i]]$fraction } aux<-reshapefield(lon,lat,t(as.matrix(pval1)))$out pvaloutput<-aux[,,1,drop=T] aux1<-reshapefield(lon,lat,t(as.matrix(frac1)))$out frac<-aux1[,,1,drop=T] th<-reshapefield(lon,lat,mthresh)$out output<-array(NA,c(2,dim(array3d)[1],dim(array3d)[2])) output[1,,]<-reshapefield(lon,lat,t(as.matrix(out1)))$out output[2,,]<-reshapefield(lon,lat,t(as.matrix(out2)))$out invisible(list(pvaloutput=pvaloutput,output=output,frac=frac,th=th)) } xparetotvtcov<-function(lon,lat,array3d,span,p=0.9,index,covariates=NULL,sigl=NULL,shl=NULL,siglink=identity,shlink=identity,nonmissing=0.5) { # Fit Generalized Pareto distribution with time-varying threshold at each grid # point for a given three-dimensional array of montly data with first two space # dimensions (e.g. longitude and latitude) and third time dimension. Allows # linear modelling of the paramters. # # Description: # # Returns fitted parameters. # # Usage: # # xparetotvtcov(lon,lat,array3d,span,p,index,covariates,sigl,shl,siglink,shlink,nonmissing) # # Input: # # lon: vector with p longitude values # # lat: vector with q latitude values # # array3d: a three-dimensional array of monthly data with p longitude points # and q latitude points as the first two dimensions and n as the # third time dimension # # span: fraction of the total number of points 'n' to be used to compute # the long term mean # # p: a value between 0 and 1 that informs the percentage of point # that will be left below the time-varying threshold # # index: index vector of length n composed of logicals (TRUE and FALSE) # defining the months to be considered when computing the # time-varying threshold # #covariates: Matrix of covariates for generalized linear modelling of # the parameters. The number of rows should be 'n' (i.e. the same as # the time dimension of 'array3d'. # # sigl, shl: Numeric vectors of integers, giving the columns of 'ydat' # that contain covariates for generalized linear modelling of # the scale and shape parameters repectively . # # siglink, shlink: Inverse link functions for generalized linear # modelling of the scale and shape parameters repectively. # #nonmissing: Only grid points with fraction given by 'nonmissing' # (between 0 and 1) of non-missing values are used to estimate the # Generalized Pareto distribution parameters. Default is 0.5. # # Output: # # $output: n' x p x q array with estimates parameters, where n' is # the total number of estimated parameters # # Authors: # # Caio Coelho 28 Feb 2006 # Chris Ferro # reshape three-dimensional array into a data matrix with # time as first dimension and space as sencond dimension data <- reshapefield(lon,lat,array3d)$out # compute percentage of non-missing values at each grid point aux <- apply(data[index,],2,function(x)sum(!is.na(x))/(length(x))) # identify grid points with less than 50% missing values indexgrid <- (1:length(aux))[aux >= nonmissing] outaux<-matrix(ncol=dim(data)[2],nrow=(length(sigl)+length(shl)+2)) # Estimate Generalized Pareto distribution parameters (scale and shape) for # exceedances above the threshold p estim.par <- function(y,p,index,covariates,sigl,shl,siglink,shlink) { threshold<-tvt(y,span,index,1-p)$threshold yy<-y[index] covs<-as.matrix(covariates[index,]) index1 <- !(is.na(yy) | apply(t(apply(covs,1,is.na)),1,any)) params<-mygpd.fit(yy[index1],threshold[index1],ydat=as.matrix(covs[index1,]),sigl=sigl,shl=shl,siglink=siglink,shlink=shlink,show=FALSE)$mle list(params=params) } out<-apply(data[,indexgrid],2,estim.par,p,index,covariates,sigl,shl,siglink,shlink) for(j in 1:(length(sigl)+length(shl)+2)){ for (i in 1:length(indexgrid)){ outaux[j,indexgrid[i]]<-out[[i]]$params[j] } } output<-aperm(reshapefield(lon,lat,outaux)$out,c(3,1,2)) invisible(list(output=output)) } zoom <- function(x, y, data, t=1,..., same = FALSE, mark = "+") { # Zoom in on either a two-dimensional map (matrix) with p longitude and q # latitude space dimensions or a three-dimensional array (field) with # p longitude and q latitude as the first two space dimensions and n as the # third time dimension, by selecting two corners that define a rectangular # region of the space and a subset of time slices # # Description: # # Returns a subset of a two-dimensional map (matrix) or a # subset of a three-dimensional array (field) and re-plot data # of the first time slice of a rectangular region selected by # clicking on an image. # # Usage: # # zoom(x, y, data, t=1, ..., same = FALSE, mark = "+") # # Arguments: # # x: Vector of longitudes # # y: Vector of latitudes # # data: two-dimensional map (matrix) or three-dimensional array of data. # Space dimensions of `data' must be `length(x)' by `length(y)' # respectively. # Three-dimensional array must have the following dimention order: # first dimension for longitude, second dimension for latitude, # and third dimension for time # # t: Subset of times slices to be extracted. Default is 1 (i.e. extract # only the first time slice # # ...: Additional arguments passed to `image' # # same: Logical: if `FALSE' (the default), zoomed images are re-plotted # in a new window; if `TRUE', each image is plotted in the same # window without cleaning the longitude and latitude coordinates # previously plotted # # mark: A vector of numbers or strings with which selected cells are # labelled. If shorter than the number of grid cells, values are # recycled; if `NULL', selections are not marked # # Details: # # Clicking on two grid cells in an image plot created with grid # lines at locations `x' and `y' causes the image of the first time slice # to be re-plotted for the rectangular region thus defined. # # Output: # # A list with three components is returned invisibly: # # x: Horizontal (lon) coordinates of grid lines used to plot the final # zoomed image # # y: Vertical (lat) coordinates of grid lines used to plot the final # zoomed image # # data: Matrix or three-dimensional array containing the selected data # # t: Extracted times slices # # See Also: # # `identify', `image' # # Author: # # Chris Ferro 27 October 2005 # Dag Johan Steinskog # Caio Coelho # # Examples: # # # Zoom in three-dimensional array # x <- seq(-20, 20, 5) # y <- seq(30, 60, 5) # dim <- c(length(x), length(y), 100) # data <- array(rnorm(prod(dim)), dim) # output <- zoom(x, y,data,t=c(1:3)) # output <- zoom(x, y,data,t=c(2,4,6)) # output <- zoom(x, y,data,t=c(10,20,30,40)) # # # Zoom in matrix # matr <- data[,,1] # output <- zoom(x, y,matr) if (length(dim(data))==2) { xnew <- x ynew <- y image(x, y, data) if (min(x)<0){map("world",add=TRUE)} else{map("world2",add=TRUE)} cat("Left-click to select two corners of zoom region; right-click to finish.\n") repeat { grid <- expand.grid(y = ynew, x = xnew)[, 2:1] if(is.null(mark)) lab <- 1:nrow(grid) else lab <- rep(mark, length.out = nrow(grid)) pts <- identify(grid, labels = lab, n = 2, plot = !is.null(mark), offset = 0) if(length(pts) < 2) break loc <- grid[pts, ] xnew <- x[x >= min(loc$x) & x <= max(loc$x)] ynew <- y[y >= min(loc$y) & y <= max(loc$y)] datanew <- data[match(xnew, x), match(ynew, y), drop = FALSE] if(!same) x11() par(mfg = par("mfg")) image(xnew, ynew, datanew, ...) if (min(x)<0){map("world",add=TRUE)} else{map("world2",add=TRUE)} } } if(length(dim(data))==3) { z <- data[,,1] xnew <- x ynew <- y tnew <- t image(x, y, data[,,1]) if (min(x)<0){map("world",add=TRUE)} else{map("world2",add=TRUE)} cat("Left-click to select two corners of zoom region; right-click to finish.\n") repeat { grid <- expand.grid(y = ynew, x = xnew)[, 2:1] if(is.null(mark)) lab <- 1:nrow(grid) else lab <- rep(mark, length.out = nrow(grid)) pts <- identify(grid, labels = lab, n = 2, plot = !is.null(mark), offset = 0) if(length(pts) < 2) break loc <- grid[pts, ] xnew <- x[x >= min(loc$x) & x <= max(loc$x)] ynew <- y[y >= min(loc$y) & y <= max(loc$y)] znew <- z[match(xnew, x), match(ynew, y), drop = FALSE] datanew <- data[match(xnew, x), match(ynew, y), t, drop = FALSE] if(!same) x11() par(mfg = par("mfg")) image(xnew, ynew, znew, ...) if (min(x)<0){map("world",add=TRUE)} else{map("world2",add=TRUE)} } } invisible(list(x = xnew, y = ynew, data = datanew, t = t)) }