GDAL-Warping MOLA Topography

In response to a post on the Matplotlib-users list I almost automatically responded that, of course, we can use GDAL for that; GDAL can read anything. Being Friday the 13th, I decided that maybe it should actually be tried first, so I downloaded the coarse 4 pixels per degree MOLA (Mars Orbiter Laser Altimeter) topography raster from http://pds-geosciences.wustl.edu/missions/mgs/megdr.html. This is not a format directly supported by GDAL, but is a raw binary file with a header descriptive enough that I can wrap it with a GDAL virtual dataset:

<VRTDataset rasterXSize="1440" rasterYSize="720">
  <GeoTransform>-180.0,0.25,0.0,90.0,0.0,-0.25</GeoTransform>
  <VRTRasterBand dataType="Int16" band="1" subClass="VRTRawRasterBand">
    <SourceFilename relativetoVRT="1">megt90n000cb.img</SourceFilename>
    <ImageOffset>0</ImageOffset>
    <PixelOffset>2</PixelOffset>
    <LineOffset>2880</LineOffset>
    <ByteOrder>MSB</ByteOrder>
  </VRTRasterBand>
 </VRTDataset>

and then warp it to an orthographic projection using the gdalwarp utility:

$ gdalwarp -tr 21875 21875 -te -7000000 -7000000 7000000 7000000 \
  -wo SOURCE_EXTRA=100 -wo SAMPLE_GRID=YES \
  -s_srs "+proj=latlong +ellps=sphere" \
  -t_srs "+proj=ortho +lon_0=0.0 +lat_0=0.0 +ellps=sphere" \
  -dstnodata -32700 topo.vrt topo_ortho.tif

The result (converted to a JPEG):

The functionality of gdalwarp is available through GDAL's gdal.py module as well, which means that we should be able to hook matplotlib's pylab up to GDAL, and thereby to almost any raster data source.