Écrire plusieurs couches dans GeoPackage en utilisant writeOGR () dans R?

11

Je suis en train d'écrire plusieurs couches à la même GeoPackage dans R, mais je reçois une erreur, Creation of output file failed. J'ai essayé de rechercher de la documentation sur la lecture et l'écriture dans des fichiers .gpkg avec RGDAL, par exemple pour déterminer si writeOGR()réellement prend en charge plusieurs couches, avec peu de succès. Est-ce même possible, si oui, comment le faire? Exemple de fonctionnement minimal:

library(sp)
library(maptools)
library(rgdal)

data(wrld_simpl)

norway <- wrld_simpl[wrld_simpl$NAME == "Norway", ]
sweden <- wrld_simpl[wrld_simpl$NAME == "Sweden", ]

file <- tempfile("scandinavia", fileext = c(".gpkg"))

writeOGR(norway, dsn = file, layer = "norway", driver = "GPKG")
writeOGR(sweden, dsn = file, layer = "sweden", driver = "GPKG")

ogrListLayers(file)

Il y a apparemment une ogr2ogr commande shell qui fait l'affaire (hat tip mdsumner ), que je peux encapsuler dans une fonction R. Cependant, ce serait bien si writeOGR () et / ou st_write () dans le sfpaquet avaient cela intégré. Je pense que cela dépend des GDAL layer_options, mais il ne semble pas y avoir d'option de type ajout pour GPKG dans GDAL .


Je pourrais écrire une fonction wrapper simple st_write()mais un support natif sfou rgdalserait mieux.

eivindhammers
la source
Pas possible afaik. Essayez avec sf, que je serai heureux d'explorer également - il est plus facile de réparer que rgdal pour une chose
mdsumner
1
@mdsumner st_write () dans sf donne le même résultat. Je pense que l'absence d'une option d'ajout dans les options de création de couches de GDAL est la source du problème pour writeOGR () et st_write ().
eivindhammers

Réponses:

10

Vous pouvez le faire en utilisant le appenddrapeau sur sf::st_write():

library(sf)

nc     <- st_read(system.file("shape/nc.shp", package="sf"))
storms <- st_read(system.file("shape/storms_xyz.shp", package="sf"))

st_write(nc,     "nc.gpkg", "nc")
st_write(storms, "nc.gpkg", "storms", append = TRUE)

st_layers("nc.gpkg")
## Driver: GPKG 
## Available layers:
##   layer_name  geometry_type features fields
## 1         nc  Multi Polygon      100     14
## 2     storms 3D Line String       71      0
jsta
la source