geographic_crs_name missing from the CF version of a CRS

[mraspaud]

Code Sample, a copy-pastable example if possible

from pyproj import CRS                                                                                                                
crs = CRS.from_proj4('+proj=omerc +alpha=8.998112717187938 +lat_0=0.0 +lonc=13.809602948622212 +gamma=0 +ellps=sphere')               
print(crs.to_cf())                                                                                                                           

Problem description

When using the latest cfchecker (4.0.0) on a netcdf file containing the output of crs.to_cf, we get an error:

------------------
Checking variable: oblique_mercator
------------------
INFO: (5.6): CF checker currently does not verify the syntax of the crs_wkt attribute which must conform to the CRS WKT specification
ERROR: (5.6): reference_ellipsoid_name, prime_meridian_name, horizontal_datum_name and geographic_crs_name must all be definied if any one is defined

Indeed, the output from to_cf in this case is:

{'crs_wkt': 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unknown based on Normal Sphere (r=6370997) ellipsoid",ELLIPSOID["Normal Sphere (r=6370997)",6370997,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Hotine Oblique Mercator (variant B)",ID["EPSG",9815]],PARAMETER["Latitude of projection centre",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8811]],PARAMETER["Longitude of projection centre",13.8096029486222,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8812]],PARAMETER["Azimuth of initial line",8.99811271718794,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8813]],PARAMETER["Angle from Rectified to Skew Grid",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8814]],PARAMETER["Scale factor on initial line",1,SCALEUNIT["unity",1],ID["EPSG",8815]],PARAMETER["Easting at projection centre",0,LENGTHUNIT["metre",1],ID["EPSG",8816]],PARAMETER["Northing at projection centre",0,LENGTHUNIT["metre",1],ID["EPSG",8817]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
 'semi_major_axis': 6370997.0,
 'semi_minor_axis': 6370997.0,
 'inverse_flattening': 0.0,
 'reference_ellipsoid_name': 'Normal Sphere (r=6370997)',
 'longitude_of_prime_meridian': 0.0,
 'prime_meridian_name': 'Greenwich',
 'horizontal_datum_name': 'Unknown based on Normal Sphere (r=6370997) ellipsoid',
 'grid_mapping_name': 'oblique_mercator',
 'latitude_of_projection_origin': 0.0,
 'longitude_of_projection_origin': 13.809602948622212,
 'azimuth_of_central_line': 8.998112717187938,
 'scale_factor_at_projection_origin': 1.0,
 'false_easting': 0.0,
 'false_northing': 0.0}

As you can see, geographic_crs_name is missing. As the CF document says:

Corresponds to a OGC WKT GEOGCS node name

and I don't see one in the wkt attached.

So in conclusion, with this projection definition, it looks like to me that the to_cf method isn't generating a CF compliant grid mapping.

Expected Output

I would expect geographic_crs_name to be present in the keys of the returned dict.

Environment Information

pyproj info:
pyproj: 2.6.0
PROJ: 6.3.1
data dir: /home/a001673/miniconda3/envs/py3/lib/python3.8/site-packages/pyproj/proj_dir/share/proj

System:
python: 3.8.2 | packaged by conda-forge | (default, Mar 23 2020, 18:16:37) [GCC 7.3.0]
executable: /home/a001673/miniconda3/envs/py3/bin/python
machine: Linux-4.18.0-147.5.1.el8_1.x86_64-x86_64-with-glibc2.10

Python deps:
pip: 20.0.2
setuptools: 46.1.3.post20200325
Cython: 0.29.14

Installation method

• conda

Conda environment information (if you installed with conda):

Environment (conda list):

$ conda list | grep -E "proj|aenum"
proj                      6.3.1                hc80f0dc_1    conda-forge
pyproj                    2.6.0            py38h95a477d_0    conda-forge

Details about conda and system ( conda info ):

$ conda info
     active environment : py3
    active env location : /home/a001673/miniconda3/envs/py3
            shell level : 2
       user config file : /home/a001673/.condarc
 populated config files : /home/a001673/.condarc
          conda version : 4.8.3
    conda-build version : not installed
         python version : 3.7.4.final.0
       virtual packages : __glibc=2.28
       base environment : /home/a001673/miniconda3  (writable)
           channel URLs : https://conda.anaconda.org/conda-forge/linux-64
                          https://conda.anaconda.org/conda-forge/noarch
                          https://repo.anaconda.com/pkgs/main/linux-64
                          https://repo.anaconda.com/pkgs/main/noarch
                          https://repo.anaconda.com/pkgs/r/linux-64
                          https://repo.anaconda.com/pkgs/r/noarch
          package cache : /home/a001673/miniconda3/pkgs
                          /home/a001673/.conda/pkgs
       envs directories : /home/a001673/miniconda3/envs
                          /home/a001673/.conda/envs
               platform : linux-64
             user-agent : conda/4.8.3 requests/2.23.0 CPython/3.7.4 Linux/4.18.0-147.5.1.el8_1.x86_64 rhel/8.1 glibc/2.28
                UID:GID : 62310:2000
             netrc file : None
           offline mode : False

[snowman2]

The code is here:

pyproj/pyproj/crs/crs.py

Lines 583 to 653 in 57eeaf5

	unknown_names = ("unknown", "undefined")
	cf_dict = {"crs_wkt": self.to_wkt(wkt_version)} # type: Dict[str, Any]
	# handle bound CRS
	if (
	self.is_bound
	and self.coordinate_operation
	and self.coordinate_operation.towgs84
	):
	sub_cf = self.source_crs.to_cf(errcheck=errcheck) # type: ignore
	sub_cf.pop("crs_wkt")
	cf_dict.update(sub_cf)
	cf_dict["towgs84"] = self.coordinate_operation.towgs84
	return cf_dict
	# handle compound CRS
	elif self.sub_crs_list:
	for sub_crs in self.sub_crs_list:
	sub_cf = sub_crs.to_cf(errcheck=errcheck)
	sub_cf.pop("crs_wkt")
	cf_dict.update(sub_cf)
	return cf_dict
	# handle vertical CRS
	elif self.is_vertical:
	vert_json = self.to_json_dict()
	if "geoid_model" in vert_json:
	cf_dict["geoid_name"] = vert_json["geoid_model"]["name"]
	if self.datum and self.datum.name not in unknown_names:
	cf_dict["geopotential_datum_name"] = self.datum.name
	return cf_dict
	# write out datum parameters
	if self.ellipsoid:
	cf_dict.update(
	semi_major_axis=self.ellipsoid.semi_major_metre,
	semi_minor_axis=self.ellipsoid.semi_minor_metre,
	inverse_flattening=self.ellipsoid.inverse_flattening,
	)
	if self.ellipsoid.name not in unknown_names:
	cf_dict["reference_ellipsoid_name"] = self.ellipsoid.name
	if self.prime_meridian:
	cf_dict["longitude_of_prime_meridian"] = self.prime_meridian.longitude
	if self.prime_meridian.name not in unknown_names:
	cf_dict["prime_meridian_name"] = self.prime_meridian.name
	# handle geographic CRS
	if self.geodetic_crs and self.geodetic_crs.name not in unknown_names:
	cf_dict["geographic_crs_name"] = self.geodetic_crs.name
	if self.is_geographic:
	if self.coordinate_operation:
	cf_dict.update(
	_INVERSE_GEOGRAPHIC_GRID_MAPPING_NAME_MAP[
	self.coordinate_operation.method_name.lower()
	](self.coordinate_operation)
	)
	if self.datum and self.datum.name not in unknown_names:
	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 and self.datum.name not in unknown_names:
	cf_dict["horizontal_datum_name"] = self.datum.name
	coordinate_operation = None
	if not self.is_bound and self.is_projected:
	coordinate_operation = self.coordinate_operation
	if self.name not in unknown_names:
	cf_dict["projected_crs_name"] = self.name

It would be a simple fix to remove the check for a name being in unknown_names. IIRC, I think I decided to do this because of the logic to create a CRS using CF uses the name if the other parameters are not supplied and the name is given:

pyproj/pyproj/crs/crs.py

Lines 728 to 731 in 57eeaf5

	elif geographic_crs_name:
	geographic_crs = CRS(geographic_crs_name)
	else:
	geographic_crs = GeographicCRS()

Anyway, suggestions & thoughts are definitely welcome 👍

Also, as far as not being CF compliant, also see #401

[snowman2]

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

Hmmm, what about if it isn't projected and you have the geographic_crs_name?

[snowman2]

Hmmm, what about if it isn't projected and you have the geographic_crs_name?

Re-read it, looks like projected_crs_name is not required in that case.

[mraspaud]

Ok, so how about:

 # handle geographic CRS 
 if self.geodetic_crs: 
     cf_dict["geographic_crs_name"] = self.geodetic_crs.name 

and

 elif geographic_crs_name and geographic_crs_name not in unknown_names: 
     geographic_crs = CRS(geographic_crs_name) 
 else: 
     geographic_crs = GeographicCRS()

?

[snowman2]

That sounds like a good plan to me.
Also, asking for clarification here:
cf-convention/cf-conventions#255

[mraspaud]

I'll bake a small PR soon then, thanks !