Question: Geostationary projection and lat_0

[snowman2]

Based on: pyproj4/pyproj#515

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 even though it persists in the PROJ 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 as latitude_of_projection_origin:
http://cfconventions.org/cf-conventions/cf-conventions#_geostationary_projection

And it currently works in the WKT 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'

But, it gets lost in the coordinate operation:

>>> crs.coordinate_operation.to_proj4()
'+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=geos +sweep=x +lon_0=0 +h=1 +x_0=0 +y_0=0 +ellps=WGS84'

I am definitely not an expert on this, so any insight into the matter would be helpful.

Thanks!

[snowman2]

It does not seem to exist in the definition here either: https://spatialreference.org/ref/sr-org/geostationary-satellite-view/

[rouault]

(spatialreference.org should not be considered at all. It has not been updated in > 10 years)

The WKT definition here is based on what GDAL 2 did, and it ignored the Latitude of origin (the netCDF driver ignores the latitude_of_projection_origin CF parameters).
In PROJ world, lat_0 can always been used, even if not mentioned in the doc, as it is a generic mechanism fir all projections.
I believe that omitting lat_0 for geos makes sense: this projection is supposed to be used by a geostationary satellite, and you can only put a geostationary satellite at the equator, and assuming it points to nadir, you have thus lat_0 = 0.
Are there use cases with latitude_of_projection_origin != 0 ?

(the fact that lat_0=1 is preserved in the PROJ.4 CRS string is due to the magic you've probably experienced a bit lately when hacking in those areas. The fact that it is lost when using the Conversion object is because the original PROJ.4 string is a hidden property of the CRS object, not the Conversion one)

[snowman2]

@rouault, thanks for the response! That makes sense and helps me towards a sensible default of 0 when the latitude of origin is not included.

The fact that it is lost when using the Conversion object is because the original PROJ.4 string is a hidden property of the CRS object, not the Conversion one

Interesting, so it gets dropped in the coordinate operation, even though it is a valid PROJ string option?

>>> from pyproj import CRS
>>> crs = CRS("+proj=geos +sweep=x +lat_0=0 +lon_0=0 +h=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs")
>>> crs.to_proj4()
'+proj=geos +sweep=x +lat_0=0 +lon_0=0 +h=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs'
>>> crs.coordinate_operation.to_proj4()
'+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=geos +sweep=x +lon_0=0 +h=1 +x_0=0 +y_0=0 +ellps=WGS84'

Are there use cases with latitude_of_projection_origin != 0 ?

@djhoese, do you have an example where this is the case?

The WKT definition here is based on what GDAL 2 did

Would it make sense to add support for it in to be explicit?

Another option, if it is included, would it be to create a PROJ-based coordinate operation:

./projinfo '+proj=geos +sweep=x +lat_0=0 +lon_0=0 +h=1 +x_0=0 +y_0=0'
PROJ string:
+proj=geos +sweep=x +lat_0=0 +lon_0=0 +h=1 +x_0=0 +y_0=0

WKT2:2019 string:
CONVERSION["PROJ-based coordinate operation",
    METHOD["PROJ-based operation method: +proj=geos +sweep=x +lat_0=0 +lon_0=0 +h=1 +x_0=0 +y_0=0"]]

[djhoese]

Are there use cases with latitude_of_projection_origin != 0 ?

@djhoese, do you have an example where this is the case?

No. I think this is one of those where it should never not be 0 until someone goes and does it. I don't think it is needed but if the math in the projection code is there to handle it it could be nice for future-proofing.

[kbevers]

In PROJ world, lat_0 can always been used, even if not mentioned in the doc, as it is a generic mechanism fir all projections.

The first part is true, the second isn't really. It is correct that PJ->phi0 is always initialized when creating a new projection but many of the projections don't use them. The docs should reflect that, so that +lat_0 is only mentioned when it actually affects the projection.

There may be a tiny bug hidden here though, since setting +lat_0=1 the parameter is marked as in-use internally (in PJ->paralist) although the actual projection isn't using it. That may be why it is displayed by projinfo. I can demonstrate this by using proj:

$ proj -v +proj=geos +h=35785831.0 +a=6378169.0 +b=6356583.8 +lat_0=1 +dummy
#Geostationary Satellite View
#	Azi, Sph&Ell
#	h=
# +proj=geos +h=35785831.0 +a=6378169.0 +b=6356583.8 +lat_0=1
#--- following specified but NOT used
# +dummy

Ideally both +dummy and +lat_0=1 should show up in the "NOT used" category.

[rouault]

PJ->phi0 is always initialized when creating a new projection but many of the projections don't use them

Ah you're right of course. I was confused with lam0 ...

[kbevers]

I was confused with lam0 ...

The same goes for PJ->lam0 although that is much rarer. Calcofi is an example of a projection that ignores the +lon_0 parameter.

[snowman2]

Insightful comment here: cf-convention/cf-conventions#246 (comment)

Most geostationary products are produced by satellites that are, in fact, geosynchronous - the satellite doesn't stay consistently above a single point on the Earth, as this would expend a lot of fuel, but rather "wobbles" in a figure-eight about its nominal point. This wobble can actually take the satellite several degrees above and below the 0° latitude line.

[snowman2]

Thanks for all the insights! Going to close now.