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

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;

 

 

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

Playing with new tools and old standards: GeoJSON, Leaflet, CartoDB across platforms

Leaflet, CartoDB, GeoJSON, and cross platform web map deployment. First some introductions:

“Leaflet is a modern, lightweight open-source JavaScript library for interactive maps for desktop and mobile web browsers, developed by CloudMade to form the core of its next generation JavaScript API. Weighting just about 21kb of gzipped JS code, it still has all the features you will ever need for you web mapping needs while providing a fast, smooth, pleasant user experience.”

Not a bad description. There’s more maturity and flexibility to OpenLayers, but still much fun to be had with Leaflet. And as both are Open Source and competing in similar markets, they inform each other, and it’s fun to watch the friendly rivalry between the communities developing the libraries.

In an earlier post, I demonstrated the use of Leaflet in conjunction with CartoDB. What’s cartodb? From their web site:

“CartoDB allows you to map data & develop location aware applications quickly and easily. With plans starting from free, take CartoDB for a test drive today!”

CartoDB is a really slick hosted service (cloud) implementation of PostGIS, similar in some respects to fusion tables functionality. Don’t like the cloud it’s hosted on, or the associated prices? Well, it’s open source, so you can put it on your private cloud too, or set up physical linux machine to host it.

Mix these two things with 3 hours of insomnia the other night, and you get this:

Nice cross platform implementation!  These are pure vector deployed maps– served through Leaflet from CartoDB as GeoJSON.  Here shown deployed to a Blackberry, iPhone, Android, and Kindle Fire (also Android), which a buddy and I pulled up as we sat in a bar tonight.  We also took a picture of this with a windows laptop also running the page for demonstrable non-webkit deployment, but this was the better (albeit fuzzy) picture.  Ah, the wonders of Open Standards.  Next time, we’ll see if we can get a Windows Mobile phone in there, along with Mac, Windows, and Linux desktops.

Still, for now GeoServer and OpenLayers are my babies for deploying GeoSpatial data, but there are a few different ways to deploy through Open Standards using Open Source software, and it’s fun to draw a little from the other geospatial tribes.

Building simple clients for MapFish — Underlying Infrastructure

In order to build simple clients for the MapFish print service, we have to understand what the protocols are that are invoked and how they function.  To do this we can read the MapFish Print Module Doc, and then modify and vamp from there.  While I was going to joke that this would be an excellent cure for insomnia, the joke fell apart when I actually read the Print Module Doc, and it was clear, concise, readable, and frankly at least as interesting as most of my blog entries (no comments here please).

But, I did that after the fact.  In other words, I did it the hard way– just like you do when you bring home that really cool electronic toy, and play with it for a couple hours before (maybe, if ever) picking up the manual.  So for me, instead of reading the excellent documentation, I sniffed the protocol using a GeoExt interface to a GeoServer/MapFish combo with Firebug.  (quick aside– if you write web stuff and are new to it, then you should know you should use Firebug to write it better and test it on the fly– or if you’re a cool kid with horn-rimmed glasses and an ironic t-shirt, do it in Google Chrome’s Javascript Console, it makes no real difference as far as outcomes, just a difference in style).

In sniffing the protocol, I saw that my request for a PDF was a POST request to the server, with a JSON object as the request.  In short, the interface converts my form information into a bunch of text (a javascript object) which it pushes to the server.  The server location in this case is: http://localhost:8080/geoserver/pdf/create.json.  The text it’s pushing is a file which reads something like this:

{
	"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", "cuyahoga_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!"
		}
	]
}

I’ll call your attention to the last little bit of code in the object:

"pages" : [{
		"center" : [2160649.7795275, 597547.8687664],
		"scale" : 6000,
		"rotation" : 0,
		"mapTitle" : "Title Here Please!"
	}

]

I was proud of myself for recognizing (through the haze of a guy who’s modified a lot of javascript, but never learned it proper-like) that this is a javascript array with just one object.  Which means, we can make it an array with more than one object.  Eureka!  multi-page pdfs with just 5-6 more lines of code:

"pages" : [{
		"center" : [2160649.7795275, 597547.8687664],
		"scale" : 6000,
		"rotation" : 0,
		"mapTitle" : "Title Here Please!"
	}, {
		"center" : [2216902.0734907, 596701.84251968],
		"scale" : 1800,
		"rotation" : 0,
		"mapTitle" : "Title Here Please!"
	}
]

Now, just to figure out how to test this out without building a web page to do it.  curl is our friend here, just a few extra flags for telling the server what we are doing with the json file (thanks to this post):

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

And now for all the json fit to print:

{
	"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", "cuyahoga_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!"
		}, {
			"center" : [2216902.0734907, 596701.84251968],
			"scale" : 1800,
			"rotation" : 0,
			"mapTitle" : "Title Here Please!"
		}
	]
}

*Updated with better formatted sourcecode, thanks to notepad++’s JSMin plugin

Ditching GeoExt– building simple clients for MapFish

I’ve been enamored with the GeoExt interface for grabbing MapFish based print services since I first saw it. It’s a slick little interface, and can even been extended for multi-page print layouts pretty easily, ala http://tinyurl.com/mapfishmultipageprint. But as I’ve started to give thought not to what an organization full of professionals needs but what a public interface should looks (and probably those interfaces for professional organizations as well, only they tend to be more tolerant of poor design…), I’ve begun to realize that there are some clever ways we can bypass the GeoExt interface for generating the print documents at all. The actual request for a document is a simple post with a JSON object that has certain properties. We can construct this object all sorts of ways.

So, in the interest of making the snake eat its tail, my objective over the next few posts is to create an entry in a PostGIS database that has a view that automatically parses the “best” multi-page print for the geometry, feeds that back like a good API to the client script, which then requests a pdf print based on those criteria through the GeoServer Mapfish extension. Are we clear as a computer program, usually running over the Internet, that allows multiple users to participate in virtual-reality role-playing games? Good.

PostGIS Cartographic Effects– Cartoonify Nearly Coincident Lines

In my previous post, a long 24-hours ago, I proposed some automatic modification of line for cartographic reasons. I had some flaws in my code. The points were over-rotated by 45 degrees. Can you spot why? Tip: it’s a basic trigonometric mistake. Here’s the corrected code (though there may be a better way):

CREATE TABLE movedRoad AS
SELECT
    ST_Translate(b.the_geom,
        0.5 * ST_Distance(b.the_geom, ST_Union(a.the_geom)) * (cos(ST_Azimuth(b.the_geom, ST_ClosestPoint(ST_Union(a.the_geom), b.the_geom))) - 3.14159/4),
        0.5 * ST_Distance(b.the_geom, ST_Union(a.the_geom)) * (sin(ST_Azimuth(b.the_geom, ST_ClosestPoint(ST_Union(a.the_geom), b.the_geom))) - 3.14159/4) )
            AS the_geom
    FROM road2_points b, road1 a
        GROUP BY a.id, b.the_geom
    ;

An alternate approach is to only move those points that are too close, ala:

CREATE TABLE movedRoadAgain AS
	SELECT
        ST_Translate(b.the_geom,
	        (100 - ST_Distance(b.the_geom, ST_Union(a.the_geom))) * (cos(ST_Azimuth(b.the_geom, ST_ClosestPoint(ST_Union(a.the_geom), b.the_geom))) - 3.14159/4) ,
	        (100 - ST_Distance(b.the_geom, ST_Union(a.the_geom)))* (sin(ST_Azimuth(b.the_geom, ST_ClosestPoint(ST_Union(a.the_geom), b.the_geom))) - 3.14159/4) ) AS the_geom
	    FROM road2_points b, road1 a
	        GROUP BY a.id, b.the_geom
		HAVING ST_Distance(b.the_geom, ST_Union(a.the_geom)) < 100

	UNION ALL

	SELECT
		b.the_geom AS the_geom
	    FROM bridle_points b, apt a
	        GROUP BY a.id, b.the_geom
		HAVING ST_Distance(b.the_geom, ST_Union(a.the_geom)) >= 100
	    ;

PostGIS Cartographic Effects– Cartoonify Nearly Coincident Lines

I’m still working on this query, but I thought I’d post what I’ve done so far. My intent is to produce scale-dependent exaggeration of the distances between quasi-parallel lines. The reason for this is so that lines such as street lines which are nearly coincident at a particular viewing scale can be spread from each other, much in the same way great cartography lies a little to display a relatively correct map. The nice thing about digital cartography is that as we zoom in, we can make the lie a smaller and smaller, until it is just barely a fib .

My thought was to start with the simplest case– move one line away from another (stationary) line. First step is to find the closest point on the stationary line, determine the distance and angle to that point, and then translate the point outward by some scale factor along that angle. Don’t take my code as gospel, though, I think there’s some deeper logical error, but see what you think:

ST_ShortestLine gives us the shortest line from a given point to the road from which we are trying to offset. This is just to help visualize the problem:

CREATE TABLE shortLine AS

SELECT
	ST_ShortestLine(b.the_geom, ST_Union(a.the_geom)) AS the_geom
    FROM road2_points b, road1 a
    GROUP BY a.id, b.the_geom
    ;

Now, we try to translate the point outward at the same angle, by half again the distance to the closest point:

CREATE TABLE movedRoad AS
SELECT
	ST_Translate(b.the_geom,
		0.5 * ST_Distance(b.the_geom, ST_Union(a.the_geom)) * cos(ST_Azimuth(b.the_geom, ST_ClosestPoint(ST_Union(a.the_geom), b.the_geom))),
		0.5 * ST_Distance(b.the_geom, ST_Union(a.the_geom)) * sin(ST_Azimuth(b.the_geom, ST_ClosestPoint(ST_Union(a.the_geom), b.the_geom))) )
			AS the_geom
    FROM road2_points b, road1 a
        GROUP BY a.id, b.the_geom
    ;

Resulting in something looking like this (red dashed line is the road we are moving away from, the dotted blue are the points making up the line we are moving, and the dotted brown are the newly transformed points, with ST_ShortestLine also visualized in thin brown):

/>