Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

CRS Issue #35459

Closed
tveinot opened this issue Mar 30, 2020 · 13 comments
Closed

CRS Issue #35459

tveinot opened this issue Mar 30, 2020 · 13 comments
Labels
Bug Either a bug report, or a bug fix. Let's hope for the latter! Feedback Waiting on the submitter for answers Projections/Transformations Related to coordinate reference systems or coordinate transformation

Comments

@tveinot
Copy link

tveinot commented Mar 30, 2020

Describe the bug
Working with EPSG 2954 all exports drop to depreciated CRS 2292

How to Reproduce
Export or add any dataset using 2954 and QGIS will assign it 2292. If using a fielgeodatabase QGIS will assign 2954 but any export to shapefile, or any other spatial format I tried, QGIS will write the CRS as 2292

QGIS and OS versions
3.10.4-A Coruña

The attached zip includes a filegeodatabase and a shapefile the shapefile is an export of the same featureclass in the filegeodatabase. You wil notice that the shapefile reads as 2292 and the fieldgeodatabase reads as 2954. When I exported I selected 2954 as the EPSG to use.

This issue extends to everything I use in QGIS in 2954 from geoprocessing outputs to adding layers; QGIS keeps reverting them to 2292.
Tanks for the help
Tyler

GIS.zip

@tveinot tveinot added the Bug Either a bug report, or a bug fix. Let's hope for the latter! label Mar 30, 2020
@gioman gioman added the Projections/Transformations Related to coordinate reference systems or coordinate transformation label Mar 30, 2020
@gioman
Copy link
Contributor

gioman commented Mar 30, 2020

@tveinot what operating system?

@gioman gioman added the Feedback Waiting on the submitter for answers label Mar 30, 2020
@tveinot
Copy link
Author

tveinot commented Mar 30, 2020

Windows 10

@gioman
Copy link
Contributor

gioman commented Mar 30, 2020

This is what ogrinfo (from GDAL 3 installation that comes with QGIS 3.10/12) says about your filegeodatabase layer:

C:\Users\qgis\Desktop\GIS>ogrinfo -so Test.gdb point
INFO: Open of `Test.gdb'
      using driver `OpenFileGDB' successful.

Layer name: point
Geometry: 3D Measured Point
Feature Count: 1
Extent: (389632.467800, 689808.842500) - (389632.467800, 689808.842500)
Layer SRS WKT:
PROJCRS["NAD83(CSRS) / Prince Edward Isl. Stereographic (NAD83)",
    BASEGEOGCRS["NAD83(CSRS)",
        DATUM["NAD83 Canadian Spatial Reference System",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4617]],
    CONVERSION["Prince Edward Isl. Stereographic (NAD83)",
        METHOD["Oblique Stereographic",
            ID["EPSG",9809]],
        PARAMETER["Latitude of natural origin",47.25,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-63,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.999912,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",400000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",800000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (E(X))",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (N(Y))",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["Canada - Prince Edward Island"],
        BBOX[45.9,-64.49,47.09,-61.9]],
    ID["EPSG",2954]]
Data axis to CRS axis mapping: 1,2
FID Column = OBJECTID
Geometry Column = SHAPE

Then you can convert it to SHP and GPKG using ogr2ogr:

C:\Users\qgis\Desktop\GIS>ogr2ogr -f "ESRI shapefile" -t_srs EPSG:2954 as_shape.shp Test.gdb point

C:\Users\qgis\Desktop\GIS>ogr2ogr -f "GPKG" -t_srs EPSG:2954 as_gpkg.gpkg Test.gdb point

and check again the results with ogrinfo:

C:\Users\qgis\Desktop\GIS>ogrinfo -so as_shape.shp as_shape
INFO: Open of `as_shape.shp'
      using driver `ESRI Shapefile' successful.

Layer name: as_shape
Metadata:
  DBF_DATE_LAST_UPDATE=2020-03-30
Geometry: 3D Measured Point
Feature Count: 1
Extent: (389632.467800, 689808.842500) - (389632.467800, 689808.842500)
Layer SRS WKT:
PROJCRS["NAD83(CSRS98) / Prince Edward Isl. Stereographic (NAD83)",
    BASEGEOGCRS["NAD83(CSRS98)",
        DATUM["NAD83 Canadian Spatial Reference System",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4140]],
    CONVERSION["Prince Edward Isl. Stereographic (NAD83)",
        METHOD["Oblique Stereographic",
            ID["EPSG",9809]],
        PARAMETER["Latitude of natural origin",47.25,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-63,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.999912,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",400000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",800000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (E(X))",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (N(Y))",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["Canada - Prince Edward Island"],
        BBOX[45.9,-64.49,47.09,-61.9]],
    ID["EPSG",2292]]
Data axis to CRS axis mapping: 1,2
FID: Integer64 (11.0)
C:\Users\qgis\Desktop\GIS>ogrinfo -so as_gpkg.gpkg point
INFO: Open of `as_gpkg.gpkg'
      using driver `GPKG' successful.

Layer name: point
Geometry: 3D Measured Point
Feature Count: 1
Extent: (389632.467800, 689808.842500) - (389632.467800, 689808.842500)
Layer SRS WKT:
PROJCRS["NAD83(CSRS) / Prince Edward Isl. Stereographic (NAD83)",
    BASEGEOGCRS["NAD83(CSRS)",
        DATUM["NAD83 Canadian Spatial Reference System",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4617]],
    CONVERSION["Prince Edward Isl. Stereographic (NAD83)",
        METHOD["Oblique Stereographic",
            ID["EPSG",9809]],
        PARAMETER["Latitude of natural origin",47.25,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-63,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.999912,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",400000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",800000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (E(X))",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (N(Y))",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["Canada - Prince Edward Island"],
        BBOX[45.9,-64.49,47.09,-61.9]],
    ID["EPSG",2954]]
Data axis to CRS axis mapping: 1,2
FID Column = OBJECTID
Geometry Column = SHAPE

As you can see the CRS definitions differ in the SHP and GPKG this should be the reason of what you are observing, anyway this seems to be a GDAL issue, not a QGIS one.

@roya0045
Copy link
Contributor

@rouault seems that this site states otherwise https://spatialreference.org/ref/epsg/2292/ is it the site that is wrong or could the deprecated status have been inverted by mistake somewhere?

@tveinot
Copy link
Author

tveinot commented Mar 30, 2020

The only difference between EPSG 2292 and EPSG 2954 is the Epoc, 2292 is Epco 98 where as 2954 is Epco 2010. 2292 is suppose to be depreciated and 2954 is the one that is the one we use on the island.
EDIT: OK, according to your link 2954 is from 2002, NR Can says 2954 is using CSRS 2010 and that 2292 is using CSRS 98. So now I am just confused on what is different.
Either way I am absolutely positive that 2954 is not depreciated and 2292 is depreciated.
I appologize for all the confusion.

@jmckenna
Copy link
Contributor

jmckenna commented Mar 30, 2020

2292 was deprecated in the 1990s, and is no longer included in the current EPSG registry. The PROJ (formerly PROJ.4) values of 2292 and 2954 are the exact same. QGIS guesses at the projection, incorrectly picking the deprecated 2292, instead of the correct 2954.

@rouault
Copy link
Contributor

rouault commented Mar 30, 2020

@rouault seems that this site states otherwise https://spatialreference.org/ref/epsg/2292/ is it the site that is wrong or could the deprecated status have been inverted by mistake somewhere?

Generally spatialreference.org should not be considered as authoritative, but here it is indeed true. I wrongly inverted the codes in my previous answer (deleted now to avoid any confusion). So EPSG:2292 has been deprecated by EPSG:2954. And indeed, there's an issue with PROJ. It wrongly identifies PROJCS["NAD_1983_CSRS_Prince_Edward_Island",GEOGCS["GCS_North_American_1983_CSRS",DATUM["D_North_American_1983_CSRS",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Double_Stereographic"],PARAMETER["False_Easting",400000.0],PARAMETER["False_Northing",800000.0],PARAMETER["Central_Meridian",-63.0],PARAMETER["Scale_Factor",0.999912],PARAMETER["Latitude_Of_Origin",47.25],UNIT["Meter",1.0]] to deprecated EPSG:2292, instead of non-deprecated EPSG:2954 . I've filed the issue upstream as OSGeo/PROJ#2116

@jmckenna
Copy link
Contributor

jmckenna commented Mar 30, 2020

@roya0045 best source is to go right to the source: go here and click the 'search by code' tab, to read of the 2292/2954 history https://www.epsg-registry.org

@tveinot
Copy link
Author

tveinot commented Mar 30, 2020

Since this went out as an issue with PROJ I have andother question; So since moving to PEI I have been using 2954 in both Esri and QGIS. Esri says it is a "Double Stereographic Projection" but QGIS says it is a "Stereographic Projection" other sites like epsg.io/2954 say it is "Stereographic" yet The Government of PEI http://www.gov.pe.ca/gis/index.php3?number=77865&lang=J calls it a "Stereographic Double Projection" University of Prince Eward Island calls it a "Double Stereographic Projection" and NR Can doesn't seem to specify. The locals and the surveyors that are using 2954 (we still have a stubbron group using NAD27 Imperial lol), they say it is indeed a Double Stereographic Projection. Should this also be addressed? How do we find out who is right?

@jmckenna
Copy link
Contributor

@tveinot maybe your comment about PROJ should be pasted into the new PROJ ticket. (be sure also to at least click the notification button on that ticket). Also, your question could be a good discussion topic to raise on the dedicated PROJ forum, which exists for this purpose (you ask good questions ha!). https://lists.osgeo.org/mailman/listinfo/proj I personally don't know the answer.

@rouault
Copy link
Contributor

rouault commented Mar 30, 2020

Esri says it is a "Double Stereographic Projection" but QGIS says it is a "Stereographic Projection"

Discussed here: http://geotiff.maptools.org/proj_list/oblique_stereographic.html

@tveinot
Copy link
Author

tveinot commented Mar 30, 2020

@tveinot maybe your comment about PROJ should be pasted into the new PROJ ticket. (be sure also to at least click the notification button on that ticket). Also, your question could be a good discussion topic to raise on the dedicated PROJ forum, which exists for this purpose (you ask good questions ha!). https://lists.osgeo.org/mailman/listinfo/proj I personally don't know the answer.

Thank you, I did as you suggested. To my understanding GeoNB (New Brunswick) is the Geomatics Authority PE has used to establish their CRS. I would say that they are the authority but again I am not too sure.

@gioman
Copy link
Contributor

gioman commented Mar 30, 2020

Closing as upstream issue.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Bug Either a bug report, or a bug fix. Let's hope for the latter! Feedback Waiting on the submitter for answers Projections/Transformations Related to coordinate reference systems or coordinate transformation
Projects
None yet
Development

No branches or pull requests

5 participants