|
Ge 151a - Lab 2Spring 2003Altimetry and Topographic AnalysisDuring this lab you will work through:
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 lab 1 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. You need to use secure shell if you want to remotely log in e.g. from the UNIX prompt use ssh -X ge151@tharsis.gps.caltech.edu Secure shell will automatically set the display settings to work (note X is uppercase) but forgets to source the .cshrc file. You have to do this manually withsource ~/.cshrc once you're logged in.Introduction and Background(top)The Mars Orbiter Laser Altimeter (MOLA) operates on the Mars Global Surveyor (MGS) spacecraft. It measures topography by pulsing a laser at the planet's surface and recording the length of time it takes 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(top)Startup IDL in the directory '/home/ge151/lab2/' (no need for your own directory in this lab) by typing 'idl -32' at the UNIX prompt. 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 that have been left in the Mars room. 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(top)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 -32' to get back in. 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 km2 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) 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. Make a sketch from this data of the bounding scarp and the locations of the breaks. Then compare your sketch to the image of Olympus Mons that comes with this lab. 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. Shane PS. Getting Help(top)We will have a lab help session,in addition you can come and see me whenever you like at 156 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. GPS homepage - Home - General Info - Schedule - Assignments - Reading - Lectures - Tutorials - Labs - Brainteasers - Projects - References |