The Problem We have to filter out the roads and ditches without removing streams that cross roads or follow them closely. I’m going to use PostGIS to find the intersection of the streams lines data with a buffered roads polygon. If the intersected line is less than 50% of the length of the stream line, […]
Posts Tagged ‘Postgre’
Posted by smathermather on June 17, 2016
Posted by smathermather on July 6, 2013
A little more refinement to the SQL for building masking fill for hand-placed text on maps.
CREATE TABLE use_area_mask AS -- We'll use the "WITH" Common Table Expression (CTE) here WITH exploded AS ( -- we want each individual letter to get it's own mask, so we Union and Dump them to break them out SELECT (ST_Dump(ST_Union(geom))).geom FROM use_area_labels ) -- now we can create a temporary table that is a 5-unit buffer of the convex hull ,buffer_cvx AS ( SELECT ST_Buffer(ST_ConvexHull(geom), 5) AS geom FROM exploded ) -- finally, we Dump these out to their own records SELECT (ST_Dump(ST_Union(geom))).geom FROM buffer_cvx
Posted by smathermather on July 2, 2013
Posted by smathermather on June 29, 2013
Recall from previous post, http://wp.me/phWh4-qm I played a bit with using convex polygons and buffers to control noisy backgrounds behind text, ala:
I was tired of doing this in the GUI in QGIS, so I wrote some simple code to do it in PostGIS.
CREATE TABLE cm_label_cvx AS WITH buffer_cvx AS ( SELECT ST_Buffer(ST_ConvexHull(geom), 250) AS geom FROM cm_label_med) SELECT ST_Union(geom) FROM buffer_cvx;
Posted by smathermather on November 30, 2012
Just discovered ST_Node, a function in PostGIS just for noding geometries. I used to use ST_Union for this, a potentially memory intensive operation. Now I can use ST_Node.
For usage, see the PostGIS docs (requires 2.0 and GEOS 3.3.2 for bug free use).
Back to coding. That is all.
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… .
Posted by smathermather on May 17, 2012
First we create a function to do our angle calculations, then use select to loop through all the records and do the calculations.
Within our function, first we find our first (in this case) five nearest streets using our distance operator , based on bounding boxes, then we further order by ST_Distance, with a LIMIT 1 to get just the closest street. By the way, this is about twice as fast as the last approach.
And our new code:
CREATE OR REPLACE FUNCTION angle_to_street (geometry) RETURNS double precision AS $$ WITH index_query as (SELECT ST_Distance($1,road.geom) as dist, degrees(ST_Azimuth($1, ST_ClosestPoint(road.geom, $1))) as azimuth FROM street_centerlines As road ORDER BY $1 <#> road.geom limit 5) SELECT azimuth FROM index_query ORDER BY dist LIMIT 1; $$ LANGUAGE SQL; DROP TABLE IF EXISTS address_points_rot; CREATE TABLE address_points_rot AS SELECT addr.*, angle_to_street(addr.geom) FROM address_points addr;
I have to credit Alexandre Neto with much of the code for this solution.
Posted by smathermather on May 12, 2012
My favorite cartography book is Eduard Imhof’s Cartographic Relief Presentation. A few years back I picked this book up (translated to English) from ESRI press for $75 if memory serves me. Now it can be gotten for much cheaper.
Imhof spends a lot of time on feature simplification and separation, a problem which keeps me up at night. For example, if you have a lot of overlapping and/or touching contours, what are the local distortions and simplifications that can be done to enhance map readability?
This problem applies to other overlapping features, such as a road following a river, where at a given scale one might not distinguish between the two, so we exaggerate their differences, in this case, while maintaining their basic topology.
I’ve played with this problem before, but now I’m considering a new approach. I pinged Martin Davis of JTS fame for some advice. He suggested using force directed layout to solve this problem, much as he does in JTS Test Builder for magnifying/investigating topology.
Posted by smathermather on May 10, 2012
In a previous post (I feel like I say that a lot), I wrote about rotating address points to match nearby roads in replicate the effect of USGS quads that represented small buildings with little squares that followed the nearby road alignment.
The function was effective:
ALTER TABLE jstein.address_points_all ADD COLUMN azimuth integer; UPDATE jstein.address_points_all addr SET azimuth = degrees(ST_Azimuth(addr.the_geom, ST_ClosestPoint(road.the_geom, addr.the_geom))) FROM jstein.street_centerlines road;
but deadly slow when applied to all 500,000 address points. And so we iterate. First, I’ll show you our reasonably fast, but elegantly written solution. Later I’ll have a post on our very fast, but somewhat hackerish approach.
A reminder of our target aesthetic:
And our new code:
CREATE TABLE address_points_rot AS WITH index_query AS ( SELECT addr.*, ST_Distance(addr.the_geom, road.the_geom) AS distance, degrees(ST_Azimuth(addr.the_geom, ST_ClosestPoint(road.the_geom, addr.the_geom))) as azimuth FROM address_pointssub as addr, street_centerlines as road WHERE ST_DWithin(addr.the_geom, road.the_geom, 2000) = true ORDER BY addr.the_geom <-> road.the_geom ) SELECT DISTINCT ON(cartodb_id) * FROM index_query ORDER BY gid, distance
In the interest of full disclosure, I stole the structure of this from Paul Ramsey’s post on K Nearest Neighbor searches in PostGIS. AFAIK, you need Postgresql 9.1 and Postgis 2.0 or later to do this… . For 500,000 points this took 16 minutes. Narrowing that search distance would help a lot. I suspect there is a much better way to structure this using true KNN for each point, but my SQL-fu failed me here… . Suggestions welcome.
Posted by smathermather on May 10, 2012
For one of our applications, we need 8-point compass roses placed at each of our points, as well as a circle 40 meters in diameter as well as one 140 meters in diameter. We did a bit of work with this a while back in GeoServer using SLDs. Now we’d like to refine it, and implementing this in an SLD is beyond my skills.
So, we move to the back end. The following is a little function to construct the roses for us, which we’ll then display though GeoServer. Just an FYI, this is running on a PostGIS 1.3 database, and I didn’t (at the time anyway) grok the use of ST_RotateZ, so there are some explicit position calls, rather than just creating the crosshairs and rotating it:
CREATE OR REPLACE FUNCTION pie ( geom Geometry, filling numeric, crust numeric ) RETURNS Geometry AS $$ DECLARE slices Geometry; BEGIN --slices = ST_MakeLine(geom, ST_Translate(geom, dough, pan)); slices = ST_Union( ST_Collect( ST_Collect( ST_Collect( ST_MakeLine(ST_Translate(geom, filling, 0), ST_Translate(geom, crust, 0)) ,ST_MakeLine(ST_Translate(geom, 0, filling), ST_Translate(geom, 0, pi()/4 * crust)) ) ,ST_Collect( ST_MakeLine(ST_Translate(geom, -filling, 0), ST_Translate(geom, pi()/4 * -crust, 0)) ,ST_MakeLine(ST_Translate(geom, 0, -filling), ST_Translate(geom, 0, pi()/4 * -crust)) ) ) ,ST_Collect( ST_Collect( ST_MakeLine(ST_Translate(geom, filling * pi()/4, filling * pi()/4), ST_Translate(geom, crust * pi()/4, crust * pi()/4)) ,ST_MakeLine(ST_Translate(geom, -filling * pi()/4, filling * pi()/4), ST_Translate(geom, -crust * pi()/4, crust * pi()/4)) ) ,ST_Collect( ST_MakeLine(ST_Translate(geom, filling * pi()/4, -filling * pi()/4), ST_Translate(geom, crust * pi()/4, -crust * pi()/4)) ,ST_MakeLine(ST_Translate(geom, -filling * pi()/4, -filling * pi()/4), ST_Translate(geom, -crust * pi()/4, -crust * pi()/4)) ) ) ) ,ST_Collect( ST_ExteriorRing(ST_Buffer(geom, filling, 32)) ,ST_ExteriorRing(ST_Buffer(geom, crust, 32)) ) ) ; RETURN slices; END; $$ LANGUAGE plpgsql;
So now to test it with a single point:
DROP TABLE IF EXISTS test_pie; CREATE TABLE test_pie AS SELECT pie(ST_MakePoint(0,0), 40, 140);
And view in uDig, or application of choice: