Makes me smile when...

Makes me smile when TAs mention that my maps were the best in class.  Remote sensing was a fun course.  I learned a lot and would love to use those skills along with GIS (perhaps coupled with Python) to extract valuable information.  Here I go into the job market; I've graduated!

More useful posts to come!

Graduating

I'm graduating very soon. Next week, I have finals.  This is why I have not posted anything recently.  I have lots of neat tricks to write about once I have the time.  Technology loves geology.

Topographic Profile From Basemap Contours, not a DEM

There are times when USGS DEMs do not match topography on basemaps.  I have a basemap that does not appear to be derived from any presently available USGS products.  DEMs are often interpolated from 7.5' quads, by using air photo stereopairs, SRTM, and so on.  The topography on the basemap I have is much more detailed than any DEM or USGS 7.5' topo I have come across.  So it was probably drafted by someone else.  Thus I set out to extract a topographic profile line using the basemap and Arc.  This is a pretty simple process, and it relies on geostatistics to interpolate elevation in between captured points.  The following instructions assumes you are familiar with principals of GIS.

1. Plot your cross section line
2. Capture elevation points on each contour that crosses your cross section line
1. It is best to turn on edge-snapping if you're using Arc Desktop 10
2. Make sure you capture the start and end points of your cross section line; you may have to interpolate elevation yourself.  In fact, there may times along the cross section line where it would make sense to estimate elevation, such as at bottoms of streams or tops of ridges.
3. Krig your elevation points or use an interpolation method that you're comfortable with 
4. Use Spatial Analyst or 3D Analyst (or open source equivalents) to generate a profile line along the resultant raster
5. The more elevation points you capture, the better.
6. Export the profile data to tabular format and use as desired (see blog post on Illustrator and profile lines)
7. Rejoice in avoiding drafting a profile line by hand!

I ended up with very reasonable results in a much shorter time than by hand.  Now my topographic profile matches the basemap and it is is much more accomodating of the geology I mapped!  The lesson here is to compare USGS DEMs to your basemap before trusting DEM derived products.  Attention to detail is always key.

Coloring Cross Sections: Adobe Illustrator Paint Groups

A while ago I posted a procedure on how to create a topographic cross section in Adobe Illutrator without vertical exaggeration ("Geologic Cross Sections and Adobe Illustrator").  The procedure used tabular elevation and distance data and Adobe Illustrator.  In that post, I did not explain what to do next.  So here's some more guidance.  Screen captures or a video may come later.

1. Export your completed geologic map WITH the cross section line(s) from Arc or your chosen GIS to Illustrator-friendly format at 300 DPI.  
2. Embed or link the AI geologic map file into your cross section Illustrator file.
3. Rotate and align your newly imported map and a cross section line to the bottom of the profile's distance (X) axis (or at the bottom of the lowest elevation on your profile).  This is where you find out if your profile is at map scale!
4. Create a new layer, call it "Cross Section Edges"
5. Select the layer and project lines of appropriate weight (contact vs. fault) from the intersection of the X axis and the map's contacts or faults up through the profile line. 
6. Begin drawing your beds using strike and dip information, remembering to use apparent dips when necessary.
7. Be sure that all your lines end at another path.  For example, a contact line might end at the boundary of your profile chart or at a fault.  Do not leave gaps open between any paths.
8. When you're done drawing all your bedding and faults, select all of your cross section art and turn it into a Live Paint Group
9. Now Live Paint all the faces and edges!  See Adobe's help for more information.
10. If you need to use patterns (say, from the FGDC Geologic Map Symbol website) and a unit color, I suggest coloring first, then copy the entire paint group, paste into place, and paint using a pattern.  Hopefully your patterns have transparent backgrounds.
11. I also suggest creating new layers for every component of your cross section such as legend, title, and other information.
12. Use the Lock layer tool and smart guides to your benefit!

This Quarter: Remote Sensing

This quarter at UC Davis, I'm in a remote sensing class.  The lab instructs on use of ENVI software and principals of remote sensing.  Yes, ENVI is powerful.  No, it does not have a modern interface.  I keep thinking I'm using some throw-back MOTIF graphical user interface.  ENVI has mile-long menus.  This is the only software I've used in recent history that has no obvious organization, and seemingly, very similar submenu categories.  But alas, it appears to be one of the standard software packages for remote sensing, and so I will learn it well.  Too bad ITT charges $200 for an academic license and their website is a nightmare from yesteryear.  ESRI, at least, has an awesome program to encourage use of its software: students in GIS courses taught by in-the-know Professors typically get free one year licenses.  What's up, ITT?

It is unfortunate that  out of thirty students, I appear to be the only geologist taking the course.  I am so very confused about that.  In fact, it is also rare that geology students take GIS courses.  So my job prospects look pretty bright, if this is the rule and not the exception.

Strike and Dip: GIS and LIDAR Data

For the last two quarters, I have had to extract strike and dips directly from LIDAR point cloud data using LidarViewer.  It was awful. Today I had some time to think about how one might do this in a GIS.  I also implemented it.

Essentially, LidarViewer fits a plane to select points in the cloud.  The problem is, it is difficult to repeat selections and extract location information for the fitted strike and dip (plane attitude).  Today, I devised a method (model) in ArcMap that uses an elevation grid, masks, and zonal statistics to find the mean dip direction and mean dip of bedding.  Of course, the dip is nothing more than the slope of topography.  Keeping that in mind, only features that are large enough relative to the scale of the LIDAR data can be measured.

I first created a polygon mask of small areas that surround previous strike and dip spatial data.  Each polygon has a unique ID.  The shapes were defined by referring to slope, aspect, and orthophoto rasters. I outlined areas with generally constant slope and aspect*.

I then derived aspects and the slopes of 0.5 m interpolated LIDAR DEM of Rainbow Basin, CA.  Using the mask, I extracted all the raster values from both derived rasters.  The polygon mask can then be used to generate a table of zonal raster statistics for both aspect and slope within each polygon.  This provides the following statistics: Minimum, Maximum, Range, Mean, Sum, Area, and Standard Deviation.  A few of these can be thrown out, such as area and sum.

So I generated two new tables without that extraneous information.  I then joined the tables together and joined the resultant table to the polygon mask using the unique IDs.  Finally, I created a point feature class by taking the centroids of each polygon within the mask.

The resultant point feature class has three important attributes: mean dip, mean dip direction (the aspect is the geographic azimuth of maximum slope), and the standard deviation.  The standard deviation is a useful measure of dispersion and may indicate an adjustment of polygons may be needed.  One can do this by referring back to the slope, aspect, and orthophoto rasters.

Here's my result.  Red are data collected using LidarViewer.  Black are strike and dips created using my model.

Polygons with blue cast are GIS sampled regions.  Notice two symbols coincide within a few degrees within strike.  A few are off by ~10 degrees in strike.  The dips are close enough, given the uncertainty in either method of extracting attitudes!  LIDAR data from OpenTopography.org, Blackwater Region, NSF EarthScope.

Here's my likely overcomplicated model:

If you want an export of this model, leave a comment below, I might just release it to the public!

I think it is fair to say, that given large enough scale (in the cartographic sense) elevation data, strike and dips can be extracted all at once very easily using a mask with multiple polygons.  The trick is comparing slope, aspect maps, and orthophotos to find areas that may have measurable attitudes.  But nothing will beat field measurements!

* I must confess that only today did I realize the value of aspect for remote mapping of geology!

Comparing Lithologic Contacts

One important thing that I have learned over the last two quarters is that geologists will rarely map contacts exactly the same.  To make things worse, if you can map contacts using remote sensing methods, the contacts may also be quite different.  Sometimes contacts are obvious since they follow specific topographic features.  Other times, they're difficult to track with only float giving some indication of presence.  Because of the variability in mapping, some geologists may want to easily compare their contacts with previous efforts, other geologists, or remotely mapped contacts.  Here is one way to quantify differences between contacts.

This method assumes that you have already input two sets of contacts in the field area into a GIS as polylines.  This method should work on any GIS, including R, but some steps are specific to Arc.  In this example, Contact A and Contact B are a set of contacts for the same related lithological area.  They are being compared using the shortest distance between contacts.

The euclidian distance raster is the background.  Central to the screen shot are the two sets of contacts being compared.  The black is Contact A, the blue is Contact B.  The point features, used for sampling the euclidian raster, should make it clear which contacts are being compared.  The points have an interval of 20 m.

1. Select only Contact A
2. Generate a Euclidian Distance raster for Contact A.  Be thoughtful about your cell size (I have used 0.5 m) and extent of the new raster.
3. If using Arc, enter Editor mode and create a point layer.  Ensure that point layer and the contact layer are in editing mode.
4. Select only Contact B, merge the contact lines (perhaps in a new feature class).
5. Generate an interval, your choice, of points that coincide with the merged Contact B.  In ArcMap 9.3., you can use the "Divide" tool under the Editor menu (... maybe).  In ArcMap 10,you can use the "Construct Points" tool under the Editor menu. Save and stop editing.
6. In ArcMap, use Spatial Analyst's Sample tool.  Use the Euclidian distance raster as your input raster and your line coinciding points as your point features ("Input Location Raster or Point Features").  Define your output table's location.  The default resampling technique is fine.  See above figure.
7. Press OK and you will now have a table with the shortest distance from the points along Contact B to the lines of Contact A
8. Now import the distance data into your favorite statistics software and see how your data compares (or just use Arc's built in statistics tools).

Output table with sampled distances between Contact A and Contact B.