Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Archive for August, 2012

Reprocessing imagery with GDAL and BASH — and then cleaning up afterward

Posted by smathermather on August 13, 2012

I show a simple example today of using gdal tools to change from raster to polygons with Bash.  That’s right.  Still no pythonista here.  I have poked ever so gently at the surface of ruby and rails recently, while watching another coder re-tune his brain to a new set of principles of least surprise.

By the way, for my regular readers, please be patient with me over the next few quarters.  I am endeavoring to undertake a much larger writing project for the next few months, which will leave little time for this blog.  But, hopefully at the end, it will be worth everyone’s while.


#!/bin/bash

for name in $( ls ????_???_gt.tif);
do

basename1=`echo $name | awk -F'.' '{print $1}'`              # Split off non-extension portion of name

if [ -f $basename1".shp" ]
then
echo $basename1".shp exists"
else
echo "Converting " $basename1".shp"

gdal_polygonize.py $name -f "ESRI Shapefile" $basename1".shp"

if [ $? -lt 1 ]
then
echo "Conversion a success."
echo "Removing original file."
rm $name
else
echo "mehmeh.  try again."
fi

fi

done

Posted in GDAL | Tagged: , | Leave a Comment »

What is the center line of a complex polygon? Routing Stream and Rivers Part Two

Posted by smathermather on August 10, 2012

The series

and

one additional example:

Posted in Analysis, Database, PostGIS, PostgreSQL, SQL | Tagged: , , , | Leave a Comment »

What is the center line of a complex polygon? Routing Stream and Rivers

Posted by smathermather on August 9, 2012

Continuing my posts on the centerline of a complex polygon (you can work your way backward from here), here’s it’s application to a riverine system (legend: light blue river, black centerline derived from routing through the skeleton of voronoi polygons from a densified stream bank set + the skipped skeleton bits in pink– read through the series if you if you don’t grok that explaination):

Posted in Analysis, Database, PostGIS, PostgreSQL, SQL | Tagged: , , , | 8 Comments »

Proper (ab)use of a database, contour interpolation using #postgresql #postgis #arcgis

Posted by smathermather on August 6, 2012

Anyone who has been following along at home knows I don’t think much like a DBA.  Sometimes that’s good; mostly it’s probably bad.  In this post, I hope it will be interesting.

The problem of the day is how to take engineering contours derived from breaklines, a lidar point cloud, and all the lot, and do a good job interpolating that to a DEM.  This is an interesting problem, as it’s another one of those problems with so many sub-optimal solutions.  In part, it’s an impossible problem– going from a surface model to contours is inherently a loss of information, so going the other way is difficult to do believably, both with respect to believable smoothness and the requisite feature retention.  This is a little like the pugilistic pie puppets problem: ‘flaky vs. tender’ of making good pie crust.  Doing one well seems to preclude doing the other well.

Now, Chris Hormann, of PovRay fame has a solution, there are solutions in ArcGIS, and in GRASS, and plenty of other attempts.  Here’s a decent list, which I haven’t checked for timeliness:  http://vterrain.org/Elevation/contour.html.  Since we have a ArcGIS license, it handles the dataset I have OK, the results checked out OK, this is what we chose to use for now.  I tried GRASS, but it was too slow for the amount of data being processed.  The command in ArcGIS ends up looking something like this:


TopoToRaster_sa "N2260635.shp elevation Contour;N2260640.shp elevation Contour;N2260645.shp elevation Contour;N2265633.shp elevation Contour;N2265640.shp elevation Contour;N2265645.shp elevation Contour;N2270635.shp elevation Contour;N2270640.shp elevation Contour;N2270645.shp elevation Contour;N2275650.shp elevation Contour" "C:\Temp\TopoToR_N2269.tif" 2 "2265031 645070 2270065 639945" 20 # # ENFORCE CONTOUR 40 # 1 0 # # # # # #

And thus we go from this:

To this:

Now to do it 623 more times… .  And in walks PostgreSQL for some abuse as a spatial engine with string manipulation to boot.

Each tile of contours I run through this approach needs the adjacent 8 tiles of contours loaded in order to avoid dreaded edge effects:

Except that sometimes, there aren’t 8 tiles available, like on the edge of large water bodies:

So, we can’t just use a heuristic that loads the 8 adjacent tiles– we need to enumerate these tiles.  An ST_Touches should do the trick here:

SELECT p.tile, a.tile AS tiles
FROM cuy_contour_tiles AS p, cuy_contour_tiles As a
WHERE ST_Touches(p.geom, a.geom);

But, that does not arrange our table very nicely:

I want something that looks like this:

N2225655, N2225660, N2220650 …

N2225660, N2225665, N220665 … etc.

Thank goodness for Postgres’ String_agg, and moreover for Leo Hsu and Regina Obe’s blog.

Adapting their entry, we can reorganize those entries like this:

SELECT
STRING_AGG(a.tile, ',' ORDER BY a.tile)
As tiles
FROM cuy_contour_tiles AS p, cuy_contour_tiles As a
WHERE ST_Touches(p.geom, a.geom)
GROUP BY p.tile
ORDER BY p.tile;

Ah, now that’s better:

And if we want, we can use an arbitrary delimiter, and build our entire command in one swell foop:

SELECT p.tile,
'TopoToRaster_sa "'
|| 'W:\contours\cuy_contours\shapefiles\' || p.tile || '.shp elevation Contour;'
|| 'W:\contours\cuy_contours\shapefiles\'
|| STRING_AGG(a.tile, '.shp elevation Contour;W:\contours\cuy_contours\shapefiles\' ORDER BY a.tile)
|| '.shp elevation Contour;" '
|| '"v:\svm\dem\cuyahoga\' || p.tile || '.tif" 2'
|| ' "' || replace(replace(replace(ST_Extent(p.geom)::text, 'BOX(', ''), ')', ''), ',', ' ') || '"'
|| ' 20 # # ENFORCE CONTOUR 40 # 1 0 # # # # '
As tiles
FROM cuy_contour_tiles AS p, cuy_contour_tiles As a
WHERE ST_Touches(p.geom, a.geom)
GROUP BY p.tile
ORDER BY p.tile

Ok.  I have 623 more tiles to do.  Until next time… .

Posted in Analysis, ArcGIS, Database, PostGIS, PostgreSQL, SQL | Tagged: , , , , , , , , , , | Leave a Comment »

What is the center line of a complex polygon? Routing as the solution

Posted by smathermather on August 2, 2012

I’ve had a series of posts (including this one) on finding the center line of a complex polygon– especially with an interest in finding the center line of streams and river bodies.  We have to give credit to my colleague Tom Kraft for solving the polygon centerline problem.  Tom struck upon the very elegant solution– derive the medial axis, and then use a few end points to route between to extract just the important features.  It’s not fully automated, but the amount of user input is so nominal as to be as close to automatic as I care about.  So far, we have this implemented in a mix of ArcGIS and PostGIS, but I can easily picture a pretty simple web interface with a PostGIS back end for extracting the centerlines of polygons in a nominally interactive way.  As a teaser, a non-stream example:

Posted in Analysis, Database, PostGIS, PostgreSQL, SQL | Tagged: , , , | 1 Comment »