inconsistent behavior of vect() when layers (in multilayer source) contain no records/geometry

[brownag]

I think I have noticed a small regression in behavior in terra::vect() when reading from a geopackage containing multiple vector layers. In the case where I encountered this I have geopackages that have a standard set of tables, which can, depending on the spatial extent, be empty (just columns, no rows)

See the comparisons below. I compare terra 1.7.73 and 1.7.71 (CRAN) with a prior version from the archive that seems to work differently (1.7.23). Note that the call terra::vect(f, "lux2") in recent versions silently returns the lux1 layer, rather than erroring with Error: [vect] cannot read this geometry type: Unknown (any). Passing proxy=TRUE gives the expected error. If "lux2" layer has more than 0 rows (i.e. has geometry) then everything works as expected.

terra 1.7.73 (same as CRAN version 1.7.71)

library(terra)
#> terra 1.7.73
x <- sf::st_as_sf(vect(system.file("ex", "lux.shp", package = "terra")))
f <- "lux.gpkg"
unlink(f)

sf::st_write(x, f, "lux1")
#> Writing layer `lux1' to data source `lux.gpkg' using driver `GPKG'
#> Writing 12 features with 6 fields and geometry type Polygon.

sf::st_write(x[0,], f, "lux2", append = TRUE)
#> Updating layer `lux2' to data source `lux.gpkg' using driver `GPKG'
#> Writing 0 features with 6 fields and geometry type Unknown (any).

 # as expected
terra::vect(f, "lux1")
#>  class       : SpatVector 
#>  geometry    : polygons 
#>  dimensions  : 12, 6  (geometries, attributes)
#>  extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
#>  source      : lux.gpkg (lux1)
#>  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#>  names       :  ID_1   NAME_1  ID_2   NAME_2  AREA   POP
#>  type        : <num>    <chr> <num>    <chr> <num> <int>
#>  values      :     1 Diekirch     1 Clervaux   312 18081
#>                    1 Diekirch     2 Diekirch   218 32543
#>                    1 Diekirch     3  Redange   259 18664

# unexpected result: silently returns lux1 layer
terra::vect(f, "lux2") 
#>  class       : SpatVector 
#>  geometry    : polygons 
#>  dimensions  : 12, 6  (geometries, attributes)
#>  extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
#>  source      : lux.gpkg (lux1)
#>  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#>  names       :  ID_1   NAME_1  ID_2   NAME_2  AREA   POP
#>  type        : <num>    <chr> <num>    <chr> <num> <int>
#>  values      :     1 Diekirch     1 Clervaux   312 18081
#>                    1 Diekirch     2 Diekirch   218 32543
#>                    1 Diekirch     3  Redange   259 18664

 # proxy=TRUE on lux2 gives expected error due to unknown geometry
terra::vect(f, "lux2", proxy = TRUE)
#> Error: [vect] cannot read this geometry type: Unknown (any)

# lux2 table does exist (it is 0 rows)
con <- DBI::dbConnect(RSQLite::SQLite(), f)
RSQLite::dbGetQuery(con, "SELECT * FROM lux2")
#> [1] fid    geom   ID_1   NAME_1 ID_2   NAME_2 AREA   POP   
#> <0 rows> (or 0-length row.names)
DBI::dbDisconnect(con)

terra 1.7.23

library(terra)
#> terra 1.7.23
x <- sf::st_as_sf(vect(system.file("ex", "lux.shp", package = "terra")))
f <- "lux.gpkg"
unlink(f)

sf::st_write(x, f, "lux1")
#> Writing layer `lux1' to data source `lux.gpkg' using driver `GPKG'
#> Writing 12 features with 6 fields and geometry type Polygon.

sf::st_write(x[0,], f, "lux2", append = TRUE)
#> Updating layer `lux2' to data source `lux.gpkg' using driver `GPKG'
#> Writing 0 features with 6 fields and geometry type Unknown (any).

# as expected
terra::vect(f, "lux1") 
#>  class       : SpatVector 
#>  geometry    : polygons 
#>  dimensions  : 12, 6  (geometries, attributes)
#>  extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
#>  source      : lux.gpkg (lux1)
#>  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#>  names       :  ID_1   NAME_1  ID_2   NAME_2  AREA   POP
#>  type        : <num>    <chr> <num>    <chr> <num> <int>
#>  values      :     1 Diekirch     1 Clervaux   312 18081
#>                    1 Diekirch     2 Diekirch   218 32543
#>                    1 Diekirch     3  Redange   259 18664

# when lux2 has no rows, terra cannot read the geometry type
terra::vect(f, "lux2") 
#> Error: [vect] cannot read this geometry type: Unknown (any)

# same error due to unknown geometry
terra::vect(f, "lux2", proxy = TRUE) 
#> Error: [vect] cannot read this geometry type: Unknown (any)

I think in an ideal solution terra would be able to read (and write) an unknown/empty table as a 0-row SpatVector.
But it would be OK to just get the original "cannot read this geometry type: Unknown (any)" error back, instead of silently returning "lux1" (which is a harder to detect and work around).

For what its worth, sf does appear to be able to write and read these types of empty/Unknown/any geometry table. Not sure if there is something specific about what sf is doing that terra doesn't support or is not accounted for--I could investigate that a bit closer.

sf::st_read(f, "lux2")
#> Reading layer `lux2' from data source 
#>   `/tmp/RtmpcQlvxK/reprex-1f1581c12684c-awake-topi/lux.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 0 features and 6 fields
#> Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
#> Geodetic CRS:  WGS 84

[rhijmans]

Hi, thanks, I think this commit fixes it.

[brownag]

Thanks! This seems closer to what I am hoping for. Now I get:

library(terra)
#> terra 1.7.73
x <- sf::st_as_sf(vect(system.file("ex", "lux.shp", package = "terra")))
f <- "lux.gpkg"
unlink(f)

sf::st_write(x, f, "lux1")
#> Writing layer `lux1' to data source `lux.gpkg' using driver `GPKG'
#> Writing 12 features with 6 fields and geometry type Polygon.

sf::st_write(sf::st_as_sf(x)[0,], f, "lux2", append = TRUE)
#> Updating layer `lux2' to data source `lux.gpkg' using driver `GPKG'
#> Writing 0 features with 6 fields and geometry type Unknown (any).

# as expected
terra::vect(f, "lux1") 
#>  class       : SpatVector 
#>  geometry    : polygons 
#>  dimensions  : 12, 6  (geometries, attributes)
#>  extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
#>  source      : lux.gpkg (lux1)
#>  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#>  names       :  ID_1   NAME_1  ID_2   NAME_2  AREA   POP
#>  type        : <num>    <chr> <num>    <chr> <num> <int>
#>  values      :     1 Diekirch     1 Clervaux   312 18081
#>                    1 Diekirch     2 Diekirch   218 32543
#>                    1 Diekirch     3  Redange   259 18664

# returns 0 row SpatVector, but should have 0 geometries, 6 attributes
terra::vect(f, "lux2")
#>  class       : SpatVector 
#>  geometry    : none 
#>  dimensions  : 0, 0  (geometries, attributes)
#>  extent      : 0, 0, 0, 0  (xmin, xmax, ymin, ymax)
#>  source      : lux.gpkg (lux2)
#>  coord. ref. : lon/lat WGS 84 (EPSG:4326)

# error due to unknown geometry
terra::vect(f, "lux2", proxy = TRUE) 
#> Error: [vect] cannot read this geometry type: Unknown (any)

There are two things I see here

• now the proxy=FALSE result returns a 0 geometry, 0 attribute SpatVector. Ideally I would want it to return the equivalent of x[0,] which has 6 attributes.
• proxy=TRUE still has an error reading the Unknown geometry; is there a way to make this work like proxy=FALSE? or have proxy=FALSE produce the same error?

---

I might be getting greedy here, so feel free to ignore this... but it would be awesome if I didn't need to use sf at all to create the geopackages in question (for my workflow where I came across this, and for this simple example)... There are a few things that stand in the way of that as I see it:

library(terra)
#> terra 1.7.73
x <- vect(system.file("ex", "lux.shp", package = "terra"))
f <- "lux.gpkg"
unlink(f)

terra::writeVector(x, f, layer = "lux1")
#> Warning in x@ptr$write(filename, layer, filetype, insert[1], overwrite[1], :
#> GDAL Message 6: dataset lux.gpkg does not support layer creation option
#> ENCODING

# my library info RE: warning above
terra::gdal(lib = TRUE)
#>     gdal     proj     geos 
#>  "3.8.1"  "9.3.1" "3.12.0"

# 0 geometries, 6 attributes
x[0,]
#>  class       : SpatVector 
#>  geometry    : none 
#>  dimensions  : 0, 6  (geometries, attributes)
#>  extent      : 0, 0, 0, 0  (xmin, xmax, ymin, ymax)
#>  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#>  names       :  ID_1 NAME_1  ID_2 NAME_2  AREA   POP
#>  type        : <num>  <chr> <num>  <chr> <num> <int>

# gives warning and writes nothing 
terra::writeVector(x[0,], f, layer = "lux2")
#> Warning: [writeVector] nothing to write

# as expected
terra::vect(f, "lux1") 
#>  class       : SpatVector 
#>  geometry    : polygons 
#>  dimensions  : 12, 6  (geometries, attributes)
#>  extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
#>  source      : lux.gpkg (lux1)
#>  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#>  names       :  ID_1   NAME_1  ID_2   NAME_2  AREA   POP
#>  type        : <num>    <chr> <num>    <chr> <num> <int>
#>  values      :     1 Diekirch     1 Clervaux   312 18081
#>                    1 Diekirch     2 Diekirch   218 32543
#>                    1 Diekirch     3  Redange   259 18664

# error because appending lux2 above fails
terra::vect(f, "lux2")
#> Error: [vect] lux2 is not a valid layer name
#> Choose one of: lux1

[rhijmans]

Yes, sorry, that addressed only a bit of it. Getting closer:

library(terra)
x <- vect(system.file("ex", "lux.shp", package = "terra"))
f <- "lux.gpkg"
unlink(f)
terra::writeVector(x, f, layer="lux1")
terra::writeVector(x[0,], f, layer = "lux2", insert=TRUE)
terra::vect(f, "lux2") 
# class       : SpatVector 
# geometry    : none 
# dimensions  : 0, 6  (geometries, attributes)
# extent      : 0, 0, 0, 0  (xmin, xmax, ymin, ymax)
# source      : lux.gpkg (lux2)
# coord. ref. : lon/lat WGS 84 (EPSG:4326) 
# names       :  ID_1 NAME_1  ID_2 NAME_2  AREA   POP
# type        : <num>  <chr> <num>  <chr> <num> <int>

Still need to fix the proxy opening:

terra::vect(f, "lux2", proxy=T) 
#Error: [vect] cannot read this geometry type: Unknown (any)

[rhijmans]

Perhaps the issue has been addressed. I now get:

terra::vect(f, "lux2", proxy=T) 
# class       : SpatVectorProxy
# geometry    : none 
# dimensions  : 0, 6  (geometries, attributes)
# extent      : 0, 0, 0, 0  (xmin, xmax, ymin, ymax)
# source      : lux.gpkg (lux2)
# layer       : lux2 
# coord. ref. : lon/lat WGS 84 (EPSG:4326) 
# names       :  ID_1 NAME_1  ID_2 NAME_2  AREA   POP
# type        : <num>  <chr> <num>  <chr> <num> <int>

[brownag]

Yes, I think so! Thanks so much Robert. Now working in my minimal example, and in my original problem workflow

library(terra)
#> terra 1.7.73
x <- vect(system.file("ex", "lux.shp", package = "terra"))
f <- "lux.gpkg"
unlink(f)

writeVector(x, f, layer = "lux1", options = "")

writeVector(x[0,], f, layer = "lux2", insert = TRUE, overwrite = TRUE, options = "")

# 12 geometries, 6 attributes
vect(f, "lux1") 
#>  class       : SpatVector 
#>  geometry    : polygons 
#>  dimensions  : 12, 6  (geometries, attributes)
#>  extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
#>  source      : lux.gpkg (lux1)
#>  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#>  names       :  ID_1   NAME_1  ID_2   NAME_2  AREA   POP
#>  type        : <num>    <chr> <num>    <chr> <num> <int>
#>  values      :     1 Diekirch     1 Clervaux   312 18081
#>                    1 Diekirch     2 Diekirch   218 32543
#>                    1 Diekirch     3  Redange   259 18664

# 0 geometries, 6 attributes
vect(f, "lux2")
#>  class       : SpatVector 
#>  geometry    : none 
#>  dimensions  : 0, 6  (geometries, attributes)
#>  extent      : 0, 0, 0, 0  (xmin, xmax, ymin, ymax)
#>  source      : lux.gpkg (lux2)
#>  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#>  names       :  ID_1 NAME_1  ID_2 NAME_2  AREA   POP
#>  type        : <num>  <chr> <num>  <chr> <num> <int>

# 0 geometries, 6 attributes
vect(f, "lux2", proxy = TRUE) 
#>  class       : SpatVectorProxy
#>  geometry    : none 
#>  dimensions  : 0, 6  (geometries, attributes)
#>  extent      : 0, 0, 0, 0  (xmin, xmax, ymin, ymax)
#>  source      : lux.gpkg (lux2)
#>  layer       : lux2 
#>  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#>  names       :  ID_1 NAME_1  ID_2 NAME_2  AREA   POP
#>  type        : <num>  <chr> <num>  <chr> <num> <int>