How do I write out a GeoTIFF with a coordinate reference system for Mars (e.g., IAU_Mars_2000_Sphere)?

I am trying to write out a digital elevation model (DEM) as a geoTIFF, with a coordinate reference system for Mars. For Earth, I am able to do:

outfile = 'synthetic_dem.tif';
geotiffwrite(outfile, int16(Z), Rgeo, 'CoordRefSysCode', 4326);  % EPSG:4326 = WGS84 (Earth)

But I am not sure how to do this for Mars, since, as far as I can tell, there is no CoordRefSysCode for Mars. Something like IAU_Mars_2000_Sphere would be ideal ( https://spatialreference.org/ref/esri/104971/wkt.html ). I am new to these systems, so I may be mixing projections and coordinate reference systems.

Thank you for your help!

[Walter Roberson]

I have no idea how to do this, but when I look around a bit, I find references that indicate that it can be done. I find specific mention of GeoTIFF format being used for Mars rover data.

[Aditya Khuller]

Thanks, Walter. I may have managed to write out the TIF using the GC MARS 2000 ellipsoid geographic coordinate reference system: https://spatialreference.org/ref/esri/104905/#

gcrs = geocrs(104905,"Authority","ESRI");
Rgeo.GeographicCRS = gcrs;
outfile = 'synthetic_dem.tif';
geotiffwrite(outfile, int16(Z), Rgeo);  

But for some reason, when I try and read this file in using readgeoraster, the geographic coordinate reference system shows up as the default WGS84...

[Z,R] = readgeoraster('synthetic_dem.tif');

Answer 1

[Moritz Tenthoff]

I had the exact same problem this week but for Mercury instead of Mars.

geotiffwrite defaults to 'EPSG:4326' if you don't provide the CoordRefSysCode Name-Value argument. Unfortunately, you cannot provide ESRI or IGNF codes like you can for projcrs or geocrs (Although I found that IAU_2015 codes would be more useful here for celestial bodies). EPSG codes are only for Earth.

That leaves two options:

1. Define the coordinate reference system manually via the GeoKeyDirectoryTag argument. That means manually building a struct with all required keys. I've found this resource for the Key IDs, but it's pretty incomprehensible.
2. Use an external program to fill in the correct information. This is what I chose after I didn't find a way to make geotiffwrite do it. I used gdal (you can execute gdal_translate with system from within MATLAB, which is nice). Inputs are the geotiff as saved by MATLAB, the WKT of your reference system and a world file specifiying your grid parameters.

The code could look something like this:

gcrs = geocrs(104905,"Authority","ESRI");
Rgeo.GeographicCRS = gcrs;
outfile = 'synthetic_dem.tif';
geotiffwrite(outfile, int16(Z), Rgeo); % You are sure you want to convert to integer? Geotiff can save floating point values.

wkt = wktstring(gcrs); % get the well-known text for your geocrs
fid = fopen('wkt.txt','w'); % write to file (or manually enter in gdal cli)
fprintf(fid,wkt);
fclose(fid);
worldfilewrite(Rgeo,replace(outfile,'.tif','.tfw')); % write world file with the same filename (this is how gdal will find it)

% Run gdal_translate from MATLAB (Windows)
gdalCommand = sprintf("gdal_translate -a_srs %s -of GTiff %s %s",'wkt.txt',outfile,'synthetic_dem_fixed.tif');
[status, cmdout] = system(gdalCommand);

% Reload
[Z,R] = readgeoraster('synthetic_dem_fixed.tif'); % R.GeographicCRS should now have the correct object

Note: If you are on Windows and have installed gdal via conda you will have to activate the conda environment first in the command.

[Jacob Halbrooks]

Here is some code that can help with trying the first option:

crs = geocrs(104971,"Authority","ESRI");

tags.GTModelTypeGeoKey = 2;
tags.GTRasterTypeGeoKey = 1;
tags.GeogAngularUnitsGeoKey =  9102;
tags.GeogCitationGeoKey =  convertStringsToChars("GCS Name = " + crs.Name + "|Datum = " + crs.Datum + "|Ellipsoid = " + crs.Spheroid.Name + "|Primem = Reference_Meridian|");
tags.GeogGeodeticDatumGeoKey = 32767;
tags.GeogEllipsoidGeoKey = 32767;
tags.GeogSemiMinorAxisGeoKey = crs.Spheroid.SemiminorAxis;
tags.GeogSemiMajorAxisGeoKey = crs.Spheroid.SemimajorAxis;
tags.GeogPrimeMeridianGeoKey = 32767; 
tags.GeographicTypeGeoKey = 32767; 


filename = "synthetic_dem.tif";
geotiffwrite(filename, Z, Rgeo, "GeoKeyDirectoryTag", tags)

[Walter Roberson]

@Jacob Halbrooks

The 32767 are almost certainly "magic numbers", that probably mean something like "other" or "extended attributes".

But is the 9102 a magic number along the same lines, or does it vary according to the planet being used?

[Jacob Halbrooks]

The number 9102 means degrees, as indicated in the angular units section of the GeoKey field names page here.

[Walter Roberson]

Ah, thank you.

[Aditya Khuller]

Thank you so much, Walter, Moritz, and Jacob! I think that, using option 1 for the Mars 2000 ellipsoid seems to work for my use case:

crs = geocrs(104905,"Authority","ESRI"); % GCS_Mars_2000

tags = struct();

% Model & raster
tags.GTModelTypeGeoKey      = 2;      % 2 = Geographic
tags.GTRasterTypeGeoKey     = 1;      % 1 = PixelIsArea

% Geographic CRS: user-defined
tags.GeographicTypeGeoKey   = 32767;  % user-defined
tags.GeodeticDatumGeoKey    = 32767;  % user-defined datum
tags.GeogEllipsoidGeoKey    = 32767;  % user-defined ellipsoid

% Units & prime meridian
tags.GeogAngularUnitsGeoKey = 9102;   % 9102 = degree
tags.PrimeMeridianLongGeoKey = 0;     % 0° (Reference_Meridian)

% Ellipsoid parameters from Mars_2000_IAU_IAG
tags.GeogSemiMajorAxisGeoKey = crs.Spheroid.SemimajorAxis;      
tags.GeogInvFlatteningGeoKey = crs.Spheroid.InverseFlattening;  

% Citation
tags.GeogCitationGeoKey = convertStringsToChars( ...
    "GCS Name = " + crs.Name + ...
    "|Datum = " + crs.Datum + ...
    "|Ellipsoid = " + crs.Spheroid.Name + ...
    "|Primem = Reference_Meridian|" );

outfile = 'synthetic_DEM.tif';
geotiffwrite(outfile, single(Z), Rgeo, "GeoKeyDirectoryTag", tags);