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 :)

2 comments:

  1. I should note that there is a plugin for QGIS that does exactly this (albeit, more elegantly) called Profile From Line. However, it does not output a tab delimited data file. But QGIS lets you copy rows of attributes and it is formatted in a way that spreadsheets can parse into columns and rows (tab delimited).

    So porting my GDAL/OGR script to QGIS would be a waste of time :-)

    ReplyDelete
  2. Hello, Could you post your profile script to take a look at it, it would be nice to have a command line tool to do this kind of work, as opposed to have to use a gui like QGIS to make use of his tools, and by the way, QGIS Rocks!, Greeting from Mexico!

    ReplyDelete