1

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

3
  • Does the shapefile use the same WKT string or is it an old-style WKT string? We're (Esri) updating our WKT2 format handling in ArcGIS Pro 2.7. Note: I work for Esri.
    – mkennedy
    Nov 20, 2020 at 3:08
  • 1
    Could you add the contents of the corresponding row in the gpkg_spatial_ref_sys table of the GeoPackage?
    – user30184
    Nov 20, 2020 at 8:19
  • @mkennedy, I'm writing out to both geopackage and shapefiles in sf which now uses WKT2. @user30184 I'm new to geopackages but I've edited to my question above to show what the gpkg_spatial_ref_sys table currently looks like when writing out the sf object in R to a geopackage in ESRI:102003 CRS.
    – mweber
    Nov 20, 2020 at 16:59

1 Answer 1

0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.