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

horizontal_datum_name missing from the CF version of a CRS, breaking conformity with CF ≥ 1.7 #1251

Closed
gerritholl opened this issue Mar 13, 2023 · 4 comments · Fixed by #1255
Labels
Milestone

Comments

@gerritholl
Copy link
Contributor

gerritholl commented Mar 13, 2023

Code Sample, a copy-pastable example if possible

import pyproj
pyproj.CRS(4326).to_cf()

Problem description

According to the documentation, to_cf should return a description conforming to CF version ≥1.8. However, the resulting dictionary does not contain any key for "horizontal_datum_name". According to the CF Convenions, Appendix F, version 1.7 or greater:

reference_ellipsoid_name, prime_meridian_name, horizontal_datum_name and geographic_crs_name must be all defined if any one is defined.

The resulting dictionary does not contain horizontal_datum_name:

{'crs_wkt': 'GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]',
 'semi_major_axis': 6378137.0,
 'semi_minor_axis': 6356752.314245179,
 'inverse_flattening': 298.257223563,
 'reference_ellipsoid_name': 'WGS 84',
 'longitude_of_prime_meridian': 0.0,
 'prime_meridian_name': 'Greenwich',
 'geographic_crs_name': 'WGS 84',
 'grid_mapping_name': 'latitude_longitude'}

Therefore, the description is not conforming with CF 1.8 or greater, despite the documentation claiming so. However, it is conform CF 1.6.

Expected Output

I expect a dictionary that includes a definition for horizontal_datum_name.

Environment Information

pyproj info:
    pyproj: 3.4.1
      PROJ: 9.1.1
  data dir: /data/gholl/mambaforge/envs/py311/share/proj
user_data_dir: /home/gholl/.local/share/proj
PROJ DATA (recommended version): 1.12
PROJ Database: 1.2
EPSG Database: v10.076 [2022-08-31]
ESRI Database: ArcGIS Pro 3.0 [2022-07-09]
IGNF Database: 3.1.0 [2019-05-24]

System:
    python: 3.11.0 | packaged by conda-forge | (main, Jan 14 2023, 12:27:40) [GCC 11.3.0]
executable: /data/gholl/mambaforge/envs/py311/bin/python
   machine: Linux-5.3.18-150300.59.76-default-x86_64-with-glibc2.31

Python deps:
   certifi: 2022.12.7
    Cython: 0.29.33
setuptools: 67.4.0
       pip: 23.0.1

Installation method

  • mamba

Conda environment information (if you installed with conda):


Environment (mamba list):
$ mamba list proj
# packages in environment at /data/gholl/mambaforge/envs/py311:
#
# Name                    Version                   Build  Channel
proj                      9.1.1                h8ffa02c_2    conda-forge
pyproj                    3.4.1           py311h945b3ca_1    conda-forge

Details about mamba and system ( mamba info ):
$ mamba info


                  __    __    __    __
                 /  \  /  \  /  \  /  \
                /    \/    \/    \/    \
███████████████/  /██/  /██/  /██/  /████████████████████████
              /  / \   / \   / \   / \  \____
             /  /   \_/   \_/   \_/   \    o \__,
            / _/                       \_____/  `
            |/
        ███╗   ███╗ █████╗ ███╗   ███╗██████╗  █████╗
        ████╗ ████║██╔══██╗████╗ ████║██╔══██╗██╔══██╗
        ██╔████╔██║███████║██╔████╔██║██████╔╝███████║
        ██║╚██╔╝██║██╔══██║██║╚██╔╝██║██╔══██╗██╔══██║
        ██║ ╚═╝ ██║██║  ██║██║ ╚═╝ ██║██████╔╝██║  ██║
        ╚═╝     ╚═╝╚═╝  ╚═╝╚═╝     ╚═╝╚═════╝ ╚═╝  ╚═╝

        mamba (1.0.0) supported by @QuantStack

        GitHub:  https://github.com/mamba-org/mamba
        Twitter: https://twitter.com/QuantStack

█████████████████████████████████████████████████████████████


     active environment : py311
    active env location : /data/gholl/mambaforge/envs/py311
            shell level : 2
       user config file : /home/gholl/.condarc
 populated config files : /data/gholl/mambaforge/.condarc
                          /home/gholl/.condarc
          conda version : 22.9.0
    conda-build version : not installed
         python version : 3.9.13.final.0
       virtual packages : __cuda=11.4=0
                          __linux=5.3.18=0
                          __glibc=2.31=0
                          __unix=0=0
                          __archspec=1=x86_64
       base environment : /data/gholl/mambaforge  (writable)
      conda av data dir : /data/gholl/mambaforge/etc/conda
  conda av metadata url : None
           channel URLs : https://conda.anaconda.org/conda-forge/linux-64
                          https://conda.anaconda.org/conda-forge/noarch
                          https://conda.anaconda.org/torch/linux-64
                          https://conda.anaconda.org/torch/noarch
          package cache : /data/gholl/mambaforge/pkgs
                          /home/gholl/.conda/pkgs
       envs directories : /data/gholl/mambaforge/envs
                          /home/gholl/.conda/envs
               platform : linux-64
             user-agent : conda/22.9.0 requests/2.28.1 CPython/3.9.13 Linux/5.3.18-150300.59.76-default opensuse-leap/15.3 glibc/2.31
                UID:GID : 16209:16217
             netrc file : None
           offline mode : False
@gerritholl gerritholl added the bug label Mar 13, 2023
@snowman2
Copy link
Member

pyproj/pyproj/crs/crs.py

Lines 714 to 722 in 8109fc4

if self.datum:
cf_dict["horizontal_datum_name"] = self.datum.name
else:
cf_dict["grid_mapping_name"] = "latitude_longitude"
return cf_dict
# handle projected CRS
if self.is_projected and self.datum:
cf_dict["horizontal_datum_name"] = self.datum.name

@djhoese
Copy link
Contributor

djhoese commented Mar 13, 2023

So just above that snippet it looks like crs.coordinate_operation is None so it never enters this if statement. The coordinate_operation is None because:

pyproj/pyproj/_crs.pyx

Lines 2563 to 2567 in 8109fc4

if not (
self.is_bound or self.is_derived
):
self._coordinate_operation = False
return None

where is_bound and is_derived are both False. I'm not sure I understand that since EPSG 4326 definitely has bounds...right?

@snowman2
Copy link
Member

I'm not sure I understand that since EPSG 4326 definitely has bounds...right?

is_bound tells you whether the CRS is a BoundCRS or not. Not necessarily if it has bounds.

@snowman2
Copy link
Member

There is definitely a chance that something is off in the logic here, the code there is just a pointer to where to start digging.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants