Clarify the projection coordinates of the geostationary projection (CF vs PROJ)

[TomLav]

Dear community,

This issue (re?)starts a discussion about the geostationary projection.

The topic I want to discuss was shortly started in another CF-convention issue led by @erget (cf-convention/cf-conventions#258) but it is better to extract the topic here, and let Issue "258" continue its life-cycle.

Here I would like to discuss the mismatch between -on the one hand- what CF stores as x/y coordinate axis for the geostationary projection (the scanning angles, with unit radians, and standard_name projection_x|y_angular_coordinate) and -on the other hand- what the popular projection package PROJ (proj.org) uses as input/output to its implementation of the Geostationary projection (https://proj.org/operations/projections/geos.html) (the distance from origin of the projection, with unit meters).

NB: This discussion does not aim at re-opening the discussion that CF stores the scanning angles in radians as coordinate axes. This choice was acted / confirmed with this issue.

Once it is accepted that CF stores the scanning angles as radians for this projection, I feel some clarifications are required in the documentations so that non-expert users do not fall in a trap.

Indeed, to my knowledge, the Geostationary projection is the first (and still only) projection supported by CF through grid_mapping that introduces a mis-match between what we store in the CF files as coordinate axis, and what the popular projection tool PROJ uses. Unless we make a clear note, non-expert users (me!) will fall in the trap and use PROJ as with the other projections, which would lead to erroneous results.

For this projection, the chain of coordinate transformations from a (lon,lat) location to a CF (x,y) position is as follow:

1. Transform (lon,lat) to "(X,Y) as meters from projection origin";
2. Transform the "(X,Y) as meters from projection origin" to "(x,y) as scanning angle in radians";

The first step is swiftly handled by PROJ:
echo "0 81" | proj +proj=geos +h=35785831.0 +lon_0=0
0.00 5416145.26

The second typically happens outside PROJ, and is a simple scaling:

scanning_angle (radians) = meters_from_origin (meters) / h (m)

Where h is one of the configuration parameters of the Geostationary projection (h is the viewing point (satellite position) height above the Earth).

To finalize this opening message I list several paths to actions:

1. Do nothing (after all, the experts know what this is about, and are used to perform these conversions).
2. Clarify the CF documentation and make it explicit that there is a mismatch in definition, and that the users/applications should take action.
3. Work with the PROJ team to further clarify their documentation of the projection (https://proj.org/operations/projections/geos.html)
4. Clarify in both CF and PROJ documentations.

My personal view would be to go with 2).

Also, note that Issue "258" already introduces a clarification to the CF document, that should clear some confusion in the current CF-1.9 draft (cf-convention/cf-conventions@bc93215).

[erget]

Dear @TomLav thanks for summarising this issue here so clearly - to add context we'd had a chat about this earlier today on another medium, and it was useful for me to understand more precisely how PROJ works - namely by projecting coordinates between coordinate reference systems (CRS).

I believe the following poitns are also relevant to this discussion:

1. PROJ uses angular coordinates only in the case of lat/lon CRS; all other coordinates are given as linear distances from an origin along well-defined axes. Thus the geostationary projection is a special case, as are other similar projections, such as orthographic. Therefore the behaviour of the PROJ software is consistent and I see no objection against it.
2. The data encoding in CF absolutely makes sense, as this saves all users from having to manage potentially very large arrays of numbers taht would not compress well (which would be the case if we were explicitly storing linear distances from origin in angular coordinates).
3. The documentation both for CF and PROJ is correct, but both are coming from "different worlds" with different needs (CF: provide coordinates for all points in a data set; PROJ: transform sets of individual coordinates from one CRS to another).

In most cases I would say that CF shouldn't include additional text explaining the behaviour of a specific software, but I think that at the end of the day, PROJ is by far the most popular library for converting between CRS that I know of. We also have several other references to PROJ throughout Appendix F, which describes grid mappings, already.

One other note: This change, if we applied it to the geostationary projection, should also be applied to any other projections used in CF that aren't aligned with PROJ. Off the top of my head I believe these are

1. Geostationary projection
2. Orthographic
3. Polar stereographic
4. Stereographic
5. Vertical perspective

It might be worth combining this with a more general cleanup section. For the sake of not having too much work I would advocate applying this only to section, although I will be bringing up some technical housekeeping at the upcoming community meeting.

So in short I would support clarifying the CF documentation IFF we agree to do it consistently across the appendix.

[TomLav]

Hei Daniel,

Yes, it was good talking f2f this morning and clarify where the potential issue is with the geostationary projection in CF.

I agree with all your points (above) except when you bring in the other projections (Orthographic, Polar Stereographic,...). To me there are no issues with these other projections since CF stores exactly what proj uses as input/output (namely the meters from origin along the projection axes). A token of this is that they all specify projection_x|y_coordinate as standard_name for their coordinate axes (and the definition of projection_x|y_coordinate matches what proj uses).

Until proven differently, I still think Geostationary is a special case of CF support for grid_mapping which is why we probably should be a bit (not much) verbose in its documentation.

[erget]

Hi @TomLav , I suspect that the other grid mappings also have the same issue with the standard name as the geostationary projection... But not being an expert on those I would need to consult with some people who are more familiar with them in order to confirm that.

[TomLav]

Status: I propose to wait until PR "Update geostationary projection to allow clean description of newer generation satellites #258" is finalized and merged before re-visiting if some extra clarification is needed of not.