Cartographically Pleasing Vegetation Boundaries– with PovRay, PostGIS, and LiDAR

Just pretty pictures today.

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

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($1, road.geom))) 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.

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

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

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

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.

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

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:

Building simple clients for MapFish — cURL as a client

I have two previous posts on using MapFish (in this case, the GeoServer version) to allow for printing to hi-resolution PDF maps from the browser.  Here we use a command-line browser (cURL) to post our json to the MapFish service in order to retrieve our PDF.

I did not keep any notes from before on making json posts to the MapFish server as a means by which to test any manual configuration of the json file, so I had to rediscover this approach from pages like this.

The “@” sign below is so that curl knows I’m feeding it a file instead of the actual json to post:

curl -i -H "Accept: application/json" -H "Content-Type: application/json" -X POST --data @mapfish_landscape.json http://localhost:8080/geoserver/pdf/create.json

Kicking the tires of PostGIS 2.0 — Testing ST_MakeValid

The feature in PostGIS 2.0 that excited me most was not topology support, raster support, or 3D functions.  Ok, raster was near the top of my list.  But what I was really excited by was the ST_MakeValid function.  Sad, isn’t it?  Lack of vision probably– excited to try to solve recurring technical snafus in a computationally inexpensive way, rather than being more excited by the shiny new toys the PostGIS team has brought us for applying to new problems.

Never shy to try the impossible before trying the practical, I threw a really nasty 4 million acre, 3 meter vectorized landscape position dataset a colleague of mine put together.  It’s a nice dendritic, complicated mess, with polygon validity issues:

Always one for doing tests, I compared the speed with which ST_MakeValid fixed the polygons vs. Horst Duester’s cleanGeometry pgsql function, which I posted about earlier.  Sadly ST_MakeValid was much slower, and completed with error.  On the other hand, the cleanGeometry function dropped polygons, although I haven’t had time to look into why.  I don’t like my functions dropping data… .

Having posted to the PostGIS forum, Martin Davis poked a bit at the data and concluded that my problem was “self-touching polygons”, and suggested I try buffering the data a zero distance.  Ahh– it makes sense that this is the issue– it’s polygonized raster data in shapefile form.  Self-touching polygons are valid in the ESRI world but not in the simple features world (see Paul Ramsey’s presentation on PostGIS for Power Users from Foss4G 2011 for more on this).

So, ST_Buffer(the_geom, 0), and we’re good to go in half the time of the cleanGeometry function, without losing data.

Stay tuned for an ST_MakeValid wrapper which will do this trick for self-touching and self-intersecting polygons before trying to fix with ST_MakeValid.  Code below.

psql -U postgres -d test -f "C:\Program Files\PostgreSQL\9.1\share\contrib\postgis-2.0\legacy.sql"
CREATE TABLE tpi_clean AS
 SELECT gid, id, gridcode, "class name", cleanGeometry(geom) AS the_geom
 FROM tpi;

CREATE TABLE tpi_valid AS
 SELECT gid, id, gridcode, "class name", ST_MakeValid(geom) AS the_geom
 FROM tpi;

CREATE TABLE facepalm AS
  SELECT gid, id, gridcode, "class name", ST_Buffer(geom, 0) AS the_geom
  FROM tpi;

 

 

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

Quick and fun post tonight.  Remember in USGS quads all the little building footprints that represented civilization?  We (me and my colleague John Stein) were contemplating how to pull off something similar with address points.  Here was our first attempt:

It looks ok, but may be a little crude to be considered cartography (click on it to see it bigger– you’ll see those buildings don’t look quite right… ).  So, what we wanted to do of course was to rotate these buildings relative to the road network.  With ST_ClosestPoint in PostGIS (1.5 and beyond) this is easier than advertised:

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;

and Walla! (or voila, or something):

Ahh, now that, with some additional work, should result in some very nice cartography.  The buildings now follow the angle of the road (so long as we always use a square… ).  USGS quads are dead.  Long live USGS quads!

Building simple clients for MapFish — Beginnings of a PL/pgSQL function

I’ve had a couple of other posts (1 and 2 and 3 and) on simple clients for MapFish.  I like the client server infrastructure for MapFish– with the client end of things built up in GeoExt, it makes for a really elegant combo.  But I’d like articulate my vision for simple clients for MapFish a little further.  One thing that seems quite feasible is to embed the JSON for the MapFish requests in a PostgreSQL table.  Why there and not just within our client?  Well, we can use PostGIS to construct really clever multi-page prints if we want to, build into PostGIS the logic to decide the orientation, number of pages, scale, and other information needed to decide how best to print this object, and we can access that JSON through a GetFeatureInfo Request through any WMS compliant server (e.g. GeoServer).  In this way, we can use the GetFeatureInfo bubble as a place where we have links (enhanced with a little javascript) to post the JSON to our MapFish service and return a PDF.

Any object we want exposed through our interface could have a link associated with it that generates a pdf map of that object.  Let’s start with the functionality we want in our PostgreSQL function and figure out what it needs to generate the JSON we want. Here’s what we want our JSON to look like, at least for a very simple example:

{
	"units" : "ft",
	"srs" : "EPSG:3734",
	"layout" : "1) LETTER 8.5x11 Portrait",
	"dpi" : 300,
	"serviceParams" : {
		"locale" : "en_US"
	},
	"resourcesUrl" : "http://maps/geoserver/www/printing",
	"layersMerging" : true,
	"preferredIntervalFractions" : [0.1, 0.2, 0.4],
	"metaTitle" : "Title Here Please! GIS Print",
	"metaAuthor" : "Title Here Please!",
	"metaSubject" : "Title Here Please! GIS Print",
	"metaKeywords" : "",
	"outputFilename" : "cm_gis",
	"legends" : [],
	"layers" : [{
			"baseURL" : "http://maps/geoserver/wms?",
			"opacity" : 1,
			"singleTile" : false,
			"type" : "WMS",
			"layers" : ["cuy_bridge_decks", "planet_osm_line_outside_cuy_map", "cuy_roads_poly", "cuy_street_centerlines", "reservation_bounds_solid"],
			"format" : "image/png",
			"styles" : [""],
			"customParams" : {
				"TILED" : "false",
				"TRANSPARENT" : true
			}
		}
	],
	"pages" : [{
			"center" : [2160649.7795275, 597547.8687664],
			"scale" : 6000,
			"rotation" : 0,
			"mapTitle" : "Title Here Please!"
		}
	]
}

As a starting point, we can split this into two sections, the global parameters, i.e. everything except “pages” (pages is what we want postgis to calculate for us).  In the most generic sense, we would want to pass all of the parameters in the global section to the function, plus the geometry of the object over which we want to print the extent, plus the actual print size of the printable area for the desired layout have it return the json, with a population of pages section done by a little PostGIS magic. PL/pgSQL to come… .

Nice post on installing pgRouting on Ubuntu

Just a quick link to a post I found on installing pgRouting on Ubuntu:

http://obsessivecoder.com/2010/02/01/installing-postgresql-8-4-postgis-1-4-1-and-pgrouting-1-0-3-on-ubuntu-9-10-karmic-koala/