Martian Cartography Options

[karolp]

I have recently freaked out a little bit about Martian maps of all sorts. And finally I was astonished with those highly detailed beauties that I list below. Nonetheless. some of them have huge inconsistencies (crater names) easily noticed when we compare the surroundings of Gusev crater. Enjoy:

http://www.ralphaeschliman.com/
http://planetologia.elte.hu/1cikkeke.phtml...arsmapinte.html
http://pubs.usgs.gov/imap/i2782/

[Andreas Plesch]

https://pds-imaging.jpl.nasa.gov/data/msl/M...CES_PDS_SIS.PDF

section 3.9.2 has a discussion on easting which I think provides a start of an explanation:

" These “easting” meters at a given latitude are related to true meters at the equator by the simple formula:

map_meters = true_meters / cos(φ)

where φ is the planetocentric latitude. "

[ hey, one can use html entities in the forum ♥ https://www.w3schools.com/html/html_symbols.asp ]

It looks like the provided meters are true meters (at the equator) and the projected easting has map meters since:

4590298.31635 = 4354494.086 / cos(18.4446271 degrees)

This is getting closer to 4590877.824, with still ca. 580m of a difference.

The gdal computed easting assumes a sphere with a radius of 3396190 m. Perhaps the reported true meter easting takes into account a local radius because these are also listed in the geojson, in the radius field.

The gdal computed easting and northing is consistent with the simple formulas given in section 4.4.2.4:

" planetocentric_latitude: Latitude of the point, measured using a planetocentric system.
Planetocentric coordinates are measured as angles from the center of the planet. Latitude (φ_pc) is
computed from northing (x) using the formula:
φ_pc = x / Re • 180 / PI
where Re is the ellipsoid radius, or 3396190 meters.

longitude: Longitude of the point. Longitude is computed from easting (y) using the formula:
θ = y / Re • 180 / PI
"

[ This is good and so simple that I can probably directly plot the traverse geojson on my equirectangular basemap ].

Since the northing matches, Re is 3396190 m. In addition to the correction to true meters from map meters easting, something else seems to be done to the easting. But what ?

[markril]

4590298.31635 = 4354494.086 / cos(18.4446271 degrees)

This is getting closer to 4590877.824, with still ca. 580m of a difference.

Solving for the latitude using the target easting we get:

18.46629998 = acos(4354494.086 / 4590877.824) degrees

Curiously, this number shows up as the center latitude of an orthographic projection in this metadata file. The relevant snippet is this:

CODE

<mapproj>
<mapprojn>Orthographic</mapprojn>
<equirect>
<stdparll>18.4663</stdparll>
<longcm>77.4298</longcm>
<feast>0.0</feast>
<fnorth>0.0</fnorth>
</equirect>
</mapproj>

The TIFF file this refers to is described on this page as the Derived JPL Surface Operations Mosaic:

Derived JPL Surface Operations Mosaic: This mosaic was finalized at JPL specifically for surface operations and mapping. It use the same individual images provide here but optimized for visual purposes using slightly differing seamlines and a 120 pixel blend across the input images. It uses an orthographic map projection centered on the landing site (77.4298 Lon, 18.4663 Lat). Download GeoTIFF (3.2 GB LZW compressed)

Thanks for bringing this up because it looks likes this may have solved a problem I've been working on, too.

Mark

[Andreas Plesch]

...
The TIFF file this refers to is described on this page as the Derived JPL Surface Operations Mosaic:

Glad we arrived at the same conclusion (at the same time !). I am using the TIFF you linked to as my base map. Confusingly, the geotiff embedded projection does have 0 as the Latitude of the 1st standard parallel, somewhat contrary to the description which mentions 'centering' on the landing ellipse:

Unknown CRS: PROJCRS["Equirectangular Mars 2000 Sphere IAU",BASEGEOGCRS["D_Mars_2000_Sphere",DATUM["Mars_2000_(Sphere)",ELLIPSOID["Mars_2000_Sphere_IAU",3396190,0,LENGTHUNIT["metre",1]],ID["ESRI",106971]],PRIMEM["Reference_Meridian",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CONVERSION["Equidistant Cylindrical",METHOD["Equidistant Cylindrical",ID["EPSG",1028]],PARAMETER["Latitude of 1st standard parallel",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["northing",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]] - Projected

And I do think the image data is supposed to be used with that 0 latitude projection because the ISIS .lbl and PDS .lbl also do not specify a standard parallel away from the equator.

It would be possible to reproject (or gdalwarp) the mosaic to the off-equator equirectangular projection and that would yield a less distorted map more appropriate for the region. In effect, it would mean squeezing the x dimension by about 5.15%. But I do not think I am going to do that because it introduces some data loss compared to the existing mosaic.

[markril]

Glad we arrived at the same conclusion (at the same time !). I am using the TIFF you linked to as my base map. Confusingly, the geotiff embedded projection does have 0 as the Latitude of the 1st standard parallel, somewhat contrary to the description which mentions 'centering' on the landing ellipse:

There are two images referenced on that page. The primary image and subject of that page is the equirectangular projection prepared by the USGS:

JEZ_hirise_soc_006_orthoMosaic_25cm_Eqc_latTs0_lon0_first.tif

I agree, the metadata seems to indicate that the 0.25 meters/pixel applies to the equator:

CODE

  Group = Mapping
    TargetName           = Mars
    EquatorialRadius     = 3396190.0 <meters>
    PolarRadius          = 3396190.0 <meters>
    LatitudeType         = Planetocentric
    LongitudeDirection   = PositiveEast
    LongitudeDomain      = 180
    MinimumLatitude      = 18.3067994496865651
    MinimumLongitude     = 77.2229330769714153
    MaximumLatitude      = 18.6693150068777456
    MaximumLongitude     = 77.5839640209350847
    ProjectionName       = Equirectangular
    CenterLongitude      = 0.0
    CenterLatitude       = 0.0
    CenterLatitudeRadius = 3396190.0
    UpperLeftCornerX     = 4577366.0
    UpperLeftCornerY     = 1106618.0
    PixelResolution      = 0.25 <meters/pixel>
    Scale                = 237098.790093224874 <pixels/degree>
  End_Group

The following relations hold:

UpperLeftCornerX = MinimumLongitude * Scale * PixelResolution
UpperLeftCornerY = MaximumLatitude * Scale * PixelResolution
Scale = EquatorialRadius (or PolarRadius or CenterLatitudeRadius) * 2 * Pi / PixelResolution / 360

So, UpperLeftCornerX appears to be an easting based on the equator and not what the JPL operations team is using in the M20_waypoints.json feed. At our latitude of interest (~18.4663 degrees), the horizontal resolution becomes 0.25*cos(18.4663) = 0.237 meters/pixel. So the pixels aren't quite square. It seems desirable to use easting and northing "meters" that are close to equal where you'll be driving your rover so that's where this image comes in:

JEZ_hirise_soc_007_orthoMosaic_25cm_Ortho_blend120.tif

An orthographic projection centered on the target landing site (77.4298 Lon, 18.4663 Lat) where "meters" are nearly equal around the landing site. I suppose they chose orthographic projection (over equirectangular) because any distortion in "meters" as you move from the center of projection will apply to both easting and northing in a symmetrical manner instead of just one dimension as you would have with an equirectangular projection.

The second image seems to be superior from what I've seen, so in order to use it for Google Earth tiles, I'll have to resample it with something like bicubic interpolation to obtain an equirectangular projection with center latitude of 18.4663 degrees.

It would be possible to reproject (or gdalwarp) the mosaic to the off-equator equirectangular projection and that would yield a less distorted map more appropriate for the region. In effect, it would mean squeezing the x dimension by about 5.15%. But I do not think I am going to do that because it introduces some data loss compared to the existing mosaic

.
I'm sure there's already been a little squeezing (and stretching) here and there in order to get these projections from the original HiRISE images, so what's a little more?

Mark