Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Posts Tagged ‘Postgre’

Using PostGIS for Hydrologic Modeling (reblog)

Posted by smathermather on June 17, 2016

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, […]

via Filtering Roads from Extracted Streams Data — GeoKota

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

Cartographic tricks and tips– making text readable, PostGIS edition, take 3

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 in PostGIS | Tagged: , , | Leave a Comment »

Cartographic tricks and tips– making text readable, PostGIS edition part 2

Posted by smathermather on July 2, 2013

I have a couple of posts number 1 and number 2 about masking for text.  One more from the actual implementation:

forest

cave

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

Cartographic tricks and tips– making text readable, PostGIS edition

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:

image010

golf_orig ——–>golf_final

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 in PostGIS | Tagged: , , | 1 Comment »

Using ST_Node in PostGIS (where once I used ST_Union)

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 in PostGIS | Tagged: , | 2 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 »

Cartography and USGS — Fake Building Footprints in PostGIS now with distance operator (part 2)

Posted by smathermather on May 17, 2012

In a previous couple of posts (this one, and this one), we dealt with point rotations, first with bounding box searches, and then with nominal use of operators.

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 in Analysis, Database, Database Optimization, PostGIS, PostgreSQL, SQL | Tagged: , , , , , , | Leave a Comment »

Going deeper into web cartography: future=past? (and Swiss cartographic genius)

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.

Not created in GeoServer (nor TileMill– from Imhof’s book)

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?

Cartographic representation choices with overlapping and tightly packed contours, also from Cartographic Relief Presentation

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.

I don’t know yet where or how (if) I’ll implement this, but it’s an interesting an potentially solvable problem.  Now where to work on this in the stack– in PostGIS, because I know it best, in the renderer (GeoServer) where I’ll be swimming in Java of a level I don’t have a chance, or in the browser/client, where there are already some examples written in Javascript… .

Posted in Analysis, Cartography, Database, PostGIS, PostgreSQL, SQL, Trail Curation, Trails | Tagged: , , , , , , , | Leave a Comment »

Cartography and USGS — Fake Building Footprints in PostGIS now with distance operator

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 in Analysis, Cartography, Database, PostGIS, PostgreSQL, SQL | Tagged: , , , , , | 3 Comments »

PostGIS for Dessert: Sketching shapes in #PostGIS– compass roses revisited

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:

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