From plotKML 0.4-1 is now possible to plot results of 'krigeST'
(http://www.inside-r.org/node/176602) directly in Google Earth. This is
an example:
library(plotKML)
library(gstat)
library(sp)
library(spacetime)
sumMetricVgm <- vgmST("sumMetric",
space=vgm( 4.4, "Lin", 196.6, 3),
time =vgm( 2.2, "Lin", 1.1, 2),
joint=vgm(34.6, "Exp", 136.6, 12),
stAni=51.7)
data(air)
rr <- rural[,"2005-06-01/2005-06-03"]
rr <- as(rr,"STSDF")
x1 <- seq(from=6,to=15,by=1)
x2 <- seq(from=48,to=55,by=1)
DE_gridded <- SpatialPoints(cbind(rep(x1,length(x2)),
rep(x2,each=length(x1))),
proj4string=CRS(proj4string(rr@sp)))
gridded(DE_gridded) <- TRUE
DE_pred <- STF(sp=as(DE_gridded,"SpatialPoints"), time=rr@time)
DE_kriged <- krigeST(PM10~1, data=rr, newdata=DE_pred,
modelList=sumMetricVgm)
gridded(DE_kriged@sp) <- TRUE
stplot(DE_kriged)
## plot in Google Earth:
png.width = DE_kriged@sp@g...@cells.dim[1]*20
png.height = DE_kriged@sp@g...@cells.dim[2]*20
z.lim = range(DE_kriged@data, na.rm=TRUE)
plotKML(DE_kriged, png.width=png.width,
png.height=png.height, z.lim=z.lim)
## add observations points:
rr.sti <- as(rr, "STIDF")
plotKML(rr.sti, z.lim=z.lim)
the output looks like this:
http://plotkml.r-forge.r-project.org/DE_kriged.kmz
Note that I had to set the PNG width and height 20 time larger because
otherwise plots in Google Earth are very fuzzy.
T. Hengl
https://www.vcard.wur.nl/Views/Profile/View.aspx?id=37263
On 1-1-2014 2:55, Hodgess, Erin wrote:
Hello again! Happy New Year!
Here is the solution for this particular situation. Note: thanks to many
people for the help. Anyhow, I started with looking at page 28 in the gstat
vignette for krigeST. If you work through that, you obtain an object called
DE_kriged, which is an STFDF.
I copied some code from the following website:
(https://r-forge.r-project.org/scm/viewvc.php/pkg/R/coerce.R?view=markup&root=spacetime
<https://r-forge.r-project.org/scm/viewvc.php/pkg/R/coerce.R?view=markup&root=spacetime&pathrev=25>)
to get the as.STIDF.STFDF lines:
# STFDF -> STIDF
as.STIDF.STFDF = function(from) {
as(as(from, "STSDF"), "STIDF")
}
setAs("STFDF", "STIDF", as.STIDF.STFDF)
Now:
DE1 <- as.STIDF.STFDF(DE_kriged)
DE1 is an STIDF object, which is good, but the @sp section is a Spatial Pixels
object.
Next I copied over the function kml_layer.STIDF to new_STIDF.R and added the
following section (see the note)
new.STIDF <-
function (obj, dtime = "", ...)
{
if (all(dtime == 0)) {
TimeSpan.begin = format(time(obj@time), "%Y-%m-%dT%H:%M:%SZ")
TimeSpan.end = TimeSpan.begin
}
else {
if (length(obj@time) > 1 & !nzchar(dtime)) {
print(obj@time)
period <- periodicity(obj@time)
dtime <- period$frequency
}
TimeSpan.begin <-
format(as.POSIXct(unclass(as.POSIXct(time(obj@time))) -
dtime/2, origin = "1970-01-01"), "%Y-%m-%dT%H:%M:%SZ")
TimeSpan.end <- format(as.POSIXct(unclass(as.POSIXct(time(obj@time))) +
dtime/2, origin = "1970-01-01"), "%Y-%m-%dT%H:%M:%SZ")
}
if (class(obj@sp) == "SpatialPoints" | class(obj@sp) ==
"SpatialPointsDataFrame") {
sp <- SpatialPointsDataFrame(obj@sp, obj@data)
kml_layer.SpatialPoints(obj = sp, TimeSpan.begin = TimeSpan.begin,
TimeSpan.end = TimeSpan.end, ...)
}
else {
if (class(obj@sp) == "SpatialPolygons" | class(obj@sp) ==
"SpatialPolygonsDataFrame") {
sp <- SpatialPolygonsDataFrame(obj@sp, obj@data)
kml_layer.SpatialPolygons(obj = sp, TimeSpan.begin =
TimeSpan.begin,
TimeSpan.end = TimeSpan.end, ...)
}
else {
if (class(obj@sp) == "SpatialLines" | class(obj@sp) ==
"SpatialLinesDataFrame") {
sp <- SpatialLinesDataFrame(obj@sp, obj@data)
kml_layer.SpatialLines(obj = sp, TimeSpan.begin =
TimeSpan.begin,
TimeSpan.end = TimeSpan.end, ...)
}
###########################################################
# New for Spatial Pixels #
##########################################################
else {
if (class(obj@sp) == "SpatialPixels" | class(obj@sp) ==
"SpatialPixelsDataFrame") {
sp <- SpatialPixelsDataFrame(obj@sp, obj@data)
kml_layer.SpatialPoints(obj = sp, TimeSpan.begin =
TimeSpan.begin,
TimeSpan.end = TimeSpan.end, ...)
}
else {
stop("The STIDF object does not extend SpatialPoints*,
SpatialLines* or SpatialPolygons*")
}
}
}
}
}
Finally, I ran this:
library(plotKML)
plotKML version 0.4-0 (2013-11-15)
URL: http://plotkml.r-forge.r-project.org/
Warning message:
replacing previous import by ‘zoo::as.zoo’ when loading ‘gstat’
library(gstat)
Loading required package: sp
kml_open("stuff2.kml")
KML file opened for writing...
new.STIDF(DE1,dtime=24*3600,colour=var1.pred)
Parsing to KML...
kml_close("stuff2.kml")
Closing stuff2.kml
All is well. Google Earth lets you run the kml file and show the time change,
as it should.
Hope this might help someone!
Thanks,
Erin
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
