layer options for PostgreSQL driver in st_write are ignored

[aazaff]

Here is an example below where I read a geojson into R using sf::st_read and then write it to a postgres/postgis database. The PostgreSQL connection is made through RPostgreSQL and DBI pacakges.

Driver = dbDriver("PostgreSQL") # Establish database driver
Ecos = dbConnect(Driver, dbname = "ecos", host = "localhost", port = 1234) # I anonymised the db credentials, but you get the idea.
Ecoregions = sf::st_read("https://services.arcgis.com/F7DSX1DSNSiWmOqh/arcgis/rest/services/Marine_Ecoregions_Of_the_World_(MEOW)/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson")
sf::st_write(dsn=Ecos,obj=Ecoregions,c("plates","meow"),layer_options=c("GEOMETRY_NAME=geom","LAUNDER=true"))

I try both the geometry name and launder layer options (which are valid options of the driver: https://gdal.org/drivers/vector/pg.html), but neither selection sticks. The command go through fine in R with no warnings or error and the table does post to my postgres database, but the fields are not laundered and the geometry field name is not changed to geom.

I suppose it is possible that I am passing the layer options to to st_write incorrectly, but then I would expect to get an error or warning of some kind?

[aazaff]

Sorry, I should add that I am R 4.0.4, psql 11.11, and postgis 2.5

[edzer]

@etiennebr ?

[etiennebr]

Hi @aazaff, this is indeed tricky because sf supports two way to connect to a database: gdal and sf "native". The first one uses the gdal driver options, while the other uses a R connection and DBI options. Here's how I would suggest to use it. Personally I prefer to use an R connection object, because it make all the DBI functions available, but the gdal drivers works too and does not require to install any driver (i.e. RPostgres).

library(DBI)
library(sf) 
#> Warning: replacing previous import 'vctrs::data_frame' by 'tibble::data_frame'
#> when loading 'dplyr'
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
library(tidyverse)

# I'm using a local docker container 
# docker run -p 25432:5432 -e POSTGRES_USER=$USER -e POSTGRES_DBNAME=$USER -e POSTGRES_PASS=$USER -t kartoza/postgis:latest
con <- dbConnect(
    RPostgres::Postgres(),
    host = "localhost",
    port = 25432,
    dbname = "gis",
    user = "docker",
    password = "docker"
)

# 1. GDAL
ecoregions <- sf::st_read("https://services.arcgis.com/F7DSX1DSNSiWmOqh/arcgis/rest/services/Marine_Ecoregions_Of_the_World_(MEOW)/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson")
#> Reading layer `OGRGeoJSON' from data source `https://services.arcgis.com/F7DSX1DSNSiWmOqh/arcgis/rest/services/Marine_Ecoregions_Of_the_World_(MEOW)/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson' using driver `GeoJSON'
#> Simple feature collection with 232 features and 12 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -180 ymin: -89 xmax: 180 ymax: 86.9194
#> geographic CRS: WGS 84

st_layers("PG:host=localhost port=25432 dbname=gis user=docker password=docker")
#> Error in max(sapply(x$geomtype, length)): invalid 'type' (list) of argument
st_write(ecoregions, 
         "PG:host=localhost port=25432 dbname=gis user=docker password=docker", 
         layer = "ecoregions_gdal", 
         layer_options=c("GEOMETRY_NAME=geom","LAUNDER=true")
)
#> Writing layer `ecoregions_gdal' to data source `PG:host=localhost port=25432 dbname=gis user=docker password=docker' using driver `PostgreSQL'
#> options:        GEOMETRY_NAME=geom LAUNDER=true 
#> Writing 232 features with 12 fields and geometry type Multi Polygon.
st_read(con, "ecoregions_gdal") %>% 
    glimpse()
#> Rows: 232
#> Columns: 14
#> $ ogc_fid       <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
#> $ fid           <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
#> $ eco_code      <int> 20192, 20053, 20072, 20194, 20228, 20109, 20110, 20223,…
#> $ ecoregion     <chr> "Agulhas Bank", "Aleutian Islands", "Amazonia", "Amster…
#> $ prov_code     <int> 51, 10, 13, 52, 61, 24, 24, 60, 19, 45, 32, 32, 62, 3, …
#> $ province      <chr> "Agulhas", "Cold Temperate Northeast Pacific", "North B…
#> $ rlm_code      <int> 10, 3, 4, 10, 12, 5, 5, 12, 5, 9, 6, 6, 12, 2, 1, 4, 2,…
#> $ realm         <chr> "Temperate Southern Africa", "Temperate Northern Pacifi…
#> $ alt_code      <int> 189, 50, 64, 191, 209, 105, 106, 215, 92, 172, 132, 133…
#> $ eco_code_x    <int> 192, 53, 72, 194, 228, 109, 110, 223, 90, 178, 140, 139…
#> $ lat_zone      <chr> "Temperate", "Temperate", "Tropical", "Temperate", "Pol…
#> $ shape__area   <dbl> 8.898574e+11, 3.242661e+12, 1.545699e+12, 6.442118e+11,…
#> $ shape__length <dbl> 3950757, 10409934, 5193985, 2924452, 25670711, 3677656,…
#> $ geom          <MULTIPOLYGON [°]> MULTIPOLYGON (((28.35993 -3..., MULTIPOLYG…

# 2. `sf` "native
ecoregions <- ecoregions %>% 
    rename_all(tolower) %>% 
    rename_all(~gsub("(_)+", "\\1", .))

sf::st_write(dsn = con, obj = ecoregions)
#> Note: method with signature 'DBIObject#sf' chosen for function 'dbDataType',
#>  target signature 'PqConnection#sf'.
#>  "PqConnection#ANY" would also be valid
dbListTables(con)
#> [1] "geography_columns" "geometry_columns"  "spatial_ref_sys"  
#> [4] "ecoregions"        "topology"          "layer"            
#> [7] "raster_columns"    "raster_overviews"  "ecoregions_gdal"
st_read(con, "ecoregions") %>% 
    glimpse()
#> Rows: 232
#> Columns: 13
#> $ fid          <int> 1, 2, 10, 3, 4, 5, 6, 7, 8, 9, 11, 12, 13, 14, 106, 15, …
#> $ eco_code     <int> 20192, 20053, 20178, 20072, 20194, 20228, 20109, 20110, …
#> $ ecoregion    <chr> "Agulhas Bank", "Aleutian Islands", "Araucanian", "Amazo…
#> $ prov_code    <int> 51, 10, 45, 13, 52, 61, 24, 24, 60, 19, 32, 32, 62, 3, 2…
#> $ province     <chr> "Agulhas", "Cold Temperate Northeast Pacific", "Warm Tem…
#> $ rlm_code     <int> 10, 3, 9, 4, 10, 12, 5, 5, 12, 5, 6, 6, 12, 2, 6, 1, 4, …
#> $ realm        <chr> "Temperate Southern Africa", "Temperate Northern Pacific…
#> $ alt_code     <int> 189, 50, 172, 64, 191, 209, 105, 106, 215, 92, 132, 133,…
#> $ eco_code_x   <int> 192, 53, 178, 72, 194, 228, 109, 110, 223, 90, 140, 139,…
#> $ lat_zone     <chr> "Temperate", "Temperate", "Temperate", "Tropical", "Temp…
#> $ shape_area   <dbl> 8.898574e+11, 3.242661e+12, 8.885664e+11, 1.545699e+12, …
#> $ shape_length <dbl> 3950757, 10409934, 3819454, 5193985, 2924452, 25670711, …
#> $ geometry     <MULTIPOLYGON [°]> MULTIPOLYGON (((28.35993 -3..., MULTIPOLYGO…

^{Created on 2021-06-10 by the reprex package (v0.3.0)}

[aazaff]

Thanks! I will try and reproduce first thing tomorrow and then close the issue if I have no further questions.

[aazaff]

I tried and verified option 1 (gdal) and it worked perfectly.

As an additional note, if you want to pass a schema then you new need to do layer="schema.table" as opposed to the layer=c("schema","table")used in my original formulation.