R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("ProbForecastGOP-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('ProbForecastGOP') Loading required package: fields fields is loaded use help(fields) for an overview of this library Loading required package: RandomFields > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "Emp.variog" > > ### * Emp.variog > > flush(stderr()); flush(stdout()) > > ### Name: Emp.variog > ### Title: Empirical variogram of forecast errors averaged over time > ### Aliases: Emp.variog > ### Keywords: file > > ### ** Examples > > ## Loading data > data(slp) > day <- slp$date.obs > id <- slp$id.stat > coord1 <- slp$lon.stat > coord2 <- slp$lat.stat > obs <- slp$obs > forecast <- slp$forecast > > ## Computing variogram > ## No specified cutpoints, no specified maximum distance > ## Default number of bins > variogram <- Emp.variog(day=day,obs=obs,forecast=forecast,id=id, + coord1=coord1,coord2=coord2,cut.points=NULL,max.dist=NULL,nbins=NULL) > ## Plotting variogram > plot(variogram$bin.midpoints,variogram$empir.variog,xlab="Distance", + ylab="Semi-variance",main="Empirical variogram") > > ## Computing variogram > ## Specified cutpoints, specified maximum distance > ## Unspecified number of bins > variogram <- + Emp.variog(day=day,obs=obs,forecast=forecast,id=id,coord1=coord1, + coord2=coord2,cut.points=seq(0,1000,by=5),max.dist=800,nbins=NULL) > ## Plotting variogram > plot(variogram$bin.midpoints,variogram$empir.variog,xlab="Distance", + ylab="Semi-variance",main="Empirical variogram") > > ## Don't show: ## The function is currently defined as > function(day,obs,forecast,id,coord1,coord2,cut.points=NULL,max.dist=NULL,nbins=300){ + # default values + if(missing(cut.points)) + cut.points <- NULL + if(missing(max.dist)) + max.dist <- NULL + if(missing(nbins)) + nbins <- NULL + #INPUT CHECK + # Here we check if the input is right. + l.day <- length(day) + l.obs <- length(obs) + l.for <- length(forecast) + l.id <- length(id) + l.coord1 <- length(coord1) + l.coord2 <- length(coord2) + if(sum((c(l.day,l.obs,l.for,l.id,l.coord1,l.coord2)/l.day)==rep(1,6))!=6){ + stop("Error in the dimension of the data") + } + if(sum(is.numeric(obs)==rep("TRUE",l.obs))=2 & (sum(is.numeric(cut.points)==rep("TRUE",l.cuts))=2 & (sum(is.numeric(cut.points)==rep("TRUE",l.cuts))==l.cuts)){ + if(sum(order(cut.points)==seq(1:l.cuts)) < l.cuts){ + stop("Cut points should be in increasing order") + } + if(sum(cut.points >= 0) < l.cuts){ + stop("Cut points should be non-negative numbers") + } + if(length(cut.points)!=length(unique(cut.points))){ + stop("The vector with cut points should not contain repeated entries") + } + } + + ## check on the max.dist + l.mdist <- length(max.dist) + if(l.mdist > 1){ + print("max.dist is a numeric field, not a vector")} + if(l.mdist==1){ + if(is.numeric(max.dist)==FALSE){ + stop("max.dist is a numeric field") + } + if(max.dist < 0){ + stop("max.dist should be a positive number") + } + } + ## check on the number of bins + l.nbins <- length(nbins) + if(l.nbins==0 & l.cuts==0){ + nbins <- 300 + } + if(l.nbins==1 & l.cuts >=2){ + nbins=NULL + } + l.nbins <- length(nbins) + if(l.nbins >1){ + stop("nbins should be an integer: not a vector") + } + + if(l.nbins==1){ + if(ceiling(nbins)!=nbins){ + stop("Invalid input: the number of bins should be a positive integer") + } + if(ceiling(nbins)==nbins & nbins < 0){ + stop("Invalid input: the number of bins should be a positive integer") + } + } + + + # ESTIMATION of THE EMPIRICAL VARIOGRAM + # here we order the data in ascending date order + day.o <- order(day) + coord1 <- coord1[day.o] + coord2 <- coord2[day.o] + obs <- obs[day.o] + id <- id[day.o] + forecast <- forecast[day.o] + # the first step in the gop software is to calculate the residuals + gop.mod <- lm(obs~forecast) + gop.res <- gop.mod$res + gop.var <- var(gop.res) + # the second step is to determine the empirical variograms: if the vector with the cutpoints + # is not specified, we determine the cutpoints by looking at the day with the median number + # of observations, we calculate the cutpoints so that the number + # of bins is equal to the one specified and each bin contains approx the same number of pairs. If the vector with the + # cutpoints is specified, then we just use that vector of cutpoints. + + if(length(cut.points)!=0 & length(max.dist)!=0){ + cut.points <- cut.points[cut.points <= max.dist] + } + + if(length(cut.points)==0){ + # all this part is to determine the day with the median number of observations + obs.day <- table(day) + obs.day.vec <- matrix(obs.day,nrow=length(unique(day)),ncol=1) + obs.median <- round(apply(obs.day.vec,2,median),0) + diff.median <- abs(obs.day-obs.median) + unique.day <- unique(day) + day.index <- min(unique.day[diff.median==min(diff.median)]) + # this part is to determine the distance among stations that have been recorded in the + # day with the median number of observations and to determine the cutpoints. + obs.day.index <- obs[day==day.index] + id.day.index <- id[day==day.index] + coord1.day.index <- coord1[day==day.index] + coord2.day.index <- coord2[day==day.index] + dist.day.index <- calc.dist(coord1.day.index,coord2.day.index,id.day.index) + dist.day.index <- dist.day.index[lower.tri(dist.day.index)] + dist.day.index <- dist.day.index[dist.day.index > 0] + ord.dist.day <- sort(dist.day.index) + + # this is in case we want to consider only those pairs of stations where the distance + # is smaller than the maximum distance allowed among locations. + if(length(max.dist)!= 0){ + ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] + } + + if(length(max.dist)==0){ + max.dist <- quantile(calc.dist(ord.dist.day),.9) + ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] + } + + l.dist.day <- length(ord.dist.day) + l.bins <- floor(l.dist.day/nbins) + for(i in 1:(nbins-1)){ + cut.points[i] <- ord.dist.day[i*l.bins] + } + } + + # this part is to calculate the empirical variogram + unique.day <- unique(day) + l.day <- length(unique.day) + n.stations <- length(unique(id)) + n.cuts <- length(cut.points) + W <- rep(0,(n.cuts-1)) + W.obs <- rep(0,(n.cuts-1)) + for(i in 1:l.day){ + distance.day <- NULL + difference.day <- NULL + new.dist.day <- NULL + s.new.dist.day <- NULL + o.new.dist.day <- NULL + new.diff.day <- NULL + s.new.diff.day <- NULL + gop.res.day <- NULL + id.day <- NULL + coord1.day <- NULL + coord2.day <- NULL + gop.res.day <- gop.res[day==unique.day[i]] + id.day <- id[day==unique.day[i]] + coord1.day <- coord1[day==unique.day[i]] + coord2.day <- coord2[day==unique.day[i]] + distance.day <- calc.dist(coord1.day,coord2.day,id.day) + new.dist.day <- distance.day[lower.tri(distance.day)] + s.new.dist.day <- sort(new.dist.day) + o.new.dist.day <- order(new.dist.day) + difference.day <- calc.difference(gop.res.day) + new.diff.day <- difference.day[lower.tri(difference.day)] + s.new.diff.day <- new.diff.day[o.new.dist.day] + W.new <- rep(0,(n.cuts-1)) + W.obs.new <- rep(0,(n.cuts-1)) + + for(s in 1:(n.cuts-1)){ + low.bound <- cut.points[s] + upp.bound <- cut.points[s+1] + v.dist <- NULL + v.dist1 <- NULL + v.diff1 <- NULL + index.v.dist <- NULL + index.v.dist1 <- NULL + + if(s < (n.cuts-1)){ + v.dist <- s.new.dist.day[s.new.dist.day >= low.bound & s.new.dist.day < upp.bound] + index.v.dist <- seq(1:length(s.new.dist.day))[s.new.dist.day >= low.bound & s.new.dist.day < upp.bound] + v.dist1 <- v.dist[v.dist!=0] + index.v.dist1 <- index.v.dist[v.dist!=0] + + if(length(v.dist1) >=1){ + v.diff1 <- s.new.diff.day[index.v.dist1] + W.new[s] <- sum(v.diff1) + W.obs.new[s] <- length(v.dist1)} + if(length(v.dist1) ==0){ + W.new[s] <- 0 + W.obs.new[s] <- 0} + } + if(s==(n.cuts-1)){ + v.dist <- s.new.dist.day[s.new.dist.day >= low.bound & s.new.dist.day <= upp.bound] + index.v.dist <- seq(1:length(s.new.dist.day))[s.new.dist.day >= low.bound & s.new.dist.day <= upp.bound] + v.dist1 <- v.dist[v.dist!=0] + index.v.dist1 <- index.v.dist[v.dist!=0] + + if(length(v.dist1) >=1){ + v.diff1 <- s.new.diff.day[index.v.dist1] + W.new[s] <- sum(v.diff1) + W.obs.new[s] <- length(v.dist1)} + if(length(v.dist1) ==0){ + W.new[s] <- 0 + W.obs.new[s] <- 0} + } + } + W <- W+W.new + W.obs <- W.obs + W.obs.new + } + avg.variog <- round(W/(2*W.obs),2) + n.h <- W.obs + x.vec <- NULL + for(i in 1:(n.cuts-1)){ + x.vec[i] <- (cut.points[i]+cut.points[i+1])/2 + } + fin.avg.variog <- c(0,avg.variog) + fin.x.vec <- c(0,x.vec) + fin.n.h <- c(n.stations,n.h) + B <- list(mar.var=round(gop.var,3),bin.midpoints=fin.x.vec,number.pairs=fin.n.h,empir.variog=fin.avg.variog) + return(B) + } function (day, obs, forecast, id, coord1, coord2, cut.points = NULL, max.dist = NULL, nbins = 300) { if (missing(cut.points)) cut.points <- NULL if (missing(max.dist)) max.dist <- NULL if (missing(nbins)) nbins <- NULL l.day <- length(day) l.obs <- length(obs) l.for <- length(forecast) l.id <- length(id) l.coord1 <- length(coord1) l.coord2 <- length(coord2) if (sum((c(l.day, l.obs, l.for, l.id, l.coord1, l.coord2)/l.day) == rep(1, 6)) != 6) { stop("Error in the dimension of the data") } if (sum(is.numeric(obs) == rep("TRUE", l.obs)) < l.obs) { stop("observations should be a numeric vector") } if (sum(ceiling(day) == day) < l.day) { stop("day should be an integer") } if (sum(is.numeric(coord1) == rep("TRUE", l.coord1)) < l.coord1 | sum(is.numeric(coord2) == rep("TRUE", l.coord2)) < l.coord2) { stop("coord1 and coord2 should be numeric vectors") } l.cuts <- length(cut.points) if (l.cuts == 1) { stop("cut.points should be a numeric vector") } if (l.cuts >= 2 & (sum(is.numeric(cut.points) == rep("TRUE", l.cuts)) < l.cuts)) { stop("cut.points should be a numeric vector") } if (l.cuts >= 2 & (sum(is.numeric(cut.points) == rep("TRUE", l.cuts)) == l.cuts)) { if (sum(order(cut.points) == seq(1:l.cuts)) < l.cuts) { stop("Cut points should be in increasing order") } if (sum(cut.points >= 0) < l.cuts) { stop("Cut points should be non-negative numbers") } if (length(cut.points) != length(unique(cut.points))) { stop("The vector with cut points should not contain repeated entries") } } l.mdist <- length(max.dist) if (l.mdist > 1) { print("max.dist is a numeric field, not a vector") } if (l.mdist == 1) { if (is.numeric(max.dist) == FALSE) { stop("max.dist is a numeric field") } if (max.dist < 0) { stop("max.dist should be a positive number") } } l.nbins <- length(nbins) if (l.nbins == 0 & l.cuts == 0) { nbins <- 300 } if (l.nbins == 1 & l.cuts >= 2) { nbins = NULL } l.nbins <- length(nbins) if (l.nbins > 1) { stop("nbins should be an integer: not a vector") } if (l.nbins == 1) { if (ceiling(nbins) != nbins) { stop("Invalid input: the number of bins should be a positive integer") } if (ceiling(nbins) == nbins & nbins < 0) { stop("Invalid input: the number of bins should be a positive integer") } } day.o <- order(day) coord1 <- coord1[day.o] coord2 <- coord2[day.o] obs <- obs[day.o] id <- id[day.o] forecast <- forecast[day.o] gop.mod <- lm(obs ~ forecast) gop.res <- gop.mod$res gop.var <- var(gop.res) if (length(cut.points) != 0 & length(max.dist) != 0) { cut.points <- cut.points[cut.points <= max.dist] } if (length(cut.points) == 0) { obs.day <- table(day) obs.day.vec <- matrix(obs.day, nrow = length(unique(day)), ncol = 1) obs.median <- round(apply(obs.day.vec, 2, median), 0) diff.median <- abs(obs.day - obs.median) unique.day <- unique(day) day.index <- min(unique.day[diff.median == min(diff.median)]) obs.day.index <- obs[day == day.index] id.day.index <- id[day == day.index] coord1.day.index <- coord1[day == day.index] coord2.day.index <- coord2[day == day.index] dist.day.index <- calc.dist(coord1.day.index, coord2.day.index, id.day.index) dist.day.index <- dist.day.index[lower.tri(dist.day.index)] dist.day.index <- dist.day.index[dist.day.index > 0] ord.dist.day <- sort(dist.day.index) if (length(max.dist) != 0) { ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] } if (length(max.dist) == 0) { max.dist <- quantile(calc.dist(ord.dist.day), 0.9) ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] } l.dist.day <- length(ord.dist.day) l.bins <- floor(l.dist.day/nbins) for (i in 1:(nbins - 1)) { cut.points[i] <- ord.dist.day[i * l.bins] } } unique.day <- unique(day) l.day <- length(unique.day) n.stations <- length(unique(id)) n.cuts <- length(cut.points) W <- rep(0, (n.cuts - 1)) W.obs <- rep(0, (n.cuts - 1)) for (i in 1:l.day) { distance.day <- NULL difference.day <- NULL new.dist.day <- NULL s.new.dist.day <- NULL o.new.dist.day <- NULL new.diff.day <- NULL s.new.diff.day <- NULL gop.res.day <- NULL id.day <- NULL coord1.day <- NULL coord2.day <- NULL gop.res.day <- gop.res[day == unique.day[i]] id.day <- id[day == unique.day[i]] coord1.day <- coord1[day == unique.day[i]] coord2.day <- coord2[day == unique.day[i]] distance.day <- calc.dist(coord1.day, coord2.day, id.day) new.dist.day <- distance.day[lower.tri(distance.day)] s.new.dist.day <- sort(new.dist.day) o.new.dist.day <- order(new.dist.day) difference.day <- calc.difference(gop.res.day) new.diff.day <- difference.day[lower.tri(difference.day)] s.new.diff.day <- new.diff.day[o.new.dist.day] W.new <- rep(0, (n.cuts - 1)) W.obs.new <- rep(0, (n.cuts - 1)) for (s in 1:(n.cuts - 1)) { low.bound <- cut.points[s] upp.bound <- cut.points[s + 1] v.dist <- NULL v.dist1 <- NULL v.diff1 <- NULL index.v.dist <- NULL index.v.dist1 <- NULL if (s < (n.cuts - 1)) { v.dist <- s.new.dist.day[s.new.dist.day >= low.bound & s.new.dist.day < upp.bound] index.v.dist <- seq(1:length(s.new.dist.day))[s.new.dist.day >= low.bound & s.new.dist.day < upp.bound] v.dist1 <- v.dist[v.dist != 0] index.v.dist1 <- index.v.dist[v.dist != 0] if (length(v.dist1) >= 1) { v.diff1 <- s.new.diff.day[index.v.dist1] W.new[s] <- sum(v.diff1) W.obs.new[s] <- length(v.dist1) } if (length(v.dist1) == 0) { W.new[s] <- 0 W.obs.new[s] <- 0 } } if (s == (n.cuts - 1)) { v.dist <- s.new.dist.day[s.new.dist.day >= low.bound & s.new.dist.day <= upp.bound] index.v.dist <- seq(1:length(s.new.dist.day))[s.new.dist.day >= low.bound & s.new.dist.day <= upp.bound] v.dist1 <- v.dist[v.dist != 0] index.v.dist1 <- index.v.dist[v.dist != 0] if (length(v.dist1) >= 1) { v.diff1 <- s.new.diff.day[index.v.dist1] W.new[s] <- sum(v.diff1) W.obs.new[s] <- length(v.dist1) } if (length(v.dist1) == 0) { W.new[s] <- 0 W.obs.new[s] <- 0 } } } W <- W + W.new W.obs <- W.obs + W.obs.new } avg.variog <- round(W/(2 * W.obs), 2) n.h <- W.obs x.vec <- NULL for (i in 1:(n.cuts - 1)) { x.vec[i] <- (cut.points[i] + cut.points[i + 1])/2 } fin.avg.variog <- c(0, avg.variog) fin.x.vec <- c(0, x.vec) fin.n.h <- c(n.stations, n.h) B <- list(mar.var = round(gop.var, 3), bin.midpoints = fin.x.vec, number.pairs = fin.n.h, empir.variog = fin.avg.variog) return(B) } > ## End Don't show > > > > cleanEx(); ..nameEx <- "EmpDir.variog" > > ### * EmpDir.variog > > flush(stderr()); flush(stdout()) > > ### Name: EmpDir.variog > ### Title: Directional empirical variogram of forecast errors averaged over > ### time > ### Aliases: EmpDir.variog > ### Keywords: file > > ### ** Examples > > ## Loading data > data(slp) > day <- slp$date.obs > id <- slp$id.stat > coord1 <- slp$lon.stat > coord2 <- slp$lat.stat > obs <- slp$obs > forecast <- slp$forecast > > ## Computing directional variogram > ## No specified cutpoints, no specified maximum distance > ## No specified tolerance angles and default number of bins > dir.variog <- EmpDir.variog(day,obs,forecast,id,coord1,coord2, + tol.angle1=NULL,tol.angle2=NULL,cut.points=NULL,max.dist=NULL, + nbins=NULL,type='E') > ## Plotting directional variogram > plot(dir.variog$bin.midpoints,dir.variog$dir.variog,xlab="Distance", + ylab="Semi-variance",main="Empirical Directional variogram") > > ## Computing directional variogram > ## Specified cutpoints, specified maximum distance > ## Specified tolerance angles and unspecified number of bins > dir.variog <- + EmpDir.variog(day,obs,forecast,id,coord1,coord2,tol.angle1=30, + tol.angle2=150,cut.points=seq(0,1000,by=5),max.dist=800,nbins=NULL, + type='N') > ## Plotting directional variogram > plot(dir.variog$bin.midpoints,dir.variog$dir.variog,xlab="Distance",ylab="Semi-variance",main="Empirical Directional variogram") > > ## Don't show: ## The function is currently defined as > function(day,obs,forecast,id,coord1,coord2,tol.angle1=45,tol.angle2=135,cut.points=NULL,max.dist=NULL,nbins=300,type){ + + # default values + if(missing(cut.points)) + cut.points <- NULL + if(missing(max.dist)) + max.dist <- NULL + if(missing(nbins)) + nbins <- NULL + if(missing(type)){ + stop("Invalid input for type: a type should be provided") + } + if(missing(tol.angle1)) + tol.angle1 <- 45 + if(missing(tol.angle2)) + tol.angle2 <- 135 + + #INPUT CHECK + # Here we check if the input is right. + l.day <- length(day) + l.obs <- length(obs) + l.for <- length(forecast) + l.id <- length(id) + l.coord1 <- length(coord1) + l.coord2 <- length(coord2) + + if(sum((c(l.day,l.obs,l.for,l.id,l.coord1,l.coord2)/l.day)==rep(1,6))!=6){ + stop("Mismatch in dimensions in the data") + } + if(sum(is.numeric(obs)==rep("TRUE",l.obs)) 1){ + stop("tol.angle1 should be a number") + } + + if(l.tol.angle1==1 & is.numeric(tol.angle1)==FALSE){ + stop("tol.angle1 should be a number between 0 and 360") + } + + if(tol.angle1 <0 | tol.angle1 > 360){ + stop("tol.angle1 should be a number between 0 and 360") + } + + l.tol.angle2 <- length(tol.angle2) + + if(l.tol.angle2==0){ + tol.angle2 <- 135} + + if(l.tol.angle2 >1){ + stop("tol.angle2 should be a number") + } + + if(l.tol.angle2==1 & is.numeric(tol.angle2)==FALSE){ + stop("tol.angle1 should be a number between 0 and 360") + } + if(tol.angle2 <0 | tol.angle2 > 360){ + stop("tol.angle1 should be a number between 0 and 360") + } + + if(tol.angle2 <= tol.angle1){ + stop("tol.angle2 should be greater than tol.angle1") + } + ## here we check the cutpoints vector + l.cuts <- length(cut.points) + if(l.cuts==1){ + stop("cut.points should be a numeric vector") + } + + if(l.cuts>=2 & (sum(is.numeric(cut.points)==rep("TRUE",l.cuts))=2 & (sum(is.numeric(cut.points)==rep("TRUE",l.cuts))==l.cuts)){ + if(sum(order(cut.points)==seq(1:l.cuts))= 0) < l.cuts){ + stop("Cut points should be non-negative numbers") + } + if(length(cut.points)!=length(unique(cut.points))){ + stop("The vector with cut points should not contain repeated entries") + } + } + + + ## check on the max.dist + l.mdist <- length(max.dist) + if(l.mdist > 1){ + stop("max.dist is a numeric field, not a vector")} + if(l.mdist==1){ + if(is.numeric(max.dist)==FALSE){ + stop("max.dist is a numeric field") + } + if(max.dist < 0){ + stop("max.dist should be a positive number") + } + } + + ## check on the number of bins + l.nbins <- length(nbins) + if(l.nbins==0 & l.cuts==0){ + nbins <- 300 + } + if(l.nbins==1 & l.cuts >=2){ + nbins=NULL + } + l.nbins <- length(nbins) + if(l.nbins >1){ + stop("nbins should be an integer: not a vector") + } + if(l.nbins==1){ + if(floor(nbins)!=nbins){ + stop("Invalid input: the number of bins should be a positive integer") + } + if(ceiling(nbins)==nbins & nbins < 0){ + stop("Invalid input: the number of bins should be a positive integer") + } + } + + + ## check on the type ## + l.type <- length(type) + if(l.type==0){ + stop("type should be entered") + } + if(type!="E" & type!="N"){ + stop("type can only be equal to E or to N") + } + # ESTIMATION of THE EMPIRICAL VARIOGRAM + # here we order the data in ascending date order + day.o <- order(day) + coord1 <- coord1[day.o] + coord2 <- coord2[day.o] + obs <- obs[day.o] + id <- id[day.o] + forecast <- forecast[day.o] + # here we calculate the residuals + gop.mod <- lm(obs~forecast) + gop.res <- gop.mod$res + gop.var <- var(gop.res) + tol.angle.rad1 <- tol.angle1/57.29577951 + tol.angle.rad2 <- tol.angle2/57.29577951 + + # if the vector with the cutpoints + # is not specified, we determine the cutpoints by looking at the day with the median number + # of observations, we calculate the cutpoints so that the number + # of bins is equal to the one specified and each bin contains approx the same number of pairs. If the vector with the + # if the vector of cutpoints is specified, then we just use that vector of cutpoints. + + if(length(cut.points)!=0 & length(max.dist)!=0){ + cut.points <- cut.points[cut.points <= max.dist] + } + if(length(cut.points)==0){ + # all this part is to determine the day with the median number of observations + obs.day <- table(day) + obs.day.vec <- matrix(obs.day,nrow=length(unique(day)),ncol=1) + obs.median <- round(apply(obs.day.vec,2,median),0) + diff.median <- abs(obs.day-obs.median) + unique.day <- unique(day) + day.index <- min(unique.day[diff.median==min(diff.median)]) + + # this part is to determine the directional distance among stations that have been recorded in the + # day with the median number of observations and to determine the cutpoints. + obs.day.index <- obs[day==day.index] + id.day.index <- id[day==day.index] + coord1.day.index <- coord1[day==day.index] + coord2.day.index <- coord2[day==day.index] + dist.day.index <- calc.dist.dir(coord1.day.index,coord2.day.index,id.day.index,tol.angle.rad1,tol.angle.rad2,type) + dist.day.index <- dist.day.index[lower.tri(dist.day.index)] + dist.day.index <- dist.day.index[dist.day.index > 0] + ord.dist.day <- sort(dist.day.index) + + # this is in case we want to consider only those pairs of stations where the distance + # is smaller than the maximum distance allowed among locations. + if(length(max.dist)!= 0){ + ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] + } + + if(length(max.dist)==0){ + max.dist <- quantile(calc.dist(ord.dist.day),.9) + ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] + } + + l.dist.day <- length(ord.dist.day) + l.bins <- floor(l.dist.day/nbins) + for(i in 1:(nbins-1)){ + cut.points[i] <- ord.dist.day[i*l.bins] + } + } + + # this part is to calculate the empirical directional variogram + unique.day <- unique(day) + l.day <- length(unique.day) + n.stations <- length(unique(id)) + n.cuts <- length(cut.points) + W <- rep(0,(n.cuts-1)) + W.obs <- rep(0,(n.cuts-1)) + + for(i in 1:l.day){ + distance.day <- NULL + difference.day <- NULL + new.dist.day <- NULL + s.new.dist.day <- NULL + o.new.dist.day <- NULL + new.diff.day <- NULL + new.diff.day <- NULL + s.new.diff.day <- NULL + gop.res.day <- NULL + id.day <- NULL + coord1.day <- NULL + coord2.day <- NULL + gop.res.day <- gop.res[day==unique.day[i]] + id.day <- id[day==unique.day[i]] + coord1.day <- coord1[day==unique.day[i]] + coord2.day <- coord2[day==unique.day[i]] + distance.day <- + calc.dist.dir(coord1.day,coord2.day,id.day,tol.angle.rad1,tol.angle.rad2,type) + new.dist.day <- distance.day[lower.tri(distance.day)] + s.new.dist.day <- sort(new.dist.day) + o.new.dist.day <- order(new.dist.day) + difference.day <- calc.difference(gop.res.day) + new.diff.day <- difference.day[lower.tri(difference.day)] + s.new.diff.day <- new.diff.day[o.new.dist.day] + W.new <- rep(0,(n.cuts-1)) + W.obs.new <- rep(0,(n.cuts-1)) + + for(s in 1:(n.cuts-1)){ + low.bound <- cut.points[s] + upp.bound <- cut.points[s+1] + v.dist <- NULL + v.dist1 <- NULL + v.diff1 <- NULL + index.v.dist <- NULL + index.v.dist1 <- NULL + + if(s < (n.cuts-1)){ + v.dist <- s.new.dist.day[s.new.dist.day >= low.bound & s.new.dist.day < upp.bound] + index.v.dist <- seq(1:length(s.new.dist.day))[s.new.dist.day >= low.bound & s.new.dist.day < upp.bound] + v.dist1 <- v.dist[v.dist!=0] + index.v.dist1 <- index.v.dist[v.dist!=0] + + if(length(v.dist1) >=1){ + v.diff1 <- s.new.diff.day[index.v.dist1] + W.new[s] <- sum(v.diff1) + W.obs.new[s] <- length(v.dist1)} + if(length(v.dist1) ==0){ + W.new[s] <- 0 + W.obs.new[s] <- 0} + } + if(s==(n.cuts-1)){ + v.dist <- s.new.dist.day[s.new.dist.day >= low.bound & s.new.dist.day <= upp.bound] + index.v.dist <- seq(1:length(s.new.dist.day))[s.new.dist.day >= low.bound & s.new.dist.day <= upp.bound] + v.dist1 <- v.dist[v.dist!=0] + index.v.dist1 <- index.v.dist[v.dist!=0] + + if(length(v.dist1) >=1){ + v.diff1 <- s.new.diff.day[index.v.dist1] + W.new[s] <- sum(v.diff1) + W.obs.new[s] <- length(v.dist1)} + if(length(v.dist1) ==0){ + W.new[s] <- 0 + W.obs.new[s] <- 0} + } + if(s==(n.cuts-1)){ + v.dist <- s.new.dist.day[s.new.dist.day >= low.bound & s.new.dist.day <= upp.bound] + index.v.dist <- seq(1:length(s.new.dist.day))[s.new.dist.day >= low.bound & s.new.dist.day <= upp.bound] + v.dist1 <- v.dist[v.dist!=0] + index.v.dist1 <- index.v.dist[v.dist!=0] + + if(length(v.dist1) >=1){ + v.diff1 <- s.new.diff.day[index.v.dist1] + W.new[s] <- sum(v.diff1) + W.obs.new[s] <- length(v.dist1)} + if(length(v.dist1) ==0){ + W.new[s] <- 0 + W.obs.new[s] <- 0} + } + } + W <- W+W.new + W.obs <- W.obs + W.obs.new + } + avg.variog <- round(W/(2*W.obs),2) + n.h <- W.obs + x.vec <- NULL + for(i in 1:(n.cuts-1)){ + x.vec[i] <- (cut.points[i]+cut.points[i+1])/2 + } + fin.avg.variog <- c(0,avg.variog) + fin.x.vec <- c(0,x.vec) + fin.n.h <- c(n.stations,n.h) + B <- list(bin.midpoints=fin.x.vec,number.pairs=fin.n.h,dir.variog=fin.avg.variog) + return(B) + } function (day, obs, forecast, id, coord1, coord2, tol.angle1 = 45, tol.angle2 = 135, cut.points = NULL, max.dist = NULL, nbins = 300, type) { if (missing(cut.points)) cut.points <- NULL if (missing(max.dist)) max.dist <- NULL if (missing(nbins)) nbins <- NULL if (missing(type)) { stop("Invalid input for type: a type should be provided") } if (missing(tol.angle1)) tol.angle1 <- 45 if (missing(tol.angle2)) tol.angle2 <- 135 l.day <- length(day) l.obs <- length(obs) l.for <- length(forecast) l.id <- length(id) l.coord1 <- length(coord1) l.coord2 <- length(coord2) if (sum((c(l.day, l.obs, l.for, l.id, l.coord1, l.coord2)/l.day) == rep(1, 6)) != 6) { stop("Mismatch in dimensions in the data") } if (sum(is.numeric(obs) == rep("TRUE", l.obs)) < l.obs) { stop("obs should be a numeric vector") } if (sum(ceiling(day) == day) < l.day) { stop("day should be a vector containing integers") } if (sum(is.numeric(coord1) == rep("TRUE", l.coord1)) < l.coord1 | sum(is.numeric(coord2) == rep("TRUE", l.coord2)) < l.coord2) { stop("coord1 and coord2 should be numeric vectors") } l.tol.angle1 <- length(tol.angle1) if (l.tol.angle1 == 0) { tol.angle1 <- 45 } if (l.tol.angle1 > 1) { stop("tol.angle1 should be a number") } if (l.tol.angle1 == 1 & is.numeric(tol.angle1) == FALSE) { stop("tol.angle1 should be a number between 0 and 360") } if (tol.angle1 < 0 | tol.angle1 > 360) { stop("tol.angle1 should be a number between 0 and 360") } l.tol.angle2 <- length(tol.angle2) if (l.tol.angle2 == 0) { tol.angle2 <- 135 } if (l.tol.angle2 > 1) { stop("tol.angle2 should be a number") } if (l.tol.angle2 == 1 & is.numeric(tol.angle2) == FALSE) { stop("tol.angle1 should be a number between 0 and 360") } if (tol.angle2 < 0 | tol.angle2 > 360) { stop("tol.angle1 should be a number between 0 and 360") } if (tol.angle2 <= tol.angle1) { stop("tol.angle2 should be greater than tol.angle1") } l.cuts <- length(cut.points) if (l.cuts == 1) { stop("cut.points should be a numeric vector") } if (l.cuts >= 2 & (sum(is.numeric(cut.points) == rep("TRUE", l.cuts)) < l.cuts)) { stop("cut.points should be a numeric vector") } if (l.cuts >= 2 & (sum(is.numeric(cut.points) == rep("TRUE", l.cuts)) == l.cuts)) { if (sum(order(cut.points) == seq(1:l.cuts)) < l.cuts) { stop("Cut points should be in increasing order") } if (sum(cut.points >= 0) < l.cuts) { stop("Cut points should be non-negative numbers") } if (length(cut.points) != length(unique(cut.points))) { stop("The vector with cut points should not contain repeated entries") } } l.mdist <- length(max.dist) if (l.mdist > 1) { stop("max.dist is a numeric field, not a vector") } if (l.mdist == 1) { if (is.numeric(max.dist) == FALSE) { stop("max.dist is a numeric field") } if (max.dist < 0) { stop("max.dist should be a positive number") } } l.nbins <- length(nbins) if (l.nbins == 0 & l.cuts == 0) { nbins <- 300 } if (l.nbins == 1 & l.cuts >= 2) { nbins = NULL } l.nbins <- length(nbins) if (l.nbins > 1) { stop("nbins should be an integer: not a vector") } if (l.nbins == 1) { if (floor(nbins) != nbins) { stop("Invalid input: the number of bins should be a positive integer") } if (ceiling(nbins) == nbins & nbins < 0) { stop("Invalid input: the number of bins should be a positive integer") } } l.type <- length(type) if (l.type == 0) { stop("type should be entered") } if (type != "E" & type != "N") { stop("type can only be equal to E or to N") } day.o <- order(day) coord1 <- coord1[day.o] coord2 <- coord2[day.o] obs <- obs[day.o] id <- id[day.o] forecast <- forecast[day.o] gop.mod <- lm(obs ~ forecast) gop.res <- gop.mod$res gop.var <- var(gop.res) tol.angle.rad1 <- tol.angle1/57.29577951 tol.angle.rad2 <- tol.angle2/57.29577951 if (length(cut.points) != 0 & length(max.dist) != 0) { cut.points <- cut.points[cut.points <= max.dist] } if (length(cut.points) == 0) { obs.day <- table(day) obs.day.vec <- matrix(obs.day, nrow = length(unique(day)), ncol = 1) obs.median <- round(apply(obs.day.vec, 2, median), 0) diff.median <- abs(obs.day - obs.median) unique.day <- unique(day) day.index <- min(unique.day[diff.median == min(diff.median)]) obs.day.index <- obs[day == day.index] id.day.index <- id[day == day.index] coord1.day.index <- coord1[day == day.index] coord2.day.index <- coord2[day == day.index] dist.day.index <- calc.dist.dir(coord1.day.index, coord2.day.index, id.day.index, tol.angle.rad1, tol.angle.rad2, type) dist.day.index <- dist.day.index[lower.tri(dist.day.index)] dist.day.index <- dist.day.index[dist.day.index > 0] ord.dist.day <- sort(dist.day.index) if (length(max.dist) != 0) { ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] } if (length(max.dist) == 0) { max.dist <- quantile(calc.dist(ord.dist.day), 0.9) ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] } l.dist.day <- length(ord.dist.day) l.bins <- floor(l.dist.day/nbins) for (i in 1:(nbins - 1)) { cut.points[i] <- ord.dist.day[i * l.bins] } } unique.day <- unique(day) l.day <- length(unique.day) n.stations <- length(unique(id)) n.cuts <- length(cut.points) W <- rep(0, (n.cuts - 1)) W.obs <- rep(0, (n.cuts - 1)) for (i in 1:l.day) { distance.day <- NULL difference.day <- NULL new.dist.day <- NULL s.new.dist.day <- NULL o.new.dist.day <- NULL new.diff.day <- NULL new.diff.day <- NULL s.new.diff.day <- NULL gop.res.day <- NULL id.day <- NULL coord1.day <- NULL coord2.day <- NULL gop.res.day <- gop.res[day == unique.day[i]] id.day <- id[day == unique.day[i]] coord1.day <- coord1[day == unique.day[i]] coord2.day <- coord2[day == unique.day[i]] distance.day <- calc.dist.dir(coord1.day, coord2.day, id.day, tol.angle.rad1, tol.angle.rad2, type) new.dist.day <- distance.day[lower.tri(distance.day)] s.new.dist.day <- sort(new.dist.day) o.new.dist.day <- order(new.dist.day) difference.day <- calc.difference(gop.res.day) new.diff.day <- difference.day[lower.tri(difference.day)] s.new.diff.day <- new.diff.day[o.new.dist.day] W.new <- rep(0, (n.cuts - 1)) W.obs.new <- rep(0, (n.cuts - 1)) for (s in 1:(n.cuts - 1)) { low.bound <- cut.points[s] upp.bound <- cut.points[s + 1] v.dist <- NULL v.dist1 <- NULL v.diff1 <- NULL index.v.dist <- NULL index.v.dist1 <- NULL if (s < (n.cuts - 1)) { v.dist <- s.new.dist.day[s.new.dist.day >= low.bound & s.new.dist.day < upp.bound] index.v.dist <- seq(1:length(s.new.dist.day))[s.new.dist.day >= low.bound & s.new.dist.day < upp.bound] v.dist1 <- v.dist[v.dist != 0] index.v.dist1 <- index.v.dist[v.dist != 0] if (length(v.dist1) >= 1) { v.diff1 <- s.new.diff.day[index.v.dist1] W.new[s] <- sum(v.diff1) W.obs.new[s] <- length(v.dist1) } if (length(v.dist1) == 0) { W.new[s] <- 0 W.obs.new[s] <- 0 } } if (s == (n.cuts - 1)) { v.dist <- s.new.dist.day[s.new.dist.day >= low.bound & s.new.dist.day <= upp.bound] index.v.dist <- seq(1:length(s.new.dist.day))[s.new.dist.day >= low.bound & s.new.dist.day <= upp.bound] v.dist1 <- v.dist[v.dist != 0] index.v.dist1 <- index.v.dist[v.dist != 0] if (length(v.dist1) >= 1) { v.diff1 <- s.new.diff.day[index.v.dist1] W.new[s] <- sum(v.diff1) W.obs.new[s] <- length(v.dist1) } if (length(v.dist1) == 0) { W.new[s] <- 0 W.obs.new[s] <- 0 } } if (s == (n.cuts - 1)) { v.dist <- s.new.dist.day[s.new.dist.day >= low.bound & s.new.dist.day <= upp.bound] index.v.dist <- seq(1:length(s.new.dist.day))[s.new.dist.day >= low.bound & s.new.dist.day <= upp.bound] v.dist1 <- v.dist[v.dist != 0] index.v.dist1 <- index.v.dist[v.dist != 0] if (length(v.dist1) >= 1) { v.diff1 <- s.new.diff.day[index.v.dist1] W.new[s] <- sum(v.diff1) W.obs.new[s] <- length(v.dist1) } if (length(v.dist1) == 0) { W.new[s] <- 0 W.obs.new[s] <- 0 } } } W <- W + W.new W.obs <- W.obs + W.obs.new } avg.variog <- round(W/(2 * W.obs), 2) n.h <- W.obs x.vec <- NULL for (i in 1:(n.cuts - 1)) { x.vec[i] <- (cut.points[i] + cut.points[i + 1])/2 } fin.avg.variog <- c(0, avg.variog) fin.x.vec <- c(0, x.vec) fin.n.h <- c(n.stations, n.h) B <- list(bin.midpoints = fin.x.vec, number.pairs = fin.n.h, dir.variog = fin.avg.variog) return(B) } > ## End Don't show > > > > cleanEx(); ..nameEx <- "Field.sim" > > ### * Field.sim > > flush(stderr()); flush(stdout()) > > ### Name: Field.sim > ### Title: Simulation of weather random field > ### Aliases: Field.sim > ### Keywords: file > > ### ** Examples > > ## Loading data > library(fields) > library(RandomFields) > data(slp) > data(gridlong) > data(gridlat) > data(forecast.grid) > day <- slp$date.obs > id <- slp$id.stat > coord1 <- slp$lon.stat > coord2 <- slp$lat.stat > obs <- slp$obs > forecast <- slp$forecast > coord1.grid <- gridlong$gridded.lon > coord2.grid <- gridlat$gridded.lat > forecast.grid <- forecast.grid$gridded.forecast > > ## Computing the empirical variogram > variogram <- Emp.variog(day=day,obs=obs,forecast=forecast,id=id, + coord1=coord1,coord2=coord2,cut.points=NULL,max.dist=NULL,nbins=NULL) > ## Estimating parameters > param.est <- Variog.fit(emp.variog=variogram,variog.model="exponential", + max.dist.fit=NULL,init.val=NULL,fix.nugget=FALSE) > > ## Simulating realizations of the weather random field > simul <- + Field.sim(obs=obs,forecast=forecast,coord1.grid=coord1.grid, + coord2.grid=coord2.grid,forecast.grid=forecast.grid,variog.model="exponential", + param.est=c(param.est$nugget,param.est$variance,param.est$range),n.sim=4, + n.displ=4,qt.displ=c(10,50,90)) Hit to see next plot: Hit to see next plot: ## A different > example function (topic, package = NULL, lib.loc = NULL, local = FALSE, echo = TRUE, verbose = getOption("verbose"), setRNG = FALSE, prompt.echo = paste(abbreviate(topic, 6), "> ", sep = "")) { topic <- substitute(topic) if (!is.character(topic)) topic <- deparse(topic)[1] INDICES <- .find.package(package, lib.loc, verbose = verbose) file <- index.search(topic, INDICES, "AnIndex", "R-ex") if (file == "") { warning(gettextf("no help file found for '%s'", topic), domain = NA) return(invisible()) } packagePath <- dirname(dirname(file)) if (length(file) > 1) { packagePath <- packagePath[1] warning(gettextf("more than one help file found: using package '%s'", basename(packagePath)), domain = NA) file <- file[1] } pkg <- basename(packagePath) lib <- dirname(packagePath) zfile <- zip.file.extract(file, "Rex.zip") if (zfile != file) on.exit(unlink(zfile)) if (!file.exists(zfile)) { warning(gettextf("'%s' has a help file but no examples file", topic), domain = NA) return(invisible()) } if (pkg != "base") library(pkg, lib = lib, character.only = TRUE) if (!is.logical(setRNG) || setRNG) { if ((exists(".Random.seed", envir = .GlobalEnv))) { oldSeed <- get(".Random.seed", envir = .GlobalEnv) on.exit(assign(".Random.seed", oldSeed, envir = .GlobalEnv)) } else { oldRNG <- RNGkind() on.exit(RNGkind(oldRNG[1], oldRNG[2])) } if (is.logical(setRNG)) { RNGkind("default", "default") set.seed(1) } else eval(setRNG) } encoding <- if (length(enc <- localeToCharset()) > 1) c(enc[-length(enc)], "latin1") else "" zz <- readLines(zfile, n = 1) if (length(grep("^### Encoding: ", zz)) > 0 && !identical(Sys.getlocale("LC_CTYPE"), "C")) encoding <- substring(zz, 15) source(zfile, local, echo = echo, prompt.echo = prompt.echo, verbose = verbose, max.deparse.length = 250, encoding = encoding) } > simul <- + Field.sim(obs=obs,forecast=forecast,coord1.grid=coord1.grid, + coord2.grid=coord2.grid,forecast.grid=forecast.grid,variog.model="gencauchy", + param.est=c(0,1,300,0.5,1.5),n.sim=4,n.displ=1,qt.displ=c(5,10,90,95)) Hit to see next plot: > ## Don't show: > ## The function is currently defined as > function(obs,forecast,coord1.grid,coord2.grid,forecast.grid,variog.model="exponential",param.est,n.sim=99,n.displ=4,qt.displ=c(10,50,90)){ + + ## Input check + # set to default values + if(missing(variog.model)) + variog.model <- "exponential" + if(missing(n.sim)) + n.sim <- 99 + if(missing(n.displ)) + n.displ <- 4 + if(missing(qt.displ)) + qt.displ <- c(10,50,90) + # check on the observations and the forecast + l.obs <- length(obs) + l.for <- length(forecast) + if(l.obs==0){ + stop("obs is a vector") + } + if(l.for==0){ + stop("forecast is a vector") + } + if(sum((c(l.obs,l.for)/l.obs)==rep(1,2))!=2){ + stop("Mismatch in dimensions in the data") + } + if(sum(is.numeric(obs)==rep("TRUE",l.obs)) 5){ + stop("Invalid input for param.est") + } + + if(l.parest== 3){ + if(sum(is.numeric(param.est)==rep(TRUE,l.parest))<3){ + stop("param.est is a numeric vector") + } + if(variog.model=="whittlematern"){ + stop("param.est should be of length 4") + } + if(variog.model=="gencauchy"){ + stop("param.est should be of length 5") + } + if(sum(param.est >= 0) < 3){ + stop("The entries of param.est should be non-negative numbers") + } + if(variog.model=="exponential"){ + case.parest <- 1 + } + if(variog.model=="spherical"){ + case.parest <- 1 + } + if(variog.model=="gauss"){ + case.parest <- 1 + } + } + + + if(l.parest==4){ + if(sum(is.numeric(param.est)==rep(TRUE,l.parest))<4){ + stop("param.est is a numeric vector") + } + if(variog.model=="gencauchy"){ + stop("param.est should be of length 5") + } + if(variog.model=="exponential"){ + stop("param.est should be of length 3") + } + if(variog.model=="gauss"){ + stop("param.est should be of length 3") + } + if(variog.model=="spherical"){ + stop("param.est should be of length 3") + } + if(sum(param.est >= 0) < 4){ + stop("The entries of param.est should be non-negative numbers") + } + if(variog.model=="whittlematern"){ + case.parest <- 1 + } + } + + if(l.parest==5){ + if(sum(is.numeric(param.est)==rep(TRUE,l.parest))<5){ + stop("param.est is a numeric vector") + } + if(variog.model=="whittlematern"){ + stop("param.est should be of length 4") + } + if(variog.model=="exponential"){ + stop("param.est should be of length 3") + } + if(variog.model=="gauss"){ + stop("param.est should be of length 3") + } + if(variog.model=="spherical"){ + stop("param.est should be of length 3") + } + if(sum(param.est >= 0) < 5){ + stop("The entries of param.est should be non-negative numbers") + } + if(variog.model=="gencauchy"){ + if(param.est[4] > 2){ + stop("a is out of the parameter space") + } + case.parest <- 1 + } + } + + if(case.parest==0){ + stop("Invalid input for param.est") + } + + ## check on the number of simulated random fields + + l.sim <- length(n.sim) + if(l.sim==0){ + n.sim <- 99} + l.sim <- length(n.sim) + + if(length(n.sim)> 1){ + stop("n.sim is a numeric field: not a vector") + } + + if(length(n.sim)==1){ + if(is.numeric(n.sim)==FALSE){ + stop("The number of simulated fields should be a positive integer") + } + if(ceiling(n.sim)!=n.sim){ + stop("The number of simulated fields should be a positive integer") + } + if(n.sim <= 0){ + stop("The number of simulated fields should be a positive integer") + } + } + + ## check on the n.displ + l.displ <- length(n.displ) + if(l.displ==0){ + n.displ <- 4} + + l.displ <- length(n.displ) + if(l.displ > 1){ + stop("n.displ is a numeric field: not a vector") + } + + if(l.displ==1){ + if(is.numeric(n.displ)!=TRUE){ + stop("n.displ is a numeric field") + } + if(ceiling(n.displ)!=n.displ){ + stop("n.displ should be an integer number") + } + if(n.displ < 0){ + stop("n.displ should be an integer non-negative number") + } + if(n.displ > n.sim){ + stop("n.displ should be less or equal to n.sim") + } + } + + ## check on the qt.displ + + l.qt <- length(qt.displ) + if(l.qt==0){ + qt.displ <- c(10,50,90)} + l.qt <- length(qt.displ) + + if(l.qt >0){ + if(sum(is.numeric(qt.displ)==rep("TRUE",l.qt)) 100){ + stop("The elements of qt.displ should be numbers between 0 and 100") + } + } + } + + simul <- sim.field(variog.model,param.est,coord1.grid,coord2.grid,n.sim) + + # here we add to the 0-mean simulated random fields the bias-corrected gridded forecasts. + + gop.mod <- lm(obs~forecast) + gop.coeff <- round(gop.mod$coeff,3) + + sim.out <- (gop.mod$coeff[1]+gop.mod$coeff[2]*forecast.grid)+ simul + + sim.out.1 <- array(0,c(65,65,n.sim)) + for(i in 1:n.sim){ + sim.out.1[,,i] <- engrid(coord1.grid,coord2.grid,sim.out[,i]) + } + + # here we determine the percentiles of the random fields + l.qtdispl <- length(qt.displ) + if(l.qtdispl==1 & qt.displ==0){ + quant.out <- 0 + qt.out.1 <- NULL} + else{ + quant.out <- matrix(0,ncol=l.qtdispl,nrow=l.coord1grid) + qt.displ <- qt.displ/100 + for(i in 1:l.qtdispl){ + quant.out[,i] <- (gop.mod$coeff[1]+gop.mod$coeff[2]*forecast.grid) + + (qnorm(qt.displ[i],0,1,TRUE,FALSE)*sqrt(param.est[1]+param.est[2])) + } + qt.out.1 <- array(0,c(65,65,l.qtdispl)) + for(i in 1:l.qtdispl){ + qt.out.1[,,i] <- engrid(coord1.grid,coord2.grid,quant.out[,i]) + } + } + + # here we return the output + if(variog.model=="whittlematern" | variog.model=="gencauchy"){ + output <- list(model=variog.model,nugget=param.est[1],variance=param.est[2], + range=param.est[3],additional.par=param.est[-seq(1:3)], + sim.fields=round(sim.out.1,4),pct.fields=round(qt.out.1,4))} + if(variog.model=="exponential" | variog.model=="gauss" | variog.model=="spherical"){ + output <- list(model=variog.model,nugget=param.est[1],variance=param.est[2], + range=param.est[3],sim.fields=round(sim.out.1,4),pct.fields=round(qt.out.1,4))} + lims <- c(min(sim.out.1[,,1:n.displ],na.rm=TRUE),max(sim.out.1[,,1:n.displ],na.rm=TRUE)) + + n.displ.4 <- ceiling(n.displ/4) + if(n.displ==1){ + x.lim <- c(min(coord1.grid,na.rm=TRUE),max(coord1.grid,na.rm=TRUE)) + y.lim <- c(min(coord2.grid,na.rm=TRUE),max(coord2.grid,na.rm=TRUE)) + par(ask=TRUE) + ens.plot(sim.out.1[,,1],lims,x.lim,y.lim,"Ensemble member 1") + } + + if(n.displ==2){ + par(mfrow=c(2,2),ask=TRUE) + plotens(coord1.grid,coord2.grid,sim.out.1[,,1:n.displ],n.displ,lims,"Ensemble member",0) + } + + if(n.displ >2){ + if(n.displ.4==1){ + par(mfrow=c(2,2),ask=TRUE) + plotens(coord1.grid,coord2.grid,sim.out.1[,,1:n.displ],n.displ,lims,"Ensemble member",0) + } + + if(n.displ.4 >1){ + for(i in 1:n.displ.4){ + n.pages <- i-1 + n.pages.4 <- 4*n.pages + if(i!= n.displ.4){ + par(mfrow=c(2,2),ask=TRUE) + first.col <- (((i-1)*4)+1) + last.col <- 4*i + plotens(coord1.grid,coord2.grid,sim.out.1[,,first.col:last.col],4,lims,"Ensemble member",n.pages.4)} + if(i==n.displ.4){ + par(mfrow=c(2,2),ask=TRUE) + first.col <- (4*(i-1)+1) + last.col <- n.displ + plotens(coord1.grid,coord2.grid,sim.out.1[,,first.col:last.col],((last.col-first.col)+1),lims,"Ensemble member",n.pages.4)} + } + } + } + + if(length(qt.out.1)!=0){ + lims <- c(min(qt.out.1[,,1:l.qtdispl],na.rm=TRUE),max(qt.out.1[,,1:l.qtdispl],na.rm=TRUE)) + l.qtdispl.4 <- ceiling(l.qtdispl/4) + + if(l.qtdispl==1){ + x.lim <- c(min(coord1.grid,na.rm=TRUE),max(coord1.grid,na.rm=TRUE)) + y.lim <- c(min(coord2.grid,na.rm=TRUE),max(coord2.grid,na.rm=TRUE)) + title <- paste(qt.displ*100,"-th Percentile") + par(ask=TRUE) + ens.plot(qt.out.1[,,1],lims,x.lim,y.lim,title) + } + + if(l.qtdispl==2){ + par(mfrow=c(2,2),ask=TRUE) + plotens.qt(coord1.grid,coord2.grid,qt.out.1[,,1:l.qtdispl],l.qtdispl,lims,qt.displ) + } + + if(l.qtdispl >2){ + if(l.qtdispl.4==1){ + par(mfrow=c(2,2),ask=TRUE) + plotens.qt(coord1.grid,coord2.grid,qt.out.1[,,1:l.qtdispl],l.qtdispl,lims,qt.displ) + } + + if(l.qtdispl.4 >1){ + for(i in 1:l.qtdispl.4){ + if(i!= l.qtdispl.4){ + par(mfrow=c(2,2),ask=TRUE) + first.col <- (((i-1)*4)+1) + last.col <- 4*i + plotens.qt(coord1.grid,coord2.grid,qt.out.1[,,first.col:last.col],4,lims,qt.displ[first.col:last.col])} + + if(i==l.qtdispl.4){ + par(mfrow=c(2,2),ask=TRUE) + first.col <- (4*(i-1)+1) + last.col <- l.qtdispl + plotens.qt(coord1.grid,coord2.grid,qt.out.1[,,first.col:last.col],(last.col-first.col)+1,lims,qt.displ[first.col:last.col])} + } + } + } + } + return(output) + } function (obs, forecast, coord1.grid, coord2.grid, forecast.grid, variog.model = "exponential", param.est, n.sim = 99, n.displ = 4, qt.displ = c(10, 50, 90)) { if (missing(variog.model)) variog.model <- "exponential" if (missing(n.sim)) n.sim <- 99 if (missing(n.displ)) n.displ <- 4 if (missing(qt.displ)) qt.displ <- c(10, 50, 90) l.obs <- length(obs) l.for <- length(forecast) if (l.obs == 0) { stop("obs is a vector") } if (l.for == 0) { stop("forecast is a vector") } if (sum((c(l.obs, l.for)/l.obs) == rep(1, 2)) != 2) { stop("Mismatch in dimensions in the data") } if (sum(is.numeric(obs) == rep("TRUE", l.obs)) < l.obs) { stop("obs should be a numeric vector") } if (sum(is.numeric(forecast) == rep("TRUE", l.for)) < l.for) { stop("forecasts should be a numeric vector") } l.coord1grid <- length(coord1.grid) l.coord2grid <- length(coord2.grid) l.forgrid <- length(forecast.grid) if (sum((c(l.coord1grid, l.coord2grid, l.forgrid)/l.coord1grid) == rep(1, 3)) != 3) { stop("Mismatch in dimensions in the data") } if (sum(is.numeric(coord1.grid) == rep(TRUE, l.coord1grid)) < l.coord1grid) { stop("Invalid input for coord1.grid") } if (sum(is.numeric(coord2.grid) == rep(TRUE, l.coord2grid)) < l.coord2grid) { stop("Invalid input for coord2.grid") } if (sum(is.numeric(forecast.grid) == rep(TRUE, l.forgrid)) < l.forgrid) { stop("Invalid input for forecast.grid") } l.variog <- length(variog.model) case.var <- 0 if (l.variog == 0) { variog.model <- "exponential" case.var <- 1 } if (l.variog == 1) { if (is.character(variog.model) == "TRUE") { if (variog.model == "exponential") { case.var <- 1 } if (variog.model == "spherical") { case.var <- 1 } if (variog.model == "whittlematern" | variog.model == "matern") { case.var <- 1 } if (variog.model == "gencauchy") { case.var <- 1 } if (variog.model == "gauss") { case.var <- 1 } } } if (case.var == 0) { stop("Incorrect variogram model specification") } l.parest <- length(param.est) case.parest <- 0 if (l.parest < 3) { stop("Invalid input for param.est") } if (l.parest > 5) { stop("Invalid input for param.est") } if (l.parest == 3) { if (sum(is.numeric(param.est) == rep(TRUE, l.parest)) < 3) { stop("param.est is a numeric vector") } if (variog.model == "whittlematern") { stop("param.est should be of length 4") } if (variog.model == "gencauchy") { stop("param.est should be of length 5") } if (sum(param.est >= 0) < 3) { stop("The entries of param.est should be non-negative numbers") } if (variog.model == "exponential") { case.parest <- 1 } if (variog.model == "spherical") { case.parest <- 1 } if (variog.model == "gauss") { case.parest <- 1 } } if (l.parest == 4) { if (sum(is.numeric(param.est) == rep(TRUE, l.parest)) < 4) { stop("param.est is a numeric vector") } if (variog.model == "gencauchy") { stop("param.est should be of length 5") } if (variog.model == "exponential") { stop("param.est should be of length 3") } if (variog.model == "gauss") { stop("param.est should be of length 3") } if (variog.model == "spherical") { stop("param.est should be of length 3") } if (sum(param.est >= 0) < 4) { stop("The entries of param.est should be non-negative numbers") } if (variog.model == "whittlematern") { case.parest <- 1 } } if (l.parest == 5) { if (sum(is.numeric(param.est) == rep(TRUE, l.parest)) < 5) { stop("param.est is a numeric vector") } if (variog.model == "whittlematern") { stop("param.est should be of length 4") } if (variog.model == "exponential") { stop("param.est should be of length 3") } if (variog.model == "gauss") { stop("param.est should be of length 3") } if (variog.model == "spherical") { stop("param.est should be of length 3") } if (sum(param.est >= 0) < 5) { stop("The entries of param.est should be non-negative numbers") } if (variog.model == "gencauchy") { if (param.est[4] > 2) { stop("a is out of the parameter space") } case.parest <- 1 } } if (case.parest == 0) { stop("Invalid input for param.est") } l.sim <- length(n.sim) if (l.sim == 0) { n.sim <- 99 } l.sim <- length(n.sim) if (length(n.sim) > 1) { stop("n.sim is a numeric field: not a vector") } if (length(n.sim) == 1) { if (is.numeric(n.sim) == FALSE) { stop("The number of simulated fields should be a positive integer") } if (ceiling(n.sim) != n.sim) { stop("The number of simulated fields should be a positive integer") } if (n.sim <= 0) { stop("The number of simulated fields should be a positive integer") } } l.displ <- length(n.displ) if (l.displ == 0) { n.displ <- 4 } l.displ <- length(n.displ) if (l.displ > 1) { stop("n.displ is a numeric field: not a vector") } if (l.displ == 1) { if (is.numeric(n.displ) != TRUE) { stop("n.displ is a numeric field") } if (ceiling(n.displ) != n.displ) { stop("n.displ should be an integer number") } if (n.displ < 0) { stop("n.displ should be an integer non-negative number") } if (n.displ > n.sim) { stop("n.displ should be less or equal to n.sim") } } l.qt <- length(qt.displ) if (l.qt == 0) { qt.displ <- c(10, 50, 90) } l.qt <- length(qt.displ) if (l.qt > 0) { if (sum(is.numeric(qt.displ) == rep("TRUE", l.qt)) < l.qt) { stop("qt.displ is a numeric vector") } if (sum(ceiling(qt.displ) == rep(qt.displ, l.qt)) < l.qt) { stop("The elements of qt.displ should be integer numbers") } for (i in 1:l.qt) { if (qt.displ[i] < 0 | qt.displ[i] > 100) { stop("The elements of qt.displ should be numbers between 0 and 100") } } } simul <- sim.field(variog.model, param.est, coord1.grid, coord2.grid, n.sim) gop.mod <- lm(obs ~ forecast) gop.coeff <- round(gop.mod$coeff, 3) sim.out <- (gop.mod$coeff[1] + gop.mod$coeff[2] * forecast.grid) + simul sim.out.1 <- array(0, c(65, 65, n.sim)) for (i in 1:n.sim) { sim.out.1[, , i] <- engrid(coord1.grid, coord2.grid, sim.out[, i]) } l.qtdispl <- length(qt.displ) if (l.qtdispl == 1 & qt.displ == 0) { quant.out <- 0 qt.out.1 <- NULL } else { quant.out <- matrix(0, ncol = l.qtdispl, nrow = l.coord1grid) qt.displ <- qt.displ/100 for (i in 1:l.qtdispl) { quant.out[, i] <- (gop.mod$coeff[1] + gop.mod$coeff[2] * forecast.grid) + (qnorm(qt.displ[i], 0, 1, TRUE, FALSE) * sqrt(param.est[1] + param.est[2])) } qt.out.1 <- array(0, c(65, 65, l.qtdispl)) for (i in 1:l.qtdispl) { qt.out.1[, , i] <- engrid(coord1.grid, coord2.grid, quant.out[, i]) } } if (variog.model == "whittlematern" | variog.model == "gencauchy") { output <- list(model = variog.model, nugget = param.est[1], variance = param.est[2], range = param.est[3], additional.par = param.est[-seq(1:3)], sim.fields = round(sim.out.1, 4), pct.fields = round(qt.out.1, 4)) } if (variog.model == "exponential" | variog.model == "gauss" | variog.model == "spherical") { output <- list(model = variog.model, nugget = param.est[1], variance = param.est[2], range = param.est[3], sim.fields = round(sim.out.1, 4), pct.fields = round(qt.out.1, 4)) } lims <- c(min(sim.out.1[, , 1:n.displ], na.rm = TRUE), max(sim.out.1[, , 1:n.displ], na.rm = TRUE)) n.displ.4 <- ceiling(n.displ/4) if (n.displ == 1) { x.lim <- c(min(coord1.grid, na.rm = TRUE), max(coord1.grid, na.rm = TRUE)) y.lim <- c(min(coord2.grid, na.rm = TRUE), max(coord2.grid, na.rm = TRUE)) par(ask = TRUE) ens.plot(sim.out.1[, , 1], lims, x.lim, y.lim, "Ensemble member 1") } if (n.displ == 2) { par(mfrow = c(2, 2), ask = TRUE) plotens(coord1.grid, coord2.grid, sim.out.1[, , 1:n.displ], n.displ, lims, "Ensemble member", 0) } if (n.displ > 2) { if (n.displ.4 == 1) { par(mfrow = c(2, 2), ask = TRUE) plotens(coord1.grid, coord2.grid, sim.out.1[, , 1:n.displ], n.displ, lims, "Ensemble member", 0) } if (n.displ.4 > 1) { for (i in 1:n.displ.4) { n.pages <- i - 1 n.pages.4 <- 4 * n.pages if (i != n.displ.4) { par(mfrow = c(2, 2), ask = TRUE) first.col <- (((i - 1) * 4) + 1) last.col <- 4 * i plotens(coord1.grid, coord2.grid, sim.out.1[, , first.col:last.col], 4, lims, "Ensemble member", n.pages.4) } if (i == n.displ.4) { par(mfrow = c(2, 2), ask = TRUE) first.col <- (4 * (i - 1) + 1) last.col <- n.displ plotens(coord1.grid, coord2.grid, sim.out.1[, , first.col:last.col], ((last.col - first.col) + 1), lims, "Ensemble member", n.pages.4) } } } } if (length(qt.out.1) != 0) { lims <- c(min(qt.out.1[, , 1:l.qtdispl], na.rm = TRUE), max(qt.out.1[, , 1:l.qtdispl], na.rm = TRUE)) l.qtdispl.4 <- ceiling(l.qtdispl/4) if (l.qtdispl == 1) { x.lim <- c(min(coord1.grid, na.rm = TRUE), max(coord1.grid, na.rm = TRUE)) y.lim <- c(min(coord2.grid, na.rm = TRUE), max(coord2.grid, na.rm = TRUE)) title <- paste(qt.displ * 100, "-th Percentile") par(ask = TRUE) ens.plot(qt.out.1[, , 1], lims, x.lim, y.lim, title) } if (l.qtdispl == 2) { par(mfrow = c(2, 2), ask = TRUE) plotens.qt(coord1.grid, coord2.grid, qt.out.1[, , 1:l.qtdispl], l.qtdispl, lims, qt.displ) } if (l.qtdispl > 2) { if (l.qtdispl.4 == 1) { par(mfrow = c(2, 2), ask = TRUE) plotens.qt(coord1.grid, coord2.grid, qt.out.1[, , 1:l.qtdispl], l.qtdispl, lims, qt.displ) } if (l.qtdispl.4 > 1) { for (i in 1:l.qtdispl.4) { if (i != l.qtdispl.4) { par(mfrow = c(2, 2), ask = TRUE) first.col <- (((i - 1) * 4) + 1) last.col <- 4 * i plotens.qt(coord1.grid, coord2.grid, qt.out.1[, , first.col:last.col], 4, lims, qt.displ[first.col:last.col]) } if (i == l.qtdispl.4) { par(mfrow = c(2, 2), ask = TRUE) first.col <- (4 * (i - 1) + 1) last.col <- l.qtdispl plotens.qt(coord1.grid, coord2.grid, qt.out.1[, , first.col:last.col], (last.col - first.col) + 1, lims, qt.displ[first.col:last.col]) } } } } } return(output) } > ## End Don't show > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "ProbForecastGOP" > > ### * ProbForecastGOP > > flush(stderr()); flush(stdout()) > > ### Name: ProbForecastGOP > ### Title: Probabilistic Weather Forecasts using the GOP method > ### Aliases: ProbForecastGOP > ### Keywords: file > > ### ** Examples > > ## Loading data > library(fields) > library(RandomFields) > data(slp) > day <-slp$date.obs > id <- slp$id.stat > coord1 <- slp$lon.stat > coord2 <- slp$lat.stat > obs <- slp$obs > forecast <-slp$forecast > > data(gridlong) > coord1.grid <- gridlong$gridded.lon > > data(gridlat) > coord2.grid <- gridlat$gridded.lat > > data(forecast.grid) > forecast.grid <- forecast.grid$gridded.forecast > > ## Specified cutpoints, default values for all the other fields. > ## Only empirical variogram computation > empirical.variog <- ProbForecastGOP(day=day,obs=obs,forecast=forecast, + id=id,coord1=coord1,coord2=coord2,cut.points=seq(0,1000,by=2), + max.dist=NULL,nbins=NULL,variog.model="exponential",max.dist.fit=NULL, + init.val=NULL,fix.nugget=FALSE,coord1.grid=coord1.grid, + coord2.grid=coord2.grid,forecast.grid=forecast.grid,n.sim=10, + out="VARIOG",n.displ=4,qt.displ=c(5,50,95)) > > ## Unspecified cutpoints. > ## Fit of a parametric variogram to an empirical variogram. > fitted.variog <- ProbForecastGOP(day=day,obs=obs,forecast=forecast,id=id, + coord1=coord1,coord2=coord2,cut.points=NULL,max.dist=NULL,nbins=NULL, + variog.model="exponential",max.dist.fit=NULL,init.val=NULL,fix.nugget=FALSE, + coord1.grid=coord1.grid,coord2.grid=coord2.grid,forecast.grid=forecast.grid, + n.sim=10,out="FIT",n.displ=4,qt.displ=c(5,50,95)) > > ## Unspecified cutpoints. > ## Whole routine. > simulation.fields <- + ProbForecastGOP(day=day,obs=obs,forecast=forecast,id=id,coord1=coord1, + coord2=coord2,cut.points=NULL,max.dist=NULL,nbins=NULL,variog.model=NULL, + max.dist.fit=NULL,init.val=NULL,fix.nugget=FALSE,coord1.grid=coord1.grid, + coord2.grid=coord2.grid,forecast.grid=forecast.grid,n.sim=4,out="SIM", + n.displ=4,qt.displ=c(5,50,95)) Hit to see next plot: Hit to see next plot: ## Don't show: > > ## The function is currently defined as > function(day,obs,forecast,id,coord1,coord2,cut.points=NULL,max.dist=NULL,nbins=300,variog.model="exponential",max.dist.fit=NULL, + init.val=NULL,fix.nugget=FALSE,coord1.grid,coord2.grid,forecast.grid,n.sim=99,out="SIM",n.displ=4,qt.displ=c(10,50,90)){ + + #################### INPUT CHECK + # default values + if(missing(cut.points)) + cut.points <- NULL + if(missing(max.dist)) + max.dist <- NULL + if(missing(nbins)) + nbins <- 300 + if(missing(variog.model)) + variog.model <- "exponential" + if(missing(max.dist.fit)) + max.dist.fit <- NULL + if(missing(init.val)) + init.val <- NULL + if(missing(fix.nugget)) + fix.nugget <- "FALSE" + if(missing(n.sim)) + n.sim <- 99 + if(missing(out)) + out <- "SIM" + if(missing(n.displ)) + n.displ <- 4 + if(missing(qt.displ)) + qt.displ <- c(10,50,90) + + # Here we check if the input is right. + + l.day <- length(day) + l.obs <- length(obs) + l.for <- length(forecast) + l.id <- length(id) + l.coord1 <- length(coord1) + l.coord2 <- length(coord2) + + if(sum((c(l.day,l.obs,l.for,l.id,l.coord1,l.coord2)/l.day)==rep(1,6))!=6){ + stop("Mismatch in dimensions in the data") + } + + if(sum(is.numeric(obs)==rep("TRUE",l.obs))=2 & (sum(is.numeric(cut.points)==rep("TRUE",l.cuts))=2 & (sum(is.numeric(cut.points)==rep("TRUE",l.cuts))==l.cuts)){ + if(sum(order(cut.points)==seq(1:l.cuts)) < l.cuts){ + stop("Cut points should be in increasing order") + } + if(sum(cut.points >= 0) < l.cuts){ + stop("Cut points should be non-negative numbers") + } + if(length(cut.points)!=length(unique(cut.points))){ + stop("The vector with cut points should not contain repeated entries") + } + } + + + ## check on the max.dist + + l.mdist <- length(max.dist) + + if(l.mdist > 1){ + stop("max.dist is a numeric field, not a vector")} + + if(l.mdist==1){ + if(is.numeric(max.dist)==FALSE){ + stop("max.dist is a numeric field") + } + if(max.dist < 0){ + stop("max.dist should be a positive number") + } + } + + ## check on the number of bins + l.nbins <- length(nbins) + if(l.nbins==0 & l.cuts==0){ + nbins <- 300 + } + + if(l.nbins==1 & l.cuts >=2){ + nbins=NULL + } + + l.nbins <- length(nbins) + + if(l.nbins >1){ + stop("nbins is a should be an integer: not a vector") + } + + if(l.nbins==1){ + if(ceiling(nbins)!=nbins){ + stop("Invalid input: the number of bins should be a positive integer") + } + if(ceiling(nbins)==nbins & nbins < 0){ + stop("Invalid input: the number of bins should be a positive integer") + } + } + + + ## check on the variog.model + + l.variog <- length(variog.model) + case.var <- 0 + if(l.variog==0){ + variog.model <- "exponential" + case.var <- 1 + } + + if(l.variog==1){ + if(is.character(variog.model)=="TRUE"){ + if(variog.model=="exponential"){ + case.var <- 1} + if(variog.model=="spherical"){ + case.var <- 1} + if(variog.model=="whittlematern" | variog.model=="matern"){ + case.var <- 1} + if(variog.model=="gencauchy"){ + case.var <- 1} + if(variog.model=="gauss"){ + case.var <- 1} + } + } + + if(case.var==0){ + stop("Incorrect variogram model specification") + } + + # check on the max.dist.fit and the initial values + + l.max.dist.fit <- length(max.dist.fit) + + if(l.max.dist.fit > 1){ + stop("max.dist.fit should be a numeric field: not a vector") + } + + if(l.max.dist.fit==1 & is.numeric(max.dist.fit)==FALSE){ + stop("max.dist.fit should be a numeric field") + } + + if(l.max.dist.fit==1 & is.numeric(max.dist.fit)==TRUE){ + if(max.dist.fit < 0){ + stop("max.dist.fit should be a positive number") + } + if(max.dist.fit > max.dist){ + stop("max.dist.fit should be less or equal than max.dist") + } + } + + + ## Input check on the initial values for the parameters + + l.init.val <- length(init.val) + + if(l.init.val > 0 & l.init.val <3 ){ + stop("init.val should be equal to NULL or to a vector of at least length 3") + } + + if(l.init.val ==3 & (sum(is.numeric(init.val)==rep("TRUE",3)) < l.init.val)){ + stop("The initial values should be numeric entries") + } + + if(l.init.val > 3 & (sum(is.numeric(init.val)==rep("TRUE",3))== l.init.val)){ + if(l.init.val==4 & variog.model!="matern"){ + stop("Incorrect number of initial values") + } + if(l.init.val==5 & variog.model!="gencauchy"){ + stop("Incorrect number of initial values") + } + if(l.init.val > 5){ + stop("Incorrect number of initial values") + } + if(sum(init.val>0)0){ + if(fix.nugget[2]!=init.val[1]){ + fix.nugget[2] <- init.val[1] + } + } + } + + if(l.nug >2){ + stop("fix.nugget is either a character field or a field of length 2") + } + + ## check on the latitude and forecast grid + + l.coord1grid <- length(coord1.grid) + l.coord2grid <- length(coord2.grid) + l.forgrid <- length(forecast.grid) + + if(sum((c(l.coord1grid,l.coord2grid,l.forgrid)/l.coord1grid)==rep(1,3))!=3){ + stop("Mismatch in dimensions in the data") + } + + if(sum(is.numeric(coord1.grid)==rep(TRUE,l.coord1grid)) < l.coord1grid){ + stop("Invalid input for coord1.grid") + } + + if(sum(is.numeric(coord2.grid)==rep(TRUE,l.coord2grid)) < l.coord2grid){ + stop("Invalid input for coord2.grid") + } + + if(sum(is.numeric(forecast.grid)==rep(TRUE,l.forgrid)) < l.forgrid){ + stop("Invalid input for forecast.grid") + } + ## check on the number of simulated random fields + + l.sim <- length(n.sim) + if(l.sim==0){ + n.sim <- 99} + l.sim <- length(n.sim) + + if(length(n.sim)> 1){ + stop("n.sim is a numeric field: not a vector") + } + + if(length(n.sim)==1){ + if(is.numeric(n.sim)==FALSE){ + stop("The number of simulated fields should be a positive integer") + } + if(ceiling(n.sim)!=n.sim){ + stop("The number of simulated fields should be a positive integer") + } + if(n.sim < 0){ + stop("The number of simulated fields should be a positive integer") + } + } + + ## check on the out + + l.out <- length(out) + case.out <- 0 + + if(l.out==0){ + out <- "SIM" + case.out <- 1 + } + + if(l.out > 1){ + stop("out is a character field: not a vector") + } + + if(l.out==1){ + if(is.character(out)!=TRUE){ + stop("out is a character field") + } + if(out=="VARIOG"){ case.out <- 1} + if(out=="FIT"){ case.out <- 1} + if(out=="SIM"){ case.out <- 1} + } + + if(case.out==0){ + stop("out is a character field, equal to either VARIOG, FIT or SIM") + } + + ## check on the n.displ + l.displ <- length(n.displ) + if(l.displ==0){ + n.displ <- 4} + + l.displ <- length(n.displ) + if(l.displ > 1){ + stop("n.displ is a numeric field: not a vector") + } + + if(l.displ==1){ + if(is.numeric(n.displ)!=TRUE){ + stop("n.displ is a numeric field") + } + if(ceiling(n.displ)!=n.displ){ + stop("n.displ should be an integer number") + } + if(n.displ < 0){ + stop("n.displ should be an integer non-negative number") + } + if(n.displ > n.sim){ + stop("n.displ should be less or equal to n.sim") + } + } + + ## check on the qt.displ + + l.qt <- length(qt.displ) + if(l.qt==0){ + qt.displ <- c(10,50,90)} + l.qt <- length(qt.displ) + + if(l.qt >0){ + if(sum(is.numeric(qt.displ)==rep("TRUE",l.qt)) 100){ + stop("The elements of qt.displ should be numbers between 0 and 100") + } + } + } + + + ### ESTIMATION of THE EMPIRICAL VARIOGRAM + # here we order the data in ascending date order + day.o <- order(day) + coord1 <- coord1[day.o] + coord2 <- coord2[day.o] + obs <- obs[day.o] + id <- id[day.o] + forecast <- forecast[day.o] + + # the first step in the gop software is to calculate the residuals + + gop.mod <- lm(obs~forecast) + gop.res <- gop.mod$res + res.var <- var(gop.res) + gop.coeff <- round(gop.mod$coeff,3) + gop.sse <- (sum(gop.res^2))/(length(gop.res)-2) + for.dev <- (length(gop.res)-1)*var(forecast) + gop.se2 <- sqrt(gop.sse/for.dev) + for.ssq <- sum(forecast^2) + gop.se1 <- sqrt((gop.sse*for.ssq)/(length(gop.res)*(length(gop.res)-1)*var(forecast))) + gop.se <- c(gop.se1,gop.se2) + gop.sumstat <- round(summary(gop.res),3) + + # the second step is to determine the empirical variograms: if the vector with the cutpoints + # is not specified, we determine the cutpoints by looking at the day with the median average + # of observations, we calculate the cutpoints so that the number + # of bins is equal to the one specified and each bin contains approx the same number of pairs. If the vector with the + # cutpoints is specified, then we just use that vector of cutpoints. + + if(length(cut.points)!=0 & length(max.dist)!=0){ + cut.points <- cut.points[cut.points <= max.dist] + } + + if(length(cut.points)==0){ + # all this part is to determine the day with the median number of observations + obs.day <- table(day) + obs.day.vec <- matrix(obs.day,nrow=length(unique(day)),ncol=1) + obs.median <- round(apply(obs.day.vec,2,median),0) + diff.median <- abs(obs.day-obs.median) + unique.day <- unique(day) + day.index <- min(unique.day[diff.median==min(diff.median)]) + + # this part is to determine the distance among stations that have been recorded in the + # day with the median number of observations and to determine the cutpoints. + + obs.day.index <- obs[day==day.index] + id.day.index <- id[day==day.index] + coord1.day.index <- coord1[day==day.index] + coord2.day.index <- coord2[day==day.index] + dist.day.index <- calc.dist(coord1.day.index,coord2.day.index,id.day.index) + dist.day.index <- dist.day.index[lower.tri(dist.day.index)] + dist.day.index <- dist.day.index[dist.day.index > 0] + ord.dist.day <- sort(dist.day.index) + + + # this is in case we want to consider only those pairs of stations where the distance + # is smaller than the maximum distance allowed among locations. + if(length(max.dist)!= 0){ + ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] + } + + if(length(max.dist)==0){ + max.dist <- quantile(ord.dist.day,.9) + ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] + } + + l.dist.day <- length(ord.dist.day) + l.bins <- floor(l.dist.day/nbins) + for(i in 1:(nbins-1)){ + cut.points[i] <- ord.dist.day[i*l.bins] + } + + } + + + + # this part is to calculate the empirical variogram + emp.variog <- avg.variog(day,coord1,coord2,id,gop.res,cut.points,max.dist,nbins) + + if(out=="VARIOG"){ + ProbForecast.VARIOG <- + list(bias.coeff=round(gop.coeff,3),se.bias.coeff=round(gop.se,3),res.var=round(res.var,3),bin.midpoints=emp.variog$bin.midpoints, + number.pairs=emp.variog$number.pairs,empir.variog=emp.variog$empir.variog) + plot(emp.variog$bin.midpoints,emp.variog$empir.variog,xlab="Distance in km",ylab="Semi-variance",main="Empirical variogram") + return(ProbForecast.VARIOG) + stop} + + + + #### FIT + # in this part we fit a parametric model to the empirical variogram and we estimate the + # parameters using the nlm function of R and a weighted least square approach. The initial + # values for the parameters are calculated by fitting a linear model to the empirical + # variogram and by taking the estimates for the coefficient. Note that the linear model + # fitted depends on the parametric model adopted for the variogram. + + init.var <- var(gop.res) + emp.variog.mod <- list(bin.midpoints=emp.variog$bin.midpoints,number.pairs=emp.variog$number.pairs,empir.variog=emp.variog$empir.variog) + param.est <- round(model.fit(variog.model,emp.variog.mod,max.dist.fit,init.val,init.var,fix.nugget),3) + + if(variog.model=="matern"){variog.model <- "whittlematern"} + if(out=="FIT"){ + ifelse((variog.model=="whittlematern" | variog.model=="gencauchy"), + ProbForecast.FIT <- + list(bias.coeff=round(gop.coeff,3),se.bias.coeff=round(gop.se,3),res.var=round(res.var,3),bin.midpoints=emp.variog$bin.midpoints, + number.pairs=emp.variog$number.pairs,empir.variog=emp.variog$empir.variog,model=variog.model,nugget=param.est[1],variance=param.est[2], + range=param.est[3],additional.par=param.est[-seq(1:3)]), + ProbForecast.FIT <- + list(bias.coeff=round(gop.coeff,3),se.bias.coeff=round(gop.se,3),res.var=round(res.var,3),bin.midpoints=emp.variog$bin.midpoints, + number.pairs=emp.variog$number.pairs,empir.variog=emp.variog$empir.variog,model=variog.model,nugget=param.est[1],variance=param.est[2], + range=param.est[3])) + + plot(emp.variog$bin.midpoints,emp.variog$empir.variog,xlab="Distance in km",ylab="Semivariance") + lines(emp.variog$bin.midpoints,linesmodel(emp.variog$bin.midpoints,variog.model,param.est),lwd=3,col="red") + return(ProbForecast.FIT) + stop} + + + ##### SIM + + # in this part we simulate realizations from the residual random field. Note to be able to + # run this part of the function, the package Random Field must have been downloaded and installed. + # this part is to simulate n.sim random fields with the specified covariance structure and + # with the estimated parameters. Note these random fields have all 0 mean. + + simul <- sim.field(variog.model,param.est,coord1.grid,coord2.grid,n.sim) + + # here we add to the 0-mean simulated random fields the bias-corrected gridded forecasts. + sim.out <- (gop.mod$coeff[1]+gop.mod$coeff[2]*forecast.grid)+ simul + + sim.out.1 <- array(0,c(65,65,n.sim)) + for(i in 1:n.sim){ + sim.out.1[,,i] <- engrid(coord1.grid,coord2.grid,sim.out[,i]) + } + + # here we determine the percentiles of the random fields + l.qtdispl <- length(qt.displ) + if(l.qtdispl==1 & qt.displ[1]==0){ + quant.out <- 0 + qt.out.1 <- NULL + } + else{ + quant.out <- matrix(0,ncol=l.qtdispl,nrow=l.coord1grid) + qt.displ <- qt.displ/100 + for(i in 1:l.qtdispl){ + quant.out[,i] <- (gop.mod$coeff[1]+gop.mod$coeff[2]*forecast.grid) + + (qnorm(qt.displ[i],0,1,TRUE,FALSE)*sqrt(param.est[1]+param.est[2])) + } + qt.out.1 <- array(0,c(65,65,l.qtdispl)) + for(i in 1:l.qtdispl){ + qt.out.1[,,i] <- engrid(coord1.grid,coord2.grid,quant.out[,i]) + } + } + + # here we return the output + if(out=="SIM"){ + if(variog.model=="whittlematern" | variog.model=="gencauchy"){ + ProbForecast.SIM <- + list(bias.coeff=round(gop.coeff,3),se.bias.coeff=round(gop.se,3),res.var=round(res.var,3),bin.midpoints=emp.variog$bin.midpoints, + number.pairs=emp.variog$number.pairs,empir.variog=emp.variog$empir.variog,model=variog.model,nugget=param.est[1],variance=param.est[2], + range=param.est[3],additional.par=param.est[-seq(1:3)],sim.fields=round(sim.out.1,4),pct.fields=round(qt.out.1,4))} + + if(variog.model!="whittlematern" & variog.model!="gencauchy"){ + ProbForecast.SIM <- + list(bias.coeff=round(gop.coeff,3),se.bias.coeff=round(gop.se,3),res.var=round(res.var,3),bin.midpoints=emp.variog$bin.midpoints, + number.pairs=emp.variog$number.pairs,empir.variog=emp.variog$empir.variog,model=variog.model,nugget=param.est[1],variance=param.est[2], + range=param.est[3],sim.fields=round(sim.out.1,4),pct.fields=round(qt.out.1,4))} + + plot(emp.variog$bin.midpoints,emp.variog$empir.variog,xlab="Distance in km",ylab="Semivariance") + lines(emp.variog$bin.midpoints,linesmodel(emp.variog$bin.midpoints,variog.model,param.est),lwd=3,col="red") + + lims <- c(min(sim.out.1[,,1:n.displ],na.rm=TRUE),max(sim.out.1[,,1:n.displ],na.rm=TRUE)) + + n.displ.4 <- ceiling(n.displ/4) + + if(n.displ==1){ + x.lim <- c(min(coord1.grid,na.rm=TRUE),max(coord1.grid,na.rm=TRUE)) + y.lim <- c(min(coord2.grid,na.rm=TRUE),max(coord2.grid,na.rm=TRUE)) + par(ask=TRUE) + ens.plot(sim.out.1[,,1],lims,x.lim,y.lim,"Ensemble member 1") + } + if(n.displ==2){ + par(mfrow=c(2,2),ask=TRUE) + plotens(coord1.grid,coord2.grid,sim.out.1[,,1:n.displ],n.displ,lims,"Ensemble member",0) + } + + if(n.displ > 2){ + if(n.displ.4==1){ + par(mfrow=c(2,2),ask=TRUE) + plotens(coord1.grid,coord2.grid,sim.out.1[,,1:n.displ],n.displ,lims,"Ensemble member",0) + } + + if(n.displ.4 >1){ + for(i in 1:n.displ.4){ + n.pages <- i-1 + n.pages.4 <- 4*n.pages + if(i!= n.displ.4){ + par(mfrow=c(2,2),ask=TRUE) + first.col <- (((i-1)*4)+1) + last.col <- 4*i + plotens(coord1.grid,coord2.grid,sim.out.1[,,first.col:last.col],4,lims,"Ensemble member",n.pages.4)} + if(i==n.displ.4){ + par(mfrow=c(2,2),ask=TRUE) + first.col <- (4*(i-1)+1) + last.col <- n.displ + plotens(coord1.grid,coord2.grid,sim.out.1[,,first.col:last.col],((last.col-first.col)+1),lims,"Ensemble member",n.pages.4)} + } + } + } + + if(length(qt.out.1)!=0){ + lims <- c(min(qt.out.1[,,1:l.qtdispl],na.rm=TRUE),max(qt.out.1[,,1:l.qtdispl],na.rm=TRUE)) + l.qtdispl.4 <- ceiling(l.qtdispl/4) + + if(l.qtdispl==1){ + x.lim <- c(min(coord1.grid,na.rm=TRUE),max(coord1.grid,na.rm=TRUE)) + y.lim <- c(min(coord2.grid,na.rm=TRUE),max(coord2.grid,na.rm=TRUE)) + title <- paste(qt.displ*100,"-th Percentile") + par(ask=TRUE) + ens.plot(qt.out.1[,,1],lims,x.lim,y.lim,title) + } + + if(l.qtdispl==2){ + par(mfrow=c(2,2),ask=TRUE) + plotens.qt(coord1.grid,coord2.grid,qt.out.1[,,1:l.qtdispl],l.qtdispl,lims,qt.displ) + } + + if(l.qtdispl > 2){ + if(l.qtdispl.4==1){ + par(mfrow=c(2,2),ask=TRUE) + plotens.qt(coord1.grid,coord2.grid,qt.out.1[,,1:l.qtdispl],l.qtdispl,lims,qt.displ) + } + + if(l.qtdispl.4 >1){ + for(i in 1:l.qtdispl.4){ + if(i!= l.qtdispl.4){ + par(mfrow=c(2,2),ask=TRUE) + first.col <- (((i-1)*4)+1) + last.col <- 4*i + plotens.qt(coord1.grid,coord2.grid,qt.out.1[,,first.col:last.col],4,lims,qt.displ[first.col:last.col])} + if(i==l.qtdispl.4){ + par(mfrow=c(2,2),ask=TRUE) + first.col <- (4*(i-1)+1) + last.col <- l.qtdispl + plotens.qt(coord1.grid,coord2.grid,qt.out.1[,,first.col:last.col],(last.col-first.col)+1,lims,qt.displ[first.col:last.col])} + } + } + } + } + return(ProbForecast.SIM) + stop + } + } function (day, obs, forecast, id, coord1, coord2, cut.points = NULL, max.dist = NULL, nbins = 300, variog.model = "exponential", max.dist.fit = NULL, init.val = NULL, fix.nugget = FALSE, coord1.grid, coord2.grid, forecast.grid, n.sim = 99, out = "SIM", n.displ = 4, qt.displ = c(10, 50, 90)) { if (missing(cut.points)) cut.points <- NULL if (missing(max.dist)) max.dist <- NULL if (missing(nbins)) nbins <- 300 if (missing(variog.model)) variog.model <- "exponential" if (missing(max.dist.fit)) max.dist.fit <- NULL if (missing(init.val)) init.val <- NULL if (missing(fix.nugget)) fix.nugget <- "FALSE" if (missing(n.sim)) n.sim <- 99 if (missing(out)) out <- "SIM" if (missing(n.displ)) n.displ <- 4 if (missing(qt.displ)) qt.displ <- c(10, 50, 90) l.day <- length(day) l.obs <- length(obs) l.for <- length(forecast) l.id <- length(id) l.coord1 <- length(coord1) l.coord2 <- length(coord2) if (sum((c(l.day, l.obs, l.for, l.id, l.coord1, l.coord2)/l.day) == rep(1, 6)) != 6) { stop("Mismatch in dimensions in the data") } if (sum(is.numeric(obs) == rep("TRUE", l.obs)) < l.obs) { stop("obs should be a numeric vector") } if (sum(is.numeric(forecast) == rep("TRUE", l.for)) < l.for) { stop("forecasts should be a numeric vector") } if (sum(ceiling(day) == day) < l.day) { stop("day should be a vector containing integers") } if (sum(is.numeric(coord1) == rep("TRUE", l.coord1)) < l.coord1 | sum(is.numeric(coord2) == rep("TRUE", l.coord2)) < l.coord2) { stop("coord1 and coord2 should be numeric vectors") } l.cuts <- length(cut.points) if (l.cuts == 1) { stop("cut.points should be a numeric vector") } if (l.cuts >= 2 & (sum(is.numeric(cut.points) == rep("TRUE", l.cuts)) < l.cuts)) { stop("cut.points should be a numeric vector") } if (l.cuts >= 2 & (sum(is.numeric(cut.points) == rep("TRUE", l.cuts)) == l.cuts)) { if (sum(order(cut.points) == seq(1:l.cuts)) < l.cuts) { stop("Cut points should be in increasing order") } if (sum(cut.points >= 0) < l.cuts) { stop("Cut points should be non-negative numbers") } if (length(cut.points) != length(unique(cut.points))) { stop("The vector with cut points should not contain repeated entries") } } l.mdist <- length(max.dist) if (l.mdist > 1) { stop("max.dist is a numeric field, not a vector") } if (l.mdist == 1) { if (is.numeric(max.dist) == FALSE) { stop("max.dist is a numeric field") } if (max.dist < 0) { stop("max.dist should be a positive number") } } l.nbins <- length(nbins) if (l.nbins == 0 & l.cuts == 0) { nbins <- 300 } if (l.nbins == 1 & l.cuts >= 2) { nbins = NULL } l.nbins <- length(nbins) if (l.nbins > 1) { stop("nbins is a should be an integer: not a vector") } if (l.nbins == 1) { if (ceiling(nbins) != nbins) { stop("Invalid input: the number of bins should be a positive integer") } if (ceiling(nbins) == nbins & nbins < 0) { stop("Invalid input: the number of bins should be a positive integer") } } l.variog <- length(variog.model) case.var <- 0 if (l.variog == 0) { variog.model <- "exponential" case.var <- 1 } if (l.variog == 1) { if (is.character(variog.model) == "TRUE") { if (variog.model == "exponential") { case.var <- 1 } if (variog.model == "spherical") { case.var <- 1 } if (variog.model == "whittlematern" | variog.model == "matern") { case.var <- 1 } if (variog.model == "gencauchy") { case.var <- 1 } if (variog.model == "gauss") { case.var <- 1 } } } if (case.var == 0) { stop("Incorrect variogram model specification") } l.max.dist.fit <- length(max.dist.fit) if (l.max.dist.fit > 1) { stop("max.dist.fit should be a numeric field: not a vector") } if (l.max.dist.fit == 1 & is.numeric(max.dist.fit) == FALSE) { stop("max.dist.fit should be a numeric field") } if (l.max.dist.fit == 1 & is.numeric(max.dist.fit) == TRUE) { if (max.dist.fit < 0) { stop("max.dist.fit should be a positive number") } if (max.dist.fit > max.dist) { stop("max.dist.fit should be less or equal than max.dist") } } l.init.val <- length(init.val) if (l.init.val > 0 & l.init.val < 3) { stop("init.val should be equal to NULL or to a vector of at least length 3") } if (l.init.val == 3 & (sum(is.numeric(init.val) == rep("TRUE", 3)) < l.init.val)) { stop("The initial values should be numeric entries") } if (l.init.val > 3 & (sum(is.numeric(init.val) == rep("TRUE", 3)) == l.init.val)) { if (l.init.val == 4 & variog.model != "matern") { stop("Incorrect number of initial values") } if (l.init.val == 5 & variog.model != "gencauchy") { stop("Incorrect number of initial values") } if (l.init.val > 5) { stop("Incorrect number of initial values") } if (sum(init.val > 0) < l.init.val) { stop("Initial values should be positive numbers") } } l.nug <- length(fix.nugget) if (l.nug == 0) { fix.nugget <- "FALSE" } l.nug <- length(fix.nugget) if (l.nug == 1) { if (fix.nugget != TRUE & fix.nugget != FALSE) { stop("fix.nugget should be either equal to TRUE or FALSE") } } if (l.nug == 2) { if (fix.nugget[1] != TRUE) { stop("Invalid input for the fix.nugget field") } if (is.numeric(fix.nugget[2]) == FALSE) { stop("The second entry of the fix.nugget field should be a numeric field") } if (fix.nugget[2] < 0) { stop("The second entry of the fix.nugget field should be a non-negative number") } if (l.init.val > 0) { if (fix.nugget[2] != init.val[1]) { fix.nugget[2] <- init.val[1] } } } if (l.nug > 2) { stop("fix.nugget is either a character field or a field of length 2") } l.coord1grid <- length(coord1.grid) l.coord2grid <- length(coord2.grid) l.forgrid <- length(forecast.grid) if (sum((c(l.coord1grid, l.coord2grid, l.forgrid)/l.coord1grid) == rep(1, 3)) != 3) { stop("Mismatch in dimensions in the data") } if (sum(is.numeric(coord1.grid) == rep(TRUE, l.coord1grid)) < l.coord1grid) { stop("Invalid input for coord1.grid") } if (sum(is.numeric(coord2.grid) == rep(TRUE, l.coord2grid)) < l.coord2grid) { stop("Invalid input for coord2.grid") } if (sum(is.numeric(forecast.grid) == rep(TRUE, l.forgrid)) < l.forgrid) { stop("Invalid input for forecast.grid") } l.sim <- length(n.sim) if (l.sim == 0) { n.sim <- 99 } l.sim <- length(n.sim) if (length(n.sim) > 1) { stop("n.sim is a numeric field: not a vector") } if (length(n.sim) == 1) { if (is.numeric(n.sim) == FALSE) { stop("The number of simulated fields should be a positive integer") } if (ceiling(n.sim) != n.sim) { stop("The number of simulated fields should be a positive integer") } if (n.sim < 0) { stop("The number of simulated fields should be a positive integer") } } l.out <- length(out) case.out <- 0 if (l.out == 0) { out <- "SIM" case.out <- 1 } if (l.out > 1) { stop("out is a character field: not a vector") } if (l.out == 1) { if (is.character(out) != TRUE) { stop("out is a character field") } if (out == "VARIOG") { case.out <- 1 } if (out == "FIT") { case.out <- 1 } if (out == "SIM") { case.out <- 1 } } if (case.out == 0) { stop("out is a character field, equal to either VARIOG, FIT or SIM") } l.displ <- length(n.displ) if (l.displ == 0) { n.displ <- 4 } l.displ <- length(n.displ) if (l.displ > 1) { stop("n.displ is a numeric field: not a vector") } if (l.displ == 1) { if (is.numeric(n.displ) != TRUE) { stop("n.displ is a numeric field") } if (ceiling(n.displ) != n.displ) { stop("n.displ should be an integer number") } if (n.displ < 0) { stop("n.displ should be an integer non-negative number") } if (n.displ > n.sim) { stop("n.displ should be less or equal to n.sim") } } l.qt <- length(qt.displ) if (l.qt == 0) { qt.displ <- c(10, 50, 90) } l.qt <- length(qt.displ) if (l.qt > 0) { if (sum(is.numeric(qt.displ) == rep("TRUE", l.qt)) < l.qt) { stop("qt.displ is a numeric vector") } if (sum(ceiling(qt.displ) == rep(qt.displ, l.qt)) < l.qt) { stop("The elements of qt.displ should be integer numbers") } for (i in 1:l.qt) { if (qt.displ[i] < 0 | qt.displ[i] > 100) { stop("The elements of qt.displ should be numbers between 0 and 100") } } } day.o <- order(day) coord1 <- coord1[day.o] coord2 <- coord2[day.o] obs <- obs[day.o] id <- id[day.o] forecast <- forecast[day.o] gop.mod <- lm(obs ~ forecast) gop.res <- gop.mod$res res.var <- var(gop.res) gop.coeff <- round(gop.mod$coeff, 3) gop.sse <- (sum(gop.res^2))/(length(gop.res) - 2) for.dev <- (length(gop.res) - 1) * var(forecast) gop.se2 <- sqrt(gop.sse/for.dev) for.ssq <- sum(forecast^2) gop.se1 <- sqrt((gop.sse * for.ssq)/(length(gop.res) * (length(gop.res) - 1) * var(forecast))) gop.se <- c(gop.se1, gop.se2) gop.sumstat <- round(summary(gop.res), 3) if (length(cut.points) != 0 & length(max.dist) != 0) { cut.points <- cut.points[cut.points <= max.dist] } if (length(cut.points) == 0) { obs.day <- table(day) obs.day.vec <- matrix(obs.day, nrow = length(unique(day)), ncol = 1) obs.median <- round(apply(obs.day.vec, 2, median), 0) diff.median <- abs(obs.day - obs.median) unique.day <- unique(day) day.index <- min(unique.day[diff.median == min(diff.median)]) obs.day.index <- obs[day == day.index] id.day.index <- id[day == day.index] coord1.day.index <- coord1[day == day.index] coord2.day.index <- coord2[day == day.index] dist.day.index <- calc.dist(coord1.day.index, coord2.day.index, id.day.index) dist.day.index <- dist.day.index[lower.tri(dist.day.index)] dist.day.index <- dist.day.index[dist.day.index > 0] ord.dist.day <- sort(dist.day.index) if (length(max.dist) != 0) { ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] } if (length(max.dist) == 0) { max.dist <- quantile(ord.dist.day, 0.9) ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] } l.dist.day <- length(ord.dist.day) l.bins <- floor(l.dist.day/nbins) for (i in 1:(nbins - 1)) { cut.points[i] <- ord.dist.day[i * l.bins] } } emp.variog <- avg.variog(day, coord1, coord2, id, gop.res, cut.points, max.dist, nbins) if (out == "VARIOG") { ProbForecast.VARIOG <- list(bias.coeff = round(gop.coeff, 3), se.bias.coeff = round(gop.se, 3), res.var = round(res.var, 3), bin.midpoints = emp.variog$bin.midpoints, number.pairs = emp.variog$number.pairs, empir.variog = emp.variog$empir.variog) plot(emp.variog$bin.midpoints, emp.variog$empir.variog, xlab = "Distance in km", ylab = "Semi-variance", main = "Empirical variogram") return(ProbForecast.VARIOG) stop } init.var <- var(gop.res) emp.variog.mod <- list(bin.midpoints = emp.variog$bin.midpoints, number.pairs = emp.variog$number.pairs, empir.variog = emp.variog$empir.variog) param.est <- round(model.fit(variog.model, emp.variog.mod, max.dist.fit, init.val, init.var, fix.nugget), 3) if (variog.model == "matern") { variog.model <- "whittlematern" } if (out == "FIT") { ifelse((variog.model == "whittlematern" | variog.model == "gencauchy"), ProbForecast.FIT <- list(bias.coeff = round(gop.coeff, 3), se.bias.coeff = round(gop.se, 3), res.var = round(res.var, 3), bin.midpoints = emp.variog$bin.midpoints, number.pairs = emp.variog$number.pairs, empir.variog = emp.variog$empir.variog, model = variog.model, nugget = param.est[1], variance = param.est[2], range = param.est[3], additional.par = param.est[-seq(1:3)]), ProbForecast.FIT <- list(bias.coeff = round(gop.coeff, 3), se.bias.coeff = round(gop.se, 3), res.var = round(res.var, 3), bin.midpoints = emp.variog$bin.midpoints, number.pairs = emp.variog$number.pairs, empir.variog = emp.variog$empir.variog, model = variog.model, nugget = param.est[1], variance = param.est[2], range = param.est[3])) plot(emp.variog$bin.midpoints, emp.variog$empir.variog, xlab = "Distance in km", ylab = "Semivariance") lines(emp.variog$bin.midpoints, linesmodel(emp.variog$bin.midpoints, variog.model, param.est), lwd = 3, col = "red") return(ProbForecast.FIT) stop } simul <- sim.field(variog.model, param.est, coord1.grid, coord2.grid, n.sim) sim.out <- (gop.mod$coeff[1] + gop.mod$coeff[2] * forecast.grid) + simul sim.out.1 <- array(0, c(65, 65, n.sim)) for (i in 1:n.sim) { sim.out.1[, , i] <- engrid(coord1.grid, coord2.grid, sim.out[, i]) } l.qtdispl <- length(qt.displ) if (l.qtdispl == 1 & qt.displ[1] == 0) { quant.out <- 0 qt.out.1 <- NULL } else { quant.out <- matrix(0, ncol = l.qtdispl, nrow = l.coord1grid) qt.displ <- qt.displ/100 for (i in 1:l.qtdispl) { quant.out[, i] <- (gop.mod$coeff[1] + gop.mod$coeff[2] * forecast.grid) + (qnorm(qt.displ[i], 0, 1, TRUE, FALSE) * sqrt(param.est[1] + param.est[2])) } qt.out.1 <- array(0, c(65, 65, l.qtdispl)) for (i in 1:l.qtdispl) { qt.out.1[, , i] <- engrid(coord1.grid, coord2.grid, quant.out[, i]) } } if (out == "SIM") { if (variog.model == "whittlematern" | variog.model == "gencauchy") { ProbForecast.SIM <- list(bias.coeff = round(gop.coeff, 3), se.bias.coeff = round(gop.se, 3), res.var = round(res.var, 3), bin.midpoints = emp.variog$bin.midpoints, number.pairs = emp.variog$number.pairs, empir.variog = emp.variog$empir.variog, model = variog.model, nugget = param.est[1], variance = param.est[2], range = param.est[3], additional.par = param.est[-seq(1:3)], sim.fields = round(sim.out.1, 4), pct.fields = round(qt.out.1, 4)) } if (variog.model != "whittlematern" & variog.model != "gencauchy") { ProbForecast.SIM <- list(bias.coeff = round(gop.coeff, 3), se.bias.coeff = round(gop.se, 3), res.var = round(res.var, 3), bin.midpoints = emp.variog$bin.midpoints, number.pairs = emp.variog$number.pairs, empir.variog = emp.variog$empir.variog, model = variog.model, nugget = param.est[1], variance = param.est[2], range = param.est[3], sim.fields = round(sim.out.1, 4), pct.fields = round(qt.out.1, 4)) } plot(emp.variog$bin.midpoints, emp.variog$empir.variog, xlab = "Distance in km", ylab = "Semivariance") lines(emp.variog$bin.midpoints, linesmodel(emp.variog$bin.midpoints, variog.model, param.est), lwd = 3, col = "red") lims <- c(min(sim.out.1[, , 1:n.displ], na.rm = TRUE), max(sim.out.1[, , 1:n.displ], na.rm = TRUE)) n.displ.4 <- ceiling(n.displ/4) if (n.displ == 1) { x.lim <- c(min(coord1.grid, na.rm = TRUE), max(coord1.grid, na.rm = TRUE)) y.lim <- c(min(coord2.grid, na.rm = TRUE), max(coord2.grid, na.rm = TRUE)) par(ask = TRUE) ens.plot(sim.out.1[, , 1], lims, x.lim, y.lim, "Ensemble member 1") } if (n.displ == 2) { par(mfrow = c(2, 2), ask = TRUE) plotens(coord1.grid, coord2.grid, sim.out.1[, , 1:n.displ], n.displ, lims, "Ensemble member", 0) } if (n.displ > 2) { if (n.displ.4 == 1) { par(mfrow = c(2, 2), ask = TRUE) plotens(coord1.grid, coord2.grid, sim.out.1[, , 1:n.displ], n.displ, lims, "Ensemble member", 0) } if (n.displ.4 > 1) { for (i in 1:n.displ.4) { n.pages <- i - 1 n.pages.4 <- 4 * n.pages if (i != n.displ.4) { par(mfrow = c(2, 2), ask = TRUE) first.col <- (((i - 1) * 4) + 1) last.col <- 4 * i plotens(coord1.grid, coord2.grid, sim.out.1[, , first.col:last.col], 4, lims, "Ensemble member", n.pages.4) } if (i == n.displ.4) { par(mfrow = c(2, 2), ask = TRUE) first.col <- (4 * (i - 1) + 1) last.col <- n.displ plotens(coord1.grid, coord2.grid, sim.out.1[, , first.col:last.col], ((last.col - first.col) + 1), lims, "Ensemble member", n.pages.4) } } } } if (length(qt.out.1) != 0) { lims <- c(min(qt.out.1[, , 1:l.qtdispl], na.rm = TRUE), max(qt.out.1[, , 1:l.qtdispl], na.rm = TRUE)) l.qtdispl.4 <- ceiling(l.qtdispl/4) if (l.qtdispl == 1) { x.lim <- c(min(coord1.grid, na.rm = TRUE), max(coord1.grid, na.rm = TRUE)) y.lim <- c(min(coord2.grid, na.rm = TRUE), max(coord2.grid, na.rm = TRUE)) title <- paste(qt.displ * 100, "-th Percentile") par(ask = TRUE) ens.plot(qt.out.1[, , 1], lims, x.lim, y.lim, title) } if (l.qtdispl == 2) { par(mfrow = c(2, 2), ask = TRUE) plotens.qt(coord1.grid, coord2.grid, qt.out.1[, , 1:l.qtdispl], l.qtdispl, lims, qt.displ) } if (l.qtdispl > 2) { if (l.qtdispl.4 == 1) { par(mfrow = c(2, 2), ask = TRUE) plotens.qt(coord1.grid, coord2.grid, qt.out.1[, , 1:l.qtdispl], l.qtdispl, lims, qt.displ) } if (l.qtdispl.4 > 1) { for (i in 1:l.qtdispl.4) { if (i != l.qtdispl.4) { par(mfrow = c(2, 2), ask = TRUE) first.col <- (((i - 1) * 4) + 1) last.col <- 4 * i plotens.qt(coord1.grid, coord2.grid, qt.out.1[, , first.col:last.col], 4, lims, qt.displ[first.col:last.col]) } if (i == l.qtdispl.4) { par(mfrow = c(2, 2), ask = TRUE) first.col <- (4 * (i - 1) + 1) last.col <- l.qtdispl plotens.qt(coord1.grid, coord2.grid, qt.out.1[, , first.col:last.col], (last.col - first.col) + 1, lims, qt.displ[first.col:last.col]) } } } } } return(ProbForecast.SIM) stop } } > ## End Don't show > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "Variog.fit" > > ### * Variog.fit > > flush(stderr()); flush(stdout()) > > ### Name: Variog.fit > ### Title: Fitting a parametric variogram model to an empirical variogram > ### Aliases: Variog.fit > ### Keywords: file > > ### ** Examples > > ## Loading data > library(fields) > library(RandomFields) > data(slp) > day <- slp$date.obs > id <- slp$id.stat > coord1 <- slp$lon.stat > coord2 <- slp$lat.stat > obs <- slp$obs > forecast <- slp$forecast > > ## Computing the empirical variogram > variogram <- + Emp.variog(day=day,obs=obs,forecast=forecast,id=id,coord1=coord1,coord2=coord2, + cut.points=NULL,max.dist=NULL,nbins=NULL) > ## Estimating parameters > ## Without specifying initial values for the parameters > param.variog <- + Variog.fit(emp.variog=variogram,variog.model="exponential",max.dist.fit=NULL, + init.val=NULL,fix.nugget=FALSE) > ## Plotting the empirical variogram with the estimated parametric variogram superimposed > plot(variogram$bin.midpoints,variogram$empir.variog,xlab="Distance",ylab="Semi-variance") > lines(variogram$bin.midpoints,linesmodel(distance=variogram$bin.midpoints, + variog.model="exponential",param=c(param.variog$nugget, + param.variog$variance,param.variog$range))) > > ## Specifying the initial values for the parameters and keeping the nugget effect fixed > param.variog <- + Variog.fit(emp.variog=variogram,variog.model="exponential",max.dist.fit=NULL, + init.val=c(0,2,100),fix.nugget=TRUE) > ## Plotting the empirical variogram with superimposed the estimated parametric variogram > plot(variogram$bin.midpoints,variogram$empir.variog,xlab="Distance",ylab="Semi-variance") > lines(variogram$bin.midpoints,linesmodel(distance=variogram$bin.midpoints, + variog.model="exponential",param=c(param.variog$nugget, + param.variog$variance,param.variog$range))) > > ## Don't show: ## The function is currently defined as > function(emp.variog,variog.model="exponential",max.dist.fit=NULL,init.val=NULL,fix.nugget=FALSE){ + # INPUT CHECK # + ## Here there should be the check on whether emp.variog is an object output of the function emp.variog + init.var <- emp.variog$mar.var + emp.variog.mod <- list(bin.midpoints=emp.variog$bin.midpoints,number.pairs=emp.variog$number.pairs,empir.variog=emp.variog$empir.variog) + # Here we do the input check on the rest of the input + # default if(missing(max.dist.fit)) + max.dist.fit <- NULL + if(missing(fix.nugget)) + fix.nugget <- FALSE + + ## check on the variog.model + + l.variog <- length(variog.model) + case.var <- 0 + if(l.variog==0){ + variog.model <- "exponential" + case.var <- 1 + } + + if(l.variog==1){ + if(is.character(variog.model)=="TRUE"){ + if(variog.model=="exponential"){ + case.var <- 1} + if(variog.model=="spherical"){ + case.var <- 1} + if(variog.model=="whittlematern" | variog.model=="matern"){ + case.var <- 1} + if(variog.model=="gencauchy"){ + case.var <- 1} + if(variog.model=="gauss"){ + case.var <- 1} + } + } + + if(case.var==0){ + print("Incorrect variogram model specification") + stop + } + + # check on the max.dist.fit and the initial values + + l.max.dist.fit <- length(max.dist.fit) + + if(l.max.dist.fit > 1){ + stop("max.dist.fit should be numeric field: not a vector") + } + + if(l.max.dist.fit==1 & is.numeric(max.dist.fit)==FALSE){ + stop("max.dist.fit should be a numeric field") + } + + if(l.max.dist.fit==1 & is.numeric(max.dist.fit)==TRUE){ + if(max.dist.fit < 0){ + stop("max.dist.fit should be positive number") + } + if(max.dist.fit > max(emp.variog[,1])){ + stop("max.dist.fit should be less or equal than max.dist") + } + } + + + ## Input check on the initial values for the parameters + + l.init.val <- length(init.val) + + if(l.init.val > 0 & l.init.val <3 ){ + stop("init.val should be equal to NULL or to a vector of at least length 3") + } + + if(l.init.val ==3 & (sum(is.numeric(init.val)==rep("TRUE",3)) < l.init.val)){ + stop("The initial values should be numeric entries") + } + + if(l.init.val > 3 & (sum(is.numeric(init.val)==rep("TRUE",3))== l.init.val)){ + if(l.init.val==4 & variog.model!="matern"){ + stop("Incorrect number of initial values") + } + if(l.init.val==5 & variog.model!="gencauchy"){ + stop("Incorrect number of initial values") + } + if(l.init.val > 5){ + stop("Incorrect number of initial values") + } + if(init.val[1] <0){ + stop("The nugget effect cannot be negative") + } + if(sum(init.val[2:l.init.val]>0)<(l.init.val-1)){ + stop("Initial values for all the parameters, but the nugget effect, should be positive numbers") + } + } + + ## Input check for the fix.nugget field + + l.nug <- length(fix.nugget) + if(l.nug==0){ + fix.nugget <- "FALSE"} + l.nug <- length(fix.nugget) + + if(l.nug==1){ + if(fix.nugget!=TRUE & fix.nugget!=FALSE){ + stop("fix.nugget should be either equal to TRUE or FALSE") + } + } + + if(l.nug==2){ + if(fix.nugget[1]!=TRUE){ + stop("Invalid input for the fix.nugget field") + } + if(is.numeric(fix.nugget[2])==FALSE){ + stop("The second entry of the fix.nugget field should be a numeric field") + } + if(fix.nugget[2] <0){ + stop("The second entry of the fix.nugget field should be a non-negative number") + } + if(l.init.val >0){ + if(fix.nugget[2]!=init.val[1]){ + fix.nugget[2] <- init.val[1] + } + } + } + + if(l.nug >2){ + stop("fix.nugget is either a character field or a 2x1 field") + } + + param.est <- round(model.fit(variog.model,emp.variog.mod,max.dist.fit,init.val,init.var,fix.nugget),3) + + if(variog.model=="matern"){variog.model <- "whittlematern"} + ifelse((variog.model=="whittlematern" | variog.model=="gencauchy"), + output <- list(model=variog.model,nugget=param.est[1],variance=param.est[2],range=param.est[3],additional.par=param.est[-seq(1:3)]), + output <- list(model=variog.model,nugget=param.est[1],variance=param.est[2],range=param.est[3])) + return(output) + } function (emp.variog, variog.model = "exponential", max.dist.fit = NULL, init.val = NULL, fix.nugget = FALSE) { init.var <- emp.variog$mar.var emp.variog.mod <- list(bin.midpoints = emp.variog$bin.midpoints, number.pairs = emp.variog$number.pairs, empir.variog = emp.variog$empir.variog) max.dist.fit <- NULL if (missing(fix.nugget)) fix.nugget <- FALSE l.variog <- length(variog.model) case.var <- 0 if (l.variog == 0) { variog.model <- "exponential" case.var <- 1 } if (l.variog == 1) { if (is.character(variog.model) == "TRUE") { if (variog.model == "exponential") { case.var <- 1 } if (variog.model == "spherical") { case.var <- 1 } if (variog.model == "whittlematern" | variog.model == "matern") { case.var <- 1 } if (variog.model == "gencauchy") { case.var <- 1 } if (variog.model == "gauss") { case.var <- 1 } } } if (case.var == 0) { print("Incorrect variogram model specification") stop } l.max.dist.fit <- length(max.dist.fit) if (l.max.dist.fit > 1) { stop("max.dist.fit should be numeric field: not a vector") } if (l.max.dist.fit == 1 & is.numeric(max.dist.fit) == FALSE) { stop("max.dist.fit should be a numeric field") } if (l.max.dist.fit == 1 & is.numeric(max.dist.fit) == TRUE) { if (max.dist.fit < 0) { stop("max.dist.fit should be positive number") } if (max.dist.fit > max(emp.variog[, 1])) { stop("max.dist.fit should be less or equal than max.dist") } } l.init.val <- length(init.val) if (l.init.val > 0 & l.init.val < 3) { stop("init.val should be equal to NULL or to a vector of at least length 3") } if (l.init.val == 3 & (sum(is.numeric(init.val) == rep("TRUE", 3)) < l.init.val)) { stop("The initial values should be numeric entries") } if (l.init.val > 3 & (sum(is.numeric(init.val) == rep("TRUE", 3)) == l.init.val)) { if (l.init.val == 4 & variog.model != "matern") { stop("Incorrect number of initial values") } if (l.init.val == 5 & variog.model != "gencauchy") { stop("Incorrect number of initial values") } if (l.init.val > 5) { stop("Incorrect number of initial values") } if (init.val[1] < 0) { stop("The nugget effect cannot be negative") } if (sum(init.val[2:l.init.val] > 0) < (l.init.val - 1)) { stop("Initial values for all the parameters, but the nugget effect, should be positive numbers") } } l.nug <- length(fix.nugget) if (l.nug == 0) { fix.nugget <- "FALSE" } l.nug <- length(fix.nugget) if (l.nug == 1) { if (fix.nugget != TRUE & fix.nugget != FALSE) { stop("fix.nugget should be either equal to TRUE or FALSE") } } if (l.nug == 2) { if (fix.nugget[1] != TRUE) { stop("Invalid input for the fix.nugget field") } if (is.numeric(fix.nugget[2]) == FALSE) { stop("The second entry of the fix.nugget field should be a numeric field") } if (fix.nugget[2] < 0) { stop("The second entry of the fix.nugget field should be a non-negative number") } if (l.init.val > 0) { if (fix.nugget[2] != init.val[1]) { fix.nugget[2] <- init.val[1] } } } if (l.nug > 2) { stop("fix.nugget is either a character field or a 2x1 field") } param.est <- round(model.fit(variog.model, emp.variog.mod, max.dist.fit, init.val, init.var, fix.nugget), 3) if (variog.model == "matern") { variog.model <- "whittlematern" } ifelse((variog.model == "whittlematern" | variog.model == "gencauchy"), output <- list(model = variog.model, nugget = param.est[1], variance = param.est[2], range = param.est[3], additional.par = param.est[-seq(1:3)]), output <- list(model = variog.model, nugget = param.est[1], variance = param.est[2], range = param.est[3])) return(output) } > ## End Don't show > > > > cleanEx(); ..nameEx <- "avg.variog" > > ### * avg.variog > > flush(stderr()); flush(stdout()) > > ### Name: avg.variog > ### Title: Empirical variogram of a random variable averaged over time > ### Aliases: avg.variog > ### Keywords: file > > ### ** Examples > > ## Loading data > data(slp) > day <- slp$date.obs > id <- slp$id.stat > coord1 <- slp$lon.stat > coord2 <- slp$lat.stat > obs <- slp$obs > forecast <- slp$forecast > > ## Computing variogram of observed temperature > ## No specified cutpoints, no specified maximum distance > ## Default number of bins > variogram <- avg.variog(day=day, coord1=coord1,coord2=coord2,id=id,variable=obs,cut.points=NULL,max.dist=NULL,nbins=NULL) > ## Plotting variogram > plot(variogram$bin.midpoints,variogram$empir.variog,xlab="Distance", ylab="Semi-variance",main="Empirical variogram") > > ## Don't show: ## The function is currently defined as > function(day,coord1,coord2,id,variable,cut.points=NULL,max.dist=NULL,nbins=300){ + # default values + if(missing(cut.points)) + cut.points <- NULL + if(missing(max.dist)) + max.dist <- NULL + if(missing(nbins)) + nbins <- NULL + #INPUT CHECK + #Here we check if the input is right. + l.day <- length(day) + l.var <- length(variable) + l.id <- length(id) + l.coord1 <- length(coord1) + l.coord2 <- length(coord2) + if(sum((c(l.day,l.var,l.id,l.coord1,l.coord2)/l.day)==rep(1,5))!=5){ + stop("Error in the dimension of the data!") + } + if(sum(is.numeric(variable)==rep("TRUE",l.var))=2 & (sum(is.numeric(cut.points)==rep("TRUE",l.cuts))=2 & (sum(is.numeric(cut.points)==rep("TRUE",l.cuts))==l.cuts)){ + if(sum(order(cut.points)==seq(1:l.cuts)) < l.cuts){ + stop("Cut points should be in increasing order") + } + if(sum(cut.points >= 0) < l.cuts){ + stop("Cut points should be non-negative numbers") + } + if(length(cut.points)!=length(unique(cut.points))){ + stop("The vector with cut points should not contain repeated entries") + } + } + ## check on the max.dist + l.mdist <- length(max.dist) + if(l.mdist > 1){ + stop("max.dist is a numeric field, not a vector") + } + if(l.mdist==1){ + if(is.numeric(max.dist)==FALSE){ + stop("max.dist is a numeric field") + } + if(max.dist < 0){ + stop("max.dist should be a positive number") + } + } + ## check on the number of bins + l.nbins <- length(nbins) + if(l.nbins==0 & l.cuts==0){ + nbins <- 300 + } + if(l.nbins==1 & l.cuts >=2){ + nbins=NULL + } + l.nbins <- length(nbins) + if(l.nbins >1){ + stop("nbins should be an integer: not a vector") + } + if(l.nbins==1){ + if(ceiling(nbins)!=nbins){ + stop("Invalid input: the number of bins should be a positive + integer") + } + } + + # ESTIMATION of THE EMPIRICAL VARIOGRAM + # here we order the data in ascending date order + day.o <- order(day) + coord1 <- coord1[day.o] + coord2 <- coord2[day.o] + variable <- variable[day.o] + marginal.var <- var(variable) + id <- id[day.o] + + # Here we determine the vector of cutpoints. + # If the vector with the cutpoints is not specified, we determine the cutpoints by looking at the day with the median number + # of observations, we calculate the cutpoints so that the number + # of bins is equal to the one specified and each bin contains approx the same number of pairs. If the vector with the + # cutpoints is specified, then we just use that vector of cutpoints. + if(length(cut.points)!=0 & length(max.dist)!=0){ + cut.points <- cut.points[cut.points <= max.dist] + } + if(length(cut.points)==0){ + # all this part is to determine the day with the median number of observations + variable.day <- table(day) + variable.day.vec <- matrix(variable.day,nrow=length(unique(day)),ncol=1) + variable.median <- round(apply(variable.day.vec,2,median),0) + diff.median <- abs(variable.day-variable.median) + unique.day <- unique(day) + day.index <- min(unique.day[diff.median==min(diff.median)]) + # this part is to determine the distance among stations that have been recorded in the + # day with the median number of observations and to determine the cutpoints. + variable.day.index <- variable[day==day.index] + id.day.index <- id[day==day.index] + coord1.day.index <- coord1[day==day.index] + coord2.day.index <- coord2[day==day.index] + dist.day.index <- calc.dist(coord1.day.index,coord2.day.index,id.day.index) + dist.day.index <- dist.day.index[lower.tri(dist.day.index)] + ord.dist.day <- sort(dist.day.index) + + # this is in case we want to consider only those pairs of stations where the distance + # is smaller than the maximum distance allowed among locations. + if(length(max.dist)!= 0){ + ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] + } + + if(length(max.dist)==0){ + max.dist <- quantile(ord.dist.day,.9) + ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] + } + + if(length(max.dist)==0){ + max.dist <- quantile(ord.dist.day,.9) + ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] + } + + l.dist.day <- length(ord.dist.day) + + l.bins <- floor(l.dist.day/nbins) + for(i in 1:(nbins-1)){ + cut.points[i] <- ord.dist.day[i*l.bins] + } + } + + ## Here is where the actual function starts + + unique.day <- unique(day) + l.day <- length(unique.day) + + n.stations <- length(unique(id)) + n.cuts <- length(cut.points) + W <- rep(0,(n.cuts-1)) + W.obs <- rep(0,(n.cuts-1)) + for(i in 1:l.day){ + distance.day <- NULL + difference.day <- NULL + new.dist.day <- NULL + s.new.dist.day <- NULL + o.new.dist.day <- NULL + new.diff.day <- NULL + s.new.diff.day <- NULL + variable.day <- NULL + id.day <- NULL + coord1.day <- NULL + coord2.day <- NULL + variable.day <- variable[day==unique.day[i]] + id.day <- id[day==unique.day[i]] + coord1.day <- coord1[day==unique.day[i]] + coord2.day <- coord2[day==unique.day[i]] + distance.day <- calc.dist(coord1.day,coord2.day,id.day) + new.dist.day <- distance.day[lower.tri(distance.day)] + s.new.dist.day <- sort(new.dist.day) + o.new.dist.day <- order(new.dist.day) + difference.day <- calc.difference(variable.day) + new.diff.day <- difference.day[lower.tri(difference.day)] + s.new.diff.day <- new.diff.day[o.new.dist.day] + W.new <- rep(0,(n.cuts-1)) + W.obs.new <- rep(0,(n.cuts-1)) + + for(s in 1:(n.cuts-1)){ + low.bound <- cut.points[s] + upp.bound <- cut.points[s+1] + v.dist <- NULL + v.dist1 <- NULL + v.diff1 <- NULL + index.v.dist <- NULL + index.v.dist1 <- NULL + + if(s < (n.cuts-1)){ + v.dist <- s.new.dist.day[s.new.dist.day >= low.bound & s.new.dist.day < upp.bound] + index.v.dist <- seq(1:length(s.new.dist.day))[s.new.dist.day >= low.bound & s.new.dist.day < upp.bound] + v.dist1 <- v.dist[v.dist!=0] + index.v.dist1 <- index.v.dist[v.dist!=0] + + if(length(v.dist1) >=1){ + v.diff1 <- s.new.diff.day[index.v.dist1] + W.new[s] <- sum(v.diff1) + W.obs.new[s] <- length(v.dist1)} + if(length(v.dist1) ==0){ + W.new[s] <- 0 + W.obs.new[s] <- 0} + } + if(s==(n.cuts-1)){ + v.dist <- s.new.dist.day[s.new.dist.day >= low.bound & s.new.dist.day <= upp.bound] + index.v.dist <- seq(1:length(s.new.dist.day))[s.new.dist.day >= low.bound & s.new.dist.day <= upp.bound] + v.dist1 <- v.dist[v.dist!=0] + index.v.dist1 <- index.v.dist[v.dist!=0] + + if(length(v.dist1) >=1){ + v.diff1 <- s.new.diff.day[index.v.dist1] + W.new[s] <- sum(v.diff1) + W.obs.new[s] <- length(v.dist1)} + if(length(v.dist1) ==0){ + W.new[s] <- 0 + } + } + W <- W+W.new + W.obs <- W.obs + W.obs.new + } + avg.variog <- round(W/(2*W.obs),2) + n.h <- W.obs + x.vec <- NULL + for(i in 1:(n.cuts-1)){ + x.vec[i] <- (cut.points[i]+cut.points[i+1])/2 + } + fin.avg.variog <- c(0,avg.variog) + fin.x.vec <- c(0,x.vec) + fin.n.h <- c(n.stations,n.h) + B <- list(mar.var=round(marginal.var,3),bin.midpoints=fin.x.vec,number.pairs=fin.n.h,empir.variog=fin.avg.variog) + return(B) + } + } function (day, coord1, coord2, id, variable, cut.points = NULL, max.dist = NULL, nbins = 300) { if (missing(cut.points)) cut.points <- NULL if (missing(max.dist)) max.dist <- NULL if (missing(nbins)) nbins <- NULL l.day <- length(day) l.var <- length(variable) l.id <- length(id) l.coord1 <- length(coord1) l.coord2 <- length(coord2) if (sum((c(l.day, l.var, l.id, l.coord1, l.coord2)/l.day) == rep(1, 5)) != 5) { stop("Error in the dimension of the data!") } if (sum(is.numeric(variable) == rep("TRUE", l.var)) < l.var) { stop("variable should be a numeric vector") } if (sum(ceiling(day) == day) < l.day) { stop("day should be a vector containing integer") } if (sum(is.numeric(coord1) == rep("TRUE", l.coord1)) < l.coord1 | sum(is.numeric(coord2) == rep("TRUE", l.coord2)) < l.coord2) { stop("coord1 and coord2 should be numeric vectors") } l.cuts <- length(cut.points) if (l.cuts == 1) { stop("cut.points should be a numeric vector") } if (l.cuts >= 2 & (sum(is.numeric(cut.points) == rep("TRUE", l.cuts)) < l.cuts)) { stop("cut.points should be a numeric vector") } if (l.cuts >= 2 & (sum(is.numeric(cut.points) == rep("TRUE", l.cuts)) == l.cuts)) { if (sum(order(cut.points) == seq(1:l.cuts)) < l.cuts) { stop("Cut points should be in increasing order") } if (sum(cut.points >= 0) < l.cuts) { stop("Cut points should be non-negative numbers") } if (length(cut.points) != length(unique(cut.points))) { stop("The vector with cut points should not contain repeated entries") } } l.mdist <- length(max.dist) if (l.mdist > 1) { stop("max.dist is a numeric field, not a vector") } if (l.mdist == 1) { if (is.numeric(max.dist) == FALSE) { stop("max.dist is a numeric field") } if (max.dist < 0) { stop("max.dist should be a positive number") } } l.nbins <- length(nbins) if (l.nbins == 0 & l.cuts == 0) { nbins <- 300 } if (l.nbins == 1 & l.cuts >= 2) { nbins = NULL } l.nbins <- length(nbins) if (l.nbins > 1) { stop("nbins should be an integer: not a vector") } if (l.nbins == 1) { if (ceiling(nbins) != nbins) { stop("Invalid input: the number of bins should be a positive\ninteger") } } day.o <- order(day) coord1 <- coord1[day.o] coord2 <- coord2[day.o] variable <- variable[day.o] marginal.var <- var(variable) id <- id[day.o] if (length(cut.points) != 0 & length(max.dist) != 0) { cut.points <- cut.points[cut.points <= max.dist] } if (length(cut.points) == 0) { variable.day <- table(day) variable.day.vec <- matrix(variable.day, nrow = length(unique(day)), ncol = 1) variable.median <- round(apply(variable.day.vec, 2, median), 0) diff.median <- abs(variable.day - variable.median) unique.day <- unique(day) day.index <- min(unique.day[diff.median == min(diff.median)]) variable.day.index <- variable[day == day.index] id.day.index <- id[day == day.index] coord1.day.index <- coord1[day == day.index] coord2.day.index <- coord2[day == day.index] dist.day.index <- calc.dist(coord1.day.index, coord2.day.index, id.day.index) dist.day.index <- dist.day.index[lower.tri(dist.day.index)] ord.dist.day <- sort(dist.day.index) if (length(max.dist) != 0) { ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] } if (length(max.dist) == 0) { max.dist <- quantile(ord.dist.day, 0.9) ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] } if (length(max.dist) == 0) { max.dist <- quantile(ord.dist.day, 0.9) ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] } l.dist.day <- length(ord.dist.day) l.bins <- floor(l.dist.day/nbins) for (i in 1:(nbins - 1)) { cut.points[i] <- ord.dist.day[i * l.bins] } } unique.day <- unique(day) l.day <- length(unique.day) n.stations <- length(unique(id)) n.cuts <- length(cut.points) W <- rep(0, (n.cuts - 1)) W.obs <- rep(0, (n.cuts - 1)) for (i in 1:l.day) { distance.day <- NULL difference.day <- NULL new.dist.day <- NULL s.new.dist.day <- NULL o.new.dist.day <- NULL new.diff.day <- NULL s.new.diff.day <- NULL variable.day <- NULL id.day <- NULL coord1.day <- NULL coord2.day <- NULL variable.day <- variable[day == unique.day[i]] id.day <- id[day == unique.day[i]] coord1.day <- coord1[day == unique.day[i]] coord2.day <- coord2[day == unique.day[i]] distance.day <- calc.dist(coord1.day, coord2.day, id.day) new.dist.day <- distance.day[lower.tri(distance.day)] s.new.dist.day <- sort(new.dist.day) o.new.dist.day <- order(new.dist.day) difference.day <- calc.difference(variable.day) new.diff.day <- difference.day[lower.tri(difference.day)] s.new.diff.day <- new.diff.day[o.new.dist.day] W.new <- rep(0, (n.cuts - 1)) W.obs.new <- rep(0, (n.cuts - 1)) for (s in 1:(n.cuts - 1)) { low.bound <- cut.points[s] upp.bound <- cut.points[s + 1] v.dist <- NULL v.dist1 <- NULL v.diff1 <- NULL index.v.dist <- NULL index.v.dist1 <- NULL if (s < (n.cuts - 1)) { v.dist <- s.new.dist.day[s.new.dist.day >= low.bound & s.new.dist.day < upp.bound] index.v.dist <- seq(1:length(s.new.dist.day))[s.new.dist.day >= low.bound & s.new.dist.day < upp.bound] v.dist1 <- v.dist[v.dist != 0] index.v.dist1 <- index.v.dist[v.dist != 0] if (length(v.dist1) >= 1) { v.diff1 <- s.new.diff.day[index.v.dist1] W.new[s] <- sum(v.diff1) W.obs.new[s] <- length(v.dist1) } if (length(v.dist1) == 0) { W.new[s] <- 0 W.obs.new[s] <- 0 } } if (s == (n.cuts - 1)) { v.dist <- s.new.dist.day[s.new.dist.day >= low.bound & s.new.dist.day <= upp.bound] index.v.dist <- seq(1:length(s.new.dist.day))[s.new.dist.day >= low.bound & s.new.dist.day <= upp.bound] v.dist1 <- v.dist[v.dist != 0] index.v.dist1 <- index.v.dist[v.dist != 0] if (length(v.dist1) >= 1) { v.diff1 <- s.new.diff.day[index.v.dist1] W.new[s] <- sum(v.diff1) W.obs.new[s] <- length(v.dist1) } if (length(v.dist1) == 0) { W.new[s] <- 0 } } W <- W + W.new W.obs <- W.obs + W.obs.new } avg.variog <- round(W/(2 * W.obs), 2) n.h <- W.obs x.vec <- NULL for (i in 1:(n.cuts - 1)) { x.vec[i] <- (cut.points[i] + cut.points[i + 1])/2 } fin.avg.variog <- c(0, avg.variog) fin.x.vec <- c(0, x.vec) fin.n.h <- c(n.stations, n.h) B <- list(mar.var = round(marginal.var, 3), bin.midpoints = fin.x.vec, number.pairs = fin.n.h, empir.variog = fin.avg.variog) return(B) } } > ## End Don't show > > > > cleanEx(); ..nameEx <- "avg.variog.dir" > > ### * avg.variog.dir > > flush(stderr()); flush(stdout()) > > ### Name: avg.variog.dir > ### Title: Directional empirical variogram of a random variable averaged > ### over time > ### Aliases: avg.variog.dir > ### Keywords: file > > ### ** Examples > > ## Loading data > data(slp) > day <- slp$date.obs > id <- slp$id.stat > coord1 <- slp$lon.stat > coord2 <- slp$lat.stat > obs <- slp$obs > forecast <- slp$forecast > > ## Computing directional variogram of observed temperature > ## No specified cutpoints, no specified maximum distance > ## No specified tolerance angles and default number of bins > dir.variog <- avg.variog.dir(day,coord1,coord2,id,variable=obs, tol.angle1=NULL,tol.angle2=NULL,cut.points=NULL,max.dist=NULL, nbins=NULL,type='E') > ## Plotting directional variogram > plot(dir.variog$bin.midpoints,dir.variog$dir.variog,xlab="Distance", ylab="Semi-variance",main="Empirical Directional variogram") > > ## Don't show: ## The function is currently defined as > "avg.variog.dir" <- + function(day,coord1,coord2,id,variable,tol.angle1=45,tol.angle2=135,cut.points=NULL,max.dist=NULL,nbins=300,type){ + # default values + if(missing(cut.points)) + cut.points <- NULL + if(missing(max.dist)) + max.dist <- NULL + if(missing(nbins)) + nbins <- NULL + if(missing(type)){ + stop("Invalid input for type: a type should be provided") + } + if(missing(tol.angle1)) + tol.angle1 <- 45 + if(missing(tol.angle2)) + tol.angle2 <- 135 + #INPUT CHECK + # Here we check if the input is right. + l.day <- length(day) + l.var <- length(variable) + l.id <- length(id) + l.coord1 <- length(coord1) + l.coord2 <- length(coord2) + if(sum((c(l.day,l.var,l.id,l.coord1,l.coord2)/l.day)==rep(1,5))!=5){ + stop("Mismatch in dimensions in the data") + } + if(sum(is.numeric(variable)==rep("TRUE",l.var)) 1){ + stop("tol.angle1 should be a number") + } + + if(l.tol.angle1==1 & is.numeric(tol.angle1)==FALSE){ + stop("tol.angle1 should be a number between 0 and 360") + } + if(tol.angle1 <0 | tol.angle1 > 360){ + stop("tol.angle1 should be a number between 0 and 360") + } + + l.tol.angle2 <- length(tol.angle2) + + if(l.tol.angle2==0){ + tol.angle2 <- 135} + + if(l.tol.angle2 > 1){ + stop("tol.angle2 should be a number") + } + + if(l.tol.angle2==1 & is.numeric(tol.angle2)==FALSE){ + stop("tol.angle1 should be a number between 0 and 360") + } + if(tol.angle2 <0 | tol.angle2 > 360){ + stop("tol.angle1 should be a number between 0 and 360") + } + + if(tol.angle2 <= tol.angle1){ + stop("tol.angle2 should be greater than tol.angle1") + } + ## here we check the cutpoints vector + l.cuts <- length(cut.points) + if(l.cuts==1){ + stop("cut.points should be a numeric vector") + } + if(l.cuts>=2 & (sum(is.numeric(cut.points)==rep("TRUE",l.cuts))=2 & (sum(is.numeric(cut.points)==rep("TRUE",l.cuts))=2 & (sum(is.numeric(cut.points)==rep("TRUE",l.cuts))==l.cuts)){ + if(sum(order(cut.points)==seq(1:l.cuts))= 0) < l.cuts){ + stop("Cut points should be non-negative numbers") + } + if(length(cut.points)!=length(unique(cut.points))){ + stop("The vector with cut points should not contain repeated entries") + } + } + ## check on the max.dist + l.mdist <- length(max.dist) + if(l.mdist > 1){ + stop("max.dist is a numeric field, not a vector")} + if(l.mdist==1){ + if(is.numeric(max.dist)==FALSE){ + stop("max.dist is a numeric field") + } + if(max.dist < 0){ + stop("max.dist should be a positive number") + } + } + ## check on the number of bins + l.nbins <- length(nbins) + if(l.nbins==0 & l.cuts==0){ + nbins <- 300 + } + if(l.nbins==1 & l.cuts >=2){ + nbins=NULL + } + l.nbins <- length(nbins) + if(l.nbins >1){ + stop("nbins should be an integer: not a vector") + } + if(l.nbins==1){ + if(floor(nbins)!=nbins){ + stop("Invalid input: the number of bins should be a positive integer") + } + if(ceiling(nbins)==nbins & nbins < 0){ + stop("Invalid input: the number of bins should be a positive integer") + } + } + ## check on the type ## + l.type <- length(type) + if(l.type==0){ + stop("type should be entered") + } + if(type!="E" & type!="N"){ + stop("type can only be equal to E or to N") + } + + # HERE THE REAL FUNCTION STARTS + # here we order the data in ascending date order + day.o <- order(day) + coord1 <- coord1[day.o] + coord2 <- coord2[day.o] + variable <- variable[day.o] + id <- id[day.o] + tol.angle.rad1 <- tol.angle1/57.29577951 + tol.angle.rad2 <- tol.angle2/57.29577951 + + # if the vector with the cutpoints + # is not specified, we determine the cutpoints by looking at the day with the median number + # of observations, we calculate the cutpoints so that the number + # of bins is equal to the one specified and each bin contains approx the same number of pairs. + # If the vector of cutpoints is specified, then we just use that vector of cutpoints. + if(length(cut.points)!=0 & length(max.dist)!=0){ + cut.points <- cut.points[cut.points <= max.dist] + } + if(length(cut.points)==0){ + # all this part is to determine the day with the median number of observations + variable.day <- table(day) + variable.day.vec <- matrix(variable.day,nrow=length(unique(day)),ncol=1) + variable.median <- round(apply(variable.day.vec,2,median),0) + diff.median <- abs(variable.day-variable.median) + unique.day <- unique(day) + day.index <- min(unique.day[diff.median==min(diff.median)]) + # this part is to determine the directional distance among stations that have been recorded in the + # day with the median number of observations and to determine the cutpoints. + variable.day.index <- variable[day==day.index] + id.day.index <- id[day==day.index] + coord1.day.index <- coord1[day==day.index] + coord2.day.index <- coord2[day==day.index] + dist.day.index <- calc.dist.dir(coord1.day.index,coord2.day.index,id.day.index,tol.angle.rad1,tol.angle.rad2,type) + dist.day.index <- dist.day.index[lower.tri(dist.day.index)] + dist.day.index <- dist.day.index[dist.day.index > 0] + ord.dist.day <- sort(dist.day.index) + + # this is in case we want to consider only those pairs of stations where the distance + # is smaller than the maximum distance allowed among locations. + if(length(max.dist)!= 0){ + ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] + } + + if(length(max.dist)==0){ + max.dist <- quantile(ord.dist.day,.9) + ord.dist.day <- ord.dist.day[ord.dist.day < max.dist] + } + + l.dist.day <- length(ord.dist.day) + l.bins <- ceiling(l.dist.day/nbins) + for(i in 1:(nbins-1)){ + cut.points[i] <- ord.dist.day[i*l.bins] + } + + tol.angle.rad1 <- tol.angle1/57.29577951 + tol.angle.rad2 <- tol.angle2/57.29577951 + + unique.day <- unique(day) + l.day <- length(unique.day) + n.stations <- length(unique(id)) + n.cuts <- length(cut.points) + W <- rep(0,(n.cuts-1)) + W.obs <- rep(0,(n.cuts-1)) + + for(i in 1:l.day){ + distance.day <- NULL + difference.day <- NULL + new.dist.day <- NULL + s.new.dist.day <- NULL + o.new.dist.day <- NULL + new.diff.day <- NULL + s.new.diff.day <- NULL + variable.day <- NULL + id.day <- NULL + coord1.day <- NULL + coord2.day <- NULL + variable.day <- variable[day==unique.day[i]] + id.day <- id[day==unique.day[i]] + coord1.day <- coord1[day==unique.day[i]] + coord2.day <- coord2[day==unique.day[i]] + distance.day <- + calc.dist.dir(coord1.day,coord2.day,id.day,tol.angle.rad1,tol.angle.rad2,type) + new.dist.day <- distance.day[lower.tri(distance.day)] + s.new.dist.day <- sort(new.dist.day) + o.new.dist.day <- order(new.dist.day) + difference.day <- calc.difference(variable.day) + new.diff.day <- difference.day[lower.tri(difference.day)] + s.new.diff.day <- new.diff.day[o.new.dist.day] + W.new <- rep(0,(n.cuts-1)) + W.obs.new <- rep(0,(n.cuts-1)) + + for(s in 1:(n.cuts-1)){ + low.bound <- cut.points[s] + upp.bound <- cut.points[s+1] + v.dist <- NULL + v.dist1 <- NULL + v.diff1 <- NULL + index.v.dist <- NULL + index.v.dist1 <- NULL + + if(s < (n.cuts-1)){ + low.bound <- cut.points[s] + upp.bound <- cut.points[s+1] + v.dist <- NULL + v.dist1 <- NULL + v.diff1 <- NULL + index.v.dist <- NULL + index.v.dist1 <- NULL + + if(s < (n.cuts-1)){ + v.dist <- s.new.dist.day[s.new.dist.day >= low.bound & + s.new.dist.day < upp.bound] + index.v.dist <- + seq(1:length(s.new.dist.day))[s.new.dist.day >= low.bound & s.new.dist.day + < upp.bound] + v.dist1 <- v.dist[v.dist!=0] + index.v.dist1 <- index.v.dist[v.dist!=0] + + if(length(v.dist1) >=1){ + v.diff1 <- s.new.diff.day[index.v.dist1] + W.new[s] <- sum(v.diff1) + W.obs.new[s] <- length(v.dist1)} + if(length(v.dist1) ==0){ + W.new[s] <- 0 + W.obs.new[s] <- 0} + } + if(s==(n.cuts-1)){ + v.dist <- s.new.dist.day[s.new.dist.day >= low.bound & + s.new.dist.day <= upp.bound] + index.v.dist <- + seq(1:length(s.new.dist.day))[s.new.dist.day >= + low.bound & s.new.dist.day <= upp.bound] + v.dist1 <- v.dist[v.dist!=0] + index.v.dist1 <- index.v.dist[v.dist!=0] + + if(length(v.dist1) >=1){ + v.diff1 <- s.new.diff.day[index.v.dist1] + W.new[s] <- sum(v.diff1) + W.obs.new[s] <- length(v.dist1)} + if(length(v.dist1) ==0){ + W.new[s] <- 0 + W.obs.new[s] <- 0} + } + } + W <- W+W.new + W.obs <- W.obs + W.obs.new + } + avg.variog <- round(W/(2*W.obs),2) + n.h <- W.obs + x.vec <- NULL + for(i in 1:(n.cuts-1)){ + x.vec[i] <- (cut.points[i]+cut.points[i+1])/2 + } + fin.avg.variog <- c(0,avg.variog) + fin.x.vec <- c(0,x.vec) + fin.n.h <- c(n.stations,n.h) + B <- list(bin.midpoints=fin.x.vec,number.pairs=fin.n.h,dir.variog=fin.avg.variog) + return(B) + } + } + } > ## End Don't show > > > > cleanEx(); ..nameEx <- "linesmodel" > > ### * linesmodel > > flush(stderr()); flush(stdout()) > > ### Name: linesmodel > ### Title: Computation of parametric variogram model > ### Aliases: linesmodel > ### Keywords: file > > ### ** Examples > > ## Loading data > data(slp) > day <- slp$date.obs > id <- slp$id.stat > coord1 <- slp$lon.stat > coord2 <- slp$lat.stat > obs <- slp$obs > forecast <- slp$forecast > > ## Computing empirical variogram > variogram <- Emp.variog(day=day,obs=obs,forecast=forecast,id=id,coord1=coord1, + coord2=coord2,cut.points=NULL,max.dist=NULL,nbins=NULL) > > ## Estimating variogram parameters > ## Without specifying initial values for the parameters > param.variog <- + Variog.fit(emp.variog=variogram,variog.model="exponential",max.dist.fit=NULL, + init.val=NULL,fix.nugget=FALSE) > > ## Plotting the empirical variogram with the estimated parametric variogram superimposed > plot(variogram$bin.midpoints,variogram$empir.variog,xlab="Distance",ylab="Semi-variance") > lines(variogram$bin.midpoints,linesmodel(distance=variogram$bin.midpoints,variog.model="exponential",param=c(param.variog$nugget, + param.variog$variance,param.variog$range))) > > ## Don't show: > ## The function is currently defined as > function(distance,variog.model="exponential",param){ + # INPUT CHECK + if(missing(distance)){ + stop("Invalid input") + } + + if(missing(variog.model)){ + stop("Invalid input") + } + + if(missing(param)){ + stop("Invalid input") + } + + # Here we check if the vector of distances is well-defined + l.dist <- length(distance) + if(l.dist <1){ + stop("Invalid input") + } + + if(l.dist>=1 & (sum(is.numeric(distance)==rep("TRUE",l.dist)) < l.dist)){ + stop("distance should be a numeric vector") + } + + if(l.dist>=1 & (sum(is.numeric(distance)==rep("TRUE",l.dist))==l.dist)){ + if(sum(distance >=0) < l.dist){ + stop("Entries in distance should be non-negative numbers") + } + if(sum(order(distance)==seq(1:l.dist))2) | param[5] <= 0){ + stop("Invalid input for at least one of the parameters") + } + } + } + + + # This is where the function really starts + + if(variog.model=="spherical"){ + fitted.values <- (param[1]+param[2]*(1.5*(distance/param[3])-0.5*(distance/param[3])^3))} + + if(variog.model=="gauss"){ + fitted.values <- (param[1]+param[2]*(1-exp(-(distance/param[3])^2)))} + + if(variog.model=="exponential"){ + fitted.values <- (param[1]+param[2]*(1-exp(-distance/param[3])))} + + if(variog.model=="whittlematern" | variog.model=="matern"){ + fitted.values <- + (param[1]+param[2]-param[2]*(((2^(1-param[4]))/gamma(param[4]))*((distance/param[3])^(param[4]))*besselK(distance/param[3],param[4],expon.scaled=FALSE)))} + + if(variog.model=="gencauchy"){ + fitted.values <- (param[1]+param[2]*(1-(1+((distance^2)/param[3])^param[4])^(-param[5]/param[4])))} + + return(fitted.values) + } function (distance, variog.model = "exponential", param) { if (missing(distance)) { stop("Invalid input") } if (missing(variog.model)) { stop("Invalid input") } if (missing(param)) { stop("Invalid input") } l.dist <- length(distance) if (l.dist < 1) { stop("Invalid input") } if (l.dist >= 1 & (sum(is.numeric(distance) == rep("TRUE", l.dist)) < l.dist)) { stop("distance should be a numeric vector") } if (l.dist >= 1 & (sum(is.numeric(distance) == rep("TRUE", l.dist)) == l.dist)) { if (sum(distance >= 0) < l.dist) { stop("Entries in distance should be non-negative numbers") } if (sum(order(distance) == seq(1:l.dist)) < l.dist) { stop("Entries in distance should be in ascending order") } } l.variog <- length(variog.model) case.var <- 0 if (l.variog == 0) { variog.model <- "exponential" case.var <- 1 } if (l.variog == 1) { if (is.character(variog.model) == "TRUE") { if (variog.model == "exponential") { case.var <- 1 } if (variog.model == "spherical") { case.var <- 1 } if (variog.model == "whittlematern" | variog.model == "matern") { case.var <- 1 } if (variog.model == "gencauchy") { case.var <- 1 } if (variog.model == "gauss") { case.var <- 1 } } } if (variog.model == "matern") { variog.model <- "whittlematern" } if (case.var == 0) { stop("Incorrect variogram model specification") } l.par <- length(param) if (l.par < 3) { stop("Invalid input: the parameter should be a vector of length at least equal to 3") } if (variog.model == "exponential" | variog.model == "spherical" | variog.model == "gauss") { if (l.par != 3) { stop("The parameter vector should have length equal to 3") } if (l.par == 3) { if (param[1] < 0 | param[2] <= 0 | param[3] <= 0) { stop("Invalid input for at least one of the parameters") } } } if (variog.model == "whittlematern") { if (l.par != 4) { stop("The parameter vector should have length equal to 4") } if (l.par == 4) { if (param[1] < 0 | param[2] <= 0 | param[3] <= 0 | param[4] <= 0) { stop("Invalid input for at least one of the parameters") } } } if (variog.model == "gencauchy") { if (l.par != 5) { stop("The parameter vector should have length equal to 5") } if (l.par == 5) { if (param[1] < 0 | param[2] <= 0 | param[3] <= 0 | (param[4] <= 0 & param[4] > 2) | param[5] <= 0) { stop("Invalid input for at least one of the parameters") } } } if (variog.model == "spherical") { fitted.values <- (param[1] + param[2] * (1.5 * (distance/param[3]) - 0.5 * (distance/param[3])^3)) } if (variog.model == "gauss") { fitted.values <- (param[1] + param[2] * (1 - exp(-(distance/param[3])^2))) } if (variog.model == "exponential") { fitted.values <- (param[1] + param[2] * (1 - exp(-distance/param[3]))) } if (variog.model == "whittlematern" | variog.model == "matern") { fitted.values <- (param[1] + param[2] - param[2] * (((2^(1 - param[4]))/gamma(param[4])) * ((distance/param[3])^(param[4])) * besselK(distance/param[3], param[4], expon.scaled = FALSE))) } if (variog.model == "gencauchy") { fitted.values <- (param[1] + param[2] * (1 - (1 + ((distance^2)/param[3])^param[4])^(-param[5]/param[4]))) } return(fitted.values) } > ## End Don't show > > > > cleanEx(); ..nameEx <- "plotfields" > > ### * plotfields > > flush(stderr()); flush(stdout()) > > ### Name: plotfields > ### Title: Plot of weather fields > ### Aliases: plotfields > ### Keywords: file > > ### ** Examples > > ## Loading data > library(fields) > library(RandomFields) > data(slp) > data(gridlong) > data(gridlat) > data(forecast.grid) > day <- slp$date.obs > id <- slp$id.stat > coord1 <- slp$lon.stat > coord2 <- slp$lat.stat > obs <- slp$obs > forecast <- slp$forecast > coord1.grid <- gridlong$gridded.long > coord2.grid <- gridlat$gridded.lat > forecast.grid <- forecast.grid$gridded.forecast > > ## Computing the empirical variogram > variogram <- Emp.variog(day,obs,forecast,id,coord1,coord2,cut.points=NULL, + max.dist=NULL,nbins=NULL) > > ## Estimating parameters > param.est <- Variog.fit(variogram,"exponential",max.dist.fit=NULL, + init.val=NULL,fix.nugget=FALSE) > > ## Simulating realizations of the weather random field > simul <- Field.sim(obs, forecast, coord1.grid, coord2.grid, forecast.grid, + variog.model="exponential", param.est=c(param.est$nugget,param.est$variance, + param.est$range), n.sim=4, n.displ=0, qt.displ=c(10,50,90)) Hit to see next plot: > ##Plotting one of the simulated weather random fields > par(mfrow=c(1,1)) > plotfields(simul$sim.fields[,,1],x.lim=c(min(coord1.grid),max(coord1.grid)), + y.lim=c(min(coord2.grid),max(coord2.grid)),title="Simulated weather field") > > ## Plotting one of the percentiles of the weather field > par(mfrow=c(1,1)) > plotfields(simul$pct.fields[,,1],x.lim=c(min(coord1.grid),max(coord1.grid)), + y.lim=c(min(coord2.grid),max(coord2.grid)),title="10th percentile") > > ## Don't show: > ## The function is currently defined as > function(field,x.lim,y.lim,title){ + + # Input check + dim.f <- dim(field) + + if(length(dim.f)==0){ + stop("Invalid input for field") + } + + if(length(dim.f)!=2){ + stop("Invalid input for field") + } + + if(dim.f[1]!=65 | dim.f[2]!=65){ + stop("Invalid input for field") + } + + dim.f2 <- dim.f^2 + + if(sum(is.numeric(field))==0){ + stop("Invalid input for field") + } + + ## input check for the x vector + l.limx <- length(x.lim) + if(l.limx < 2){ + stop("Invalid input for x.lim") + } + + if(l.limx > 2){ + stop("Invalid input for x.lim") + } + + if(l.limx==2){ + if(sum(is.numeric(x.lim)==rep(TRUE,2))<2){ + stop("Invalid input for x.lim") + } + if(x.lim[1] >= x.lim[2]){ + stop("Invalid input for x.lim") + } + } + + ## input check on the lims for the y vector + l.limy <- length(y.lim) + if(l.limy < 2){ + stop("Invalid input for y.lim") + } + + if(l.limy > 2){ + stop("Invalid input for y.lim") + } + + if(l.limy==2){ + if(sum(is.numeric(y.lim)==rep(TRUE,2))<2){ + stop("Invalid input for y.lim") + } + if(y.lim[1] >= y.lim[2]){ + stop("Invalid input for y.lim") + } + } + + lims <- c(min(field,na.rm=TRUE),max(field,na.rm=TRUE)) + size=65 + hor.crd <- seq(x.lim[1],x.lim[2],length = size) + ver.crd <- seq(y.lim[1],y.lim[2],length = size) + n.level <- 5 + US.map <- 0 + par(ask=FALSE) + image.plot(hor.crd, ver.crd, field, xlim=c(min(hor.crd), max(hor.crd)), ylim = c(min(ver.crd),max(ver.crd)),zlim=lims,legend.width=.015,legend.shrink=.8, + main=title,xlab="",ylab="",offset = 0.02, col=rainbow(100,start=0,end=0.85)[100:1]) + + if(min(hor.crd) <= -124.7 & max(hor.crd) >= -124.7){ + US.map <- 1 + } + + if(min(hor.crd) >= -124.7 & max(hor.crd) <= -67.1){ + US.map <- 1 + } + + if(min(hor.crd) <= -67.1 & max(hor.crd) >= -67.1){ + US.map <- 1 + } + + if(min(ver.crd) <= 25.2 & max(ver.crd) >= 25.2){ + US.map <- 1 + US.map + } + + if(min(ver.crd) >= 25.2 & max(ver.crd) <= 49.4){ + US.map <- 1 + US.map + } + + if(min(ver.crd) <= 49.4 & max(ver.crd) >= 49.4){ + US.map <- 1 + US.map + } + + if(US.map==2){ + US(xlim=c(min(hor.crd), max(hor.crd)), ylim =c(min(ver.crd), max(ver.crd)),add=TRUE,col=1,lwd=2) + world(xlim=c(min(hor.crd), max(hor.crd)), ylim =c(min(ver.crd), max(ver.crd)),add=TRUE,col=1,lwd=2) + } + + if(US.map < 2){ + world(xlim=c(min(hor.crd), max(hor.crd)), ylim =c(min(ver.crd), max(ver.crd)),add=TRUE,col=1,lwd=2) + } + + } function (field, x.lim, y.lim, title) { dim.f <- dim(field) if (length(dim.f) == 0) { stop("Invalid input for field") } if (length(dim.f) != 2) { stop("Invalid input for field") } if (dim.f[1] != 65 | dim.f[2] != 65) { stop("Invalid input for field") } dim.f2 <- dim.f^2 if (sum(is.numeric(field)) == 0) { stop("Invalid input for field") } l.limx <- length(x.lim) if (l.limx < 2) { stop("Invalid input for x.lim") } if (l.limx > 2) { stop("Invalid input for x.lim") } if (l.limx == 2) { if (sum(is.numeric(x.lim) == rep(TRUE, 2)) < 2) { stop("Invalid input for x.lim") } if (x.lim[1] >= x.lim[2]) { stop("Invalid input for x.lim") } } l.limy <- length(y.lim) if (l.limy < 2) { stop("Invalid input for y.lim") } if (l.limy > 2) { stop("Invalid input for y.lim") } if (l.limy == 2) { if (sum(is.numeric(y.lim) == rep(TRUE, 2)) < 2) { stop("Invalid input for y.lim") } if (y.lim[1] >= y.lim[2]) { stop("Invalid input for y.lim") } } lims <- c(min(field, na.rm = TRUE), max(field, na.rm = TRUE)) size = 65 hor.crd <- seq(x.lim[1], x.lim[2], length = size) ver.crd <- seq(y.lim[1], y.lim[2], length = size) n.level <- 5 US.map <- 0 par(ask = FALSE) image.plot(hor.crd, ver.crd, field, xlim = c(min(hor.crd), max(hor.crd)), ylim = c(min(ver.crd), max(ver.crd)), zlim = lims, legend.width = 0.015, legend.shrink = 0.8, main = title, xlab = "", ylab = "", offset = 0.02, col = rainbow(100, start = 0, end = 0.85)[100:1]) if (min(hor.crd) <= -124.7 & max(hor.crd) >= -124.7) { US.map <- 1 } if (min(hor.crd) >= -124.7 & max(hor.crd) <= -67.1) { US.map <- 1 } if (min(hor.crd) <= -67.1 & max(hor.crd) >= -67.1) { US.map <- 1 } if (min(ver.crd) <= 25.2 & max(ver.crd) >= 25.2) { US.map <- 1 + US.map } if (min(ver.crd) >= 25.2 & max(ver.crd) <= 49.4) { US.map <- 1 + US.map } if (min(ver.crd) <= 49.4 & max(ver.crd) >= 49.4) { US.map <- 1 + US.map } if (US.map == 2) { US(xlim = c(min(hor.crd), max(hor.crd)), ylim = c(min(ver.crd), max(ver.crd)), add = TRUE, col = 1, lwd = 2) world(xlim = c(min(hor.crd), max(hor.crd)), ylim = c(min(ver.crd), max(ver.crd)), add = TRUE, col = 1, lwd = 2) } if (US.map < 2) { world(xlim = c(min(hor.crd), max(hor.crd)), ylim = c(min(ver.crd), max(ver.crd)), add = TRUE, col = 1, lwd = 2) } } > ## End Don't show > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > ### *