A Close Encounter with Mars

Submitted on Monday, 9/13/2010, at 8:19 PM

Mars at Opposition 2001The planet Mars is currently enjoying its biannual dalliance with the Earth. Marching through the constellation Cancer, the Red Planet's nearness and opposition to the Sun produce an unusually bright spessartine gem against the black velvet sky. It now hangs above the eastern horizon in the evening, where it will remain for several months.

Cydonian "Face" on Mars

Mars has been the subject of intense human scrutiny for centuries, and a battery of NASA probes have visited the planet since the dawn of the space age. The first wave began in the 1960s with the Mariner flybys and orbiters, and ended in the early 1980s with the Viking orbiters and landers. Starting in the late 1990s, a new series of sophisticated craft invaded Mars. Probably most famous are the three robotic rovers deployed to directly explore and analyze its surface (Sojourner, Spirit, and Opportunity). Of more interest for geography are the three orbiting satellites that have remotely imaged the entirety of Mars, returning unprecedented amounts of information and transforming our understanding of the planet.

Gorganum Chaos and Probable Water Gullies

The Mars Global Surveyor, operating from 1999 to 2006, carried a number of scientific instruments, including a high-resolution camera with a horizontal resolution as small as 50 cm, a laser altimeter with a vertical resolution of 30 cm, a thermal emission spectrometer, and a magnetometer, and provided evidence for relatively recent water flow!

Mars Odyssey began studying the material composition of the planet in 2002, using another thermal emission imaging system and a gamma ray spectrometer, and confirmed the existence of subsurface frozen water.

Water Content of Martian Soil

The Mars Reconnaissance Orbiter reached Mars in 2006, and includes a telescopic camera with a resolution around half a meter per pixel. Using stereoimaging, improved measurements of elevation are also being calculated.

Currently the highest resolution topographic data is only available in certain areas. There is so much data that it isn't being processed in its entirety, but only where scientists have been particularly interested!

The best global topographic data currently available comes from the Mars Global Surveyor's Mars Orbiter Laser Altimeter (MOLA), which covers the planet in 1/128th-degree (28 arc second) measurements (and four times more detailed in the polar regions). For example, this "slice" of the planet extends from the north pole on the left to the south pole on the right along the prime meridian, revealing the generally higher elevations in the southern hemisphere:

Mars Topography Along the Prime Meridian

One representation of MOLA data as a Digital Elevation Model (DEM) is in the Planetary Data System (PDS) IMG raster format, which can be viewed with the NASAView software:

http://starbeam.jpl.nasa.gov/tools/nasa-view.shtml

PDS/IMG is somewhat obscure; NASAView can save it as GIF or JPEG format, which is one way to quickly use it in other tools, but a lot of information is lost. Specifically, Martian elevation varies from -8,068 m to +21,134 m (a larger range than on Earth!), so the data is stored as 16-bit signed integers, while GIF and JPEG can only handle 8-bit unsigned integers (grayscale), forcing the output to be scaled and offset into values from 0 to 255.

Not surprisingly, the specialized image analysis program ENVI understands PDS/IMG as an external/generic format, which can be opened and resaved in another format such as ESRI's Grid for use with ArcGIS.

More generally, PDS/IMG is actually a pretty simple format, being a file of uncompressed data with a separate human-readable label file (.lbl) describing its layout. The data can be purely sequential or interleaved (if it's color/multi-band). It's therefore the same as what ESRI ArcGIS and IDL-VIS ENVI call Band Sequential (.bsq), Band Interleaved by Pixel (.bip) or Band Interleaved by Line (.bil). (Note that these formats are identical for grayscale/single-band data.) The only difference, it appears, is that ArcGIS expects a "header" file (.hdr) with a different layout of information. So only the latter file needs to be created, and the data file can simply have its file extension changed. A basic file is described below, but there is also a detailed description of the .hdr format here:

http://webhelp.esri.com/arcgisdesktop/9.3/index.cfm?topicname=BIL,_BIP,_and_BSQ_raster_files

The procedure to view the Mars data is as follows:

  1. The MOLA data is stored here:

    http://pds-geosciences.wustl.edu/missions/mgs/megdr.html

    As an example, the lowest resolution topography can be studied by clicking on the link MEG004, downloading megt90n000cb.img and megt90n000cb.lbl, and renaming the former to megt90n000cb.bsq .
  2. Open the label file with a text editor and read its description. For these files, the data is sampled in geographic coordinates, 0.25°/pixel, and its values are referenced to an estimate of the "areoid" (like sea level but without the water!). The datum used is spherical, with a radius of 3396 Km. With this information you can create an "areographic coordinate system". Note that the Mars standard is to take the Prime Meridian at its Meridian Planum (coincidence?), and to measure longitude from 0° to 360°, which puts the Prime Meridian at the far left of the typical Mars image. In ArcGIS this means assigning the prime meridian to -180°.

    The necessary projection file (e.g. megt90n000cb.prj) will therefore look like this one line:

    GEOGCS["ACS_Mars_MGS_MOLA",DATUM["D_<custom>",SPHEROID["_custom",3396000.0,0.0]],PRIMEM["Meridian Planum",-180.0],UNIT["Degree",0.017453292519943295]]
  3. Using information from the label file, the header file (e.g. megt90n000cb.hdr) can be created with a text editor, and it should look something like the first column of this table (after the first row):

    Corresponding Descriptions of Raster Characteristics
    .hdr parameter Corresponding .lbl parameter(s)
    nrows 720 lines = 720

    Number of rows of pixels in the image

    ncols 1440 line_samples = 1440

    Number of pixels per row (columns) in the image

    skipbytes 0 record_bytes = 2880
    file_records = 720
    lines = 720

    Number of bytes of header information to skip over to get to the image data; it will be padded to the end of a row, so that it will always be an integer multiple of the number of bytes per row:
    record_bytes * (file_records - lines)

    nbands 1 bands = 1

    Number of bands in each pixel (e.g. 1 for grayscale & 3 for RGB)

    nbits 16 sample_bits = 16

    Number of bits in each band (e.g. 32 bits for four-byte values)

    pixeltype unsignedint
    pixeltype signedint
    pixeltype float
    sample_type = *_unsigned_integer
    sample_type = *_signed_integer
    sample_type = *_real

    The * represents a byte order (described next).
    The default value is msb_unsigned_integer.

    byteorder I

    byteorder M
    sample_type = lsb_*
    sample_type = pc_*

    sample_type = msb_*
    sample_type = sun_*

    lsb = least significant byte first = pc = Intel => I
    msb = most significant byte first = sun = Motorola => M
    The * represents a pixel type (described previously).
    The default value is msb_unsigned_integer.

    layout BSQ
    layout BIP
    layout BIL
    band_storage_type = band_sequential
    band_storage_type = sample_interleaved
    band_storage_type = line_interleaved

    These are the same when bands = 1 .

    nodata -3.4028226550889045e+38 missing_constant = 16#FF7FFFFB#

    xdim 0.25 map_resolution = 4.0 <PIX/DEG>

    Width of an easting pixel in map units
    = 1 / map_resolution

    ydim 0.25 map_resolution = 4.0 <PIX/DEG>

    Width of a northing pixel in map units
    = 1 / map_resolution

    ulxmap 0.125
    westernmost_longitude = 0 <DEG>
    map_resolution = 4.0 <PIX/DEG>

    Map coordinate of the center of the first column
    = westernmost_longitude + 1 / map_resolution / 2

    ulymap 89.875 maximum_latitude = 90 <DEG>
    map_resolution = 4.0 <PIX/DEG>

    Map coordinate of the center of the first row
    = maximum_latitude - 1 / map_resolution / 2


    This provides the same information as in a world file, so one isn't necessary for this raster.
  4. At this point you can preview the image in ArcCatalog! Optionally you can build statistics for the image, which will establish the maximum and minimum values and reduce the necessary color range, thereby enhancing the grayscale/color representation. In ArcCatalog, double-click on the data file (e.g. megt90n000cb.bsq) to open its Properties dialog, scroll down to Statistics, click on the button Options, and select Build Statistics….
  5. Finally, using 3D Analyst, create an elevation profile of an interesting location, e.g. the gigantic extinct volcano Olympus Mons. Add the raster you just created to ArcGIS, then menu Tools => Extensions…, and turn on 3D Analyst. Then menu View => Toolbars => 3D Analyst.
    In the 3D Analyst toolbar, click on the button Interpolate Line, and create a polyline along the path where you want to see the profile, by clicking once at each vertex and ending with a double-click. Then click on the button Create Profile Graph.

    The resulting profile for Olympus Mons is below:
Olympus Mons Elevation

An alternative source of data is the Precision Experiment Data Records (PEDR), a point set (Digital Terrain Model/DTM) from which the raster DEM is derived:

http://pds-geosciences.wustl.edu/missions/mgs/pedr.html

As noted above, there are also some bits and pieces of newer, high-resolution Mars DTM available from the Mars Reconnaissance Orbiter, available here:

http://www.uahirise.org/dtm/

This data has an embedded header, which ArcGIS 9.3 is able to partially recognize, but it makes a few mistakes along the way. In particular:

  1. ArcMap understands that the data has a projection of Equirectangular (Plate Carée), but it sets the linear unit incorrectly, and it misinterprets the value described as "CENTER_LATITUDE" as the origin of latitude, rather than as the standard parallel (secant line of the projection cylinder). Creating a projection file that looks like this will fix this:

    PROJCS["Plate_Carree",GEOGCS["EQUIRECTANGULAR MARS",
    DATUM["D_MARS",SPHEROID["MARS",3393830.0,0.0]],PRIMEM["Reference_Meridian",-180.0],
    UNIT["Degree",0.0174532925199433]],PROJECTION["Plate_Carree"],
    PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],
    PARAMETER["central_meridian",341.0599975585938],
    PARAMETER["Standard_Parallel_1",20.0],UNIT["Millimeter",0.001]]


    Note that this data will vary from image to image; be sure to check not only Standard_Parallel_1 and the Central_Meridian, but also the spheroid radius.
  2. ArcMap doesn't properly interpret the NoData value, expressed as a hexadecimal value in the LBL file.  To fix this, load the BSQ (renamed from IMG) file, and then export the layer as a TIFF file; in that dialog, in the field NoData as:, assign the following value: -3.4028226550889045e+38 .
  3. Finally, Arc somehow sees a maximum value that's way too large. In the LBL file, the actual maximum is listed, and you can use that to change the symbology stretch to a proper maximum-minimum range.

The cleanest way to import this data is to use a utility like tail (standard on Unix, an add-on for Windows) to remove the label from the beginning of the data file; with properly constructed .hdr and .prj files as described above, the data is ready to use.

An alternative, more black-box approach is described at ISIS support.

Permalink