Aviso data should be georeferenced using the TOPEX/Poseidon ellipsoid, not the WGS 1984 ellipsoid

[jjr8]

A few years ago I tried to find out from Aviso what reference ellipsoid was used for the Mercator-projected merged satellite products. I could not obtain a satisfactory answer--I could not figure out how to get the person I was communicating with to understand the question. But recently I inquired again, after getting a helpful answer about the problem described in #511.

On 5 January 2012, I received email from Francoise Mertz at Aviso User Services saying "In the DUACS products (SLA, MSLA; Corssh are not concerned) , all the missions (ERS-1, ERS-2, Envisat,GFO) are adjusted on the ellipsoid of the reference mission (TP,J1 and J2) which is: revolution with equatorial radius of 6378.1363 kilometers and a flattening coefficient of 1/298.257."

Those are the parameters for the standard TOPEX/Poseidon ellipsoid, as can be seen in Table 6 of Tapley et al. (1994).

Tapley, B.D, Ries, J.C., Davis, G.W., Eanes, R.J., Schutz, B.E., Shum, C.K., Watkins, M.M., Marshall, J.A., Nerem, R.S., Putney, B.H., Klosko, S.M., Luthcke, S.B., Pavlis, D.E., Williamson, R.G. and Zelensky, N.P. (1994) Precision Orbit Determination for TOPEX/POSEIDON. Journal of Geophysical Research 99: 24383-24404.

In Aviso's netCDF files, the lower-left corner of their global grid is longitude 0, latitude -82. Note that until January 2012, Aviso's FAQ stated:

Question: In the gridded data, to which point lat/lon refer to?
Answer: The lon/lat given is the SW corner of the box, as specified in the metadata of the NetCDF file

I suspected that was wrong and Francoise confirmed my suspicion, saying "indeed, you are right: in the gridded data, the data are stored as points and not boxes that means that the lat/lon refers to the exact point of the given value (and not the SW corner as mentionned in the FAQ)." (This issue may be discussed more in #511.)

So, given that the CornerCoords property stores the x and y coordinates of the center of the lower-left cell, we simply need to transform lon 0, lat -82 into the Mercator system:

>>> import math
>>> from GeoEco.Datasets import Dataset
>>> sr1 = Dataset.ConvertSpatialReference('proj4', '+proj=latlong +a=6378136.3 +rf=298.257 +no_defs', 'obj')              # If you use WGS 1984 instead (i.e. replace "+a=6378136.3 +rf=298.257" with "+ellps=WGS84 +datum=WGS84", the results are the same
>>> sr2 = Dataset.ConvertSpatialReference('proj4', '+proj=merc +a=6378136.3 +rf=298.257 +lon_0=180.0 +no_defs', 'obj')
>>> transformer = Dataset._osr().CoordinateTransformation(sr1, sr2)
>>> xLLC, yLLC = transformer.TransformPoint(0.0, -82.0)[:2]
>>> xLLC, yLLC
(-20037506.143674389, -16925420.02271352)

We can also compute the x and y CoordIncrements:

>>> 6378136.3 * 2 * math.pi / 360 / 3                   # 1/3 degree cell size at equator (1080 columns)
37106.492858656267

MGET's Aviso code should be corrected to use the TOPEX/Poseidon ellipsoid (sr2 in the code above) and the CornerCoords and CoordIncrements I computed above.