Geodetic transformation with Mars data: 'ographic to 'ocentric latitude

Question

Introduction

I am working with a few raster datasets covering the south polar region of Mars. Some use a coordinate system with aerographic latitudes, which are equivalent in concept to geodetic latitude on Earth (measured relative to the equatorial plane and a line normal to the surface). Others use a coordinate system with aerocentric latitudes, which are like geocentric latitudes (measured relative to the surface and the center of the ellipsoid). In planetary science, these are referred to generically as planetographic and planetocentric coordinate systems.

I would like to apply a geodetic transformation to the 'ographic data so that it will properly align with the 'ocentric data (there can be an offset of several km in the polar regions). However, I have been unable to figure out how to properly define projections for the two systems to use for a transformation.

What I've been able to figure out so far

• Typical ways of defining coordinate systems (PROJ, WKT) do not seem to have a way to explicitly say whether latitude is 'ographic or 'ocentric (i.e., geodetic or geocentric).
• Defining a geocentric coordinate system (such as with a GEOCCS in WKT) means defining X, Y, and Z coordinates rather than lat/lon as I need to do.
• PROJ allows you to convert between geodetic and geocentric latitudes, but it only works on pre-defined ellipsoids for the Earth, not Mars.
  
  

Question

How can I run a geodetic transformation of raster data between 'ographic and 'ocentric latitude systems with a custom ellipsoid?

The way I envision this going is:

1. define a polar stereo projection in which latitude is 'ographic, then assign that to the raster
2. define a polar stereo projection in which latitude is 'ocentric
3. use gdalwarp or some other utility to re-project the 'ographic raster to be 'ocentric.
   
   ... but maybe something else is needed, too.

Answer 1

Assuming PROJ 6.3.1 and GDAL 3.0.4:

---

There is registered in the PROJ database the ESRI:104905 aerographic CRS:

C:\>projinfo ESRI:104905
PROJ.4 string:
+proj=longlat +a=3396190 +rf=169.894447223612 +no_defs +type=crs

WKT2:2019 string:
GEOGCRS["GCS_Mars_2000",
    DATUM["D_Mars_2000",
        ELLIPSOID["Mars_2000_IAU_IAG",3396190,169.894447223612,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Reference_Meridian",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["unknown"],
        AREA["World"],
        BBOX[-90,-180,90,180]],
    ID["ESRI",104905]]

---

To define the Cartesian 'ocentric coordinate system, we can use the +proj=geocent parameter:

C:\>projinfo "+proj=geocent +a=3396190 +rf=169.894447223612 +no_defs +type=crs"
PROJ.4 string:
+proj=geocent +a=3396190 +rf=169.894447223612 +units=m +no_defs +type=crs

WKT2:2019 string:
GEODCRS["unknown",
    DATUM["unknown",
        ELLIPSOID["unknown",3396190,169.894447223612,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]],
    PRIMEM["Reference meridian",0,
        ANGLEUNIT["degree",0.0174532925199433,
            ID["EPSG",9122]]],
    CS[Cartesian,3],
        AXIS["(X)",geocentricX,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(Y)",geocentricY,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(Z)",geocentricZ,
            ORDER[3],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

---

About the Universal Polar Stereographic projection, for the South Pole, we can use the +proj=ups +south parameters:

C:\>projinfo "+proj=ups +south +a=3396190 +rf=169.894447223612 +no_defs +type=crs"
PROJ.4 string:
+proj=ups +south +a=3396190 +rf=169.894447223612 +no_defs +type=crs

WKT2:2019 string:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["unknown",
            ELLIPSOID["unknown",3396190,169.894447223612,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Reference meridian",0,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CONVERSION["unknown",
        METHOD["PROJ ups south"],
        PARAMETER["rf",169.894447223612,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    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]]]]

About the projection from 'ocentric coordinates, it can be performed as a transformation (with a pipeline), but I don't know if we can define the CRS from PROJ.

---

To transform from 'ographic to 'ocentric coordinates, we need to define a pipeline:

+proj=pipeline +step +proj=longlat +a=3396190 +rf=169.894447223612 +step +proj=geoc +a=3396190 +rf=169.894447223612

For instance, the 45 degrees latitude:

C:\>cct +proj=pipeline +step +proj=longlat +a=3396190 +rf=169.894447223612 +step +proj=geoc +a=3396190 +rf=169.894447223612
0 45 0
  0.0000000000   44.6617680466        0.0000           inf

Take into account that 'ographic and 'ocentric coordinate systems may have the latitude as first axis order, so we need to swap the axes:

C:\>cct +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=latlong +a=3396190 +rf=169.894447223612 +step +proj=geoc +a=3396190 +rf=169.894447223612 +step +proj=axisswap +order=2,1
45 0 0
 44.6617680466    0.0000000000        0.0000           inf

We can use the pipeline to transform a raster file with the gdalwarp utility and the -ct parameter. The command could be:

gdalwarp -ct "+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=latlong +a=3396190 +rf=169.894447223612 +step +proj=geoc +a=3396190 +rf=169.894447223612 +step +proj=axisswap +order=2,1" input.tif output.tif

---

Lastly, to project a raster file defined in 'ocentric coordinates to a Universal South Pole Stereographic CRS, we need to go through a 'ographic one:

gdalwarp -ct "+proj=pipeline +step +proj=axisswap +order=2,1 +proj=geoc +inv +a=3396190 +rf=169.894447223612 +step +proj=ups +south +a=3396190 +rf=169.894447223612" input.tif output.tif

Answer 2

With huge thanks to Gabriel De Luca for getting 95% of the way there, here's the method that worked:

Run gdalinfo to get the PROJ string for the input data and use it to start a pipeline:

$ gdalinfo -proj4 ographic_input.tif
...
'+proj=stere +lat_0=-90 +lon_0=0 +k=1 +x_0=0 +y_0=0 +R=3396000 +units=m +no_defs'
...  

$ echo 0 45 0 | cct +proj=pipeline +step +proj=stere +lat_0=-90 +lon_0=0 +k=1 +x_0=0 +y_0=0 +R=3396000 +units=m +no_defs
       0.0000  16397338.5156        0.0000           inf

The +inv parameter can now be added to the source projection so that input coordinates in meters are passed as degrees to a longlat CRS. In this case the spherical radius for +R in the input projection is used to make a Mars-shaped ellipse +a=3396000 +rf=169.894447223612:

$ echo 0 16397338.5156 0 | cct +proj=pipeline +step +inv +proj=stere +lat_0=-90 +lon_0=0 +k=1 +x_0=0 +y_0=0 +R=3396000 +units=m +no_defs +step +proj=longlat +a=3396000 +rf=169.894447223612
0.0000000000   44.9999999999        0.0000           inf  

Now you can tack on +step +proj=geoc +a=3396000 +rf=169.894447223612 to go from the 'ographic coordinates on the sphere to 'ocentric on the ellipse. 45º latitude is now 44.66º:

$ echo 0 16397338.5156 0 | cct +proj=pipeline +step +inv +proj=stere +lat_0=-90 +lon_0=0 +k=1 +x_0=0 +y_0=0 +R=3396000 +units=m +no_defs +step +proj=longlat +a=3396000 +rf=169.894447223612 +step +proj=geoc +a=3396000 +rf=169.894447223612
0.0000000000   44.6617680465        0.0000           inf   

The final step is to convert back to a spherical polar stereo projection modeled after the 'ocentric data that the raster needs to align with. I added +lon_0=-0.0915 to adjust for a change in the definition of Mars's prime meridian. This pipeline is then used for the +ct parameter in gdalwarp.

$ gdalwarp -ct "+proj=pipeline +step +inv +proj=stere +lat_0=-90 +lon_0=0 +k=1 +x_0=0 +y_0=0 +R=3396000 +units=m +no_defs +step +proj=longlat +a=3396000 +rf=169.894447223612 +step +proj=geoc +a=3396000 +rf=169.894447223612 +step +proj=stere +lat_0=-90 +lat_ts=-90 +lon_0=-0.0915 +x_0=0 +y_0=0 +R=3396000 +units=m +no_defs" ographic_input.tif ocentric_output.tif  

I'm not sure why the transformation is needed in the first place given that the 'ographic source and 'ocentric target use spherical datums, but this corrected for the offset perfectly.