Lab - Processing and Interpreting Spectral Data Sets

Ge151: Spring 2011

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, or your own computer if you are on the GPS network. 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

Ask your TA for a list of Linux machines that you can use, if they haven't already given them to you. Try not to use a sun machine, as there's no guarantee these programs will work as intended.

Source the .cshrc file once logged in with:

source ~/.cshrc

I - Multispectral Image Analysis - THEMIS

During this section of the lab you will work through:

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/lab5 directory, the filename is I02207005_geo_radiance.envi.ers. This is a multispectral image taken by the THEMIS instrument, onboard the Mars Odyssey spacecraft.

Log on to the ge151 account with the password given out in class. First, you need to create a folder to save any output files you might create. Move into the lab5 directory and use the 'mkdir' command to create your own subdirectory (mkdir directoryname).

Now change to your directory by typing:

cd directoryname

1. Introduction and Background

The Mars Odyssey Thermal Emission Imaging System (THEMIS) consists of a 5-band visible imager with 18 m/pixel spatial resolution, covering a wavelength range of 0.42 to 0.86 microns, and a 10-band infrared imager (100 m/pixel) covering 9 unique spectral bands between 6.8 and 14.9 microns [Christensen et al., 2004]. (Bands 1 and 2 both cover the same wavelength range). In this lab we will use data from the infrared subsystem, over the region Nili Patera, a caldera of Syrtis Major. In Lab 4 you learned that minerals exhibit diagnostic spectral features at infrared wavelengths. Here will use that information to interpret thermal infrared data from Mars.

THEMIS data are available from the PDS at http://themis.asu.edu, in the form of calibrated radiance and brightness temperature. Note that THEMIS data products at the PDS have not been geometrically projected or band registered. Calibrated radiance images contain both surface and atmospheric contributions. The radiance measured by THEMIS, L(&lambda), is given by:

Where &epsilon(&lambda) is the surface emissivity, B[T_surf, &lambda] is blackbody radiance as a function of surface temperature (T_surf), &mu is the cosine of the emission angle, &tau(&lambda,p) is the normal opacity profile as a function of pressure, and T(p) is the atmospheric temperature profile. The integral is taken through the atmosphere from &tau=0 (opacity at the spacecraft) to &tau =&tau₀ (opacity at the surface). The first term describes absorption of surface radiance (&epsilon(&lambda)B[T_surf, &lambda]) by the atmosphere, while the second term describes radiation emitted by the atmosphere [Bandfield et al., 2004]. For this lab, the second term has already been accounted for in the radiance images, using the method described in Bandfield et al., [2004].

2. Daytime radiance and introduction to the area of interest

Start ENVI by typing:

envi_rt

in a terminal window. Open the THEMIS radiance image for this lab by clicking File &rArr Open External File &rArr IP Software &rArr ER Mapper. In the dialog window navigate the filepath to /home/ge151/lab5/I02207005_geo_radiance.envi.ers. This image has already been band registered for you using ISIS (which you encountered in the Imaging and Topography Lab). An 'Available Bands' window should pop up, showing 10 available bands. Before you start, first crop the image to make it more manageable. In the main tool bar, click Basic Tools &rArr Resize Data (Spatial/Spectral). In the 'Resize Data Input File', select the 'I02207005_geo_radiance.envi' file, and click the 'Spatial Subset' button. Change the Samples to 41:300, NS to 260, Lines to 2701:4500, and NL to 1800, and click OK, and then OK again. You will be asked to provide a filename. Do that, and then click OK. The cropped data should appear in your Available Bands window. Select Band 5 of the cropped image, click 'Gray Scale', and click 'Load Band'. Three new windows should open: 'Scroll', 'Image' and 'Zoom'.

How many surface units would you define based on what you see in the radiance image? Briefly describe them. (there's no right answer, just use morphology and temperature differences to decide)

3. Deriving surface temperature and emissivity from radiance

Variations in surface radiance (&epsilon(&lambda)B[T_surf, &lambda]) from one region to the next can be due to either surface temperature, emissivity variations (likely composition), or both (eq. 1). To extract compositional information from a radiance cube, you must first perform a temperature-emissivity separation. For each radiance spectrum, there are n data points (radiance value at each channel n), but n+1 unknowns (emissivity at n, and T_surf). Thus either an emissivity at one channel must be assumed, or a temperature must be assumed to solve the equation [e.g., Kahle and Alley, 1992]. In this lab, we will use the Emissivity Normalization technique to separate temperature and emissivity [Gillespie, 1985; Realmuto,1990]. For each pixel, a brightness temperature is calculated for each band using an emissivity value of 1.0. The maximum brightness temperature found for that spectrum is assumed to be the surface temperature (T_surf), and then Tsurf is used to solve for emissivity for the other channels. In doing this, the emissivity values will be normalized to 1.0 and to the band with the maximum brightness temperature [Christensen, 1998]. Note that surface emissivity is likely not equal to 1 across a full THEMIS band (~1 μm wide). Some researchers use a lower emissivity value (usually between 0.95 or 0.99), especially if a priori knowledge about the surface materials is available [Kahle and Alley, 1992]. In any case, it is important to remember that the derived emissivity values are relative to that of the channel set to 1.0; absolute values may be lower.

Convert 10-band radiance to surface temperature and 10-band emissivity. To do the calculation, you first have to give ENVI the wavelength centers of each THEMIS band, so that it can generate the required Planck functions. On the main tool bar, click File &rArr Edit ENVI Header. Select your cropped radiance file and click OK. In the new window that pops up, Click Edit Attributes &rArr Wavelength. Enter the following wavelength values for each band: Band 1 (6.78), 2(6.78), 3(7.93), 4(8.56), 5(9.35), 6(10.21), 7(11.04), 8(11.79), 9(12.57), 10 (14.88). For the wavelength units, select 'Micrometers' from the pull down menu, then click OK and OK. You will have to reload your Band 5 radiance image from the Available Bands list.

Now you're ready to convert to emissivity. On the main tool bar, click Spectral &rArr PreProcessing &rArr Calibration Utilities &rArr Calculate Emissivity &rArr Emissivity Normalization. Select your cropped radiance file that has the wavelength information, and click Spectral Subset. Select only Bands 3-10, click OK and OK. In the 'Calculate Emittance Parameters' window, change the Data Scale Factor to 10000. This is done to convert radiance units from Wcm^-2μm^-1sr^-1 to Wm^-2μm^-1sr^-1. Note that the assumed emissivity value defaults to 0.96 in ENVI. For this lab, and to be consistent with published analyses of THEMIS data, we will set that value to 1.0. Choose to output the temperature image, and give filenames for both the emissivity and temperature files, and click OK.

The files you have created will appear in the Available Bands window. You can open them by selecting the band you want, and clicking Display 1 &rArr New Display and then click Load Band (or you can just click 'Load Band' first, and ENVI will load the data into the image space you already have open).

Examine the surface temperature image (in the image window, click Tools &rArr Cursor Location/Value to see pixel values under your cursor). What is the temperature range observed? What are possible reasons for the observed temperature variations?

Now, examine Band 5 emissivity. How does it look different than the radiance image? Why?

In the Available Bands list, select emissivity bands 8, 7, and 5 and plot them as red, green, and blue using the 'Open RGB' button. How does the image appear? Can you identify differing spectral units on the surface? (Note: Spectral variations are often due to compositional variations, but not always, thus until full analysis is completed, we will call them spectral units).

4. Visualizing spectral variations

Emissivity images are highly correlated from channel to channel. By this I mean that areas that have low emissivity values in Band 4 will likely have low values in Band 5, Band 6, etc. This is due to the broad similarities in spectra of most rocks between 8-12 microns at low spectral resolution, and is why the 3-band color image from 2c appears somewhat pale. To visualize spectral differences that are present, the image must be transformed, or stretched.

Perform a Principal Components Analysis (PCA) to try to determine the number of end-member components of the scene. PCA is a common technique used in image analysis to reduce the dimensionality of the dataset. On the main tool bar, click Transform &rArr Principal Components &rArr Forward PC Rotation &rArr Compute New Statistics and Rotate. Select your emissivity file, click on 'Mask Options', and click 'Mask NaNs [All Bands]', then click OK. Provide filenames when prompted, then click OK.

A plot of eigenvalues for each band will pop-up. Based on the eigenvalues, how many principal components do you think are significant? (no right answers here)

Examine each of the PC bands. How many of the PCs identified represent real surface components? How many are noise? (or here)

Perform a Decorrelation Stretch (DCS) [Gillespie et al., 1986] on the emissivity bands. Decorrelation stretching is similar to Principal Component Analysis, except that after the PC transformation occurs, the data is stretched orthogonal to the PC axes, and then rotated back to the original coordinates and displayed as RGB. In PC analysis, you may choose the number of PC component images to be output, and are limited only by the number of available bands. In decorrelation stretching, however, because the data is displayed as RGB, the user must pre-select 3 bands.

On the main tool bar, select Transform &rArr Decorrelation Stretch. From the pop-up window, select Available Bands and select OK. Select EMISSIVITY Bands 8, 7, and 5 from the available bands list and click OK. Give a filename and click OK. Then load the RGB DCS image into a new display.

Now, based on examination of the DCS image, how many spectrally distinct units are present? (no right answers again for these questions...)

Note how the dune field appears in the DCS image. Which unit would you guess is the source of the dunes?

Try a decorrelation stretch using Bands 9, 6, and 4. Now that you can see additional spectral information (by looking at two separate stretches), would you answer the previous question the same way from just looking at the 8-7-5 image? Why do the two DCS images look slightly different?

The 8-7-5 and 9-6-4 band combinations were pre-selected for you simply because they provide broad spectral coverage. Choose a third band combination and try the DCS. Describe any new spatial variability in spectral properties that you find.

5. Quantify differences between spectral units.

Remember that the DCS is a type of stretch, therefore it is possible that differences observed in the stretch may be not that great in actuality. Examining the spectra from each 'spectral unit' allows you to quantify the differences observed in the DCS images. Here you will extract spectra from units you identified from the DCS images and plot them. The emissivity image from this region has been atmospherically corrected for you, in /home/ge151/lab5/I02207005_surf_emiss_cropped. Use File &rArr Open Image File and navigate the file path to this data cube. Choose any band from the 'Available Bands' list (except band 3 - that is the band most commonly set to unity in Step 2) and open it in a new display window. Link the atmospherically corrected emissivity cube with your favorite DCS RGB image. In the Image window of your DCS image click Tools &rArr Link &rArr Link Displays. Then select the displays you want to be linked and click 'Yes' for those displays, and then OK. When the images are linked, you can click around in one image, and the cursor will move to that location for all linked images.

Next, in the Image window containing your atmospherically corrected emissivity band, click Tools &rArr Profiles &rArr Z Profile (Spectrum). A plot window should open up. Now, click around in your DCS image (the small skinny window) and watch the profile window. The plot will change every time you move the crosshairs. In the Plot Window, click Options &rArr Collect Spectra. Now click around again. The spectra from each pixel you click on will remain in the plot box. Next, instead of using the mouse, try using direction arrows on the keyboard to move the crosshairs around one pixel at a time. To clear the plot, turn off the 'Collect Spectra' option and move to a new location. Using your new ENVI plotting skills, do the following three tasks:

From one of the spectral units you have identified, extract a few adjacent spectra. Plot the spectra and note the variation between them. Would you feel comfortable with interpreting surface emissivity variations from one pixel to the next?

Now extract ~100 spectra from one unit, and note the improvement in signal to noise. Save these as an ascii file by clicking (in the plot window) File &rArr Save Plot As &rArr Ascii, and entering a file output name. Do this for other spectral units you had identified from the DCS/PCA exercise. Open the spectra in a spreadsheet application. Average the spectra from each unit and plot the averages. What are the major emissivity differences between each unit?

Open the file /home/ge151/lab5/laboratory_spectra.xls with your spreadsheet application. This file contains laboratory emission spectra of various rock types. Compare the spectra with your spectral unit averages. What do you think are the broad compositional differences between these units?

II - Mineralogical Composition - TES

1. Introduction and Background

The Thermal Emission Spectrometer (TES), one of the instruments on the Mars Global Surveyor (MGS) spacecraft, has lower spatial resolution than THEMIS, but many more spectral bands [Christensen et al., 1992]. This allows for extraction of much more mineralogic detail. When used together, the TES and THEMIS data are a powerful compositional mapping tool. TES spectra from two distinct units identified in the THEMIS image have been provided for you. Continuing in ENVI, the locations from which the TES spectra were taken may be viewed by opening an ROI file. In any of your image windows, click Tools &rArr Region of Interest &rArr ROI Tool. In the ROI Tool dialog box click File &rArr Restore ROIs, and select /home/ge151/lab5/TES_spectra_locations.roi.

The TES spectra from each unti are from MGS mapping orbit 1927. These spectra contain contributions from the surface and atmosphere. Here you will perform the atmospheric correction using the method of Smith et al., [2000]. Using a library of atmospheric spectral end-members derived by Bandfield et al. [2000] and mineral spectra, you will perform a linear least squares fit to the TES spectra.

cd /home/ge151/lab5
matlab

in a terminal window to start. First, we need to read in the data files. In matlab, type

specinit

to read in the data for units 1 and 2, as well as the library of laboratory mineral spectra, and the library of atmospheric endmembers. Each of the items read in is a structure, with several variables each holding information about the data. Type the name of one of the variables (lib, atm, unit1, unit2) to see what's inside. Take a look at the unit1 spectrum with the command:

plotwb2(squeeze(unit1.xaxis),squeeze(unit1.data(1,:,:)))

This is a plot of emissivity (the ratio of the observed energy radiated by a material to that of a blackbody at the same temperature) vs. wavenumber (cm^-1). You should see that the spectrum is dominated by an absorption feature at ~700 wavenumbers. Unfortunately, this is due to atmospheric CO₂ absorption, and only gets in the way of determining the surface mineralogy. There are also some areas at the ends of the spectrum that aren't considered for mineralogical analysis. In order to proceed, we'll need to crop out these areas of the spectrum with the commands:

unit1.data=cat(3,unit1.data(:,:,9:35),unit1.data(:,:,65:110));
unit2.data=cat(3,unit2.data(:,:,9:35),unit2.data(:,:,65:110));

The deconvolution code we will use will search for a linear combination of the library and atmospheric spectra that best fits the measured spectrum. Take a look at the mineral endmembers we'll be considering here, by typing:

lib.label

Here you should see a bunch of mineral names (Quartz, Microcline, etc., followed by various letters and numbers. The letters and numbers represent the catalog from which that mineral sample was acquired (e.g. WAR=Wards's).

2. Fitting a spectrum, the hard way

u1=decon(unit1.data,lib,[200 1305],[],atm);
u2=decon(unit2.data,lib,[200 1305],[],atm);

A bunch of text will scroll by showing the results of the deconvolution. We'll end up with the variables u1 and u2, which contain various results their respective structures. Before we actually look at the results, let's try to get a more intuitive feel for the fitting process. At the matlab prompt, type:

specfit(u1,lib,atm);

This should bring up a plot window with a bunch of slider bars. The point of this routine is to give you a sense for how the automatic deconvolution works. The slider bars add different amounts of a given mineral to your fitted spectrum, while the plot shows the total sum of all the components you select, along with the current rms error between your fit and the measured spectrum. Both the raw (measured) and atmosphere-removed spectra are plotted. If your atmosphere component is set to zero, the rms is calculated for the atmosphere removed spectra. Once you start adding atmosphere, the rms switches over to the measured spectrum. Try and find the best fit for the atmospherically corrected spectrum. What is the lowest rms error you can find? (If you can beat the least squares value of .0019, you get an A!!) What are the best fit values you find for each mineral group?

Now, looking at the individual mineral group spectra with the slider bars, you should notice that each spectrum is dominated by one main absorption feature. For many of the mineral groups, this occurs at 10-11 microns, but for others, it occurs at longer wavelengths. Try to guess what transitions might be (generally) responsible for these features. (Hint: these energies are too low for electronic transitions, so we're primarily looking at vibrational modes).

3. Fitting a spectrum, the easy way

Okay, now let's look at the best fit model of the spectrum. At the matlab prompt, type:

specview2(u1,lib,atm);

You should see another window pop up, with a plot and some checkboxes for various the mineral groups. These checkboxes let you turn on and off the various mineral components (as determined by the decon function). Which are the most important minerals to the best fit spectrum? How does this compare with the fit you found by hand? What are the fewest number of components which provide a reasonably good fit to the spectrum? (no right answer for this one...)

Okay, now lets compare the unit1 and unit2 (which appeared different to THEMIS) spectra, to see how they appear to TES. We'll look at the atmosphere removed spectra, since we're interested in mineralogical differences. Type:

plotwb2(u1.xaxis,u1.rematm(:,1));
hold on
plotwb2(u2.xaxis,u2.rematm(:,1),'r'); hold off

Unit1 will be blue, while unit2 will be red. You should see that there are some clear differences. Now, lets open up the fitted solutions. Type:

[u1.groupcoeff(:,1) u2.groupcoeff(:,1)]

This will bring up the coefficients for each of mineral groups, for unit1 (left column) and unit2 (right column). What are the major differences in mineral abundances? You can open up u1.groupnames if you need to see the names of the groups again. Now let's look at the standard deviations of the derived coefiecients, to determine the significance of some of the results.Type:

[u1.groupcoeff(:,1) u1.groupstd(:,1) u2.groupcoeff(:,1) u2.groupstd(:,1)]

Now, looking at the standard error magnitudes, are there any mineral detections that are not statistically significant? Is the difference in say, carbonate between the two units significant? What about the difference in high-silica phases?

III - Interpreting the geology

Using the compositional, morphologic, and stratigraphic relationships present in the various datasets you have analyzed, what can you say about the geologic history of this area? Obviously, you can't know every single event and process that has occurred, but try to say what you can. Additional context images from MOC and THEMIS are located in the /home/ge151/lab5 directory, and can be viewed by typing display imagename in the unix terminal. They should also be linked on the class webpage for this lab.

Now that you're all finished, we can tell you that a similar analysis of this region was presented in Christensen et al. [2005]. Check it out if you are interested!

That's it, good luck!

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

Back to top

These pages are currently maintained by Oded Aharonson