st_write, GPKG and ESRI

[fendabenda]

Bug description

When I export data as a gpkg using st_write, I am unable to work with the data as expected in ArcGIS Pro (I cannot add a new field to the attribute table or modify existing values, though I can move points), and in ArcMap I cannot open the data (I get a layer drawing error). However when I export the same data as a shapefile, I am able to work with the data as expected.

To Reproduce

`install.packages("osmdata") library(osmdata)

bb <- getbb('Montreal') #get bounding box for montreal

q <- opq(bbox = bb) #constructing an overpass API query (opq)

x <- q %>% add_osm_feature(key = 'amenity', value = "hospital", value_exact = FALSE)

hospitals <- osmdata_sf(x) #Now that the query is designed, we can send it to Overpass, get the data, and store it as an sf object with osmdata_sf()

library(sf) centroid_hospitals <- hospitals$osm_polygons %>% st_centroid()

st_write(centroid_hospitals, "hospitals.gpkg", "hospitals" )

st_write (centroid_hospitals, "hospitals.shp")`

Additional Context This issue is not specific to OSM datasets. I have reproduced this problem with other datasets. I have also reviewed Issue #392 but I am unable to reproduce the OP's successful export of the geopackage using something like st_write(nc, file.path(getwd(), 'nc.gpkg'), quiet = TRUE).

[rsbivand]

Please reproduce directly in GDAL. Few developers have access to any version of ArcGIS. For example report sf::sf_extSoftVersion() and sessionInfo(), and run sf::gdal_utils("vector_translate", ...) on any shapefile that ArcGIS can use to convert to GPKG outside R and sf using only GDAL.

[fendabenda]

Hi,

I have never used GDAL before. Are these scripts I can run directly in R without installing anything? I'm just a little confused about how to proceed.

On Tue, Nov 8, 2022 at 2:18 AM Roger Bivand @.***> wrote:

Please reproduce directly in GDAL. Few developers have access to any version of ArcGIS. For example report sf::sf_extSoftVersion() and sessionInfo(), and run sf::gdal_utils("vector_translate", ...) on any shapefile that ArcGIS can use to convert to GPKG outside R and sf using only GDAL.

— Reply to this email directly, view it on GitHub https://github.com/r-spatial/sf/issues/2029#issuecomment-1306741321, or unsubscribe https://github.com/notifications/unsubscribe-auth/AO4DNZPAXKZCRDZIPRDRS3DWHH5ERANCNFSM6AAAAAARZY352M . You are receiving this because you authored the thread.Message ID: @.***>

[edzer]

Maybe a problem of ArcGIS Pro?

[fendabenda]

The issue also occurs when I attempt import the data into ArcMap -- I get a layer drawing error upon import.

On Tue, Nov 8, 2022 at 2:27 PM Edzer Pebesma @.***> wrote:

Maybe a problem of ArcGIS Pro?

— Reply to this email directly, view it on GitHub https://github.com/r-spatial/sf/issues/2029#issuecomment-1307724674, or unsubscribe https://github.com/notifications/unsubscribe-auth/AO4DNZNDSXZRFD6ARCYW4FLWHKSQHANCNFSM6AAAAAARZY352M . You are receiving this because you authored the thread.Message ID: @.***>

[fendabenda]

Oh I think I see now what you want me to do. I will run the scripts and post what I did + results.

On Tue, Nov 8, 2022 at 2:18 AM Roger Bivand @.***> wrote:

Please reproduce directly in GDAL. Few developers have access to any version of ArcGIS. For example report sf::sf_extSoftVersion() and sessionInfo(), and run sf::gdal_utils("vector_translate", ...) on any shapefile that ArcGIS can use to convert to GPKG outside R and sf using only GDAL.

— Reply to this email directly, view it on GitHub https://github.com/r-spatial/sf/issues/2029#issuecomment-1306741321, or unsubscribe https://github.com/notifications/unsubscribe-auth/AO4DNZPAXKZCRDZIPRDRS3DWHH5ERANCNFSM6AAAAAARZY352M . You are receiving this because you authored the thread.Message ID: @.***>

[fendabenda]

Let me know if I did this correctly...

` > sf::sf_extSoftVersion() GEOS GDAL proj.4 GDAL_with_GEOS "3.9.1" "3.2.1" "7.2.1" "true" USE_PROJ_H PROJ "true" "7.2.1"

sessionInfo() R version 4.0.3 (2020-10-10) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 18363)

Matrix products: default

locale: [1] LC_COLLATE=English_United States.1252 [2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252 [4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages: [1] stats graphics grDevices utils datasets methods base

other attached packages: [1] sf_1.0-7

loaded via a namespace (and not attached): [1] rnaturalearth_0.1.0 Rcpp_1.0.8.3 rstudioapi_0.13
[4] magrittr_2.0.2 units_0.8-0 tidyselect_1.1.2
[7] lattice_0.20-41 R6_2.5.1 rlang_1.0.2
[10] fansi_1.0.2 dplyr_1.0.8 tools_4.0.3
[13] grid_4.0.3 xfun_0.30 tinytex_0.37
[16] KernSmooth_2.23-17 utf8_1.2.2 cli_3.2.0
[19] e1071_1.7-9 DBI_1.1.2 ellipsis_0.3.2
[22] class_7.3-17 assertthat_0.2.1 readxl_1.4.0
[25] tibble_3.1.6 lifecycle_1.0.1 crayon_1.5.0
[28] purrr_0.3.4 vctrs_0.3.8 glue_1.6.2
[31] sp_1.4-6 proxy_0.4-26 compiler_4.0.3
[34] pillar_1.7.0 cellranger_1.1.0 generics_0.1.2
[37] classInt_0.4-3 pkgconfig_2.0.3

in_file <- read_sf("hospitals.shp")

out_file <- paste0(tempfile(), ".gpkg")

gdal_utils(

• util = "vectortranslate",
• source = in_file,
• destination = out_file, # output format must be specified for GDAL < 2.3
• options = c("-f", "GPKG")
• ) Error in CPL_gdalvectortranslate(source, destination, options, oo, doo, : Not compatible with STRSXP: [type=list]. `

[Bevann]

I was running into the same problem, based on this post https://gis.stackexchange.com/questions/250165/arcmap-error-drawing-postgis-layer-the-number-of-points-is-less-than-required it seems like a mismatch between the tolerance/precision of the gpkg and Arcgis defaults causes the issue. Even though the shapes are valid in the R environment, when Arcgis trys to draw them with a different tolerance you get invalid shapes. I found that using SF_DATA %>% lwwgeom::st_snap_to_grid(0.1) %>% st_make_valid() fixed this issue.

[rsbivand]

No, the documentation referred to in the SE post is relevant for low-precision representations used before ArcGIS 9.2 (p.19 of http://downloads.esri.com/support/whitepapers/ao_/J9655_Understanding_Coordinate_Management_Esri_Whitepaper.pdf). This says:

High-precision geodatabases store snapped feature coordinates using 53 bits, the equivalent of double-precision storage (8 bytes). The coordinate grid for a dataset is approximately 9.007 x 10^15 grid cells, as shown in figure 21. Very fine distinctions between snapped coordinate values can be maintained because the minimum distance that can be represented between grid mesh points is small. The default XY resolution for a high-precision dataset is 0.0001 meters (1/10 mm) or its equivalent in map units (footnote 8: For unknown coordinate systems, users will have to set XY resolution values appropriate to the type of data without explicitly setting the unit of measure). For example, if the coordinate system units are feet, the default value is approximately 0.0003281 feet; if the coordinate system units are decimal degrees, the default is approximately 8.983 x 10^-10 degrees. Note that the default XY resolution for a coordinate system with decimal degree units is based on the major axis of the GCS spheroid; ArcGIS uses World Geodetic System (WGS) 1984 as the default. In most cases, the default XY resolution will be adequate for most data.

My suggestion without access to Arc is that users have left stale/legacy XY resolution values in their systems. It seems that Arc has supported the double precision coordinate values written by GDAL vector drivers for at least 15 years.

[Bevann]

Yes the 0.1 metre tolerance is much more coarse than the default tolerances in Arcmap, but was sufficient for the data I was working with, but at least in my case the snapping to grid makes it so data loads properly in esri software (ArcgisPro 2.9.5 in my case). Since the goal of my use of gpkg is to share outputs with other groups that can be loaded by a range of software it is preferable to have a product that can be loaded without requiring them modify their default Argis environment settings if that is their software of choice.