Adding the geostationary projection.

[martin.raspaud]

There seems to be a need from the weather satellite people to add the "geos" projection to the list of projections in the CF metadata.

"geos" stands for geostationary. The projection is described in detail in the "CGMS 03, LRIT/HRIT Global specification" (Eumetsat) document. (see here also: ​http://remotesensing.org/geotiff/proj_list/geos.html )

A possible text for this projection is:

Geostationary projection

---

grid_mapping_name = geos

Map parameters:

latitude_of_projection_origin

longitude_of_projection_origin

perspective_point_height

false_easting

false_northing

Map coordinates:

The x (abscissa) and y (ordinate) rectangular coordinates are

identified by the standard_name attribute value projection_x_coordinate and projection_y_coordinate respectively. Notes:

Notes on using the PROJ.4 software packages for computing the

mapping may be found at ​http://www.remotesensing.org/geotiff/proj_list/geos.html . These notes assume the point of observation is directly over the equator. The projection coordinates in this projection are directly related to the scanning angle of the satellite instrument.

[jonathan]

Dear Martin

Thanks for making this proposal. Technically it looks fine to me and I support adding this to Appendix F. It would be helpful if others who are familiar with the projection could also add their support.

Cheers

Jonathan

[caron]

It seems to me the main stumbling block will be getting the x and y coordinate values correct. Your proposal is to put the scale factor into these. So perhaps we need some language to explain this. GRIB-2 defines this projection using Grid Template 3.90, described eg here:

​http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp3-90.shtml

Its worth looking at those parameters, including "apparent diameter of the earth in grid lengths". I think these are, or are part of, this mysterious scale factor needed to figure out the coordinates. The proposal should perhaps explain how one would translate GRIB-2 parameters into the equivilent CF.

(sidenote: the original Eumetsat GRIB-2 encoding had an error in it, since corrected. Simon Eliot is the expert there I think.)

[rhorne@…]

To make this projection work with GOES requires the addition of a "gimbal" parameter to account for GOES' different imager scanning design.

This gimbal parameter will take on a value of "ns_ew" for GOES and "ew_ns" for MSG. This gimbal parameter will provide the mapping algorithm the order in which to apply the coordinate transformation rotations associated with the East to West scanning angle and the North to South elevation angle. Note that I have pre-vetted this with the author (Martin) of this Trac item.)

[caron]

Im concerned that theres not enough info in this text for someone to actually figure out the transformation:

1) the geotiff page has an out of date link to the main document, so we should add the reference in CF:

​http://www.eumetsat.int/idcplg?IdcService=GET_FILE&dDocName=PDF_CGMS_03&RevisionSelectionMethod=LatestReleased

2) perspective_point_height might be from the surface of the earth, although distance from the center of the earth could also be used. this should be clarified.

3) its not clear what the reference implementation is. the closest i see in proj4 is "PJ_geos.c". I have ported the code from

​http://www.eumetsat.int/idcplg?IdcService=GET_FILE&dDocName=zip_tools_msg_nav_c&RevisionSelectionMethod=LatestReleased

from eumetsat. are these the same?

4) clarifying what the x,y coordinates are is essential. If scale_x and scale_y are to be folded into the coordinates, im unclear what they are, eg what units they are in.

[caron]

Replying to rhorne@excaliburlabs.com:

To make this projection work with GOES requires the addition of a "gimbal" parameter to account for GOES' different imager scanning design.

This gimbal parameter will take on a value of "ns_ew" for GOES and "ew_ns" for MSG. This gimbal parameter will provide the mapping algorithm the order in which to apply the coordinate transformation rotations associated with the East to West scanning angle and the North to South elevation angle. Note that I have pre-vetted this with the author (Martin) of this Trac item.)

Hi Randy:

I dont see the "gimbal" parameter in any code I have access to. Do you have code that implements this? If so, perhaps you could post it.

Regards, John

[martin.raspaud]

Replying to caron:

Im concerned that theres not enough info in this text for someone to actually figure out the transformation:

1) the geotiff page has an out of date link to the main document, so we should add the reference in CF:

​http://www.eumetsat.int/idcplg?IdcService=GET_FILE&dDocName=PDF_CGMS_03&RevisionSelectionMethod=LatestReleased

Yes, thanks for that.

2) perspective_point_height might be from the surface of the earth, although distance from the center of the earth could also be used. this should be clarified.

The geos projection in proj.4 uses the distance from the surface of the ellipsoid, while the CGMS document uses the distance from the center of the earth. I would tend to use the former though: to my understanding it fits better with "perspective_point_height".

3) its not clear what the reference implementation is. the closest i see in proj4 is "PJ_geos.c". I have ported the code from

​http://www.eumetsat.int/idcplg?IdcService=GET_FILE&dDocName=zip_tools_msg_nav_c&RevisionSelectionMethod=LatestReleased

from eumetsat. are these the same?

They are not the same, in the sense that the formulas are not literally written in the the C code. However, we wrote the CGMS formulas in a test program, and the values returned are very similar (within 10e-5 meter in precision). I didn't have time to check the math, but I guess the formulas are thus equivalent.

4) clarifying what the x,y coordinates are is essential. If scale_x and scale_y are to be folded into the coordinates, im unclear what they are, eg what units they are in.

I would use add_offset and scale factor, and specify that the x and y coordinates are in radians. Here is the latest example I posted:

dimensions:

y = 328; x = 580; reflective_bands = 5;

variables:

short reflective_band_data(y, x, reflective_bands)

:long_name = "Radiance of reflective bands of the SEVIRI instrument

(MSG)";

:units = "milliwatts/m2/steradian/cm-1"; :grid_mapping = "Projection"; :coordinates = "y x";

short x(x);

:standard_name = "projection_x_coordinate"; :add_offset = 1856 :scale_factor = -208.16554260253906 :units = "rad"

short y(y);

:standard_name = "projection_y_coordinate"; :add_offset = 1856 :scale_factor = -208.16554260253906 :units = "rad"

char Projection;

:grid_mapping_name = "geos"; :perspective_point_height = 35785831.; :semi_major_axis = 6378169.; :semi_minor_axis = 6356583.8; :latitude_of_projection_origin = 0.; :longitude_of_projection_origin = 0.;

[martin.raspaud]

Replying to caron:

Replying to rhorne@excaliburlabs.com:

To make this projection work with GOES requires the addition of a "gimbal" parameter to account for GOES' different imager scanning design.

This gimbal parameter will take on a value of "ns_ew" for GOES and "ew_ns" for MSG. This gimbal parameter will provide the mapping algorithm the order in which to apply the coordinate transformation rotations associated with the East to West scanning angle and the North to South elevation angle. Note that I have pre-vetted this with the author (Martin) of this Trac item.)

Hi Randy:

I dont see the "gimbal" parameter in any code I have access to. Do you have code that implements this? If so, perhaps you could post it.

Regards, John

Hi John,

The gimbal parameter is not implemented anywhere to my knowledge. I plan to implement it in proj.4 though.

I attach the modifications to the CGMS 03 document that Randy sent to me to support the GOES version of the projection.

[martin.raspaud]

I attach the modifications to the CGMS 03 document that Randy sent to me to support the GOES version of the projection.

Looks like the document is too big to attach here (650k). So I can send it by mail to anyone interested.

[caron]

Replying to martin.raspaud:

2) perspective_point_height might be from the surface of the earth, although distance from the center of the earth could also be used. this should be clarified.

The geos projection in proj.4 uses the distance from the surface of the ellipsoid, while the CGMS document uses the distance from the center of the earth. I would tend to use the former though: to my understanding it fits better with "perspective_point_height".

i suppose as long as

distance from the surface of the ellipsoid + earth major axis = distance from the center of the earth

then either should do.

3) its not clear what the reference implementation is. the closest i see in proj4 is "PJ_geos.c". I have ported the code from

​http://www.eumetsat.int/idcplg?IdcService=GET_FILE&dDocName=zip_tools_msg_nav_c&RevisionSelectionMethod=LatestReleased

from eumetsat. are these the same?

They are not the same, in the sense that the formulas are not literally written in the the C code. However, we wrote the CGMS formulas in a test program, and the values returned are very similar (within 10e-5 meter in precision). I didn't have time to check the math, but I guess the formulas are thus equivalent.

thanks, thats a relief.

4) clarifying what the x,y coordinates are is essential. If scale_x and scale_y are to be folded into the coordinates, im unclear what they are, eg what units they are in.

I would use add_offset and scale factor, and specify that the x and y coordinates are in radians. Here is the latest example I posted:

dimensions:

y = 328; x = 580; reflective_bands = 5;

variables:

short reflective_band_data(y, x, reflective_bands)

:long_name = "Radiance of reflective bands of the SEVIRI instrument

(MSG)";

:units = "milliwatts/m2/steradian/cm-1"; :grid_mapping = "Projection"; :coordinates = "y x";

short x(x);

:standard_name = "projection_x_coordinate"; :add_offset = 1856 :scale_factor = -208.16554260253906 :units = "rad"

short y(y);

:standard_name = "projection_y_coordinate"; :add_offset = 1856 :scale_factor = -208.16554260253906 :units = "rad"

char Projection;

:grid_mapping_name = "geos"; :perspective_point_height = 35785831.; :semi_major_axis = 6378169.; :semi_minor_axis = 6356583.8; :latitude_of_projection_origin = 0.; :longitude_of_projection_origin = 0.;

to clarify, from ​BestPractices

"The units attribute applies to unpacked values."

that is, after the scale/offset is applied. is that your intent?

[martin.raspaud]

to clarify, from ​BestPractices

"The units attribute applies to unpacked values."

that is, after the scale/offset is applied. is that your intent?

Yes.

[rhorne@…]

I mentioned this to Martin, but I figure I would post a comment regarding the example above. The scale factor used in the example "scale_factor = -208.16554260253906", does not make sense to me. The individual pixels on the x,y grid from the persepctive of a satellite in space for weather satellites prividing kilometer-scale horizontal spatial resolution will be in tens of microradians. The implication is that a plausible scale_factor would be something like "0.000028", which is a 1 km horizontal spatial resolution.

[caron]

1) just for reference:

unpacked_data_value = packed_data_value * scale_factor + add_offset

the whole scale/offset representation should probably be redacted from the conversation, as it doesnt have anything to do with the transformation.

2) The question is, what should these x,y coordinates be, and how does one generate them??

Ideally they would be something like "km on the projection plane", but then (I think) you need scaling factors that depend on the instrument. If you include the scaling factors into the coordinates, then you have (I guess) radians.

So if anyone has generated these coordinates, can you explain how to generate them? I know of at least one group (besides myself) who needs to understand this better.

[martin.raspaud]

Replying to caron:

1) just for reference:

unpacked_data_value = packed_data_value * scale_factor + add_offset

the whole scale/offset representation should probably be redacted from the conversation, as it doesnt have anything to do with the transformation.

2) The question is, what should these x,y coordinates be, and how does one generate them??

Ideally they would be something like "km on the projection plane", but then (I think) you need scaling factors that depend on the instrument. If you include the scaling factors into the coordinates, then you have (I guess) radians.

So if anyone has generated these coordinates, can you explain how to generate them? I know of at least one group (besides myself) who needs to understand this better.

As far as I am concerned, I understand this as follows: x and y, when packed, are pixel positions on the grid, ie integers between 0 and the size of the image minus one. When unpacked, using add_offset and scale_factor, we get angles, so radians for example. This is what is done in the CGMS document, and seems logical since the satellites are scanning at regular angle intervals. In this case, the scale_factor is thus the angle increment between two scanning position, and the add_offset is the pixel number of nadir.

However, I have to admit it is quite unusual to have angles as input to a projection. For example, the geos projection in the proj4 library takes as input distances on the projection plane (which in this case is actually not a plane, but a curved surface...). So one has to multiply the angles previously mentioned (in radians) by the distance to the surface of the earth (perspective_point_height) to get the input to proj4.

[rhorne@…]

Essentially, this geos projection will provide the capability to map the EW scanning and NS elevation angles from the perspective of the imaging instrument in geostationary satelllite orbit to lat/lon, and also provide the inverse function (i.e. lat/long to EW scanning and NS elevation angles).

The "units" of the x and y coordinate variables that will be natively available by the data production systems will be "radians", as Martin shows in his example. That is, the EW scanning and NS elevation angles will be expressed in "radians". As a result, that is what we on the GOES-R ground system are looking for, and very likely what Martin is looking for with MSG.

To answer your question on how the x and y coordinate variables are generated....

This is based on the characteristics of the imaging instrument and the specification of the "grid" that the level 1 products lie. For all three geostationary weather satellites in the western world, the origin of this grid is at the ideal satellite sub-point (i.e. directly under the satellite corresponding to a specific longitude at the equator.) The grid points are equidistant in terms of viewing angle from the imaging instrument in orbit. However, from a surface of the earth standpoint, these grid points cover increasing amounts of surface area as the distance from the sub-point increases. Note that this is why "radians" is the native unit on this grid.

[martin.raspaud]

the add_offset is the pixel number of nadir.

I meant to say the angle offset at nadir.

[caron]

Replying to martin.raspaud:

the add_offset is the pixel number of nadir.

I meant to say the angle offset at nadir.

Possibly that should be "the angle offset at pixel 0".

[martin.raspaud]

Here is the final proposal. Maybe a moderator could put this in the head of the ticket instead ?

1. Title

Adding the geostationary projection.

2. Moderator

3. Requirement

There seems to be a need from the weather satellite people to add the "geos" projection to the list of projections in the CF metadata.

"geos" stands for geostationary. The projection is described in detail in the "CGMS 03, LRIT/HRIT Global specification" (Eumetsat) document. (see here also: ​http://remotesensing.org/geotiff/proj_list/geos.html )

4. Initial Statement of Technical Proposal

Geostationary projection

The geostationary projection represents the view of earth as scanned by a geostationary satellite.

grid_mapping_name = geostationary

Map parameters:

latitude_of_projection_origin

longitude_of_projection_origin

perspective_point_height

false_easting

false_northing

sweep_angle_axis

fixed_angle_axis

Map coordinates:

The x (abscissa) and y (ordinate) rectangular coordinates are identified by the standard_name attribute value projection_x_coordinate and projection_y_coordinate respectively. In the case of this projection, the projection coordinates in this projection are directly related to the scanning angle of the satellite instrument, and their units are radians.

Notes:

The algorithm for computing the mapping may be found at ​http://www.eumetsat.int/idcplg?IdcService=GET_FILE&dDocName=PDF_CGMS_03&RevisionSelectionMethod=LatestReleased . This document assumes the point of observation is directly over the equator, and that the sweep_angle_axis is y.

Notes on using the PROJ.4 software packages for computing the mapping may be found at ​http://trac.osgeo.org/proj/wiki/proj%3Dgeos and ​http://remotesensing.org/geotiff/proj_list/geos.html .

The "perspective_point_height" is the distance to the surface of the ellipsoid. Adding the earth major axis gives the distance from the centre of the earth.

The "sweep_angle_axis" attribute indicates which axis the instrument sweeps. The value = "y" corresponds to the spin-stabilized Meteosat satellites, the value = "x" to the GOES-R satellite.

The "fixed_angle_axis" attribute indicates which axis the instrument is fixed. The values are opposite to "sweep_angle_axis". Only one of those two attributes are mandatory.

5. Benefits

This benefits to the satellite community, through adding the satellite projection of geostationary satellites.

Use case:

char GOES_R;

:grid_mapping_name = "geostationary"; :perspective_point_height = 35785831.; :semi_major_axis = 6378137.; :semi_minor_axis = 6356752.31414; :latitude_of_projection_origin = 0.; :longitude_of_projection_origin = -137.0; :sweep_angle_axis = "x" :coordinates = "y x" :units = "m"

double x(x) ;

:long_name = "x coordinate of projection" ; :standard_name = "projection_x_coordinate" ; :units = "rad" ;

double y(y) ;

:long_name = "y coordinate of projection" ; :standard_name = "projection_y_coordinate" ; :units = "rad" ;

6. Status Quo

<< describe any ad-hoc approaches that might achieve the same effect without changes to the standard >>

[rhorne@…]

Not sure whether the following issue should be addressed in the notes section of the projection definition ...

I think there is an implicit assumption that the projection x and y coordinates are defined such that the projection y coordinates are positive north of the satellite sub-point and the projection x coordinates are positive east of the satellite sub-point.

Should this implicit assumption be made explicit ?

If so, maybe state something like ...

The sign of the projection x and y coordinates are positive east and north of the satellite subpoint.

[jonathan]

John Caron has given his support to this proposal in a ​comment on the email list. Randy's comment above has also been ​answered on the email list. I support this proposal too, as I said earlier. According to the rules, this ticket will be accepted for the next version of CF if there are no further objections within three weeks. In the next version, Martin Raspaud should be named in the list of additional authors of the CF convention.

Jonathan

[martin.raspaud]

This doesn't seem to be included in the 1.7 draft document... How can I help to make that happen ?

BR, Martin

[rhorne@…]

Dear All:

This is a follow-up to Martin's most recent post with the hope it will solicit a response capturing the status of adding this projection to CF.

very respectfully,

randy

[rkaiser]

Folks,

GOES-16 has launched and needs this projection added in order to be CF compliant. Can we change the state of this ticket to "accepted" and get the suggested text added to v 1.7?

Rob Kaiser (follow-on to Randy Horne in this area)

[jonathan]

Dear Rob

There are many tickets that have been agreed, including this one, but unfortunately we have not yet tackled implementing them all in 1.7, although the document is now in a form which can be updated. Although it's not in the 1.7 working version yet, it will be included in 1.7, so please proceed as though it was already there.

Best wishes

Jonathan