Skip to content
New issue

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

#

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

Already on GitHub? # to your account

Jump to bottom

Proj fails with OGC WKT given in CF-convention specifications #4321

Closed
observingClouds opened this issue Nov 12, 2024 · 4 comments
Closed

Proj fails with OGC WKT given in CF-convention specifications #4321

observingClouds opened this issue Nov 12, 2024 · 4 comments
Labels
bug

Comments

@observingClouds
Copy link

I have a particular issue with a WKT CRS from the CF-conventions, particularly the example 5.12. I want to read it with pyproj, but get a pyproj.exceptions.CRSError: (Internal Proj Error: proj_create: Missing BASEGEODCRS / BASEGEOGCRS / GEOGCS node)

Upon a quick inspection of the OGC WKT specifications 18-010r11 I could not find any hint that those values are always required, so I raise this issue here first (instead of mentioning it to the CF-conventions community) in the hope that someone could comment on whether this WKT is actually invalid or not.

I previously raised this issue in pyproj but since this issue is caused by proj and not pyproj itself, this seems the more appropriate place.

Example of problem

import pyproj

ogc_wkt = """COMPOUNDCRS [
    "OSGB 1936 / British National Grid + ODN",
    PROJCRS [
        "OSGB 1936 / British National Grid",
        GEODCRS [
            "OSGB 1936",
            DATUM [
                "OSGB 1936",
                ELLIPSOID [
                    "Airy 1830",
                    6377563.396,
                    299.3249646,
                    LENGTHUNIT ["metre", 1.0]
                ],
                TOWGS84 [375, -111, 431, 0, 0, 0, 0]
            ],
            PRIMEM ["Greenwich", 0],
            UNIT ["degree", 0.0174532925199433]
        ],
        CONVERSION [
            "OSGB",
            METHOD [
                "Transverse Mercator",
                PARAMETER ["False easting", 400000, LENGTHUNIT ["metre", 1.0]],
                PARAMETER ["False northing", -100000, LENGTHUNIT ["metre", 1.0]],
                PARAMETER ["Longitude of natural origin", -2.0, ANGLEUNIT ["degree", 0.0174532925199433]],
                PARAMETER ["Latitude of natural origin", 49.0, ANGLEUNIT ["degree", 0.0174532925199433]],
                PARAMETER ["Longitude of false origin", -7.556, ANGLEUNIT ["degree", 0.0174532925199433]],
                PARAMETER ["Latitude of false origin", 49.766, ANGLEUNIT ["degree", 0.0174532925199433]],
                PARAMETER ["Scale factor at natural origin", 0.9996012717, SCALEUNIT ["Unity", 1.0]],
                AUTHORITY ["EPSG", "27700"]
            ],
            CS [Cartesian, 2],
            AXIS ["easting (X)", east],
            AXIS ["northing (Y)", north],
            LENGTHUNIT ["metre", 1.0]
        ]
    ],
    VERTCRS [
        "Newlyn",
        VDATUM ["Ordnance Datum Newlyn", 2005],
        AUTHORITY ["EPSG", "5701"],
        CS [vertical, 1],
        AXIS ["gravity-related height (H)", up],
        LENGTHUNIT ["metre", 1.0]
    ]
]"""

crs = pyproj.crs.CRS.from_wkt(ogc_wkt)

This results in:

CRSError: Invalid projection: COMPOUNDCRS ["OSGB 1936 / British National Grid + ODN", PROJCRS ["OSGB 1936 / British National Grid", GEODCRS ["OSGB 1936", DATUM ["OSGB 1936", ELLIPSOID ["Airy 1830", 6377563.396, 299.3249646, LENGTHUNIT["metre",1.0]], TOWGS84[375, -111, 431, 0, 0, 0, 0]], PRIMEM ["Greenwich", 0], UNIT ["degree", 0.0174532925199433]], CONVERSION["OSGB", METHOD["Transverse Mercator", PARAMETER["False easting", 400000, LENGTHUNIT["metre",1.0]], PARAMETER["False northing", -100000, LENGTHUNIT["metre",1.0]], PARAMETER["Longitude of natural origin", -2.0, ANGLEUNIT["degree",0.0174532925199433]], PARAMETER["Latitude of natural origin", 49.0, ANGLEUNIT["degree",0.0174532925199433]], PARAMETER["Longitude of false origin", -7.556, ANGLEUNIT["degree",0.0174532925199433]], PARAMETER["Latitude of false origin", 49.766, ANGLEUNIT["degree",0.0174532925199433]], PARAMETER["Scale factor at natural origin", 0.9996012717, SCALEUNIT["Unity",1.0]], AUTHORITY["EPSG", "27700"]] CS[Cartesian,2], AXIS["easting (X)",east], AXIS["northing (Y)",north], LENGTHUNIT["metre", 1.0]], VERTCRS ["Newlyn", VDATUM ["Ordnance Datum Newlyn", 2005], AUTHORITY ["EPSG", "5701"] CS[vertical,1], AXIS["gravity-related height (H)",up], LENGTHUNIT["metre",1.0]]]]: (Internal Proj Error: proj_create: Missing BASEGEODCRS / BASEGEOGCRS / GEOGCS node)
Full traceback 
CRSError                                  Traceback (most recent call last)
Cell In[3], line 1
----> 1 pyproj.crs.CRS.from_wkt(ogc_wkt)

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/pyproj/crs/crs.py:434, in CRS.from_wkt(cls, in_wkt_string)

432 if not is_wkt(in_wkt_string):

433     raise CRSError(f"Invalid WKT string: {in_wkt_string}")

--> 434 return cls.from_user_input(_prepare_from_string(in_wkt_string))


File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/pyproj/crs/crs.py:503, in CRS.from_user_input(cls, value, **kwargs)

501 if isinstance(value, cls):

502     return value

--> 503 return cls(value, **kwargs)


File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/pyproj/crs/crs.py:350, in CRS.init(self, projparams, **kwargs)

348     self._local.crs = projparams

349 else:

--> 350     self._local.crs = _CRS(self.srs)


File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/pyproj/_crs.pyx:2364, in pyproj._crs._CRS.init()


CRSError: Invalid projection: COMPOUNDCRS ["OSGB 1936 / British National Grid + ODN", PROJCRS ["OSGB 1936 / British National Grid", GEODCRS ["OSGB 1936", DATUM ["OSGB 1936", ELLIPSOID ["Airy 1830", 6377563.396, 299.3249646, LENGTHUNIT["metre",1.0]], TOWGS84[375, -111, 431, 0, 0, 0, 0]], PRIMEM ["Greenwich", 0], UNIT ["degree", 0.0174532925199433]], CONVERSION["OSGB", METHOD["Transverse Mercator", PARAMETER["False easting", 400000, LENGTHUNIT["metre",1.0]], PARAMETER["False northing", -100000, LENGTHUNIT["metre",1.0]], PARAMETER["Longitude of natural origin", -2.0, ANGLEUNIT["degree",0.0174532925199433]], PARAMETER["Latitude of natural origin", 49.0, ANGLEUNIT["degree",0.0174532925199433]], PARAMETER["Longitude of false origin", -7.556, ANGLEUNIT["degree",0.0174532925199433]], PARAMETER["Latitude of false origin", 49.766, ANGLEUNIT["degree",0.0174532925199433]], PARAMETER["Scale factor at natural origin", 0.9996012717, SCALEUNIT["Unity",1.0]], AUTHORITY["EPSG", "27700"]] CS[Cartesian,2], AXIS["easting (X)",east], AXIS["northing (Y)",north], LENGTHUNIT["metre", 1.0]], VERTCRS ["Newlyn", VDATUM ["Ordnance Datum Newlyn", 2005], AUTHORITY ["EPSG", "5701"] CS[vertical,1], AXIS["gravity-related height (H)",up], LENGTHUNIT["metre",1.0]]]]: (Internal Proj Error: proj_create: Missing BASEGEODCRS / BASEGEOGCRS / GEOGCS node)

Environment Information

  • PROJ version (proj): Rel. 9.5.0, September 15th, 2024
  • Operation System Information: macOS
  • pyproj: pyproj.version: '3.7.0'

Installation method
@jjimenezshaw
Copy link
Contributor

You should use projinfo to test it. It will give you some clues.
For instance, the geographic systems should be BASEGEODCRS and not GEODCRS. The CS is not part of the CONVERSION. TOWGS84 is not in WKT2, etc

@jjimenezshaw
Copy link
Contributor

jjimenezshaw commented Nov 12, 2024 •
edited
Loading

And looks like you want to use this CRS
projinfo EPSG:27700+5701 or EPSG:7405
https://spatialreference.org/ref/epsg/7405/wkt2.html

@rouault
Copy link
Member

rouault commented Nov 12, 2024 •
edited
Loading

The above WKT is invalid.

The minimum fix for it is:

COMPOUNDCRS [
    "OSGB 1936 / British National Grid + ODN",
    PROJCRS [
        "OSGB 1936 / British National Grid",
        BASEGEODCRS [
            "OSGB 1936",
            DATUM [
                "OSGB 1936",
                ELLIPSOID [
                    "Airy 1830",
                    6377563.396,
                    299.3249646,
                    LENGTHUNIT ["metre", 1.0]
                ]
            ],
            PRIMEM ["Greenwich", 0],
            UNIT ["degree", 0.0174532925199433]
        ],
        CONVERSION [
            "OSGB",
            METHOD ["Transverse Mercator"],
            PARAMETER ["False easting", 400000, LENGTHUNIT ["metre", 1.0]],
            PARAMETER ["False northing", -100000, LENGTHUNIT ["metre", 1.0]],
            PARAMETER ["Longitude of natural origin", -2.0, ANGLEUNIT ["degree", 0.0174532925199433]],
            PARAMETER ["Latitude of natural origin", 49.0, ANGLEUNIT ["degree", 0.0174532925199433]],
            PARAMETER ["Longitude of false origin", -7.556, ANGLEUNIT ["degree", 0.0174532925199433]],
            PARAMETER ["Latitude of false origin", 49.766, ANGLEUNIT ["degree", 0.0174532925199433]],
            PARAMETER ["Scale factor at natural origin", 0.9996012717, SCALEUNIT ["Unity", 1.0]]
        ],
        CS [Cartesian, 2],
        AXIS ["easting (X)", east],
        AXIS ["northing (Y)", north],
        LENGTHUNIT ["metre", 1.0],
        ID["EPSG",27700]
    ],
    VERTCRS [
        "Newlyn",
        VDATUM ["Ordnance Datum Newlyn"],
        CS [vertical, 1],
        AXIS ["gravity-related height (H)", up],
        LENGTHUNIT ["metre", 1.0],
        ID["EPSG",5701]
    ]
]

And clean WKT:2015 output is given by:

$ bin/projinfo EPSG:27700+5701 -o WKT2:2015 -q
COMPOUNDCRS["OSGB36 / British National Grid + ODN height",
    PROJCRS["OSGB36 / British National Grid",
        BASEGEODCRS["OSGB36",
            DATUM["Ordnance Survey of Great Britain 1936",
                ELLIPSOID["Airy 1830",6377563.396,299.3249646,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]]],
        CONVERSION["British National Grid",
            METHOD["Transverse Mercator",
                ID["EPSG",9807]],
            PARAMETER["Latitude of natural origin",49,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8801]],
            PARAMETER["Longitude of natural origin",-2,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8802]],
            PARAMETER["Scale factor at natural origin",0.9996012717,
                SCALEUNIT["unity",1],
                ID["EPSG",8805]],
            PARAMETER["False easting",400000,
                LENGTHUNIT["metre",1],
                ID["EPSG",8806]],
            PARAMETER["False northing",-100000,
                LENGTHUNIT["metre",1],
                ID["EPSG",8807]]],
        CS[Cartesian,2],
            AXIS["(E)",east,
                ORDER[1],
                LENGTHUNIT["metre",1]],
            AXIS["(N)",north,
                ORDER[2],
                LENGTHUNIT["metre",1]],
        SCOPE["Engineering survey, topographic mapping."],
        AREA["United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore."],
        BBOX[49.75,-9.01,61.01,2.01],
        ID["EPSG",27700]],
    VERTCRS["ODN height",
        VDATUM["Ordnance Datum Newlyn"],
        CS[vertical,1],
            AXIS["gravity-related height (H)",up,
                LENGTHUNIT["metre",1]],
        SCOPE["Geodesy, engineering survey."],
        AREA["United Kingdom (UK) - Great Britain onshore - England and Wales - mainland; Scotland - mainland and Inner Hebrides."],
        BBOX[49.93,-7.06,58.71,1.8],
        ID["EPSG",5701]]]

This should be reported to the authors of the CF conventions

@rouault rouault closed this as completed Nov 12, 2024
@observingClouds
Copy link
Author

Thanks @rouault and @jjimenezshaw for your amazingly quick and super helpful response.

@observingClouds observingClouds mentioned this issue Nov 13, 2024
Closed
# for free to join this conversation on GitHub. Already have an account? # to comment
Assignees
No one assigned
Labels
bug
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
3 participants
@rouault @jjimenezshaw @observingClouds