[R-sig-Geo] problems with plotting STFDF

[Tomislav Hengl]

 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 at sp)))
gridded(DE_gridded) <- TRUE
DE_pred <- STF(sp=as(DE_gridded,"SpatialPoints"), time=rr at time)
DE_kriged <- krigeST(PM10~1, data=rr, newdata=DE_pred,
                      modelList=sumMetricVgm)
gridded(DE_kriged at sp) <- TRUE
stplot(DE_kriged)

## plot in Google Earth:
png.width = DE_kriged at sp@grid at cells.dim[1]*20
png.height = DE_kriged at sp@grid at cells.dim[2]*20
z.lim = range(DE_kriged at 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 at time), "%Y-%m-%dT%H:%M:%SZ")
>          TimeSpan.end = TimeSpan.begin
>      }
>      else {
>          if (length(obj at time) > 1 & !nzchar(dtime)) {
>          print(obj at time)
>              period <- periodicity(obj at time)
>              dtime <- period$frequency
>          }
>          TimeSpan.begin <- format(as.POSIXct(unclass(as.POSIXct(time(obj at time))) -
>              dtime/2, origin = "1970-01-01"), "%Y-%m-%dT%H:%M:%SZ")
>          TimeSpan.end <- format(as.POSIXct(unclass(as.POSIXct(time(obj at time))) +
>              dtime/2, origin = "1970-01-01"), "%Y-%m-%dT%H:%M:%SZ")
>      }
>      if (class(obj at sp) == "SpatialPoints" | class(obj at sp) == "SpatialPointsDataFrame") {
>          sp <- SpatialPointsDataFrame(obj at sp, obj at data)
>          kml_layer.SpatialPoints(obj = sp, TimeSpan.begin = TimeSpan.begin,
>              TimeSpan.end = TimeSpan.end, ...)
>      }
>      else {
>          if (class(obj at sp) == "SpatialPolygons" | class(obj at sp) ==
>              "SpatialPolygonsDataFrame") {
>              sp <- SpatialPolygonsDataFrame(obj at sp, obj at data)
>              kml_layer.SpatialPolygons(obj = sp, TimeSpan.begin = TimeSpan.begin,
>                  TimeSpan.end = TimeSpan.end, ...)
>          }
>          else {
>              if (class(obj at sp) == "SpatialLines" | class(obj at sp) ==
>                  "SpatialLinesDataFrame") {
>                  sp <- SpatialLinesDataFrame(obj at sp, obj at data)
>                  kml_layer.SpatialLines(obj = sp, TimeSpan.begin = TimeSpan.begin,
>                    TimeSpan.end = TimeSpan.end, ...)
>              }
>
> ###########################################################
> #  New for Spatial Pixels                                 #
> ##########################################################
>          else {
>              if (class(obj at sp) == "SpatialPixels" | class(obj at sp) ==
>                  "SpatialPixelsDataFrame") {
>                  sp <- SpatialPixelsDataFrame(obj at sp, obj at 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 at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>