In Envirometrix/plotKML: Visualization of Spatial and Spatio-Temporal Objects in Google Earth

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)

plotKML

plotKML: Visualization of Spatial and Spatio-Temporal Objects in Google Earth

How to cite:

• Hengl, T., Roudier, P., Beaudette, D., & Pebesma, E. (2015). plotKML: Scientific visualization of spatio-temporal data. Journal of Statistical Software, 63(5), 1-25. https://www.jstatsoft.org/article/view/v063i05

Installing

Install development versions from github:

library(devtools)
install_github("envirometrix/plotKML")
# install_github("envirometrix/plotKML", ref = "stars")

Integration with sf

We will now present new sf methods and compare them with current sp approach. Start with POINTS

suppressPackageStartupMessages({
  library(plotKML)
  library(sp)
  library(rgdal)
  library(sf)
})
# Load data
data(eberg)
coordinates(eberg) <- ~ X + Y
proj4string(eberg) <- CRS("+init=epsg:31467")
## subset to 20 percent:
eberg <- eberg[runif(nrow(eberg)) < .1, ]
# sp methods
plotKML(eberg["CLYMHT_A"], open.kml = FALSE)
plotKML(eberg["CLYMHT_A"], colour_scale = rep("#FFFF00", 2), points_names = "", open.kml = FALSE)

# convert eberg to sf format
eberg_sf <- st_as_sf(eberg)

# and apply sf methods
plotKML(eberg_sf["CLYMHT_A"], open.kml = FALSE)
plotKML(eberg_sf["CLYMHT_A"], points_names = "", colour_scale = SAGA_pal[[2]], open.kml = FALSE) #Change palette
plotKML(eberg_sf["CLYMHT_A"], colour_scale = rep("#FFFF00", 2), points_names = "", open.kml = FALSE) # Degenerate palette

Recent version of sf could show a warning message relative to the old PROJ4string specified by get("ref_CRS", envir = plotKML.opts).

We can now compare the two implementations:

all.equal(readLines("eberg__CLYMHT_A__.kml"), readLines("eberg_sf__CLYMHT_A__.kml"))

The two string mismatches are simply caused by the different names and classes:

id_mismatches <- readLines("eberg__CLYMHT_A__.kml") != readLines("eberg_sf__CLYMHT_A__.kml")
readLines("eberg__CLYMHT_A__.kml")[id_mismatches]
readLines("eberg_sf__CLYMHT_A__.kml")[id_mismatches]

We can also use kml_layer functions:

data(eberg_grid)
gridded(eberg_grid) <- ~ x + y
proj4string(eberg_grid) <- CRS("+init=epsg:31467")
eberg_grid_sf <- st_as_sf(eberg_grid)

# sfc objects
kml_open("eberg_grids.kml")
kml_layer(st_geometry(eberg_grid_sf))
kml_close("eberg_grids.kml")

# sf objects
kml_open("eberg_grid_sf.kml")
kml_layer(eberg_grid_sf, colour = DEMSRT6, colour_scale = R_pal[["terrain_colors"]])
kml_layer(eberg_grid_sf, colour = TWISRT6, colour_scale = SAGA_pal[[1]])
kml_close("eberg_grid_sf.kml")

This is also a more complex example where we add a legend to the plot:

shape = "http://maps.google.com/mapfiles/kml/pal2/icon18.png" 
kml_file = "kml_file_point.kml"
legend_file = "kml_legend.png"
kml_open(kml_file)
kml_layer(
  eberg_sf["CLYMHT_A"], 
  colour = CLYMHT_A, 
  points_names = eberg_sf[["CLYMHT_A"]], 
  size = 1, 
  shape = shape, 
  alpha = 0.75, 
  colour_scale = SAGA_pal[[2]]
)
kml_legend.bar(eberg_sf$CLYMHT_A, legend.file = legend_file, legend.pal = SAGA_pal[[1]])
kml_screen(image.file = legend_file)
kml_close(kml_file)
kml_View(kml_file)

Then we present sf methods for MULTIPOINT objects. The idea is exactly the same, but MULTIPOINT object are converted to POINT.

eberg_sf_MULTIPOINT <- eberg_sf %>% 
  dplyr::mutate(random_ID = sample(1:4, size = dplyr::n(), replace = TRUE)) %>% 
  dplyr::group_by(random_ID) %>% 
  dplyr::summarise()
plotKML(eberg_sf_MULTIPOINT["random_ID"], open.kml = FALSE)

Then we present LINESTRING methods, starting from the sp example:

# sp
data(eberg_contours)
plotKML(eberg_contours, open.kml = FALSE)
plotKML(eberg_contours, colour = Z, altitude = Z, open.kml = FALSE)
# sf
eberg_contours_sf <- st_as_sf(eberg_contours)
plotKML(eberg_contours_sf, open.kml = FALSE)
plotKML(eberg_contours_sf, colour = Z, altitude = Z, open.kml = FALSE)

Again, the kml files are identical but for super small differences (if they are present) due to rounding errors:

all.equal(readLines("eberg_contours.kml"), readLines("eberg_contours_sf.kml"))

id_mismatches <- which(readLines("eberg_contours.kml") != readLines("eberg_contours_sf.kml"))
readLines("eberg_contours.kml")[id_mismatches[1:5]]
readLines("eberg_contours_sf.kml")[id_mismatches[1:5]]

The same ideas can be applied to MULTILINESTRING objects. The only problem is that some of the features in eberg_contous_sf are not valid according to the simple feature definition of LINESTRING. For example

eberg_contours@lines[[4]]

since it would represent a LINESTRING with only 1 point:

st_is_valid(st_geometry(eberg_contours_sf)[[4]], NA_on_exception = FALSE, reason = TRUE)

So we need to exclude these elements and create a MULTILINESTRING object,

ID_valid <- vapply(
  st_geometry(eberg_contours_sf), 
  function(x) isTRUE(st_is_valid(x)), 
  logical(1)
)
eberg_contous_sf_multi <- eberg_contours_sf %>% 
  dplyr::filter(ID_valid) %>% 
  dplyr::group_by(Z) %>% 
  dplyr::summarise()
plotKML(eberg_contous_sf_multi, colour = Z, altitude = Z, open.kml = FALSE)

We can also use kml_layer functions:

# sfc LINESTRING object
kml_open("eberg_contours_sfc.kml")
kml_layer(st_geometry(eberg_contours_sf))
kml_close("eberg_contours_sfc.kml")

Then we can work with POLYGON geometry, starting from sp example:

# sp 
set.seed(1) # I set the seed to compare the kml files
data(eberg_zones)
plotKML(eberg_zones["ZONES"], open.kml = FALSE)
## add altitude:
zmin = 230
plotKML(eberg_zones["ZONES"], altitude = zmin + runif(length(eberg_zones)) * 500, open.kml = FALSE)

# sf objects with sfc_POLYGON geometry
eberg_zones_sf <- st_as_sf(eberg_zones)
set.seed(1)
plotKML(eberg_zones_sf["ZONES"], open.kml = FALSE)
plotKML(eberg_zones_sf["ZONES"], altitude = zmin + runif(length(eberg_zones)) * 500, open.kml = FALSE)

and compare the results:

all.equal(readLines("eberg_zones__ZONES__.kml"), readLines("eberg_zones_sf__ZONES__.kml"))

The differences here (if they are present) are caused by extremely small rounding problems:

id_mismatches <- readLines("eberg_zones__ZONES__.kml") != readLines("eberg_zones_sf__ZONES__.kml")
readLines("eberg_zones__ZONES__.kml")[id_mismatches][1:5]
readLines("eberg_zones_sf__ZONES__.kml")[id_mismatches][42]

Again, the same ideas can be applied to MULTIPOLYGON objects:

eberg_zones_sf_MULTI <- eberg_zones_sf %>% 
  dplyr::group_by(ZONES) %>% 
  dplyr::summarise() %>% 
  st_cast("MULTIPOLYGON") 
plotKML(eberg_zones_sf_MULTI["ZONES"], open.kml = FALSE)

We will now present plotKML methods applied to stars objects. We extended the plotKML() function to stars objects representing Raster Data Cubes creating a wrapper to RasterLayer method:

library(stars)
# Start with a raster data cube
(g <- read_stars(system.file("external/test.grd", package = "raster")))
plotKML(g, open.kml = FALSE)

Default methods do not work if the third dimension does not represent a Date or POSIXct object:

# Test with another example: 
(L7 <- read_stars(system.file("tif/L7_ETMs.tif", package = "stars")))
plotKML(L7, open.kml = FALSE)

It fails. The problems is that there is no method definition for plotKML() function for objects of class RasterBrick. The code behind as(x, "Raster") is defined here and it says that if length(dim(x)) > 2 (i.e. there is an extra dimension other than x and y), then it creates a RasterBrick object. We can still plot the raster associated to one of the layers in the RasterBrick with a little bit of extra work:

plotKML(raster::raster(as(L7, "Raster"), layer = 1), open.kml = FALSE)

We are working on defining an extension to stars object with a temporal dimension as a wrapper around RasterBrickTimeSeries methods.

Let's focus now on Vector Data cubes:

# Vector data cube
data(air, package = "spacetime")
d = st_dimensions(station = st_as_sfc(stations), time = dates)
(aq = st_as_stars(list(PM10 = air), dimensions = d))
(agg = aggregate(aq, "months", mean, na.rm = TRUE))
plotKML(agg, open.kml = FALSE) 

Unfortunately there is a known bug in plotKML/spacetime and the following example fails:

# Spatial aggregation: 
(a = aggregate(agg, st_as_sf(DE_NUTS1), mean, na.rm = TRUE))
plotKML(a) # error

# The error is in as(obj, "STIDF"), i.e. the conversion between STFDF and STIDF

Moreover, all the examples fail if the third dimension is not a temporal object (since the methods we defined as wrappers around STFDF class):

nc = read_sf(system.file("gpkg/nc.gpkg", package="sf"))
nc.df = st_set_geometry(nc, NULL)
mat = as.matrix(nc.df[c("BIR74", "SID74", "NWBIR74", "BIR79", "SID79", "NWBIR79")])
dim(mat) = c(county = 100, var = 3, year = 2) # make it a 3-dimensional array
dimnames(mat) = list(county = nc$NAME, var = c("BIR", "SID", "NWBIR"), year = c(1974, 1979))
(nc.st = st_as_stars(pop = mat))
(nc.geom <- st_set_dimensions(nc.st, 1, st_geometry(nc)))
plotKML(nc.geom)

We can somehow fix this problem converting the stars object into sf format:

# We can use the following approach
plotKML(st_as_sf(nc.geom), open.kml = FALSE)

In some cases, it is still possible to redefine a temporal dimension (but unfortunately the example fails for the same bug as the previous case):

(nc.geom <- st_set_dimensions(nc.geom, 3, as.Date(c("1974-01-01","1979-01-01"))))
(nc.geom <- split(aperm(nc.geom, c(1,3,2))))
plotKML(nc.geom)

Envirometrix/plotKML documentation built on June 13, 2022, 11:21 p.m.