Geology Research Methodologies

Introduction

The aim of this independent study is to determine whether or not a correlation exists between shifts in the obliquity of Mars and the deposition of rhythmic sedimentary beds on the Martian surface. Previous studies (e.g. Lewis et al, 2008) have documented the deposition of interior layered deposits (ILD’s) on the central mounts of Martian craters that seem to have a rhythmic pattern in their thickness. These studies conclude that obliquity seems to be the driving factor in the periodicity. They give a general time scale for the deposition based on predicted obliquity values, but they do not seem to directly correlate the time it takes to deposit these layers with the change in thickness observed, which is the goal of this study.

Several programs are used in this study, including ArcMap, Davinci, Excel, and AnalySeries. Each Section of this document will correspond to the use of one of these programs, and will show the methodologies used to obtain my results. The aim of this document is to fully show all methodologies, with examples when needed, so these techniques can be reproduced.

I present my methodologies in a numbered list; anything stated after a number is a process I’ve done, and anything said after a letter is either additional information about this process or a more detailed explanation of what needs to be done.

Part 1: ArcMap

Elevation data and spatial data are needed to conduct this study. These data are stored in image files, of which the actual data need to be extracted. ArcMap is the preferred program to extract these data.

1. Download the data files.
   1. These files can be found at http://hirise.lpl.arizona.edu/dtm/ . For this analysis to be done, two files need to be downloaded: A Digital Terrain Map (DTM), which shows elevation data, and one of the two stereo pairs used to create the DTM, which is a grayscale image and is needed to correlate the elevation data to the area we wish to study.
   2. For this study, the DTM “Layers in Becquerel Crater” is used, along with the left observation stereo image.
2. Open a new document in ArcMap.
3. Click Add Data, and add the visual image first.
   1. The visual image starts with “PSP_” and is a .JP2 file.
4. Next, add the DTM.
   1. The DTM starts with “DTEEC_” and is a .IMG file.
5. Georeference the images.
   1. The images pertain to the same area, but do not display that way on ArcMap. This means they need to be georeferenced. This can be done by using the georeferencing toolbar in ArcMap.
   2. Use the “Add control points” button on the Georeferencing toolbar to select a point on the visual image, and match it up to the same point on the DTM. It is recommended to use a unique region of the image, as we want the images to line up as closely as possible, or else the data will be inaccurate. Use enough control points to ensure the two images match as closely as possible.
   3. The area that is to be clipped should be extremely well georeferenced-- at the expense of other areas of the map, if necessary. If another area needs to be clipped, this step can always be repeated.
6. Once the images are properly georeferenced, clip them to the area of study. It is recommended to overestimate clipping size, since the resulting image will be cropped later anyway. (note that image files that are too big may not be able to be displayed in Davinci/MS Paint. As such, balance must be found between an extent that is large enough to work with but small enough to not cause errors associated with a large file size)
   1. ArcToolbox -> Data Management Tools -> Raster -> Raster Processing -> Clip
   2. For the input, put the visual image (.JP2) first. Then set the minimum and maximum X-Y values using the coordinates given on the lower right corner of the screen. Remember to overestimate the clipping area.
   3. For the output raster dataset, be sure to save the file in a place you know where it is!
   4. Next, navigate to Clip again, but this time select the DTM (.IMG) as the input. For the output extent, choose the clipped image, which should have been added as a layer. The minimum and maximum x and y values should change to match the area of the clipped visual image. Do not change them. Save the output in a known directory and clip the DTM.
7. Extract the map from the clipped images to image files that can be opened with Davinci.
   1. Right click on the image and go to Data -> Export Data. Keep the default cell size as they respond to the resolution of the image. Save the image as a TIFF in a known directory. Do this for both the visual image and the DTM.

Using ArcMap, we have now taken the elevation and visual data and selected only the areas we want to study. The next step is to define a transect and extract the data along the transect. This is done using Davinci.

Part 2: Davinci

Davinci is a command-line program, and it is assumed that the user has some prior experience with the program. It is also assumed that Davinci’s default image viewer is MS Paint, and that a .dvrc file has already been set up. Any code to be inputted into Davinci will be written in italics . The easiest way to store information in Davinci is by the use of a structure file, and my research was done in this fashion. Therefore, my instructions will use this format to save the data as well. If at any time you wish to save data, go to step 9.

1. Open Davinci and define a new structure file, naming it whatever you wish. I will call it “a”.
   1. a={}
2. Create another structure file to load the full image into. It is easier to use multiple structures inside of the main structure to work with specific files. I will call the first one “full” since it uses the full image extent.
   1. a.full={}
3. Load the image files into the structure, giving them distinguishable names. Throughout this paper, I will use “vis” for the visual image and “topo” for the DTM.
   1. a.full.vis=load(“C:/.../’vis.tif”)
      1. Note that the information in parenthesis represents the pathway to where the tif image is stored. This part of the code will change based on where the images are stored and what the clipped images were named; the above line is an example for simplicity. (My images were stored on the M drive with a different name).
   2. a.full.topo=load(“C:/.../topo.tif”)
      1. For the write() and load() command, the information will always be different based on the name and location of the file.
      2. This code is telling Davinci to load the image file made from the DTM into the “full” struct within the “a” struct as “topo”. In subsequent steps, when I work with other images in other structure files (i.e. not named “full”), I will assume a new structure file has been created in the same manner as step 2a, but with a different name.
4. Begin by displaying the visual image and rotating it to the point where the layers look like they are vertical (going from top to bottom).
   1. display(sstretch(a.full.vis))
      1. Sstretch is needed to convert the image into a format that is able to be displayed. In subsequent steps when I mention displaying an image, using a sstretch will be assumed.
   2. b=rotate(a.full.vis,##)
      1. In this command, “##” refers to the number of degrees the image is to be rotated so the layers appear vertical. Some experimentation may be needed, and the file can be overwritten if desired. (e.g. b=rotate(b,##)) )
      2. To test to see if the layers appear vertical, Paint can be used to show the transect. The easiest way to do this is to turn the gridlines on in Paint (under view), and if the horizontal grid lines (representing your transect) go equally through all layers, the image is rotated sufficiently. (Because of the nature of the layers, it is not possible to rotate an image so each layer is entirely vertical).
      3. For this image, a rotation of 37.5 degrees was used.
   3. Save the files in a different structure file, but still in the same master structure (in my case, “a”). I will call mine a.rotated.
      1. a.rotated.vis=b
         1. In this case, “b” refers to the final rotation of the original image.
   4. Rotate the topographic image in exactly the same manner as the visual image, and save it in the same structure (e.g. as a.rotated.topo)
5. Once the images are rotated, decide on a transect.
   1. The transect should be an area that covers every layer wished to study, so it has to be on a pixel number that passes through every layer. The transect can only be one pixel in width; techniques such as averaging pixels in the y-direction are not advised because of the non-linear shape of the layers.
   2. a.transect.data=a.rotated.topo[x1:x2,y:y,]
      1. In the above code, x1 is the start point in the x-direction, and x2 is the end point. The ‘y’ value must be the same so the image is now only one pixel in the y-direction.
      2. To take brightness values, do the same thing for a.rotated.vis . If you are taking both values, give them different names (e.g. a.transect.topo and a.transect.vis)
   3. To see a graph of elevation vs. distance, use the following code: pplot(a.transect)
6. Create an x-axis.
   1. Since the spatial resolution of the DTM is 1m/pixel, an xaxis is not needed for pplot. However, it is still useful to create an x-axis to import into Excel. The create() command will be used to do this.
      1. a.transect.xaxis=create(####,1,1,start=0,step=1,format=float)
         1. In the above code, ‘####’ refers to the amount of pixels in the x-direction, which can be found by entering a.transect.data (or whatever you called your transect image) into Davinci. It will then say “####x1x1 array of float…”). #### is the number you want in your xaxis.
         2. For brightness data, the step is 0.25.
7. Prepare the topographic data to be exported into Excel.
   1. The exported data needs to be in bsq format. To do this, simply overwrite the existing file as a BSQ.
      1. a.transect.data=bsq(a.transect.data)
8. Export the data as a CSV file to be read in Excel. Use the comma as the separator to make it easier to read in Excel.
   1. write(a.transect,”C:/.../transect.csv”,type=csv,separator=”,”,header=1)
      1. The header will make the first steps easier in Excel.
9. Save the Davinci data as an HDF file, in case you need to return to the data.
   1. write(a,”C:/.../a.hdf”,type=hdf)
   2. The data can now be loaded. If you wish to overwrite the file, simply add ,force=1 between “type=hdf” and “)”

Part 3: Excel

The use of Excel is needed to further prepare the data for analysis.

1. Open Excel and load the CSV file created from Davinci. Also open a new Excel workbook.
   1. You should see two rows of data: The first row is the header, which shows whether the data is elevation/brightness or xaxis data. The second row is the data itself.
2. Find the column where the xaxis data begins. Copy the values (second row) from the beginning of the xaxis data to the end of the spreadsheet (it should be 0 to whatever number represents the last pixel in the image)
3. Once these values are copied, go to the blank Excel workbook, click on cell A1, and paste Transpose.
   1. The values should now be pasted in the first column of the workbook.
4. Now, copy the values of data from the beginning of the CSV file (again, second row) to the last value in the “data” series (the last value that has a “data[##]” value in the first row)
   1. If you exported something other than “data”, it may have a different name (e.g. brightness, topo, vis, or whatever you named the dataset in Davinci)
5. Paste these values into the second column of the Excel workbook that now contains the xaxis data.
   1. Click on cell B1 and click Transpose Paste.
   2. The spreadsheet should now have 2 columns with values that end at the same row. (for me it was 1600).
6. If you want to apply a boxcar filter to the data, go to the column next to the data (most likely column C) and the row corresponding to the first value outside the “margin”. Then use Excel to average the value in the column to the left as well as the values up to the first datapoint, and the equal number of values down. Do this all the way down, then make the “margins” the same as their original value.
   1. For this study, a 17-band filter was used. This means the first and last 8 datapoints were in the “margin” and were the same as the original values, while every other datapoint was an average of the original datapoint plus the 8 values above and below it. The 8th datapoint was the datapoint that began the filter.
7. Save this new workbook with a name other than the name you gave the csv file.
   1. It is advisable to save as both a tab-delimited text file and an Excel workbook. The tab file will be needed for AnalySeries.

Part 4: AnalySeries

The use of AnalySeries is required to perform the Blackman-Tukey analysis. It should be noted that this program is only available for Macintosh computers.

1. Open AnalySeries and load the tab-delimited text file created from Excel.
2. Name the columns so you know what they are. The first column should be distance (meters), and the second column should be brightness/elevation. The third column may be the boxcar filter of the second column..
3. From the drop-down menu, select BTukey. Choose the correct columns to do the analysis on (the data--with or without boxcar filter--and distance/xaxis)
4. Save the data and open it.
   1. Go to “Save Worksheet As” to do this.

The data can now be viewed. Graphs can be made either using the “draw” tool in AnalySeries or (preferably) in an Excel sheet containing the BTukey data.

Part 5: Obliquity Data (Excel)

The goal of this study is to match the ratio of peaks of the Blackman-Tukey power spectrum for the layers (the previous steps) with the Blackman-Tukey power spectrum for the obliquity of Mars, which will be discussed in the next two parts. The first part involves preparing the data for analysis in Excel.

1. Download the Mars orbital data.
   1. The data can be downloaded from http://vo.imcce.fr/insola/earth/online/mars/DATA/index.html . Right click the file you want to download and click “save link as”. The default file type is fine.
   2. For this study, the dataset “301003BIN_A_P001_N.ASC” was used because it was found by a previous study (Andrews-Hanna and Lewis, 2011.) to be consistent with modeling of Martian paleoclimates.
2. Start Excel, and open the downloaded file.
3. The text import wizard should appear. The type that best describes the data is Fixed Width. If it is not already selected, click that. Then click Next.
4. Set the column widths so that each column is separated. This should be already done, but if it is not, click on the spaces in between columns until each of the four columns are separate. Click next.
5. The column data format should be General. Click Finish.
6. To match the x-axis to years, multiply the x-axis (column A) by -1000.
   1. This can be accomplished a number of ways, but the easiest may be to type =-1000*(A1) in cell E1. Then, paste throughout the column (double click the green square in the lower right hand corner of the cell) and copy the entire E column and paste values into column A. Then, delete column E.
7. The data should be imported into four separate columns. Save the data as tab-delimited text (and Excel, if you wish).

Part 6: Obliquity Data (AnalySeries)

The final part of the methods is the performing of a Blackman-Tukey analysis on the obliquity data and comparing results to the layered sediment data.

1. Start Analyseries and open the Mars data.
2. Name the columns. The first is time (years), the second is eccentricity, the third is obliquity, and the fourth is Ls.
3. Run a BTukey analysis on the obliquity data with time in the x-axis.
4. Save the data and open it, either in AnalySeries (draw) or in another program such as Excel.
   1. If opening in Excel, it will recognize a tab-delimited text file. Just click “finish” on the pop-up that opens and it will correctly sort the data into columns.

The goal is to analyze the graphs of both the brightness and obliquity data in the hopes of correlating them. The value on the x-axis of the BTukey graph is equal to the reciprocal of the value (in meters or years, respectively). For instance, if there is a peak at roughly 8.5e-6, this corresponds to (1)/(8.5e-6)= every 117,000 years. There should be two peaks in the brightness spectral plot and ~7 peaks in obliquity. This study looked at the ratio of the second peak in the brightness plot divided by the first, and compared that to the ratios of the peaks in obliquity (e.g. second divided by first, third divided by second, etc.). We found a ratio of 1.75 between the two peaks, which corresponds almost perfectly with the ratio between the first two peaks in obliquity. For more information on the results and conclusions, see the poster!