faulty spatial ref when writing gpkg

[tghoward]

While this could be strictly a GDAL thing, I'm going through sf and thus it could possibly be an implementation thing here.

In short: writing some projections to geopackage create layers that are unreadable (arcgis) or in the wrong part of the world (gvSig and Qgis) in other systems. More detail:

> library(sf)
Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
Warning message:
package ‘sf’ was built under R version 3.5.3 
> 
> fn <- system.file("shape/nc.shp", package="sf")
> nc <- st_read(fn, quiet = TRUE)
> 
> gpkg <- "temp.gpkg"
> st_write(nc, dsn = gpkg, layer = "nc_LLnad27")
Updating layer `nc_LLnad27' to data source `temp.gpkg' using driver `GPKG'
Writing 100 features with 14 fields and geometry type Multi Polygon.
> 
> nc_aea <- st_transform(nc, 5071)
> st_crs(nc_aea)
Coordinate Reference System:
  EPSG: 5071 
  proj4string: "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
> 
> st_write(nc_aea, dsn = gpkg, layer = "nc_aea")
Updating layer `nc_aea' to data source `temp.gpkg' using driver `GPKG'
Writing 100 features with 14 fields and geometry type Multi Polygon.

Here's what the spatial reference table in the geopackage looks like:

Note the details of the last line (the Albers Equal Area definition):

1. srs_name is "unnamed"
2. srs_id is assigned a sequential integer value from the previous record.
3. the definition also begins with "unnamed"

Back in R, sf sees reference system properly, as noted by the st_crs(nc_aea) call.

The nc_aea layer is what is giving me problems in other systems.

Thanks in advance for any help.

[tghoward]

And trying to do a full circle (reading the gpkg layer back into r), we lose the SRID.

> nc_aea_getBack <- st_read(dsn = gpkg, layer = "nc_aea")
Reading layer `nc_aea' from data source `G:\_Projects\_NM_BLM\HSM_repo\_data\species\temp5.gpkg' using driver `GPKG'
Simple feature collection with 100 features and 14 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 1054293 ymin: 1348025 xmax: 1833499 ymax: 1689236
epsg (SRID):    NA
proj4string:    +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
> st_crs(nc_aea_getBack)
Coordinate Reference System:
  No EPSG code
  proj4string: "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
> 

[edzer]

On sf devel (soon release) I see:

library(sf)
# Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
fn <- system.file("shape/nc.shp", package="sf")
nc <- st_read(fn, quiet = TRUE)
gpkg <- "temp.gpkg"
st_write(nc, dsn = gpkg, layer = "nc_LLnad27")
# Updating layer `nc_LLnad27' to data source `temp.gpkg' using driver `GPKG'
# Updating existing layer nc_LLnad27
# Writing 100 features with 14 fields and geometry type Multi Polygon.
nc_aea <- st_transform(nc, 5071)
st_crs(nc_aea)
# Coordinate Reference System:
#   User input: EPSG:5071 
#   wkt:
# BOUNDCRS[
#     SOURCECRS[
#         PROJCRS["NAD83(HARN) / Conus Albers",
#             BASEGEOGCRS["NAD83(HARN)",
#                 DATUM["NAD83 (High Accuracy Reference Network)",
#                     ELLIPSOID["GRS 1980",6378137,298.257222101,
#                         LENGTHUNIT["metre",1]]],
#                 PRIMEM["Greenwich",0,
#                     ANGLEUNIT["degree",0.0174532925199433]],
#                 ID["EPSG",4152]],
#             CONVERSION["Conus Albers",
#                 METHOD["Albers Equal Area",
#                     ID["EPSG",9822]],
#                 PARAMETER["Latitude of false origin",23,
#                     ANGLEUNIT["degree",0.0174532925199433],
#                     ID["EPSG",8821]],
#                 PARAMETER["Longitude of false origin",-96,
#                     ANGLEUNIT["degree",0.0174532925199433],
#                     ID["EPSG",8822]],
#                 PARAMETER["Latitude of 1st standard parallel",29.5,
#                     ANGLEUNIT["degree",0.0174532925199433],
#                     ID["EPSG",8823]],
#                 PARAMETER["Latitude of 2nd standard parallel",45.5,
#                     ANGLEUNIT["degree",0.0174532925199433],
#                     ID["EPSG",8824]],
#                 PARAMETER["Easting at false origin",0,
#                     LENGTHUNIT["metre",1],
#                     ID["EPSG",8826]],
#                 PARAMETER["Northing at false origin",0,
#                     LENGTHUNIT["metre",1],
#                     ID["EPSG",8827]]],
#             CS[Cartesian,2],
#                 AXIS["easting (X)",east,
#                     ORDER[1],
#                     LENGTHUNIT["metre",1]],
#                 AXIS["northing (Y)",north,
#                     ORDER[2],
#                     LENGTHUNIT["metre",1]],
#             USAGE[
#                 SCOPE["unknown"],
#                 AREA["USA - CONUS - onshore"],
#                 BBOX[24.41,-124.79,49.38,-66.91]],
#             ID["EPSG",5071]]],
#     TARGETCRS[
#         GEOGCRS["WGS 84",
#             DATUM["World Geodetic System 1984",
#                 ELLIPSOID["WGS 84",6378137,298.257223563,
#                     LENGTHUNIT["metre",1]]],
#             PRIMEM["Greenwich",0,
#                 ANGLEUNIT["degree",0.0174532925199433]],
#             CS[ellipsoidal,2],
#                 AXIS["latitude",north,
#                     ORDER[1],
#                     ANGLEUNIT["degree",0.0174532925199433]],
#                 AXIS["longitude",east,
#                     ORDER[2],
#                     ANGLEUNIT["degree",0.0174532925199433]],
#             ID["EPSG",4326]]],
#     ABRIDGEDTRANSFORMATION["NAD83(HARN) to WGS 84 (1)",
#         VERSION["EPSG-Usa"],
#         METHOD["Geocentric translations (geog2D domain)",
#             ID["EPSG",9603]],
#         PARAMETER["X-axis translation",0,
#             ID["EPSG",8605]],
#         PARAMETER["Y-axis translation",0,
#             ID["EPSG",8606]],
#         PARAMETER["Z-axis translation",0,
#             ID["EPSG",8607]],
#         USAGE[
#             SCOPE["Approximation at the +/- 1m level assuming that NAD83(HARN) is equivalent to WGS 84."],
#             AREA["USA - HARN"],
#             BBOX[-14.59,144.58,71.4,-64.51]],
#         ID["EPSG",1580],
#         REMARK["For many purposes NAD83(HARN) can be considered to be coincident with WGS 84."]]]
#Coordinate Reference System:
#  EPSG: 5071
#  proj4string: "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
st_write(nc_aea, dsn = gpkg, layer = "nc_aea")
# Updating layer `nc_aea' to data source `temp.gpkg' using driver `GPKG'
# Updating existing layer nc_aea
# Writing 100 features with 14 fields and geometry type Multi Polygon.
nc_aea_getBack <- st_read(dsn = gpkg, layer = "nc_aea")
# Reading layer `nc_aea' from data source `/tmp/temp.gpkg' using driver `GPKG'
# Simple feature collection with 300 features and 14 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: 1054293 ymin: 1348002 xmax: 1833467 ymax: 1689212
# projected CRS:  NAD83(HARN) / Conus Albers

but this may well be an issue of GDAL/PROJ versions, or mixes of those, which will remain an issue for a while if you use OSX or windows binary CRAN installs.

[tghoward]

Thanks for the reply. If I build sf dev, will I also get the GDAL/PROJ updates? I have an updated version of GDAL installed:

G:\>gdalinfo --version
GDAL 3.0.4, released 2020/01/28

but that's not rgdal (which remained at 2.2.3 using the binary). Do I need to build rgdal from source too? Or, can I point sf to my installed version of GDAL? Sorry if this is spelled out somewhere I just can't seem to find it.
Tim

[Robinlovelace]

If I build sf dev, will I also get the GDAL/PROJ updates?

No, I think you'll need to install the upstream versions (although the version of GDAL seems pretty recent). There is discussion of upstream dependencies on the r-spatial/discuss GitHub org and I think this will get PROJ 7 in Docker: https://github.com/Nowosad/rspatial_proj7

[tghoward]

@Robinlovelace ok, thanks. Note that, yes, while I have GDAL 3.0.4 installed, my install of sf is still using GDAL 2.2.3. I also have PROJ 6.3.1 installed and sf is using 4.9.3. Can I point sf to these newer versions?

Thanks for pointing out the discuss repository. I'm going through that as well.

[tghoward]

Closing this so it doesn't linger and it looks like the issue is resolved down the pipeline. But, for me, sigh, I was hoping to finally make the leap away from shapefiles and it looks like I may need to wait just a little bit longer.

[edzer]

On windows, if you have rtools installed, you may want to experiment with installing #1275

[florisvdh]

It can be confirmed this now works properly on Linux Mint 18.1 with sf at 0.9-1 and using:

> sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H 
       "3.5.1"        "2.2.2"        "4.9.2"         "true"        "false" 

With older sf versions it seems at our institute we have made GPKG files having the problem indeed (link to gpkg-file of below screenshot):

... where all 'unnamed' CRSes should be shown as EPSG 31370. Up to now we re-setted the CRS as 31370 in R functions that read these GPKG files - actually just in order to show back the EPSG-code in the resulting sf object. We did not discover the error until my colleague @w-jan saw an odd response in ArcGIS, complaining about the International 1909 (Hayford) attribute in the GPKG (as seen in the screenshot above).

Further exploration

It seems that the hack in R of re-applying EPSG:31370 will not longer be needed for new GPKG files, as those are correctly read in as 'CRS:31370'. With recent sf, reading those old GPKG files now no longer shows a 'missing' CRS (NA) but only shows the proj4string.

It appears that the WKT parameters in the above (old) GPKG-file, apart from the GEOGCS, DATUM and SPHEROID names and EPSG identifiers, essentially match the WKT of CRS 31370, and the derived proj4strings are exactly the same! (The International 1909 (Hayford) ellipsoid is the same as International 1924).

> old_gpkg_layer <- st_read("old.gpkg"), "old_gpkg_layer", quiet = TRUE)
> st_crs(old_gpkg_layer)
Coordinate Reference System:
  No user input
  wkt:
PROJCS["unnamed",
    GEOGCS["International 1909 (Hayford)",
        DATUM["unknown",
            SPHEROID["intl",6378388,297],
            TOWGS84[-106.8686,52.2978,-103.7239,0.3366,-0.457,1.8422,-1.2747]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",51.16666723333333],
    PARAMETER["standard_parallel_2",49.8333339],
    PARAMETER["latitude_of_origin",90],
    PARAMETER["central_meridian",4.367486666666666],
    PARAMETER["false_easting",150000.013],
    PARAMETER["false_northing",5400088.438],
    UNIT["Meter",1]]

> st_crs(old_gpkg_layer)$proj4string
[1] "+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-106.8686,52.2978,-103.7239,0.3366,-0.457,1.8422,-1.2747 +units=m +no_defs "

> st_crs(31370)
Coordinate Reference System:
  User input: EPSG:31370 
  wkt:
PROJCS["Belge 1972 / Belgian Lambert 72",
    GEOGCS["Belge 1972",
        DATUM["Reseau_National_Belge_1972",
            SPHEROID["International 1924",6378388,297,
                AUTHORITY["EPSG","7022"]],
            TOWGS84[-106.8686,52.2978,-103.7239,0.3366,-0.457,1.8422,-1.2747],
            AUTHORITY["EPSG","6313"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4313"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",51.16666723333333],
    PARAMETER["standard_parallel_2",49.8333339],
    PARAMETER["latitude_of_origin",90],
    PARAMETER["central_meridian",4.367486666666666],
    PARAMETER["false_easting",150000.013],
    PARAMETER["false_northing",5400088.438],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    AUTHORITY["EPSG","31370"]]

> st_crs(31370)$proj4string
[1] "+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-106.8686,52.2978,-103.7239,0.3366,-0.457,1.8422,-1.2747 +units=m +no_defs "

The WKT stored within a newly written GPKG file is the one from st_crs(31370) :

PROJCS["Belge 1972 / Belgian Lambert 72",GEOGCS["Belge 1972",DATUM["Reseau_National_Belge_1972",SPHEROID["International 1924",6378388,297,AUTHORITY["EPSG","7022"]],TOWGS84[-106.8686,52.2978,-103.7239,0.3366,-0.457,1.8422,-1.2747],AUTHORITY["EPSG","6313"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4313"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",51.16666723333333],PARAMETER["standard_parallel_2",49.8333339],PARAMETER["latitude_of_origin",90],PARAMETER["central_meridian",4.367486666666666],PARAMETER["false_easting",150000.013],PARAMETER["false_northing",5400088.438],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","31370"]]