Thursday, December 23, 2010

Raplee Ridge in 3D

Raplee Ridge was infamous in my structural geology lab.  It was our last project due by the end of finals week.  We had to measure strike and dip using lidar and lidarviewer and create a geologic map using an aerial photo.  We did so by using color descriptions and relative age of formations.  We used ArcDesktop and utilized topology.

The dataset provided to us was of very high quality.  After a bit of processing using what I learned about GRASS today, I loaded up nvis.  Raplee Ridge is beautiful.  I didn't quite get that from the project.  In fact, looking at these views, one can very quickly get an idea what the structure was that we were trying to tease apart from the lidar data using lidarviewer!


I just received Open Source GIS: A Grass GIS Approach, 3rd edition.  Time to learn!

Tuesday, December 21, 2010

QGIS Plugin Scripting: Revisited

I learned, rather hurriedly, how to script QGIS plugins for a class.  I came away unimpressed with the class hierarchy (I got lost and certain methods I needed seemed misplaced and therefore hard to find).

But now that QGIS alone is becoming a viable replacement for other software, I'm motivated to try scripting plugins again.

So I'll be taking a look at moving over the apparent dip script and cross section script to the QGIS interface.  If all goes well, I'll actually have a nice way to distribute them via the QGIS plugin repository.

Monday, December 20, 2010

Topo Profile Script Results

For the better part of the day, I wrote a topo profile script.  It allows the use of a cross section line of any geometry in a GIS to determine a topographic profile.  I probably implemented it in all the wrong ways, but it works well!    The distances for every point do in fact add up to the total line distance... and the cell values are spot on.  Very good.

As usual, scripted using Python with all the usual suspects (GDAL, OGR) but with a new comer: Numpy.

Basically, you tell it: where to find your cross section line, DEM, and a sample interval distance.  Assuming everything is projected right, it plots points.  It uses those points to sample the raster for elevations (cell values).  It also keeps cumulative track of distances between points.  It generates a shapefile with that data.  It also produces a tab delimited text file with the same data (but arguably easier to import into various graphing tools).

So we end up with something like the following (note that the profile graph was not generated by the script, just the elevation and distance data):
Generated points on a cross section line of various trends (El Mayor rupture region).  Point interval is 20 m; line vertices, start and end points, are always sampled.

Units are in meters, no vertical exaggeration.  Was created at 1:6000 scale. Zero is sea-level.
Yes, I was lazy and didn't label my axes.

I don't recommend reading ARC GRID rasters into this program as it has no smart way of dealing with massive datasets.  I kind of recommend using GeoTiffs for that reason... and clip to your work area or cross section.

The script is a mess... I'll release it eventually.  Especially if someone asks for it.

Thanks to OpenTopography for the data that I used to test this script :)

Sunday, December 19, 2010

Topo Profile: The Next Script

I'll be scripting a topo profile extractor based on cross section lines of arbitrary geometry.  It'll have the ability to scale how many elevation points you want based on distance between points.  Eventually it might find some use in conjunction with my apparent dip script (e.g., another column in the output with the apparent dip value, so it can be plotted as a marker with distance).

It should be pretty painless.  I've downloaded (for fun) a tiny slice of point cloud data for the El Mayor-Cucapah Earthquake rupture from OpenTopography.  I'll be testing cross section lines across that slice.  It's a good set since there are areas where there's no data (I didn't fill the nulls).

I think the goal of my winter is to simply script the tools to make life easier for the next two quarters.

And have fun while doing it :-0

I should note that the El Mayor data is affiliated with UC Davis.  Cool.

Friday, December 17, 2010

Updated: Geologic Symbology for QGIS

I am currently working on creating markerline SVGs for QGIS, since 1.7 trunk now allows for creation of most symbols found on geologic maps.  So far I have the following SVGs that are directly using the FGDC specifications (from their EPS files):

  • Anticline arrows (black and magenta)
  • Syncline arrows (black and magenta)
  • Asymmetric syncline arrows (black and magenta)
  • Downthrown symbol for normal faults
  • Sawtooth symbol for thrust faults
  • Left and right lateral fault symbols
Note that these do not include lines.  Those can be easily created using the QGIS symbology interface.  Generally, lines are 0.25 and 0.375 mm (0.38 in QGIS 1.7 trunk; does not allow thousandths) wide.  See FGDC guidelines if you're interested in following their standard; but frankly, as someone dear to me would say, cartography isn't always about following rules...

Take a look and/or use them here.

QGIS + GRASS = Topology

Today I embarked on a mission to figure out a workflow in QGIS that includes topology.  After searching around, I found that GRASS is fully-topological.

After playing around with GRASS, I think I have found my workflow.  I can now create contact boundaries, create centroids within them in a separate layer (boundaries + centroids = areas/polygons), and generate appropriately labeled polygons automatically.

The trick was (1) understanding how to attribute using the QGIS interface to GRASS, and (2) realizing that the Split tool needs a double click to create a vertex that can snap.

So perhaps my days requiring Arc are numbered.  Now I just need to get QGIS Trunk to compile for its awesome new symbology features.

Oh, if you haven't tried using GRASS since before QGIS 1.7, you might want to try.  Its interface is a bit more intuitive and it is way less crash prone.

Thickness of Bedding (Stratigraphy)

I imagine a script that could take elevation data (say, a raster derived from DEM) and calculate bedding thickness might be useful.

The workflow would involve:
  1. Loading a DEM
  2. Overlaying the DEM with remote imagery that has visible bedding
  3. Drawing a polygon over bedding where you need to find the approximate thickness
  4. Process the elevation data within the polygon to find the bedding thickness
    1. This would probably require generating a slope raster, perhaps average slope over the area within the polygon, measuring the various required distances and hope for the best...
Python w/Numpy and GDAL would be utilized.

Wednesday, December 15, 2010

Measuring Azimuth of a Line in a GIS

In order to find the apparent dip of bedding, strikes must be projected into a cross section line. The azimuth of that segment of cross section line must be found.

I did this in my GIS apparent dip script by creating a geospatial straight-edge and using it as a protractor.

Using a point on the cross section line at which a projected strike line intersects as the origin, I offset to the west a known distance (in this case, 0.5 map units, or as is often the case, 0.5 meters).  From that offset, I have OGR create a line feature that goes due north for 20 or so map units.  The line also extends 20 map units due south from the offset point (we want to make sure we can intersect the cross section line no matters its trend).  The idea is the straight-edge/protractor will intersect the cross section line providing it isn't trending due north.

Once the straight-edge intersects, we have a distance from our known offset point.  It is then a simple matter of using inverse tangent to find the angle that the cross section line makes from the north.

But what happens if the due north straight-edge NEVER intersects the cross section line?  Given how small of a distance it is offset from the line to the west, we can easily assume that the two lines are essentially parallel.  Thus we call the azimuth of that segment of cross section line 000° azimuth.

All of this, of course, is automated in a Python script and it is a measure made every time a projected strike line intersects the cross section line.  It does this because cross section lines can change direction.

By the way, if you have a geologic problem that would benefit from scripting, leave a comment.  If it is something that would be useful to me, I just might write a script.  Keep in mind I am an undergraduate student, so there will be a limited number of useful ideas :-)

Tuesday, December 14, 2010

Apparent Dip Script Initial Release

I've made the Apparent Dip script much more robust.  Get the latest gzip'd tarball here with my test datasets.  It requires Python with OGR.  I'll probably license it under the GPL (haven't really thought about it).

It is very cool though still not very polished and has a few caveats (read the commented code).

I think I've solved some of the obvious and easiest problems.

Brown lines are projected strike lines.  Magenta is the cross section line.  The apparent dip in the red box has its information shown in the dialog.

adip is the calculated apparent dip.  Dipdir, strike, and tdip (true dip) are copied from the original strike and dip and are really only useful for cross referencing or checking validity of the calculated apparent dip.  xsec is the name of the cross section and xsecaz is the azimuth of the cross section approximately where the projected strike line crosses.

Monday, December 13, 2010

Using GIS as an Apparent Dip Calculator

Today I prototyped half the Python script to calculate apparent dips from strike and dip point data and a cross section line.  So far, the script takes an arbitrary cross section, buffers it, selects the strike and dips within that buffer, and then projects lines from the strike of each of those points.  Then, the script finds the intersection of each of those lines with the cross section line.

The difficulty in finding the apparent dip is that one cannot assume the cross section line maintains one trend.  Nor is there a built-in OGR function (to my knowledge) that finds the angle(s) of intersection between two lines.

I solved this problem by creating an arbitrary north-south trending line offset a known distance from the determined intersection point.  From there it is super easy to find the azimuth of the cross section line!  Note that this method will fail on cross section lines oriented exactly north-south.

Next step: making the calculated apparent dips accessible with distances along the cross section line.

Projected Strikes in Magenta, Cross Section Line in Green.  Note that rotated strike and dip symbols pivot around the point's geographic coordinate, thus the appearance of offset (the strike lines go right through the point).  Buffer was 100 meters.  Arbitrary north-south oriented lines aren't in this screen shot. Viewed in QGIS.

Sunday, December 12, 2010

Winter Break Plans

I hope that, within the next few weeks, I will have done at least two of the following:

  1. Create geologic map symbols using SVG format for QGIS (I already have a E-W oriented S/D symbol for rotation using Dip Direction; removes a step ;-))
  2. Figure out a good way to emulate Arc's topology in QGIS: contact lines to unit polygons.
  3. Write a script to calculate apparent dips using strike and dip point data within a buffered cross section line
  4. Find out what Vrui does with regards to projection of geographic coordinates
  5. Relax

Wednesday, December 8, 2010

LidarViewer and Projected Coordinates

For Structural Geology lab, we were instructed to use the open sourced LidarViewer to measure strike and dips on a LIDAR point cloud of Raplee Ridge.  LidarViewer relies on an open sourced library and GUI called Vrui.  Both are products of UC Davis.

LidarViewer uses Vrui to display the LIDAR point cloud in 3D.  Vrui uses WGS 84 as its geodetic coordinate system.

Oddly, I cannot figure out what that PCS is.  It looks like a Mercator given the magnitude of the units.  Yet, the eastings and northings always mismatches the various projections found in QGIS.  World Mercator?  Nope.  At times it seemed rather non-linear.

Does anyone know what the heck Vrui is doing to get its 2D coordinates?  It is using WGS 84, so a false easting and northing would be nice to have...

I attempted understanding the code, but Vrui, being C++ code, is so abstracted that I can't make heads or tails of what is going on (I profess my love of Python!).  Perhaps it is taking whatever the projection was of the LIDAR point cloud?  If so, then the Raplee Ridge data has a very odd projection.  At least, odd in that I can't find a good match for it.

Knowing what the PCS is would be helpful in the future: I could take strike and dips plus have the coordinates saved to a file.  I like efficiency.

Monday, December 6, 2010

GIS and Cross Sections

I had the realization that a python script could take all your strike and dip data, project the strikes out some distance from each point, intersect a cross section line, find the apparent dip using the dip attribute from the S/D feature class, and provide distances along that section line to each projected strike along with apparent dips in a table.  I do believe that a QGIS plugin that does this would be extremely handy.  Heck, it might be handier yet to allow extracting a topographic profile along that section line too.

Sounds like a good winter project.  And it would save a load of time later on.