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 '20 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 '20 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 '20 at 16:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.