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

## 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.

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
``````
• Fantastic, I think this might be it.. Rather than try to amend your proj pipeline, I used gdalwarp to get from the embedded projection (+proj=stere +lat_0=-90 +lon_0=0 +k=1 +x_0=0 +y_0=0 +R=3396000 +units=m +no_defs) to the 'ographic latlong you define (+proj=longlat +a=3396190 +rf=169.894447223612), but it throws up "ERROR 1: Too many points failed to transform, unable to compute output bounds." I assume I'm messing up with the proj syntax somewhere. What's missing?
– Sam
Jun 13, 2020 at 4:35
• If you want to use the strings as CRSes instead of conversion/transformation pipeline, must include the `+type=crs` parameter. Also, make sure that the raster data is georeferenced. Jun 13, 2020 at 12:15

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.