Monday, June 13, 2011

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!

Friday, June 3, 2011


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.

Wednesday, May 11, 2011

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.

Sunday, April 24, 2011

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!

Wednesday, April 6, 2011

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.

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.

Saturday, February 26, 2011

Geological Society of America References Style for BibTex

I hunted around long enough for a GSA BibTeX style, and I found one.  Direct link to .bst.  It appears to work well enough for my class papers.  Thanks to the Journal of Cave and Karst Studies!

I'll be keeping my copy safe.

Wednesday, February 23, 2011

Geologic Cross Sections and Adobe Illustrator

Last academic quarter, I spent time trying to understand how to easily remove vertical exaggeration (VE) from tabulated elevation and distance data in Adobe Illustrator.

Adobe Illustrator has a chart tool.  It isn't very robust or powerful.  In fact, it can be a pain to work with.

However, there is a way to manipulate its initial size such that (hopefully) VE = 0.

Perhaps you have found a better way? Here's my final solution.
  1. Import the data into a temporary scatter plot and take note of its automatic maximum vertical and horizontal tick values.  These values seem to be consistent no matter how you size the chart area.  Delete it.
  2. Take those maximum values and convert to the scale that you need to use (e.g., 1:24000) 
  3. The final two numbers should be in units of your artboard
  4. Create a new scatter plot chart but taking care to snap its size to your determined horizontal and vertical map distances using snap guides
  5. You should have VE at or near zero; verify this!
  6. Configure your chart: remove point markers, set line weights (tricky...), etc.
  7. Select your chart and ungroup it.  Yes, you do want to destroy the ability to change graph data because you'll really want to modify your new profile graph!
  8. See "Coloring Cross Sections: Adobe Illustrator Paint Groups" for what to do next.
This method hinges on the assumption that the drag box is exactly where your data will show up, with the axes coinciding or just barely exterior to it.  It may not be the best assumption but it seems to be the case.  Be sure that all the ticks are equidistant and there are no distortions on the axes.

In an ideal world, we would scale elevations and horizontal distances to artboard distances.  Then, import and plot the points so they automatically coincide with the units of the artboard WITHOUT using the chart tool.  This would remove the guesswork... is it possible?

Sunday, February 13, 2011

Measuring Lengths of Many Lines in Illustrator CS5

When balancing a cross section, you may need to know the lengths of beds between your pin lines.

Adobe Illustrator does not make this process obvious, and the built-in measure tool can measure only one line segment at a time.

To find the lengths of a group of lines (grouped or selected), turn on the Document Info dialog using the Window menu.

Select the options drop-down menu (the icon looks like justified text) on the upper right-hand corner of the Document Info dialog.  Select Objects.

Now you'll have the nitty-gritty details about your selected path(s)!

This could also be useful for measuring distances on non-georeferenced raster maps (with an absolute scale) by using a series of line paths. 

Thursday, February 10, 2011

Balanced Cross Sections Made Easy with Illustrator CS5's Smart Guides

Should you ever find yourself needing to professionally draft a balanced cross section, you may want to try Adobe Illustrator.  Under the Smart Guide preferences, you can customize the angles that lines snap to (so long as you have snapped the beginning node to another node, I think).

Sometimes Smart Guides has a difficult time snapping to angles if you start drawing a line without attaching it to an anchor. One thing that you can do is define the angle using the line dialog.  The other may be to draw a boundary line at the margin of the cross section.  Or, they just won't snap... and I have no idea why.  The Mac version of Illustrator seems rather inconsistent with its functions, at times.

But once you get past that initial point of setting up a line, Smart Guides is usually content snapping between lines that represent dip domains (angles that bisect between dip domains).  Remember, by default, smart guide activity is symbolized with a green line.

Once you've balanced a cross section, use the measure tool to make sure bed thicknesses and angles come out as expected.  It's pretty slick stuff.
30 and 60 degree angles automatically snapped between bisectors of dip domains.

So much for the protractor.

Monday, February 7, 2011

OSXStereonet 1.0.0 Released

N. Cardozo and R. Allmendinger has released a brand new version of Stereonet called OSXStereonet.  With a new name, this release is version 1.0.0.  It makes use of modern OS X interface and display frameworks.  You can now export to EPS and PDF.  In fact, you can drag and drop your stereonet plots onto the desktop.  Nice.

The software makes use of the Inspector interface, doing away with the Preferences menu.  In the Inspector you can set how the stereonet displays, tweak symbology,  and adjust Kamb contour parameters.

In addition, data now plots as you enter them within a brand new data editor interface.  Results for calculations are detailed in another panel.  You can load multiple datasets and turn them on and off at any time.

The data file format has changed slightly though older data format from the previous generation of Stereonet should import.  However, I had issues with importing trend and plunge with the plunge and quadrant format (trends didn't show up) I used for Stereonet 6.3.3.

This looks like a slick release!

Saturday, February 5, 2011

The Strain Ellipse and Adobe Illustrator

When you have stretches and orientations for a population of features, you might have a strain ellipse.  Why plot it by hand when you can let Adobe Illustrator do it for you?

Illustrator CS5's Line tool allows you to specify an angle and a length (stretch, in this case).  The line tool uses a slightly different system of orientation than in geology.  360/0 is east, 90 is north, 180 is west, and 270 is south.  So you must keep that in mind when plotting geologic orientations.

First, setup two guides that have an intersection at a convenient location on the artboard; lock your guides.  Use the line tool and snap-click on that intersection.  Define your angle (keeping in mind differences in orientation), and length (stretch).  Do this over and over at the same intersection point until you have exhausted your supply of data.

At this point you may need to copy sets of stretch lines, group them, duplicate, and horizontally or vertically reflect them to get a complete ellipse.  Ensure that their origins snap to your intersection of guides.

Now, use the Ellipse tool and setup a unit circle with radius  S=1 and snap the center of that circle to the intersected guide.  Duplicate the circle and use the shear tool to fit your data; use the scale tool if area cannot be preserved.  If all goes well, you'll end up with Lines of No Finite Extension at the intersection of the unit circle and your ellipse.  Snap drag lines to each intersection to define your fields of lengthening and shortening.

You can then use the incredibly handy Shape builder tool and text tools to clean up your ellipse for others to view.

Remember, to preserve the angular relationships of your resultant strain ellipse when scaling, use the same scaling factor for the horizontal and vertical.  It is really useful to scale so that you can draw (by hand) reasonable tangent lines to your ellipse to find ψ at any geologic orientation.  You could measure this in Illustrator, but I think that it would not be very efficient.

Sunday, January 30, 2011

Apparent Dip Script Revisited

The apparent dip script mentioned a few posts below has some issues with due north trending cross section lines.  It's calculating bogus apparent dips.  I need to revisit it at the conclusion of Winter quarter to fix bugs.  Until then, I really don't recommend relying on it.  Unless, of course, you want to make it better.  But then I'd recommend just writing the script from scratch, since I used the rapid prototyping model.  It is really ugly work.

Just a warning.

Friday, January 21, 2011

Stereonet Apps

For OS X, there is a stereonet app called Stereonet, by R. W. Allmendinger.  It is a little long in the tooth but it works great for plotting massive amounts of planes and lines.  It can also do a best cylindrical fit, Kamb contours, plot poles, get intersections, rotate your data, and a heck of a lot more.  It also exports to PICT using vectors: nice for use in Illustrator or your preferred vector editing application.  Highly recommended for geologists of all levels.

If I were a masochist (Objective C isn't my cup of tea), I'd code a Stereonet App for iOS.  It'd be beautiful on the iPad!  And I wouldn't charge anything for it.

Maybe this summer...

Thursday, January 20, 2011

Taking Data with a Consumer-Level GPS

Last year I wrote a Python script that could take a GPX file loaded with specially described tracks and waypoints and generate a series of shapefiles with appropriate geometries.

For example, you could name a track "contact" and it would close the track to create a ring and thus a polygon.  Or describe two points such as "f1,rl" (where 1 is a unique ID) to get a fault line attributed as right lateral.  Finally, it could also take waypoints with the description "sd,270,30" to generate strike and dip points.

I haven't messed with that script in a long time and probably won't start because it's often more work using a consumer level GPS rather than plotting directly onto a basemap.

But I think that using a GPS to record strike and dip data is a good idea.  Strikes and dips are often taken at contact boundaries, and with a GPS fix, you could make adjustments to your mapped contacts.  It also can save a lot of time as you don't have to digitize your S/Ds.  In addition, you'll have them in a tabular format (DBF) ready to process in any other software.

There is one caveat of course:  your GPS unit has built-in uncertainty due to government restrictions.  If your GPS provides an uncertainty value, make sure it isn't resolvable at your mapscale!  Or is at least within the uncertainty of your own mapping abilities (e.g., plotting strike and dip based on topo lines and real world topography).

That said...

You can download your waypoint data with all kinds of software.  Garmin Basecamp for Mac, gpsbabel, etc.  QGIS can directly import or download the GPX data as a layer from your GPS.  You then simply export the layer, load the resultant shapefile, and give it a projection (it should already have a datum such as WGS 84).  The shapefile has all the attribute data provided by the GPX file such as x,y coordinates, elevation, comment/description, and time and date.

On my Garmin, I put the following into my waypoint comment/description field: 330,10.

That's it.  I didn't put any other identifiers since I only took strike and dips  using the GPS.  I used the waypoint numbers as my station number in my field book and put them on my map.  I recorded the strike and dips into my field book as well as the UTM coordinates and the unit name that I took the reading on.  I also put the strike and dips on my map.

Within my GPX file, the strike and dip ends up in a field called "comment" and "description".  I used QGIS to convert the GPX waypoints into a point shapefile.  Then I used ArcGIS 10 to take care of the rest.

In ArcGIS 10, you can write a simple python function within field calculator to extract the strike or dip:

Pre-logic code:

def getStrike(val):
 sd = val.split(",")
 if sd[0].isalpha() == True: # in case you accidentally put a alphabetic character in your comment field
  return # return nothing. null.
  return int(sd[0]) # return first element in the list: strike.  Change data type to integer, as we get it as a string and that'll fail in a short int field.

Field calculator code:
getStrike( !comment!)

You could easily modify the above so that you can tell the function if you want a strike or dip and then return the appropriate value.  This splits "330,10" into a Python list of "330" and "10".  It then returns the first element of that list with sd[0].  Easy processing of GPS comments/descriptions into useful attribute data!

I should note that Arc 10 has a strike and dip symbol with the strike trending due north.  You no longer have to create a special expression to change your strike azimuth into a dip direction for symbol rotation.

Thursday, January 13, 2011

Initial thoughts on ArcGIS 10

I am fortunate enough to now have an educational license of ArcGIS 10, as I am taking a formal class on concepts and methods in GIS.

Overall, I like what I see.  But I am a little confused as to why ESRI removed the drop down menus for the Spatial Analyst and 3D analyst.  Is there a new reclassify raster tool?  A different way to generate a hillshade?  Or generate contours over an entire raster extent?  I found the old tools by customizing the toolbar, but their icons are an unsightly hammer that makes you think there's something better to be using... (and considering the revamp of Arc's raster infrastructure, there must be).

Update: See this help document on why Spatial Analyst toolbar changed.  Hint:  Integration.  Use the Toolbox and remember to set your Geoprocessing Environment!

Also, the editor has been completely revamped and it took me a few minutes to understand what was going on.

Output quality for JPEG (used in my class) has been improved considerably.  So too has loading speeds, in general.  Very nice.

One thing I do not like: ArcGlobe is very difficult to navigate around... especially compared to Google Earth.  But it is a great idea to make import of elevation data and draping imagery super easy.

Tuesday, January 4, 2011

Winter Quarter, 2011

I have begun class at UC Davis.  This means I will not be scripting unless there is a great reason to do so.  However, I may script calculations that are often repeated (such as equations of the Mohr Circle from last quarter).

In addition to the second part of structural geology, I'm taking another formal GIS course that focuses on concepts and methods.  It doesn't hurt to improve my formal coursework.

I am also taking Geologic Oceanography / Marine Geology.

For Structural Geology lab, I will try to use GRASS and QGIS whenever possible.  However, since I have access to a very modern GIS lab (all but Arc 10), I may just appease the folks at Redlands, CA this quarter.

But already I have a GRASS location for Rainbow Basin, and I've looked at the data using nvis.  The dataset includes 2009 NAIP imagery (Cal-Atlas), a LiDAR DEM product (OpenTopography!), and a USGS DRG (Cal-Atlas). We'll see how much use that location will get.  I hope to use r.profile this quarter and my ApparentDip script.