Monday, March 21, 2011

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, 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.

Tuesday, March 15, 2011


I'm in the midst of finals right now.  I plan on making a post or two about this quarter's geology related GIS work.    In general, I mostly learned a few tricks when exporting ArcMap layouts to Illustrator.

One of the very last things I did was extracting measurements of error between two sets of contacts of the same lithological unit.

Stay tuned!

Thursday, March 10, 2011

Dip Direction and Dip stdout patch for LidarViewer

If you have used LidarViewer or are about to for taking measurements of bed attitude, you might be interested in a little patch I have.  It outputs the dip direction and dip to the Terminal when taking attitude measurements in LidarViewer (using the Strike and Dip primitive).

This is important because the dip direction and dip numbers can be obscured by point cloud data and so one has to waste time manipulating the world to find the numbers.

Download patch (right click and Save As in your LidarViewer2.5 directory)

Apply by saving into the LidarViewer-2.5/ directory and using:
patch -p1 -i diprdirdip_stdout.patch

Then recompile using make.

I do not warrant this patch. Use at your own risk and remember that I really dislike C++, so I may have broken some rules and made computer scientists cry :-)

Tuesday, March 8, 2011

The Magic Wand

Illustrator's Magic Wand tool is precisely what you need when it is necessary to select only the fault lines from the contact lines, as exported from Arc.  Set the tool's stroke width and color constraints, click on a fault, and viola, all the faults are selected!  Also, it might be useful to remember that the eyedropper tool can also pick up stroke width, color, and style... not just colors.


Tuesday, March 1, 2011

Crusta (from UCD KeckCAVES) and Arc GRID Raster Projection Issues

It may be a long shot, but in case someone is out there wondering why their elevation data is offset from their drape imagery while using crusta, I have a little hint.

GDAL, to my knowledge, has of yet to learn how to find the projection of Arc GRID data.  So you need to specifically point to the prj.adf file.  Then it is smart enough to find the files that contain the image data and process acccording to the projection found in the prj.adf file that is part of the GRID data structure.  Otherwise, you'll be mislead into thinking Construo found projection information when it begins processing data... it defaults to WGS 84. A keen eye, and awareness of your data's projection, would be the tip off that something may be wrong. So, using construo:

construo -dem ./topo.globeFile ./lidargrid/prj.adf

The Crusta wiki instructs users to point to the GRID's container directory/folder, but GDAL won't find the GRID's projection and so Construo may create a globeFile with the wrong projection.

Construo uses its default projection... and transforms/registers the elevation data incorrectly.

Elevation data offset solved by using the prj.adf file within the GRID data structure.

This example is a good reason to understand what a projection is.   They can cause a lot of trouble.