Astro 3310Fall 2015
LAB #8
------
Please copy this document to the REPORT sub-directory from the expanded LAB8_Data_Package_FA15.tar.gz. Then, edit it to write your answers in all the "______" areas. When finished, create a tar.gz archive of the REPORT directory and all of its contents,then scp thefile to datafarm.astro.cornell.edu and place it in:
/data/Courses/A3310/FA15/”your netid”/LAB8/
Remember that you will only get credit for the files that you put in the REPORT sub-directory and copy to datafarm. Please make sure that youkeep a Matlab workbook with all of the commands you used to answer the questions in the lab. Feel free to comment and organize your workbook so that it will be easy for us to follow your algorithms when we execute the doe. If you generate any functions for the lab, ensure that they are also in the REPORT sub-directory and properly called from the workbook file. For you convenience, there is a template for the workbook file already in the REPORT sub-directory.
YOUR NAME: ______YourNetID: ______
For this lab, we are going to combine ArcGIS and Matlab to study the orientation and thickness of periodic bedding in the Arabia Terra region of Mars. The lab is based off of a paper entitled “Quasi-Periodic Bedding in the Sedimentary Rock Record of Mars” by Kevin Lewis. The paper was published in the journal Science in 2008 and is located in the BACKGROUND sub-directory of the lab data archive. Before getting started on the lab, please read through the paper.
Problem 1: Measuring the thickness of individual beds.
- The goal of this problem is to recreate Figure 1 in Lewis et al. 2008 (reproduced below). To get started, open a new project in ArcMap and load the image PSP_002733_1880_RED_A_01_ORTHO_crop.cub from the DATA directory of the LAB8 archive. Using ArcCatalog, create a new shapefile (store it in the REPORT sub-directory) with a “PolyLine” geometry with the geographic projection “Mars 2000”. Then, add a line to your new shapefile that defines the transect shown in Figure 1 of Lewis et al. 2008. Be sure to save your edits and click the “Stop Editing” option when you are done. Before closing ArcMap, export a Figure similar to the bottom panel of Lewis et al. 2008 Figure 1 that showing the location of your transect.
Next, open matlab and move to the DATA sub-directory of LAB8 and add the SUBROUTINES directory to your matlab path by typing “addpath(genpath(‘../.’))” from within the DATA sub-direcotry. In the SUBROUTINES directory of the lab you will find a function called “extract_dtm_values.m”. Open the function up in the matlab editor by typing “edit extract_dtm_values” and look at. This function will grab the topographic heights and image brightness values along the lines you made in ArcMap. To run function, type:
“shp = extract_dtm_values(shpfile,dtmfile,imgfile)”, where:
shpfile = the path+filename of your ArcMapshapefile
Dtmfile = DTEEC_002878_1880_002733_1880_U01_crop.cub
imgfile = PSP_002733_1880_RED_A_01_ORTHO_crop.cub
The function will return a structure variable “shp” that has, in addition to their standard content of an ArcMapshapefile structure, fields named “samples”, ”lines”, ”lat”, ”lon”, “x”, ”y”, ”z”, ”b”, and “track”. These fields correspond to the position of the line, sample, latitude, longitude, and xyz values of the resolution elements along each polyline in your shapefile. The “b” field is the IOF (ratio of received to incidence flux) pixel values from the image file along your transect and “track” is the integrated along-track distance along the transect.
Use the vectors stored in the “track” and “z” fields to generate a figure similar to the upper panel of Figure 1 in Lewis et al. 2008. Do you see the same step-stair pattern that is pointed out in the paper?
Place your two-panel figure here
(Replace the below reproduction of Lewis et al. 2008 Figure 1).
- What is the average thicknesses of the beds along your transect? Measure a few using either ArcMap or Matlab and report the average.
Hint: Bed thickness should be measured using the elevation (“z”) differences.
Answer: Average bed thickness = ______
Problem 2: Measuring the periodicity of stratified layers
- The goal of this problem will be to generate plots that are similar to Figure 2 of Lewis et al. 2008 (reproduced below). To get started, open a new project in ArcMap and load the image PSP_001546_2015_RED_A_01_ORTHO_crop.cub from the DATA directory of the LAB8 archive. Using ArcCatalog, create a new shapefile (store it in the REPORT sub-directory) with “PolyLine” geometry with the geographic projection “Mars 2000”. Then, add a series of lines to your new shapefile that traces several (~10) of the beds shown in Figure 3 of Lewis et al. 2008. After tracing the beds, add a final line that runs perpendicular to the bedding and defines a line that runs perpendicular to the beds (similar to what you did in Figure 1). Be sure to save your edits and click the “Stop Editing” option when you are done. Before closing ArcMap, export a Figure similar to the bottom panel of Lewis et al. 2008 Figure 3 that shows the location of your lines (don’t worry about trying to recreate the three-dimensional view in Panel A).
Insert your Figure here.
- Open matlab and move to the DATA sub-directory of LAB8(remember to add the SUBROUTINES directory to your matlab path by typing “addpath(genpath(‘../.’))” from within the DATA sub-directory). Use “extract_dtm_values.m”to grab the topographic heights and image brightness values along the lines you made in ArcMap. To run the function, type:
“shp = extract_dtm_values(shpfile,dtmfile,imgfile)”, where:
shpfile = the path+filename of your ArcMapshapefile
Dtmfile = DTEEC_002878_1880_002733_1880_U01_crop.cub
imgfile = PSP_002733_1880_RED_A_01_ORTHO_crop.cub
For the lines that trace bedding planes, calculate the dip of the bedding. To calculate the dip, fit the xyz coordinates for each line to a plane defined by the equation:
Z = A .* X + B.*y + C
The dip (slope of the plane) is then equal to:
dip = 180/pi*atan( (A^2+B^2)/ sqrt(A^2+B^2) )
(If you have then time and are interested, you can also try and determine the error in your dip calculations)
You can you any of the methods discussed in class or used for previous labs to find the best-fit coefficients A, B, and C.
Report your dip values below.
Answer: Measured dips of traced beds = ______
Average dip of traced beds = ______
Are your dip values similar? If there are disparate values, do you notice anything about the beds with dissimilar dips?
- Now lets looks at the line that is running perpendicular to bedding. To correct for bedding orientation, divide the z values for this transect by the cosine of the average dip calculated above (should be between 1-3 degrees). This is the stratigraphic height. Next, sort the image IOF values such that they are ordered by increasing stratigraphic height:
zs = z / cosd(average_dip);
[zs_sort, sort_ind] = sort(zs);
bs_sort = b(sort_ind);
[zs_sort, unique_ind] = unique(zs_sort);
bs_sort = bs_sort(unique_ind);
The last line above removes any repeat Z-values and ensures that the elevation vector is monotonically increasing (required for the next step).
Our goal is to take the Fourier Transform of the detrended image brightness (see Lewis 2008 Figure 2A). In order to do this, we want to create a uniformly sampled vector of the detrended brightness as a function of stratigraphic height. In order to do this, we will interpolate the vectors “zs_sort” and “bs_sort”:
zs_sort2 = min(zs_sort):0.1:max(zs_sort);
bs_sort2 = interp1(zs_sort,bs_sort,zs_sort2,’pchip’);
Make a plot of zs_sort2 vs. bs_sort2, similar to Figure 2A of Lewis et al. 2008 and paste it below.
Insert figure here.
Do you see any periodic trends in the data? If show, what is the main period in m?
- Using the “FFT” function in matlab in order to generate the Fast Fourier Transform of “bs_sort2”.Since we interpolated “bs_sort2” to have uniform 0.1 m spacing, our frequency sampling is 10 m^(-1):
Fs = 1/0.1;
L = length(bs_sort2);
NFFT = 2^(nextpow2(L));
num = numel(bs_sort2);
Y = fft( (bs_sort2’ – mean(bs_sort2)).*hann(num),NFFT)/L;
Y = 2*abs(Y(1:NFFT/2+1));
F = Fs/2*linspace(0,1,NFFT/2+1);
where Y is the single-sided amplitude spectrum of “bs_sort2” and F is the spatial frequency in units of inverse meters.
Be sure you understand what the above lines are doing and why (if not, come ask).
Make a plot of F vs. Y (similar to Figure 2B of Lewis et al. 2008) and insert it below.
Do you see any peaks in the FFT? If so, at what spatial frequency do they occur?
Note that we used a hanning window to taper the function prior to taking its FFT. Different tapering techniques may help clean up the amplitude spectrum. If you have the time, try to taper the function differently.
Insert figure here.