GOES projection and Cartopy

[blaylockbk]

A comment emailed to me about GOES projection and Cartopy from Ilan Gonzalez-Hirshfeld

I wanted to point out a detail regarding how Cartopy assigns ellipsoids when creating projection objects in the cartopy.crs module, which is relevant to the “Using other projections” section of your tutorial. GOES imagery is rectified to the GRS80 ellipsoid (https://www.goes-r.gov/users/docs/PUG-L1b-vol3.pdf). This information is in the GOES NetCDF file and Metpy passes this along to the Cartopy Globe class when you use it to create a cartopy.crs object. Following the variable names used in your workflow, we can see the details of the geos Globe as follows.

In:
vars(geos.globe)


Out:
{'datum': None,
'ellipse': None,
'semimajor_axis': 6378137.0,
'semiminor_axis': 6356752.31414,
'flattening': None,
'inverse_flattening': 298.2572221,
'towgs84': None,
'nadgrids': None}

These are the correct GRS80 parameters. It turns out, however, that when you instantiate a new cartopy.crs object in another projection, Cartopy by default uses the WGS84 ellipsoid (see Globe hyperlink above). If we take a look at the Globe class of the lc cartopy.crs variable from the tutorial we see:

In:
vars(lc.globe)

Out:
{'datum': None,
'ellipse': 'WGS84',
'semimajor_axis': None,
'semiminor_axis': None,
'flattening': None,
'inverse_flattening': None,
'towgs84': None,
'nadgrids': None}

The GRS80 and WGS84 ellipsoids are very similar and use the same semimajor axis. However, they have slightly different inverse flattening values, making the semiminor axis slightly different. We can directly examine the values Cartopy uses for the WGS84 ellipsoid and compare them with the geos.globe object:

In:
from cartopy._crs import WGS84_SEMIMAJOR_AXIS, WGS84_SEMIMINOR_AXIS
print("WGS84 semimajor axis:", WGS84_SEMIMAJOR_AXIS)
print("WGS84 semiminor axis:", WGS84_SEMIMINOR_AXIS)
print("GOES semimajor axis:", vars(geos.globe)['semimajor_axis'])
print("GOES semiminor axis:", vars(geos.globe)['semiminor_axis'])

Out:
WGS84 semimajor axis: 6378137.0
WGS84 semiminor axis: 6356752.3142
GOES semimajor axis: 6378137.0
GOES semiminor axis: 6356752.31414

We see confirmed that while the semimajor axes are identical, the semiminor axes are ever so slightly different. It is unclear to me how Cartopy handles transformations between ellipses/datums, especially since the ‘towgs84’ parameter in the geos.globe object is empty. I suspect in this case Cartopy might simply handle the coordinate transformation as if they were on the same Globe (although I am uncertain). In any case, it’s easy enough to avoid the uncertainly by simply specifying the Globe of whatever new projection object you create. Even if Cartopy is performing an ellipsoid transformation, it’s always best practice to avoid having to transform between datums/ellipses unless it is required. In the above example, we can simply add one argument to define lc as follows:

lc = ccrs.LambertConformal(
  central_longitude=-97.5,
  standard_parallels=(38.5,38.5),
  globe=geos.globe,
)

This way the new cartopy.crs object has the exact same Globe as the native geos object. Of course, with this particular data and with visualizations at this scale, the difference between using GRS80 and WGS84 is entirely unnoticeable. But the geographer in me thinks it is worthwhile to make sure datums are being handled appropriately. In other situations (like calculating very precise geodesic distances, etc.), a miniscule change in the ellipsoid axes can indeed lead to noteworthy discrepancies in downstream calculations.