Faulty spatial ref when writing geopackage from sf

Question

I am using sf version 0.9-5 in R 3.6.1 in Windows, and am transforming an sf feature to the USA_Contiguous_Albers_Equal_Area_Conic projection that has ESRI authority 102003 and no EPSG code - I can transform in R:

new <- st_transform(old, "ESRI:102003")

And I can see by plotting data it is correctly projected and the projection information appears correct:

st_crs(test)

Coordinate Reference System:
User input: ESRI:102003 
wkt:
PROJCRS["USA_Contiguous_Albers_Equal_Area_Conic",
BASEGEOGCRS["NAD83",
    DATUM["North American Datum 1983",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["Degree",0.0174532925199433]]],
CONVERSION["USA_Contiguous_Albers_Equal_Area_Conic",
    METHOD["Albers Equal Area",
        ID["EPSG",9822]],
    PARAMETER["Latitude of false origin",37.5,
        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["(E)",east,
        ORDER[1],
        LENGTHUNIT["metre",1]],
    AXIS["(N)",north,
        ORDER[2],
        LENGTHUNIT["metre",1]],
USAGE[
    SCOPE["unknown"],
    AREA["USA - CONUS - onshore"],
    BBOX[24.41,-124.79,49.38,-66.91]],
ID["ESRI",102003]]

But when I write to a geopacakge file:

st_write(new, dsn="new.gpkg", layer = "test")

The geopackage is seen as corrupt by ArcGIS Pro and as not having projection information by QGIS. If, however, I write the sf object in ESRI:102003 projection out to a shapefile, the shapefile is read just fine with it's projection information in both QGIS and ArcGIS Pro.

Is this an issue with non-EPSG code CRSs and geopackages? Is there a workaround? The gpkg_spatial_ref_sys that is written out to geopackage in sf using ESRI:102003 CRS looks like:

srs_name srs_id
1 Undefined cartesian SRS     -1
2 Undefined geographic SRS      0
3 WGS 84 geodetic   4326
4 USA_Contiguous_Albers_Equal_Area_Conic 102003

organization organization_coordsys_id
1 NONE    -1
2 NONE     0
3 EPSG     4326
4 ESRI     102003

definition
1 undefined
2 undefined

3GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
4PROJCS["USA_Contiguous_Albers_Equal_Area_Conic",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",37.5],PARAMETER["longitude_of_center",-96],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["ESRI","102003"]]

description 1 undefined cartesian coordinate reference system 2 undefined geographic coordinate reference system 3 longitude/latitude coordinates in decimal degrees on the WGS 84 spheroid 4

Answer 1

We had almost exactly the same problem when writing to GeoPackage from GeoDMS. The difference was that we used Esri:102022 projection definition. That being said, we were able to solve this issue by changing the projection definition to one that ArcGis provides in the gpkg_spatial_ref_sys field when writing to geopackage:

ogrinfo -sql "UPDATE gpkg_spatial_ref_sys SET definition='PROJCS[\"Africa_Albers_Equal_Area_Conic\",GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Albers\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",25.0],PARAMETER[\"Standard_Parallel_1\",20.0],PARAMETER[\"Standard_Parallel_2\",-23.0],PARAMETER[\"Latitude_Of_Origin\",0.0],UNIT[\"Meter\",1.0],AUTHORITY[\"Esri\",102022]]' WHERE srs_id=102022" path/to/dataset.gpkg

Detailed reason why the initial created geopackage fails can be found in our github discussion section, short version: gdal calls proj to create a OSRSpatialReference in our case from wkt. The wkt1 interpreter fails to interpret, however it is able to recover by using the Esri:102022 definition found in proj.db. This definition is not understood by ArcGis.