Unexpected geostationary CF conversion #515

djhoese · 2020-01-23T17:04:22Z

Code Sample, a copy-pastable example if possible

A "Minimal, Complete and Verifiable Example" will make it much easier for maintainers to help you:
http://matthewrocklin.com/blog/work/2018/02/28/minimal-bug-reports

from pyproj import CRS                                                                                                                                                                                       
a = 6378169.
b = 6356583.8
h = 35785831.

CRS.from_dict({'proj': 'geos', 'h': h, 'a': a, 'b': b}).to_cf()
# Output
{'crs_wkt': 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6378169,295.488065897001,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Geostationary Satellite (Sweep Y)"],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Satellite Height",35785831,LENGTHUNIT["metre",1,ID["EPSG",9001]]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],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]]]]',
 'grid_mapping_name': 'geostationary',
 'longitude_of_projection_origin': 0,
 'perspective_point_height': 35785831,
 'false_easting': 0,
 'false_northing': 0,
 'semi_major_axis': 6378169,
 'semi_minor_axis': 6356583.8,
 'unit': 'm'}

CRS.from_dict({'proj': 'geos', 'h': h, 'a': a, 'b': b, 'lat_0': 0}).to_cf()
# Output
{'crs_wkt': 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6378169,295.488065897001,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Geostationary Satellite (Sweep Y)"],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Satellite Height",35785831,LENGTHUNIT["metre",1,ID["EPSG",9001]]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],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]]]]',
 'grid_mapping_name': 'geostationary',
 'perspective_point_height': 35785831.0,
 'semi_major_axis': 6378169.0,
 'semi_minor_axis': 6356583.8,
 'latitude_of_projection_origin': 0}

Problem description

As seen in the output above, when lon_0 and lat_0 are missing in the PROJ definition, to_cf is providing a default longitude_of_projection_origin, false_easting, and false_northing, but no latitude_of_projection_origin (default of 0). When we add a lat_0 parameter the longitude_of_projection_origin, false_easting, and false_northing are removed.

Expected Output

I'm not sure. It is not completely unreasonable to me that a descriptive PROJ description has to be provided to get the related CF attributes to come out on the other end. However, there are some parameters for CF that are missing if they are not provided and have known defaults. That said, I understand with the way the CF conversion is implemented there is no way to know how to convert a parameter that doesn't exist.

Any idea if CF considers parameters like latitude_of_projection_origin optional? What about sweep_angle_axis? I didn't see a mention of them as being optional in the CF spec. Maybe they should be?

Just wanted to point this all out in case it was unexpected for you.

The text was updated successfully, but these errors were encountered:

djhoese · 2020-01-23T17:17:52Z

Similarly, if you provide the datum/ellps then the resulting CF has horizontal_datum_name and reference_ellipsoid_name, but does that mean that semi_minor_axis and semi_major_axis are optional?

snowman2 · 2020-01-23T17:46:23Z

@djhoese, the current version of from CF and to CF are definitely imperfect. It is entirely based on the output of CRS.to_dict()/PROJ strings which are not great. This means that if the parameter is not there, it is not included.

If you look at the output of the WKT in both examples, it has all of that information included. I have plans to change this in version 2.5 to be based on the WKT form in import/export as it is more complete. See: https://pyproj4.github.io/pyproj/latest/build_crs.html

djhoese · 2020-01-23T17:59:06Z

I have plans to change this in version 2.5 to be based on the WKT form

Interesting. That makes a lot of sense. As for the how to build a CRS in pyproj 2.5, it makes me a little scared. I'm sure it is more accurate and flexible to specify things one at a time and define each portion (ellipsoid, datum, conversion, etc), but it seems really complicated. I suppose as long as there are some helper functions for converting between things it shouldn't be too bad.

snowman2 · 2020-01-24T01:18:42Z

As for the how to build a CRS in pyproj 2.5, it makes me a little scared.

The good news is that it is not the only way to build a CRS. The new methods are there to give users more options when creating a CRS.

djhoese · 2020-02-12T02:41:25Z

Awesome! Thanks @snowman2!

snowman2 · 2020-02-12T03:17:02Z

Not a problem 👍. If you wouldn't mind giving it a test drive sometime, that would be helpful to iron out some bugs.

djhoese · 2020-02-12T14:31:17Z

I tried my original example from above and get an error for both cases:

In [1]: from pyproj import CRS
   ...: a = 6378169.
   ...: b = 6356583.8
   ...: h = 35785831.
   ...:
   ...: CRS.from_dict({'proj': 'geos', 'h': h, 'a': a, 'b': b}).to_cf()
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-1-9a9b176b025b> in <module>
      4 h = 35785831.
      5
----> 6 CRS.from_dict({'proj': 'geos', 'h': h, 'a': a, 'b': b}).to_cf()

~/miniconda3/envs/satpy_py37/lib/python3.7/site-packages/pyproj/crs/crs.py in to_cf(self, wkt_version, errcheck)
    655         cf_dict.update(
    656             _INVERSE_GRID_MAPPING_NAME_MAP[coordinate_operation.method_name.lower()](
--> 657                 coordinate_operation
    658             )
    659         )

~/miniconda3/envs/satpy_py37/lib/python3.7/site-packages/pyproj/crs/_cf1x8.py in _geostationary__to_cf(conversion)
    413         "sweep_angle_axis": sweep_angle_axis,
    414         "perspective_point_height": params["satellite height"],
--> 415         "latitude_of_projection_origin": params["latitude of natural origin"],
    416         "longitude_of_projection_origin": params["longitude of natural origin"],
    417         "false_easting": params["false easting"],

KeyError: 'latitude of natural origin'

In [2]: CRS.from_dict({'proj': 'geos', 'h': h, 'a': a, 'b': b, 'lat_0': 0}).to_cf()
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-2-edafb7fb94b0> in <module>
----> 1 CRS.from_dict({'proj': 'geos', 'h': h, 'a': a, 'b': b, 'lat_0': 0}).to_cf()

~/miniconda3/envs/satpy_py37/lib/python3.7/site-packages/pyproj/crs/crs.py in to_cf(self, wkt_version, errcheck)
    655         cf_dict.update(
    656             _INVERSE_GRID_MAPPING_NAME_MAP[coordinate_operation.method_name.lower()](
--> 657                 coordinate_operation
    658             )
    659         )

~/miniconda3/envs/satpy_py37/lib/python3.7/site-packages/pyproj/crs/_cf1x8.py in _geostationary__to_cf(conversion)
    413         "sweep_angle_axis": sweep_angle_axis,
    414         "perspective_point_height": params["satellite height"],
--> 415         "latitude_of_projection_origin": params["latitude of natural origin"],
    416         "longitude_of_projection_origin": params["longitude of natural origin"],
    417         "false_easting": params["false easting"],

KeyError: 'latitude of natural origin'

snowman2 · 2020-02-12T15:16:46Z

Ha - I had a geostationary test and thought that was good. I guess I should have tried your specific example 🙃. Thanks for trying it out 👍. Will add a fix later.

snowman2 · 2020-02-13T02:29:14Z

@djhoese, I don't see lat_0 in the options for the geostationary projection: https://proj.org/operations/projections/geos.html

And when I do projinfo it is not there in the WKT form:

./projinfo '+proj=geos +h=35785831.0 +a=6378169.0 +b=6356583.8 +lat_0=1 +type=crs'
PROJ.4 string:
+proj=geos +h=35785831.0 +a=6378169.0 +b=6356583.8 +lat_0=1 +type=crs

WKT2:2019 string:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["unknown",
            ELLIPSOID["unknown",6378169,295.488065897001,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["Geostationary Satellite (Sweep Y)"],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Satellite Height",35785831,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    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]]]]

However, it is mentioned for the geostationary projection for CF:
http://cfconventions.org/cf-conventions/cf-conventions#_geostationary_projection

And it currently works with how I have it setup if you use the CF form:

>>> crs = CRS.from_cf(
...     dict(
...         grid_mapping_name="geostationary",
...         perspective_point_height=1,
...         fixed_angle_axis="y",
...     )
... )
>>> print(crs.coordinate_operation.to_wkt(pretty=True))
CONVERSION["unknown",
    METHOD["Geostationary Satellite (Sweep X)"],
    PARAMETER["Satellite height",1,
        LENGTHUNIT["metre",1,
            ID["EPSG",9001]]],
    PARAMETER["Latitude of natural origin",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8801]],
    PARAMETER["Longitude of natural origin",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8802]],
    PARAMETER["False easting",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8806]],
    PARAMETER["False northing",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8807]]]
>>> crs.to_proj4()
'+proj=geos +sweep=x +lon_0=0 +h=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs'

So, I think leaving it in its current form where it will throw an exception is safer because it drops the lat_0 parameter in the WKT form when created from a PROJ string when provided and I cannot access it through coordinate_operation.params.

>>> crs = CRS(
...     {"proj": "geos", "h": 35785831.0, "a": 6378169.0, "b": 6356583.8, "lat_0": 1}
... )
>>> crs.to_proj4()
'+proj=geos +h=35785831.0 +a=6378169.0 +b=6356583.8 +lat_0=1 +type=crs'
>>> print(crs.coordinate_operation.to_wkt(pretty=True))
CONVERSION["unknown",
    METHOD["Geostationary Satellite (Sweep Y)"],
    PARAMETER["Longitude of natural origin",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8802]],
    PARAMETER["Satellite Height",35785831,
        LENGTHUNIT["metre",1,
            ID["EPSG",9001]]],
    PARAMETER["False easting",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8806]],
    PARAMETER["False northing",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8807]]]

I am going to re-open this issue and see what the PROJ devs have to say about it.

snowman2 · 2020-02-13T02:30:50Z

It also seems to be lost in the coordinate operation PROJ string:

>>> crs.coordinate_operation.to_proj4()
'+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=geos +lon_0=0 +h=35785831 +x_0=0 +y_0=0 +a=6378169 +b=6356583.8'

snowman2 · 2020-02-13T02:38:40Z

Forwarded to: OSGeo/PROJ#1930

snowman2 · 2020-02-13T13:21:44Z

@djhoese, added a fix in #539. Hopefully that should enable your testing to get a little further 😄

djhoese · 2020-02-13T14:59:34Z

Thanks! My examples above don't raise an error any more. I have a pull request in Satpy to switch our "CF" writer to use pyproj for all the projection conversions. That will be the real test for all of this functionality. @mraspaud, if you could use the pyproj master branch with converting some of your omerc projections that would be a good test.

djhoese added the bug label Jan 23, 2020

djhoese mentioned this issue Jan 23, 2020

Switch to pyproj for projection to CF NetCDF grid mapping pytroll/satpy#1038

Merged

5 tasks

snowman2 mentioned this issue Feb 12, 2020

ENH: Update from_cf/to_cf to use WKT instead of PROJ strings for internal management #534

Merged

3 tasks

snowman2 closed this as completed in #534 Feb 12, 2020

snowman2 reopened this Feb 13, 2020

snowman2 mentioned this issue Feb 13, 2020

Question: Geostationary projection and lat_0 OSGeo/PROJ#1930

Closed

snowman2 mentioned this issue Feb 13, 2020

make lat_0 optional for geostationary & add warning for omerc when gamma lost #539

Merged

3 tasks

snowman2 closed this as completed in #539 Feb 13, 2020

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Unexpected geostationary CF conversion #515

Unexpected geostationary CF conversion #515

djhoese commented Jan 23, 2020

djhoese commented Jan 23, 2020

snowman2 commented Jan 23, 2020

djhoese commented Jan 23, 2020

snowman2 commented Jan 24, 2020

djhoese commented Feb 12, 2020

snowman2 commented Feb 12, 2020

djhoese commented Feb 12, 2020 •

edited

Loading

snowman2 commented Feb 12, 2020

snowman2 commented Feb 13, 2020

snowman2 commented Feb 13, 2020

snowman2 commented Feb 13, 2020

snowman2 commented Feb 13, 2020

djhoese commented Feb 13, 2020

Unexpected geostationary CF conversion #515

Unexpected geostationary CF conversion #515

Comments

djhoese commented Jan 23, 2020

Code Sample, a copy-pastable example if possible

Problem description

Expected Output

djhoese commented Jan 23, 2020

snowman2 commented Jan 23, 2020

djhoese commented Jan 23, 2020

snowman2 commented Jan 24, 2020

djhoese commented Feb 12, 2020

snowman2 commented Feb 12, 2020

djhoese commented Feb 12, 2020 • edited Loading

snowman2 commented Feb 12, 2020

snowman2 commented Feb 13, 2020

snowman2 commented Feb 13, 2020

snowman2 commented Feb 13, 2020

snowman2 commented Feb 13, 2020

djhoese commented Feb 13, 2020

djhoese commented Feb 12, 2020 •

edited

Loading