kriging is working

jobonaf · Mar 25, 2015 · cb1c241 · cb1c241
1 parent 479a011
commit cb1c241
Show file tree

Hide file tree

Showing 8 changed files with 455 additions and 23 deletions.
diff --git a/R/kriging.R b/R/kriging.R
@@ -1,4 +1,4 @@
-kriging <- function (x.pnt, y.pnt, obs, x.grd, y.grd, 
+kriging <- function (x.pnt, y.pnt, obs,
                      model,
                      proxy.1, 
                      proxy.2=NULL, 
@@ -8,11 +8,31 @@ kriging <- function (x.pnt, y.pnt, obs, x.grd, y.grd,
                      K.pairs.min=1) {
   # remove points from the target grid
   # where model or any proxy is missing
-  idx <- which(!is.na(model$grid$z) & !is.na(proxy.1$grid$z) & !is.na(proxy.2$grid$z))
-  model$grid$z
+  clean.grids <- function(objs) {
+    for (i in 1:length(objs)) {
+      idx <- !is.na(objs[[i]]$grid$z)
+      if(1==1) {Idx <- idx} else {Idx <- Idx&idx}
+    }
+    for (i in 1:length(objs)) {
+      objs[[i]]$grid$x <- objs[[i]]$grid$x[which(Idx)]
+      objs[[i]]$grid$y <- objs[[i]]$grid$y[which(Idx)]
+      objs[[i]]$grid$z <- objs[[i]]$grid$z[which(Idx)]
+    }
+    return(objs)
+  }
+  if(is.null(proxy.2)) {
+    out <- clean.grids(list(model=model, proxy1=proxy.1))
+    model <- out$model
+    proxy.1 <- out$proxy.1
+  } else {
+    out <- clean.grids(list(model=model, proxy.1=proxy.1, proxy.2=proxy.2))
+    model <- out$model
+    proxy.1 <- out$proxy.1
+    proxy.2 <- out$proxy.2
+  }
 
   ### Trends
-  if(is.null(proxy2)){
+  if(is.null(proxy.2)){
     trend.d <- ~ boxcox(model$points$z,lambda) + proxy.1$points$z 
     trend.l <- ~ boxcox(model$grid$z,lambda)   + proxy.1$grid$z   
   } else {
@@ -48,7 +68,7 @@ kriging <- function (x.pnt, y.pnt, obs, x.grd, y.grd,
 
   # Kriging
   kkk<-krige.conv(geodata,
-                  loc=cbind(x.grd,y.grd),
+                  locations=cbind(model$grid$x, model$grid$y),
                   krige= krige.control(
                     type.krige="OK",
                     dist.epsilon= max(K.min.dist,1e-10),
@@ -61,7 +81,9 @@ kriging <- function (x.pnt, y.pnt, obs, x.grd, y.grd,
 
   # Correction
   out <- list(K     = kkk$predict, 
-              K.var = kkk$krige.var)
+              K.var = kkk$krige.var,
+              x     = model$grid$x,
+              y     = model$grid$y)
   return(out)
 }
 

diff --git a/R/prepare.R b/R/prepare.R
@@ -3,7 +3,8 @@
 
 ## CTM
 prepare.ctm <- function(ctm.daily, day,
-                        x.pnt, y.pnt, x.grd, y.grd) {
+                        x.pnt, y.pnt, x.grd, y.grd,
+                        conc.min=10^-6) {
   time <- sort(unique(as.POSIXct(as.character(ctm.daily$time),tz="Africa/Algiers")))
   if(min(diff.POSIXt(time)) < as.difftime(1,units="days")) {
     stop("Cannot deal with sub-daily data.")
@@ -25,10 +26,12 @@ prepare.ctm <- function(ctm.daily, day,
                     xp=x.grd, yp=y.grd,
                     type="grid",method="linear")
   dum <- expand.grid(ctm.grd$x, ctm.grd$y, KEEP.OUT.ATTRS=F)
+  ctm.pnt$z <- pmax(ctm.pnt$z, conc.min)
+  ctm.grd$z <- pmax(as.vector(ctm.grd$z), conc.min)
   out <- list(points=ctm.pnt,
               grid=list(x=dum[[1]],
                         y=dum[[2]],
-                        z=as.vector(ctm.grd$z)))
+                        z=ctm.grd$z))
   return(out)
 }
 
@@ -83,7 +86,8 @@ prepare.elev <- function(elev,
 }
 
 ## Observed data
-prepare.obs <- function(obs.daily, day) {
+prepare.obs <- function(obs.daily, day,
+                        conc.min=10^-6) {
   time <- sort(unique(as.POSIXct(as.character(obs.daily$Time),tz="Africa/Algiers")))
   if(min(diff.POSIXt(time)) < as.difftime(1,units="days")) {
     stop("Cannot deal with sub-daily data.")
@@ -92,6 +96,7 @@ prepare.obs <- function(obs.daily, day) {
                 subset=(format(Time,format="%Y-%m-%d")==
                           format(as.POSIXct(day),format="%Y-%m-%d")))
   out <- out[which(!is.na(out[,2])),]
+  out[,2] <- pmax(out[,2], conc.min)
   rownames(out) <- 1:nrow(out)
   return(out)
 }

diff --git a/TODO.md b/TODO.md
@@ -15,7 +15,7 @@ TODO
 * ~~read data~~
 * ~~prepare data~~
 * ~~kriging~~
-* test the kriging
+* ~~test the kriging~~
 * ~~all-inclusive function to prepare all you need to process a day~~
 * write the output of kriging
 * build source package
@@ -31,5 +31,5 @@ TODO
 * manual (*.Rd)
 * ~~demo for reading data~~
 * ~~demo for preparing data~~
-* demo for kriging
+* ~~demo for kriging~~
 * demo for cross-validation
diff --git a/demo/00Index b/demo/00Index
@@ -1,2 +1,3 @@
 daily.synthesis      How to prepare daily data
 prepare.day          Prepare the data to perform the kriging on a specific day
+kriging              Kriging
diff --git a/demo/kriging.R b/demo/kriging.R
@@ -0,0 +1,88 @@
+#' ---
+#' title: "Demo: How to perform the kriging on a specific day"
+#' author: "Giovanni Bonafe'"
+#' date: "March 25rd, 2015"
+#' ---
+
+#' ## How to perform the kriging on a specific day
+#' You can reproduce the following R code with 
+#' ```demo(kriging)``` after loading package ```pesco```. Note that
+#' some of the commands used here have been already used in ```demo(prepare.day)```.
+
+## load package
+require(pesco)
+
+## load data
+data(PM10.obs)
+data(emissions)
+data(PM10.ctm)
+data(elevation)
+
+## prepare data for the day
+myDay <- "2015-03-02"
+PM10.ctm.ave <- dailyCtm(PM10.ctm, statistic  =  "mean")
+dataDay <- prepare.day(day = myDay,
+                       obs.daily = PM10.obs,
+                       ctm.daily = PM10.ctm.ave,
+                       emis.winter = emissions$PM10.winter,
+                       emis.summer = emissions$PM10.summer,
+                       elev = elevation,
+                       verbose = TRUE)
+
+## kriging
+K <- kriging (x.pnt = dataDay$ctm.day$points$x, 
+              y.pnt = dataDay$ctm.day$points$y, 
+              obs = dataDay$obs.day$PM10, 
+              model = dataDay$ctm.day,
+              proxy.1 = dataDay$emis.day, 
+              proxy.2 = dataDay$elev.day, 
+              lambda = 0)
+
+#' Let's see the results on a map.
+
+bb <- round(quantile(c(K$K, dataDay$ctm.day$grid$z),
+                     probs=seq(0,1,length.out=11)))
+cc <- tim.colors(10)
+ll <- paste(bb[-length(bb)],bb[-1],sep="-")
+xr <- range(dataDay$ctm.day$grid$x)
+yr <- range(dataDay$ctm.day$grid$y)
+
+## CTM
+plot(dataDay$ctm.day$grid$x, dataDay$ctm.day$grid$y, 
+     main="interpolated CTM",
+     col=cc[cut(dataDay$ctm.day$grid$z,breaks=bb)], 
+     pch=19, cex=0.4, xlab="x UTM", ylab="y UTM")
+legend("bottomleft", fill=rev(cc), legend=rev(ll), 
+       bg="white", y.intersp=0.8)
+
+## observations
+plot(dataDay$ctm.day$points$x, dataDay$ctm.day$points$y, 
+     main="observations",
+     col=cc[cut(dataDay$obs.day$PM10,breaks=bb)],
+     pch=19, cex=1, xlim=xr, ylim=yr, 
+     xlab="x UTM", ylab="y UTM")
+points(dataDay$ctm.day$points$x, dataDay$ctm.day$points$y, 
+       col="black", pch=1, cex=1)
+legend("bottomleft", fill=rev(cc), legend=rev(ll), 
+       bg="white", y.intersp=0.8)
+
+## kriging: predicted values
+plot(K$x, K$y, main="kriging: predicted values",
+     col=cc[cut(K$K,breaks=bb)], 
+     pch=19, cex=0.4, xlim=xr, ylim=yr, 
+     xlab="x UTM", ylab="y UTM")
+legend("bottomleft", fill=rev(cc), legend=rev(ll), 
+       bty="n", y.intersp=0.8)
+
+## kriging: predicted variances
+bb <- round(quantile(K$K.var,
+                     probs=seq(0,1,length.out=11)))
+cc <- designer.colors(10)
+ll <- paste(bb[-length(bb)],bb[-1],sep="-")
+plot(K$x, K$y, main="kriging: predicted variances",
+     col=cc[cut(K$K.var,breaks=bb)], 
+     pch=19, cex=0.4, xlim=xr, ylim=yr, 
+     xlab="x UTM", ylab="y UTM")
+legend("bottomleft", fill=rev(cc), legend=rev(ll), 
+       bty="n", y.intersp=0.8)
+
diff --git a/demo/prepare.day.R b/demo/prepare.day.R
@@ -87,12 +87,12 @@ signif(data.frame(Obs=PM10.obs.day$PM10,
 #' The same result can be achieved with one single function.
 
 dataDay <- prepare.day(day=myDay,
-            obs.daily=PM10.obs,
-            ctm.daily=PM10.ctm.ave,
-            emis.winter=emissions$PM10.winter,
-            emis.summer=emissions$PM10.summer,
-            elev=elevation,
-            verbose=TRUE)
+                       obs.daily=PM10.obs,
+                       ctm.daily=PM10.ctm.ave,
+                       emis.winter=emissions$PM10.winter,
+                       emis.summer=emissions$PM10.summer,
+                       elev=elevation,
+                       verbose=TRUE)
 
 ## data at the stations points
 signif(data.frame(Obs=dataDay$obs.day$PM10, 

diff --git a/doc/03_PrepareDay.html b/doc/03_PrepareDay.html
@@ -315,12 +315,12 @@ <h2>How to prepare the data to perform the kriging on a specific day</h2>
 <p>The same result can be achieved with one single function.</p>
 
 <pre><code class="r">dataDay &lt;- prepare.day(day=myDay,
-            obs.daily=PM10.obs,
-            ctm.daily=PM10.ctm.ave,
-            emis.winter=emissions$PM10.winter,
-            emis.summer=emissions$PM10.summer,
-            elev=elevation,
-            verbose=TRUE)
+                       obs.daily=PM10.obs,
+                       ctm.daily=PM10.ctm.ave,
+                       emis.winter=emissions$PM10.winter,
+                       emis.summer=emissions$PM10.summer,
+                       elev=elevation,
+                       verbose=TRUE)
 </code></pre>
 
 <pre><code>## [1] &quot;Prepared observations for day 2015-03-02&quot;