Writing empty table in a Geopackage of a certain geometry type using sf

Question

I'd like to create an empty table (i.e.: no rows) in Geopackage but with a specific geometry type (i.e.: polygons).

The following code creates a simple feature collection with 0 features and 0 fields, but nowhere is the geometry type specified.

I looked to see if anything could be provided to the layer_options argument for sf::st_write but I didn't come up with anything.

file_path <- "example/mypackage.gpkg"

empty_table <- sf::st_sf(geometry = sf::st_sfc())

sf::st_write(empty_table, file_path, layer = "empty_table")

Answer 1

This is somewhat tricky; kind of a zen-like: the geometry type of layer that is not...

While I kind of fail to see the benefit of doing this (I am a man of limited imagination) it should be in principle doable at low level - utilizing direct SQL calls via DBI backend. I have used it successfully to create views - wrapping a spatial table in a kind of a star schema.

You would need to:

• connect to your geopackage as a RSQLite database; doable from R or even python if so inclined
• create an empty data table of desired structure
• create a record for your data table in gpkg_contents table, which is obligatory in a geopackage
• create for your data table a record in gpkg_geometry_columns which is again obligatory for a geopackage, and sets stuff like geometry type and CRS.

In the last step you should be able to declare your empty table to be of geometry type polygon or what not...

Answer 2

I can't see a way to do it from sf because it seems to drop geometry type for empty geometries. But you can create a table with one row and then use SQL to delete all rows.

Example, here's a simple 1-row POINT spatial data frame:

> d1
Simple feature collection with 1 feature and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 0.1055935 ymin: 0.4832159 xmax: 0.1055935 ymax: 0.4832159
CRS:           NA
  Z                    geometry
1 a POINT (0.1055935 0.4832159)

Write it to a geopackage:

> st_write(d1, "empty.gpkg", "empty")
writing GPKG: substituting LOCAL_CS["Undefined Cartesian SRS"] for missing CRS
Writing layer `empty' to data source `empty.gpkg' using driver `GPKG'
Writing 1 features with 1 fields and geometry type Point.

Using the RSQLite package, make a connection:

> con = dbConnect(SQLite(), "empty.gpkg")

Delete every row:

> dbExecute(con, "delete from empty where 1 = 1")
[1] 5

And now when we read back it has zero rows:

> st_read("empty.gpkg")
Reading layer `empty' from data source 
  `/home/rowlings/Downloads/SO/empty/empty.gpkg' using driver `GPKG'
Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: Undefined Cartesian SRS

And note sf has no clue of the geometry, although it is in the geopackage metadata table as "POINT":

> dbReadTable(con, "gpkg_geometry_columns")
  table_name column_name geometry_type_name srs_id z m
1      empty        geom              POINT 100000 0 0