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:
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;
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… .