Envirometrix / plotKML

plotKML: Visualization of Spatial and Spatio-Temporal Objects in Google Earth
29 stars 7 forks source link

plotKML

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

Build
Status R build
status CRAN\_Status\_Badge Github\_Status\_Badge

How to cite:

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")
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
#> datum Deutsches_Hauptdreiecksnetz in CRS definition
## subset to 20 percent:
eberg <- eberg[runif(nrow(eberg)) < .1, ]
# sp methods
plotKML(eberg["CLYMHT_A"], open.kml = FALSE)
#> Plotting the first variable on the list
#> KML file opened for writing...
#> Warning in proj4string(obj): CRS object has comment, which is lost in output
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs ...
#> Writing to KML...
#> Closing  eberg__CLYMHT_A__.kml
#> Object written to: eberg__CLYMHT_A__.kml
plotKML(eberg["CLYMHT_A"], colour_scale = rep("#FFFF00", 2), points_names = "", open.kml = FALSE)
#> Plotting the first variable on the list
#> KML file opened for writing...
#> Warning in proj4string(obj): CRS object has comment, which is lost in output
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs ...
#> Writing to KML...
#> Closing  eberg__CLYMHT_A__.kml
#> Object written to: eberg__CLYMHT_A__.kml

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

# and apply sf methods
plotKML(eberg_sf["CLYMHT_A"], open.kml = FALSE)
#> Plotting the first variable on the list
#> KML file opened for writing...
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs
#> Writing to KML...
#> Closing  eberg_sf__CLYMHT_A__.kml
#> Object written to: eberg_sf__CLYMHT_A__.kml
plotKML(eberg_sf["CLYMHT_A"], points_names = "", colour_scale = SAGA_pal[[2]], open.kml = FALSE) #Change palette
#> Plotting the first variable on the list
#> KML file opened for writing...
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs
#> Writing to KML...
#> Closing  eberg_sf__CLYMHT_A__.kml
#> Object written to: eberg_sf__CLYMHT_A__.kml
plotKML(eberg_sf["CLYMHT_A"], colour_scale = rep("#FFFF00", 2), points_names = "", open.kml = FALSE) # Degenerate palette
#> Plotting the first variable on the list
#> KML file opened for writing...
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs
#> Writing to KML...
#> Closing  eberg_sf__CLYMHT_A__.kml
#> Object written to: eberg_sf__CLYMHT_A__.kml

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"))
#> Warning in readLines("eberg__CLYMHT_A__.kml"): incomplete final line found on
#> 'eberg__CLYMHT_A__.kml'
#> Warning in readLines("eberg_sf__CLYMHT_A__.kml"): incomplete final line found on
#> 'eberg_sf__CLYMHT_A__.kml'
#> [1] "2 string mismatches"

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]
#> [1] "    <name>eberg__CLYMHT_A__</name>"       
#> [2] "      <name>SpatialPointsDataFrame</name>"
readLines("eberg_sf__CLYMHT_A__.kml")[id_mismatches]
#> [1] "    <name>eberg_sf__CLYMHT_A__</name>"
#> [2] "      <name>sfdata.frame</name>"

We can also use kml_layer functions:

data(eberg_grid)
gridded(eberg_grid) <- ~ x + y
proj4string(eberg_grid) <- CRS("+init=epsg:31467")
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
#> datum Deutsches_Hauptdreiecksnetz in CRS definition
eberg_grid_sf <- st_as_sf(eberg_grid)

# sfc objects
kml_open("eberg_grids.kml")
#> KML file opened for writing...
kml_layer(st_geometry(eberg_grid_sf))
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs
#> Writing to KML...
kml_close("eberg_grids.kml")
#> Closing  eberg_grids.kml

# sf objects
kml_open("eberg_grid_sf.kml")
#> KML file opened for writing...
kml_layer(eberg_grid_sf, colour = DEMSRT6, colour_scale = R_pal[["terrain_colors"]])
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs
#> Writing to KML...
kml_layer(eberg_grid_sf, colour = TWISRT6, colour_scale = SAGA_pal[[1]])
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs
#> Writing to KML...
kml_close("eberg_grid_sf.kml")
#> Closing  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()
#> `summarise()` ungrouping output (override with `.groups` argument)
plotKML(eberg_sf_MULTIPOINT["random_ID"], open.kml = FALSE)
#> Casting the input MULTIPOINT objct into POINT object.
#> Warning in st_cast.sf(obj, "POINT"): repeating attributes for all sub-geometries
#> for which they may not be constant
#> Plotting the first variable on the list
#> KML file opened for writing...
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs
#> Writing to KML...
#> Closing  eberg_sf_MULTIPOINT__random_ID__.kml
#> Object written to: eberg_sf_MULTIPOINT__random_ID__.kml

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

# sp
data(eberg_contours)
plotKML(eberg_contours, open.kml = FALSE)
#> KML file opened for writing...
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs ...
#> Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
#> to PROJ string
#> Warning in spTransform(xSP, CRSobj, ...): +init dropped in PROJ string
#> Writing to KML...
#> Closing  eberg_contours.kml
#> Object written to: eberg_contours.kml
plotKML(eberg_contours, colour = Z, altitude = Z, open.kml = FALSE)
#> KML file opened for writing...
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs ...
#> Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
#> to PROJ string

#> Warning in spTransform(xSP, CRSobj, ...): +init dropped in PROJ string
#> Writing to KML...
#> Closing  eberg_contours.kml
#> Object written to: eberg_contours.kml
# sf
eberg_contours_sf <- st_as_sf(eberg_contours)
#> Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
#> deprecated. It might return a CRS with a non-EPSG compliant axis order.
plotKML(eberg_contours_sf, open.kml = FALSE)
#> KML file opened for writing...
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs
#> Writing to KML...
#> Closing  eberg_contours_sf.kml
#> Object written to: eberg_contours_sf.kml
plotKML(eberg_contours_sf, colour = Z, altitude = Z, open.kml = FALSE)
#> KML file opened for writing...
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs
#> Writing to KML...
#> Closing  eberg_contours_sf.kml
#> Object written to: eberg_contours_sf.kml

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"))
#> [1] "2 string mismatches"

id_mismatches <- which(readLines("eberg_contours.kml") != readLines("eberg_contours_sf.kml"))
readLines("eberg_contours.kml")[id_mismatches[1:5]]
#> [1] "    <name>eberg_contours</name>"         
#> [2] "      <name>SpatialLinesDataFrame</name>"
#> [3] NA                                        
#> [4] NA                                        
#> [5] NA
readLines("eberg_contours_sf.kml")[id_mismatches[1:5]]
#> [1] "    <name>eberg_contours_sf</name>" "      <name>sfdata.frame</name>"   
#> [3] NA                                   NA                                  
#> [5] NA

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]]
#> An object of class "Lines"
#> Slot "Lines":
#> [[1]]
#> An object of class "Line"
#> Slot "coords":
#>         [,1]    [,2]
#> [1,] 3579413 5714807
#> 
#> 
#> 
#> Slot "ID":
#> [1] "3"

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)
#> Error in CPL_geos_is_valid_reason(x): Evaluation error: IllegalArgumentException: point array must contain 0 or >1 elements.

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()
#> `summarise()` ungrouping output (override with `.groups` argument)
plotKML(eberg_contous_sf_multi, colour = Z, altitude = Z, open.kml = FALSE)
#> Casting the input MULTILINESTRING objct into LINESTRING object.
#> Warning in st_cast.sf(obj, "LINESTRING"): repeating attributes for all sub-
#> geometries for which they may not be constant
#> KML file opened for writing...
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs
#> Writing to KML...
#> Closing  eberg_contous_sf_multi.kml
#> Object written to: eberg_contous_sf_multi.kml

We can also use kml_layer functions:

# sfc LINESTRING object
kml_open("eberg_contours_sfc.kml")
#> KML file opened for writing...
kml_layer(st_geometry(eberg_contours_sf))
#> Casting the input MULTILINESTRING object into LINESTRING.
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs
#> Writing to KML...
kml_close("eberg_contours_sfc.kml")
#> Closing  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)
#> Plotting the first variable on the list
#> KML file opened for writing...
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs ...
#> Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
#> to PROJ string
#> Warning in spTransform(xSP, CRSobj, ...): +init dropped in PROJ string
#> Writing to KML...
#> Closing  eberg_zones__ZONES__.kml
#> Object written to: eberg_zones__ZONES__.kml
## add altitude:
zmin = 230
plotKML(eberg_zones["ZONES"], altitude = zmin + runif(length(eberg_zones)) * 500, open.kml = FALSE)
#> Plotting the first variable on the list
#> KML file opened for writing...
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs ...
#> Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
#> to PROJ string

#> Warning in spTransform(xSP, CRSobj, ...): +init dropped in PROJ string
#> Writing to KML...
#> Closing  eberg_zones__ZONES__.kml
#> Object written to: eberg_zones__ZONES__.kml

# 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)
#> Plotting the first variable on the list
#> KML file opened for writing...
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs
#> Writing to KML...
#> Closing  eberg_zones_sf__ZONES__.kml
#> Object written to: eberg_zones_sf__ZONES__.kml
plotKML(eberg_zones_sf["ZONES"], altitude = zmin + runif(length(eberg_zones)) * 500, open.kml = FALSE)
#> Plotting the first variable on the list
#> KML file opened for writing...
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs
#> Writing to KML...
#> Closing  eberg_zones_sf__ZONES__.kml
#> Object written to: eberg_zones_sf__ZONES__.kml

and compare the results:

all.equal(readLines("eberg_zones__ZONES__.kml"), readLines("eberg_zones_sf__ZONES__.kml"))
#> [1] "2 string mismatches"

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]
#> [1] "    <name>eberg_zones__ZONES__</name>"      
#> [2] "      <name>SpatialPolygonsDataFrame</name>"
#> [3] NA                                           
#> [4] NA                                           
#> [5] NA
readLines("eberg_zones_sf__ZONES__.kml")[id_mismatches][42]
#> [1] NA

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") 
#> `summarise()` ungrouping output (override with `.groups` argument)
plotKML(eberg_zones_sf_MULTI["ZONES"], open.kml = FALSE)
#> Casting the input MULTIPOLYGON objct into POLYGON object.
#> Warning in st_cast.sf(obj, "POLYGON"): repeating attributes for all sub-
#> geometries for which they may not be constant
#> Plotting the first variable on the list
#> KML file opened for writing...
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs
#> Writing to KML...
#> Closing  eberg_zones_sf_MULTI__ZONES__.kml
#> Object written to: eberg_zones_sf_MULTI__ZONES__.kml

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)
#> Loading required package: abind
# Start with a raster data cube
(g <- read_stars(system.file("external/test.grd", package = "raster")))
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#>    test.grd      
#>  Min.   : 138.7  
#>  1st Qu.: 294.0  
#>  Median : 371.9  
#>  Mean   : 425.6  
#>  3rd Qu.: 501.0  
#>  Max.   :1736.1  
#>  NA's   :6022    
#> dimension(s):
#>   from  to offset delta                       refsys point values x/y
#> x    1  80 178400    40 +proj=sterea +lat_0=52.15...    NA   NULL [x]
#> y    1 115 334000   -40 +proj=sterea +lat_0=52.15...    NA   NULL [y]
plotKML(g, open.kml = FALSE)
#> Plotting the first variable on the list
#> KML file opened for writing...
#> Warning in proj4string(obj): CRS object has comment, which is lost in output
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs ...
#> Writing to KML...
#> Closing  obj.kml
#> Object written to: obj.kml

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")))
#> stars object with 3 dimensions and 1 attribute
#> attribute(s):
#>   L7_ETMs.tif    
#>  Min.   :  1.00  
#>  1st Qu.: 54.00  
#>  Median : 69.00  
#>  Mean   : 68.91  
#>  3rd Qu.: 86.00  
#>  Max.   :255.00  
#> dimension(s):
#>      from  to  offset delta                       refsys point values x/y
#> x       1 349  288776  28.5 UTM Zone 25, Southern Hem... FALSE   NULL [x]
#> y       1 352 9120761 -28.5 UTM Zone 25, Southern Hem... FALSE   NULL [y]
#> band    1   6      NA    NA                           NA    NA   NULL
plotKML(L7, open.kml = FALSE)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
#>  but +towgs84= values preserved
#> Error in .local(obj, ...): Select only one Raster layer or provide a Time series

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)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
#>  but +towgs84= values preserved
#> Plotting the first variable on the list
#> KML file opened for writing...
#> Warning in proj4string(obj): CRS object has comment, which is lost in output
#> Reprojecting to +proj=longlat +datum=WGS84 +no_defs ...
#> Writing to KML...
#> Closing  raster__raster(as(L7,__Raster_),_layer_=_1).kml
#> Object written to: raster__raster(as(L7,__Raster_),_layer_=_1).kml

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))
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#>      PM10        
#>  Min.   :  0.00  
#>  1st Qu.:  9.92  
#>  Median : 14.79  
#>  Mean   : 17.70  
#>  3rd Qu.: 21.99  
#>  Max.   :274.33  
#>  NA's   :157659  
#> dimension(s):
#>         from   to     offset  delta                     refsys point
#> station    1   70         NA     NA +proj=longlat +datum=WGS84  TRUE
#> time       1 4383 1998-01-01 1 days                       Date FALSE
#>                                                          values
#> station POINT (9.585911 53.67057),...,POINT (9.446661 49.24068)
#> time                                                       NULL
(agg = aggregate(aq, "months", mean, na.rm = TRUE))
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#>      PM10       
#>  Min.   : 2.99  
#>  1st Qu.:13.01  
#>  Median :16.43  
#>  Mean   :17.68  
#>  3rd Qu.:20.83  
#>  Max.   :68.55  
#>  NA's   :4926   
#> dimension(s):
#>         from  to offset delta                     refsys point
#> time       1 144     NA    NA                       Date    NA
#> station    1  70     NA    NA +proj=longlat +datum=WGS84  TRUE
#>                                                          values
#> time                                  1998-01-01,...,2009-12-01
#> station POINT (9.585911 53.67057),...,POINT (9.446661 49.24068)
plotKML(agg, open.kml = FALSE) 
#> Plotting the first variable on the list
#> KML file opened for writing...
#> Warning in proj4string(obj): CRS object has comment, which is lost in output

#> Warning in proj4string(obj): CRS object has comment, which is lost in output
#> Writing to KML...
#> Closing  obj.kml
#> Object written to: obj.kml

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))
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#>      PM10        
#>  Min.   : 6.037  
#>  1st Qu.:14.072  
#>  Median :16.704  
#>  Mean   :17.897  
#>  3rd Qu.:20.378  
#>  Max.   :59.218  
#>  NA's   :832     
#> dimension(s):
#>          from  to offset delta                     refsys point
#> geometry    1  16     NA    NA +proj=longlat +datum=WGS84 FALSE
#> time        1 144     NA    NA                       Date    NA
#>                                                                     values
#> geometry MULTIPOLYGON (((9.65046 49....,...,MULTIPOLYGON (((10.77189 51...
#> time                                             1998-01-01,...,2009-12-01
plotKML(a) # error
#> Error in from@sp[from@index[, 1], ]: SpatialPolygons selection: can't find plot order if polygons are replicated

# 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))
#> stars object with 3 dimensions and 1 attribute
#> attribute(s):
#>       pop       
#>  Min.   :    0  
#>  1st Qu.:    8  
#>  Median :  538  
#>  Mean   : 1657  
#>  3rd Qu.: 1784  
#>  Max.   :30757  
#> dimension(s):
#>        from  to offset delta refsys point              values
#> county    1 100     NA    NA     NA    NA  Ashe,...,Brunswick
#> var       1   3     NA    NA     NA    NA BIR  , SID  , NWBIR
#> year      1   2     NA    NA     NA    NA          1974, 1979
(nc.geom <- st_set_dimensions(nc.st, 1, st_geometry(nc)))
#> stars object with 3 dimensions and 1 attribute
#> attribute(s):
#>       pop       
#>  Min.   :    0  
#>  1st Qu.:    8  
#>  Median :  538  
#>  Mean   : 1657  
#>  3rd Qu.: 1784  
#>  Max.   :30757  
#> dimension(s):
#>      from  to offset delta refsys point
#> sfc     1 100     NA    NA  NAD27 FALSE
#> var     1   3     NA    NA     NA    NA
#> year    1   2     NA    NA     NA    NA
#>                                                                 values
#> sfc  MULTIPOLYGON (((-81.47276 3...,...,MULTIPOLYGON (((-78.65572 3...
#> var                                                BIR  , SID  , NWBIR
#> year                                                        1974, 1979
plotKML(nc.geom)
#> Error in SpatialPolygonsDataFrame(x, y, match.ID = match.ID, ...): Object length mismatch:
#>      x has 100 Polygons objects, but y has 200 rows

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)
#> Casting the input MULTIPOLYGON objct into POLYGON object.
#> Warning in st_cast.sf(obj, "POLYGON"): repeating attributes for all sub-
#> geometries for which they may not be constant
#> Plotting the first variable on the list
#> KML file opened for writing...
#> Writing to KML...
#> Closing  st_as_sf(nc.geom).kml
#> Object written to: st_as_sf(nc.geom).kml

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"))))
#> stars object with 3 dimensions and 1 attribute
#> attribute(s):
#>       pop       
#>  Min.   :    0  
#>  1st Qu.:    8  
#>  Median :  538  
#>  Mean   : 1657  
#>  3rd Qu.: 1784  
#>  Max.   :30757  
#> dimension(s):
#>      from  to     offset     delta refsys point
#> sfc     1 100         NA        NA  NAD27 FALSE
#> var     1   3         NA        NA     NA    NA
#> year    1   2 1974-01-01 1826 days   Date    NA
#>                                                                 values
#> sfc  MULTIPOLYGON (((-81.47276 3...,...,MULTIPOLYGON (((-78.65572 3...
#> var                                                BIR  , SID  , NWBIR
#> year                                                              NULL
(nc.geom <- split(aperm(nc.geom, c(1,3,2))))
#> stars object with 2 dimensions and 3 attributes
#> attribute(s):
#>       BIR             SID             NWBIR      
#>  Min.   :  248   Min.   : 0.000   Min.   :    1  
#>  1st Qu.: 1177   1st Qu.: 2.000   1st Qu.:  206  
#>  Median : 2265   Median : 5.000   Median :  742  
#>  Mean   : 3762   Mean   : 7.515   Mean   : 1202  
#>  3rd Qu.: 4451   3rd Qu.: 9.000   3rd Qu.: 1316  
#>  Max.   :30757   Max.   :57.000   Max.   :11631  
#> dimension(s):
#>      from  to     offset     delta refsys point
#> sfc     1 100         NA        NA  NAD27 FALSE
#> year    1   2 1974-01-01 1826 days   Date    NA
#>                                                                 values
#> sfc  MULTIPOLYGON (((-81.47276 3...,...,MULTIPOLYGON (((-78.65572 3...
#> year                                                              NULL
plotKML(nc.geom)
#> Error in from@sp[from@index[, 1], ]: SpatialPolygons selection: can't find plot order if polygons are replicated