Lab - Imaging & Altimetry

Ge151: Spring 2009

This lab requires that you use certain programs on the gps computers. You may use the computers in the lab in the basement of S. Mudd. Login with 'ge151' and use the password given out in class.

Alternatively, you can remotely log into the ge151 account from your own computer if you have one with the appropriate programs. From the UNIX prompt use:

ssh -X ge151@machinename.gps.caltech.edu

Available machines are:tharsis, iupika, hermes, europa, fizeau, xanth

Source the .cshrc file once logged in with:

source ~/.cshrc

Imaging Techniques & Analysis

During this lab you will work through:

• How to find, retrieve and process raw imaging data.
• Apply this to some area of Mars you are interested in.
• Some more advanced imaging techniques.

We're going to start with a pre-selected image to demonstrate some software packages and how they're used. The image which we'll be working with is stored in the /home/ge151/lab3 directory, the filename is m1000782.imq. This is an image taken in the south polar region of Mars in December 1999 by the Mars Observer Camera (MOC), part of the ongoing Mars Global Surveyor (MGS) mission.

Introduction and Background

This image has been prepared for release by the Planetary Data System (PDS). These people are responsible for providing a consistent format for all current and future planetary missions (they're slowly converting the older ones as well) known unsurprisingly as PDS format. PDS images have either an .img or an .imq suffix (.imq is simply a compressed version of .img). They have a text header with information specific to that image followed by the binary image data itself.

Log on to the ge151 account with the password given out in class. First, you need to create a folder to keep all of our source files and the output files you'll be creating. (This is where I'll look to grade your assignments). Move into the lab3 directory and use the 'mkdir command to create your own subdirectory (mkdir directoryname). Copy each of the files in /lab3 to your directory. (cp filename directoryname).

Now we're ready to look at some images. From your subdirectory, take a look at the image file either in a text editor or by typing 'less m1000782.imq' at the unix prompt. There are other files there but ignore them for now. The information contained in the header usually relates to things which won't change such as the time the image was taken, the gain state of the camera or the number of lines and samples in the image. Other information such as the position of the spacecraft relative to the planet, the longitude and latitude of the image or the size of the pixels in meters is subject to change as better solutions for the orbit of the spacecraft are developed. This other information is not part of the image file as it would be a huge effort to revise the whole dataset each time these new and improved numbers are derived.

MOC images are identified by an 8-digit name; the first 3 characters denote the mission phase m10 in this case means the ninth mapping phase. Each mapping phase lasts roughly one calendar month. Other prefixes exist such as ab1, sp2, fha and cal which stand for aerobraking phase 1, science phasing 2, full high gain antenna and calibration phase respectively. These other phases took place early in the mission and were pretty short. The remaining 5 characters represent the orbit number (first 3 characters) and the image number (last 2 characters) within that orbit. So in the case of our image m1000782, it is the 82^nd image taken on the 7^th orbit in mapping phase 9.

MOC is a line scan camera i.e. it has one row of detectors that sweeps across the surface of the planet as the spacecraft moves in its orbit.The resolution along this line (known as crosstrack resolution) is set by the height of the orbit i.e. if the spacecraft is far away from the planet each pixel sees more of the surface. The mapping orbit was planned to be ~400 Km high, however there is 35 Km of relief on the planet (from the summit of Olympus Mons to the bottom of the Hellas impact basin) so the crosstrack resolution can vary from place to place on the planet (by almost 10%). The downtrack resolution is determined by the groundspeed of the spacecraft (~ 3 Km s^-1) combined with the length of time the camera exposes each line. Ideally you want the spacecraft to move only a fraction of a pixel (to prevent smearing) in the downtrack direction during an exposure, however the exposure must be long enough to collect sufficient light for a high signal to noise ratio. The camera sensitivity, crosstrack resolution and mapping orbit parameters have been designed so that the pixels have roughly the same resolution in the cross and down track directions. Changes in distance to the surface and surface reflectivity mean that the pixels are not entirely square so all MOC images have some non-unity aspect ratio. Sometimes this can be quite severe and must be removed before the image makes any sense. In general it is always nicer to correct the aspect ratio when looking at images (there is something fundamentally disturbing about elliptical craters).

MOC images are always taken at the full resolution of the camera (~1.4 meters/pixel for the narrow angle) but to cope with the voluminous amounts of data the spacecraft computer intentionally degrades the resolution by summing and averaging pixels in both the cross and down track directions. This summing is not necessarily the same in each direction so in addition to the inherent aspect ratio there is sometimes an induced aspect ratio due to this differential summing.

ISIS (Integrated Software for Imagers and Spectrometers)

ISIS is a software package developed by the United States Geological Survey (USGS) for use with spacecraft data. ISIS is really a collection of stand-alone programs, which perform operations on datasets such as map-projecting an image onto a sphere. In this case we are going to use some of the programs in the ISIS package to calibrate our image.

ISIS has its own environment from which you can call these programs known as tae. Its possible to call them as stand alone programs however there is little benefit in that extra complication here. Log on to the ge151 account and make sure you are in the correct directory. Enter the tae environment by typing 'tae' at the unix prompt. Some text blurb should appear along with the cursor at the tae prompt.

There are two ways to call ISIS programs from here. Each program needs some set of parameters that you have to feed it e.g. the input file or what format to put the output in etc… If you know all the parameters you want to use in advance you can just type the program name followed by the parameters and their values and that will run the program. If (in our case) you don't know everything in advance then you can use the user-friendlier 'tutor mode'. To enter tutor mode just type t, a space and then the program name e.g. to run the program banana in tutor mode enter 't banana' at the tae prompt and press return.Tutor mode will show you all the parameters for your chosen command and allow you to fill them in individually and get help on each one.

ISIS works with its own image format called 'cubes' with a .cub extension. An ISIS cube contains header information like the PDS format does. The first step will be to convert this image into the ISIS format. This is known, as level 0 processing i.e. it isn't really processing at all but just data translation. There is a script called moclev0 that will run all the necessary ISIS programs to do this conversion.

Type

't moclev0'

Use the down arrow to cycle through the different parameters, there aren't too many since this is a pretty straightforward procedure. Use the arrow keys to move to the 'FROM' field. Enter the filename, 'm1000782.imq' and press enter. The other fields are fine as they are. Type run and press enter (make sure you are not entering run as a parameter somewhere :) ). The ISIS programs should whirl seamlessly into action. Text messages will scroll past, as these programs run letting you know what's going on.

OK, now we can take a look at this image. You should be back at the tae prompt now. ISIS has created a file called m1000782.lev0.cub; this is the raw data in ISIS format. Run the program 'qview' (no need for tutor mode this time, just type the program name at the tae prompt and press enter). Qview is a cube-viewing program, after it loads use it to open the new cube file, which you just made (File -> Open Display), just click the 'ok' button for the load options.

Pretty ugly! It's streaky, got some weird garbage lines in the middle, looks like it has an aspect ratio problem and looks pretty uninspiring in general. However you will turn this into a thing of beauty and scientific worth in a few easy steps.

You can use the magnifying glass to zoom in (left button), zoom to image resolution (middle button) and zoom out (right button). Take a moment to look over the image zooming up on the interesting stuff. The world button on the upper right will bring you back to the state you started in if you need it. The other buttons do other things not of interest to us at the moment apart from the 'MD' button (5^th from the left). This is the 'doctor' button and lets you edit the image. Before that can happen though go to the options menu and select doctor options. On the panel that comes up select the rectangle option on top and the null option on the bottom. Ok those choices and go back to the image display window. Zoom up on the garbage lines one or twice (you should still be able to see the entire line in the window, zoom out if you can't). Select the MD button and drag a rectangle around the bad lines they should disappear. Go to the file menu and exit qview, no need to save the file has been automatically updated.

Getting rid of those garbage lines was necessary to make the rest of the processing go smoothly. You should be back at the tae prompt now. We'll run the next program without tutor mode since all we need to provide it with is a filename and we know how to do that. We'll use the program 'moclev1' next. As before this calls several ISIS programs in succession, we need only tell it what the input file is.

Type

moclev1 from="m1000782.lev0.cub"

This takes the level zero file we made and converts it to a level one cube file. Level one means that the image has been radiometricaly corrected i.e. the previous DN values have been converted into meaningful I/F values (see tutorial one for the explanation of what I/F means). More text blurb will scroll past; basically listing the parameters used in the DN --> I/F conversion. When the dust has settled and you have the tae prompt back again, startup qview as before. There should be a new cube file there called m1000782.lev1.cub, open it up and take a look. It should be looking a whole lot better now; the streakiness was due to different pixels on the line array having different efficiencies. The bright ones were very efficient at recording light and the darker one were not. The calibration process has taken account of that and 'flattened' all the columns. Press the button on the top right marked with the 'paw' symbol. This turns on the information display at the bottom left, pan the mouse cursor around the image and watch the numbers in the bottom left.

This is still not the ideal situation the image is still distorted because of its aspect ratio and we have no idea of what the scale is and which way is north. We'll answer some of those concerns now. What we really need to do to answer all of them is map project the image into some coordinate system. Map projection is something we'll worry about later in the course. For now we'll work around it. You can exit qview from the File menu in the main box.

Two useful routines in ISIS that can give you some information are lev1pt and avg_sd. lev1pt must be used with a level one cube. It works for one pixel in the image, which by default is the center pixel. Try running it…

Type

lev1pt from="m1000782.lev1.cub"

Lots of information comes up such as the latitude and longitude (of the center pixel) and the pixel size in both the horizontal and vertical (sample and line) directions (in km/pixel) along with the value of the aspect ratio. Take note of the pixel sizes, that information will be useful later. Also displayed is the value of north azimuth, which gives the direction of north in the image. To find north you need to draw an imaginary line from the center to the right edge of the image and then move it clockwise around the center from there by the number of degrees in the north azimuth figure. The line should now point north.

Lets get rid of the annoying aspect ratio; the program that will do this for us is called mocaspect. There's a small complication here so lets use the tutor mode for this one. Enter tutor mode in the usual way (type

't mocaspect'

Now try the other information gathering program 'avg_sd' on this new aspect corrected cube. Type

avg_sd from="m1000782.asp.cub"

To communicate with other image processing programs some other file format other than ISIS cubes must be used. Firstly we'll output the final image we have (m1000782.asp.cub) as simple raw 8-bit data. At the moment the pixel values are the I/F numbers, which are stored as floating point values. We just want the picture for the moment so 8-bit numbers from 0 to 255 will do fine. When converting to 8-bit we must decide what the zero represents in I/F units and what the 255 represents in I/F units. All I/F values will be scaled, any above 255 will be set equal to 255 and any below 0 will be set equal to zero. So we could potentially lose a lot of information here and we certainly don't want to do that after all this effort. So we must choose a low and a high I/F value so that all the other I/F values fall between them. So its lucky indeed that we happen to know what the min and max I/F values are since you just noted them down (you did note them down….didn't you?). We're going to use the program isis2raw to get out of cube format and just get a raw picture.

Type

't isis2raw'

We're briefly going to startup IDL and convert that raw file to something a little more normal like a tiff file. Quit out of the tae environment by typing exit at the tae prompt. Type IDL at the unix prompt to start that program. You should be in the IDL environment with an IDL prompt. This will only take a few commands, type:

openr,1,'m1000782.raw'
a=assoc(1,bytarr(XXXX,YYYY))
a=a(0)
close,1
write_tiff,'m1000782.tif',a

In place of XXXX type the number of samples and in place of YYYY type the number of lines. After this you can exit IDL just by typing exit at the prompt. You now have a processed image in tiff format. This part of the lab is over, you need to transfer the tiff file to a PC or a MAC that has photoshop running for the next part of the lab. If you need help with the transfer then let me know. Make sure you have the north azimuth and sample and line resolution figures noted for the next section of the lab.

To recap the whole process:

• Run moclev0 to convert the PDS file to an ISIS cube.
• Use qview's doctor to remove any garbage lines.
• Run moclev1 to radiometrically calibrate the image.
• Run mocaspect to remove the aspect ratio.
• Run avg_sd to find the min/max pixel values.
• Run isis2raw to convert the final cube to a raw file.
• Run a small piece of IDL code to rewrite the raw image as a tiff.

Something for you to do yourself now…. In the directory there was a second .imq file m1000783.imq. This is the context image for the narrow angle image we just processed i.e. it's a small low-resolution image taken at the same time to show the regional setting. Using the exact same technique, process this image and create a tiff file. The method will be exactly the same although this time there won't be any garbage lines you need to cut out with the qview doctor. Note down the north azimuth and line and sample resolution number for the next section.

Photoshop

Full on scientific image processing programs usually lack a fast simple interface and performing simple tasks like contrast stretching, zooming and panning often take more effort than is really warranted. Adobe photoshop is a good compromise between convenience and usefulness. We're going to use it to touch up the tiff images you processed in the last section.

Firstly load in the two tiff images you've made and take a moment to look them over zoom in and out a few times. In the case of the narrow angle image (m1000782) try and stretch the contrast in the image --> adjust --> levels menu. You'll probably notice that the best stretch is a compromise between the light and dark areas.

Here's a trick to try, in the filters -> others menu select the high pass filter and a value of around 100 (play around with this and experiment with very high and very low values). The best value still preserves the image details but the large scale (low frequency) bright/dark variations have been removed, you can contrast stretch the entire image now. The filter also serves to emphasize small-scale detail, which is usually a good thing unless the image is very noisy.

Experiment around with other filters and techniques; you can always undo changes from the edit menu. For this exercise you need to produce some products.

Narrow angle: Select a subset of the m1000782 image that has some interesting features. Touch it up as best you can with different filters and stretches. Put a scale bar on the image that is some appropriate number of meters long. You have the resolution of the image so you can figure out how many pixels long a 500-meter scale bar should be. Also use the line tool with the arrowhead option to indicate on the image which way is south.

Wide Angle: As before touch up the image as best you think and include a scale bar and a south-pointing arrow. Also draw a rough outline where the narrow angle image lies in this scene.

I want you to make these as if they were figures in an article you're about to publish, include a caption on why the features you choose are worthy of note and what they mean. Include any information about the image you think is relevant in the caption.

Searching for your own images

Staying with the case of MOC images, several websites are now offering various search mechanisms. Each one has its advantages and disadvantages a quick summary of the main three follow.

http://ida.wr.usgs.gov

This site is run by the USGS. Its main benefit is that if you know the image name then you can find and download it very fast as everything is arranged by orbit number. The database runs only up until the end of the mapping mission its not as useful as it could have been.

http://www.msss.com

This is the website of the company that actually built the MOC camera. They are also interested in the science it returns and have constructed a graphical database for all the publicly released images. Its difficult to find a specific image in here but this site is great for browsing regions of the planet if you don't know exactly what you're looking for.

http://www-pdsimage.jpl.nasa.gov/PDS/public/Atlas/Atlas.html

This site has recently been finished and is provided by the PDS. It's an extension of a search engine from the previous major Mars mission (Viking). This allows searching graphically by zooming up on a map of the planet. However, the thing that makes this site so very useful is the ability to search the MOC dataset using forms. You can specify any number of image parameters (including latitude and longitude ranges) and get a list of all the images, which match your search. Preview thumbnails are available and each image can be viewed online and downloaded in any number of formats.

Use any one of these sites (I recommend the third) to search in some region you find interesting, such as the Valles Marineris canyon system or the Tharsis volcanoes. Download the PDS version of the image and process it as before. As before make up a respectable figure and caption saying what is of interest in the picture. Include things you deem relevant in the image such as scale bars, north arrows, context views, direction to the sun etc...

Advanced techniques: Photoclinometry

This section is aimed at being an introduction to the concept of photoclinometry (deriving topography from surface shading).If you've read tutorial one by now you should know how local slopes as well as albedo can affect the brightness of a particular patch of ground (if you haven't read it yet now would be a good time). Photoclinometry is difficult to do in any quantitative way but this exercise will illustrate some of the qualitative concepts associated with it.

In the lab3 directory start IDL by typing

'idl'

You need to add the path of the idl routines for this lab by typing:

!path=!path+':/home/ge151/idl_routines/'

This program uses the I/F values in the image m0201989.asp.cub (stored in the lab3 directory). Some of the program lines are reproduced below.

name = 'm0201989.asp.cub' ;Name of file to read in
res = 7.879 ;Size of a pixel in meters
sunaz = 52.42 ;Solar Azimuth
inc_ang = 88.11 ;Inclination Angle

These lines set up the variables specific to this image.

b = b-shadow

b is the variable that stores the I/F values alone the sun's line of sight. Here we remove some estimate of the I/F value, which is due solely to the atmosphere (called shadow brightness because this is the only way shadows can be illuminated)

c = acos(b * !PI * cos(inc_ang) / albedo)
c = c - inc_ang
c = res * tan(c)
for i=1,n_elements(c)-1 do c(i) = c(i-1)+c(i)

These 4 lines are the guts of the program. The first calculates the incidence angle for each pixel based on its I/F value. The second then removes the incidence angle that a flat surface would have, leaving just the local slopes for each pixel. The third line converts the slope of each pixel to the change in height of that pixel. The fourth and final line adds up all these height changes to find the actual height of each pixel.

This isn't a programming class so you don't need to learn this, only notice what important. Notice for example that we have to assume some value for the albedo and that the atmosphere can have an important role due to its ability to scatter light.

You can start the program by typing 'pclin' at the IDL prompt. The standard case will come up i.e. the program selects an interesting part of the image, guesses an albedo and assumes no atmospheric effects. You can force it to use some particular albedo by calling it like 'pclin, albedo=0.25' for an albedo of 25%. You can make some correction for the atmosphere by assuming some amount of I/F is due to scattering e.g. 'pclin, shadow=0.02' will assume that 0.02 needs to be subtracted from each I/F value as a first order atmospheric correction. You can use a third keyword 'pclin, sl=3000' to select the starting line for the section of image you want to look at. You can use all three of these keywords in any combination you like.

The solar azimuth means that only tracks that cut across the image at the right angle can be used for photoclinometry so the profiles will be quite short.

Experiment with different albedos, letting the program guess first and then trying higher and lower values. Try this for a number of areas and try to figure out what goes wrong when the albedo is far away from the real answer. Explain it as best you can in words. Do the same for the shadow illumination parameter.

On the website is a section of image along with the MOLA altimetry track (it's also the lab3 directory as m0201989.jpg). Use the UNIX command 'xv m0201989.jpg' to view it. We'll learn more about MOLA in the next section. This figure should guide you in what the true topography is (but note that the MOLA track and photoclinometry profiles are not parallel).

Altimetry and Topographic Analysis

During this lab you will work through:

• Looking at single-track data from the MOLA instrument.
• Deducing some simple things about prominent topographic features on Mars.
• Making some gridded products of MOLA data and what they are useful for.

This lab is going to make extensive use of IDL (interactive data language). You used it briefly in lab 1 for converting raw files to tiff. I've tried to write all the software in advance so that this doesn't turn into a class on computer programming. Note however that this software is not bulletproof and if its given funky data you'll get funky answers. Also the MOLA data is much more difficult to extract from its archived PDS form into useful numbers than is the MOC data you used in the first part but I'll try to make this as painless as possible.

For the second part of this lab you'll probably need to use a computer thats better than average, otherwise IDL will run out of memory. I recommend tharsis. Also most machines in GPS have undergone an operating system upgrade withing the last week or so meaning that telnet no longer works.

Introduction and Background

The Mars Orbiter Laser Altimeter (MOLA) operated on the Mars Global Surveyor (MGS) spacecraft. It measured topography by pulsing a laser at the planet's surface and recording the length of time it took for that pulse to be reflected back. It is a nadir-pointing instrument and so can only see what's directly below it. The laser fires 10 shots a second, the spacecraft has a groundspeed of 3 km s^-1 so each shot is separated from the next by about 300 meters. The laser illuminates a patch on the surface that's about 120 meters in size. As the spacecraft orbits the planet a track of elevation measurements is built up of the terrain that MGS has passed over.

In the first part of the lab we'll look at one of these tracks in particular and the features that they cross over. In the second part of the lab we will make some common grided products that you've probably seen displayed like shaded relief maps.

I. The single MOLA trace

Startup IDL in the directory '/home/ge151/lab3/' (no need for your own directory in this lab) by typing 'idl' at the UNIX prompt. Remember to modify the !path variable as stated above. We have all the MOLA data on disks in the gps computer system. The data structure of the MOLA files is long and complicated so to spare you the details and let you get straight to the data we're going to use a small program that will automatically go and get any MOLA track you ask for. This program will return 4 arrays, one of the longitude of each shot, one for the latitude, one for the topographic height and one for the time each shot occurred. All you need give it is the number of the MOLA track and the names of the arrays in which you want to return the above information. In IDL type the command:

get_mola_converted_track,11368,lat,lon,z,tm

This will retrieve track number 11368 and store the information described above in the arrays lat, lon, z, tm (latitude, longitude, height and time of shot respectively).

A quick aside on MOLA track numbers, the number is actually the orbit number + 10000. So this is the track obtained on orbit 1368. If you want to find the MOLA track that goes along with a specific MOC image for example you need only know the orbit number that the image was taken in.

In IDL you can find out information about arrays by typing 'help,ARRAY-NAME' if you leave off the array name and just type 'help' you'll get information on all the variables IDL has stored in memory. Type 'help', you can see each array has 68298 elements. You can display the values in any array with the print command but I don't recommend that unless you want to see 68298 numbers scroll past on your screen. Try things like 'print,lat(0:9)' which will show you the first 10 elements in lat (IDL indexes arrays starting from zero). Try 'print,max(lon)' or 'print,min(tm)' to see maximum and minimum values of some of the arrays. Note that the time is in awkward units of seconds after JD2000 (Julian date of 00:00 1/1/2000). Computers have no problem with numbers like this but humans generally do (at least I do). We can convert the times into something less obnoxious looking by subtracting the lowest value from each shot time. This will convert the shot times into seconds from the start of the orbit instead of seconds from some arbitrary time over a year ago. So for the first part of the lab print out all the shot times (all 68298 values) and manually subtract the time of the first shot from each one. Just kidding, IDL will do this for us with a single command, type:

tm = tm - min(tm)

This redefines the array tm to be seconds from the start of the orbit i.e. from 0 to just under 2 hours (7200 seconds). To see the elevation try the command:

plot,tm,z,ps=3

This will plot the topography recorded against the time during that orbit. Remember each second represents about 3 Km along the ground (the ps=3 part just tells IDL to use a dot for each point in the plot). This is the easiest way to view the MOLA topography, if you try to plot it against latitude for example you'll get two lines on top of each other one for each side of the planet. Take a look at the height profile, it includes both the highest and lowest parts of Mars. What is the total relief on the planet? You can use the plot command to plot longitude vs. latitude. This will show you the ground track of the spacecraft. Do this and then look at the geologic maps in 056 S. Mudd. Make a list of interesting looking features and their latitudes that should appear in the plot. We'll do all the height plots vs. time but you can relate time to latitude or longitude by plotting them against each other e.g. plot,tm,lat,ps=3 gives the latitude at each time in the track. To look at specific parts of the track in more detail we'll specify an x range to display over with the xr keyword. Since the x-axis is displaying time in all these cases our x range will represent a slice of time during the orbit. Lets try from 3700 to 4000 seconds.

plot,tm,z,ps=3,xr=[3700,4000]

You can use this x range with any of the plots you tried above also. You should see something pretty damn big in this part of the topography plot. Remember the y-axis is elevation in meters. Find the latitude and longitude by plotting time vs. lat and time vs. lon with the same x range. Check the geologic map to see what the feature is.

What are the average slopes on the flanks of this feature? It looks pretty high but the answer is you could probably hike it (assuming you had a spacesuit and enough time). Estimate the volume of the feature, assume its a cone, measure the base in seconds and convert to kilometers using the 1 s = 3 Km rule. What is the mass of the feature? If the strength of rock is about 100 mega-pascals what's the maximum height of a conical mountain on Mars?

You can go through the list of features that you saw on the map earlier and find what times in the orbit correspond their longitude and latitude. Find the part of the topography profile that crosses the feature and describe how the topography behaves there.

View the entire profile again by plotting without an xr keyword. Find the times around the large depression early in the orbit and zoom up on it. This is the Hellas impact structure. Estimate the volume ejected during the impact. It's fair to say that we're firmly in the gravity regime here. Use the scaling relationship derived in the crater brainteaser to figure out the energy of the impactor. If the impact was at 10 Kms ^-1 what was the size of the asteroid that hit (assume a density of 2700 kg m^-3)? How does that compare with the size of the largest asteroids visible in the asteroid belt today?

Now pick a random number between 10000 and 15000. Use that as the track number and retrieve the data in the same way. Look for some interesting feature and identify it in the geological maps. Describe how the topography behaves over it. Is that what you expect for a feature like that? Note try for a fairly large feature otherwise you might not find it on the geologic map.

II. Making gridded products

If you have enough data (and we do) you can create maps of the surface either showing plain topography, slopes or shaded relief (both of which can be derived from the topography). Contoured maps can also be created but they're rapidly going out of fashion even for earth geologists (info courtesy of Kerry Sieh). So we're not going to go into contouring.

First we need data and lots of it. There is an automatic tool which will return ALL the MOLA data within a specified lat/lon region, that's a lot of data and since you don't have 5 hours to watch it all download and I don't want to get hate emails from Goddard space flight center we'll skip this step and use data I downloaded earlier. We're going to try the region surrounding the Olympus Mons Shield Volcano.

It might be nice to restart IDL at this point so that all the old variables you were using are wiped out. Use 'exit' to get out and 'idl' to get back in while still in the lab3 directory. Set the !path variable.

Try the following command to get the data:

restore,'ge151_lab2_data.sav'

Take a look at the arrays, this time were looking at over a million points! You can never have enough data. We have latitude, longitude, track number and height. To make a map we have to somehow convert the lat/lon values into x and y values in some projection. We're looking at a locality pretty close to the equator try 'print,mean(lat),mean(lon)' to see our average position. So a reasonable good projection to use is sinusoidal. It's commonly used in Mars research when you're not too close to the poles. In this projection latitude lines are straight and parallel and longitude lines are straight and converge to a point at the poles. Use the command:

latlon2xy_sinusoidal_sph,lat,lon,133.0,3393400.0,x=x,y=y

This makes two new arrays called x and y, and fills them with the planar coordinates in meters needed. You now have point coordinates in standard Cartesian space, i.e. x,y & z. The sampling is pretty irregular though. Try 'plot,x,y,ps=3,/isotropic' to see the distribution of the points. The /isotropic keyword makes the x & y axis have the same scale (a good thing since all coordinates are now in meters). We're going to make a regular surface out of this jumble. We'll choose a pixel size for our output image of 1 Km since that will ensure that most pixels have some topography points within them. This routine fills in the gaps with interpolation so we don't want to leave too many gaps. Type the command:

make_surface,z,x,y,res=1000.0,xra=xra,yra=yra,output=pic

You'll get 3 new variables out of this. 'pic' is the new regular grid. Its a 2-d array of topography values where each cell in the array is 1 km² on the surface of mars. So it can be displayed like an image. 'xra' and 'yra' contain the x and y coordinates of the samples and lines in the image respectively. This process can take a few minutes so be patient. You can view the result by typing 'tvscl,pic', make sure the window is big enough for the image. IDL has many color tables try 'loadct,33' its similar to the color scheme that the MOLA team uses in their publications. To make the shaded relief and slope maps use the two following commands.

shad = shade_topo(pic,xra,yra,elev=30.0)
angl = topo_to_slopes(pic,xra,yra)

These can be viewed using the same commands e.g. 'tvscl,shad', it's probably best to set the color table back to normal for this, use 'loadct,0' to do this. You may notice that there are nasty edge effects where the data ends in the images. You can use the following few lines of code to clean up the images. It just slices sharp edges into the images and sets all the surrounding pixels to zero kind of like using a cookie cutter to throw away the ragged edges.

lim=[11.0,127.1,26.0,139.9]

slice, pic,xra,yra,lim=lim

slice,shad,xra,yra,lim=lim

slice,angl,xra,yra,lim=lim

Note in the command to create the shaded relief map there is a parameter 'elev' this is the suns elevation above the horizon. Experiment with low and high sun angles and see what the difference is. Low sun angles tend to bring out subtle details of the topography better.

Display the slope map once again 'tvscl,angl'. IDL allows you to display only pixels with certain values if you wish e.g. 'tvscl,angl<20' will only display slopes less than 20 degrees, everything higher will look saturated. 'tvscl,angl>5<20' will only display slopes with values between 5 and 20 degrees etc..

Use the command 'tvscl,angl>15<15.01' this makes everything with slopes more than 15.01 degrees to be bright white and everything with slopes less than 15 degrees to be black. This shows the bounding scarp very well. It also shows that it has broken in several places in what look like large landslides. Compare this to the image of Olympus Mons that comes with this lab (oly.jpg). Can you see the same details in the visible image?

Try the command:

plot,histogram(angl),ps=10

This will show what the distribution of slopes is, ignore the big spike at zero degrees, that's caused by all the empty space surrounding the image.

Lastly output the shaded relief map as a tiff file using the write_tiff procedure (see lab1 for details).

write_tiff,'tiffname.tif',bytscl(shad)

The bytscl command just scales the pixel values from 0-255 and so makes them suitable for use in a tif file. Import into photoshop, contrast stretch and attach a scale bar as before, remember we set one pixel equal to 1 Km. Make an estimate of the caldera's size, how many Rosebowl sized stadiums can you fit in there?

That's it, good luck.

Troy

PS. Getting Help

PS. Getting Help

We will have a lab help session, in addition you can come and see me whenever you like at 167 South Mudd. If you email in advance we can fix a time or you could just try and catch me there.

Photoshop has an extensive help section that can be accessed from the program.

IDL help manager can be started by typing '?' at the IDL prompt.

ISIS can provide help for any routine or any parameter when in tutor mode. In addition the USGS web-site at http://wwwflag.wr.usgs.gov/isis-bin/isis.cgi, documents all ISIS programs very well.

GPS homepage - Home - General Info - Schedule - Assignments - Labs - Field Trip - Reading - Projects - References - Feedback

Back to top

These pages are currently maintained by O. Aharonson