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:

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!

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/

Mapfish Play cont. Musings on PostGIS driven Mapfish requests– Code only

Code only post-- which only means it's been out here not-quite ready to post for a month.  Now I post out of shear annoyance with myself... .

SELECT '{"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":"GIS Print","metaAuthor":"","metaSubject":"GIS Print","metaKeywords":"","outputFilename":"cm_gis","legends":[{"name":"","classes":[{"name":"Reservation","icons":["http://maps/geoserver/wms?&VERSION=1.1.0&REQUEST=GetLegendGraphic&LAYER=reservation_bounds&HEIGHT=10&WIDTH=10&FORMAT=image%2Fpng&TRANSPARENT=true&LEGEND_OPTIONS=forceLabels%3Afalse&EXCEPTIONS=application%2Fvnd.ogc.se_xml&RULE=unrestricted"]},{"name":"Restricted","icons":["http://maps/geoserver/wms?&VERSION=1.1.0&REQUEST=GetLegendGraphic&LAYER=reservation_bounds&HEIGHT=10&WIDTH=10&FORMAT=image%2Fpng&TRANSPARENT=true&LEGEND_OPTIONS=forceLabels%3Afalse&EXCEPTIONS=application%2Fvnd.ogc.se_xml&RULE=restricted"]}]},{"name":"Detailed Hydro","classes":[{"name":"Ditch","icons":["http://maps/geoserver/wms?&VERSION=1.1.0&REQUEST=GetLegendGraphic&LAYER=detailed_hydro_view&HEIGHT=10&WIDTH=10&FORMAT=image%2Fpng&TRANSPARENT=true&LEGEND_OPTIONS=forceLabels%3Afalse&EXCEPTIONS=application%2Fvnd.ogc.se_xml&RULE=cm_streams_ditch"]},{"name":"Non-Stream Waterway","icons":["http://maps/geoserver/wms?&VERSION=1.1.0&REQUEST=GetLegendGraphic&LAYER=detailed_hydro_view&HEIGHT=10&WIDTH=10&FORMAT=image%2Fpng&TRANSPARENT=true&LEGEND_OPTIONS=forceLabels%3Afalse&EXCEPTIONS=application%2Fvnd.ogc.se_xml&RULE=cm_streams_ditch"]},{"name":"Stream","icons":["http://maps/geoserver/wms?&VERSION=1.1.0&REQUEST=GetLegendGraphic&LAYER=detailed_hydro_view&HEIGHT=10&WIDTH=10&FORMAT=image%2Fpng&TRANSPARENT=true&LEGEND_OPTIONS=forceLabels%3Afalse&EXCEPTIONS=application%2Fvnd.ogc.se_xml&RULE=cm_streams_stream"]},{"name":"Stream or River","icons":["http://maps/geoserver/wms?&VERSION=1.1.0&REQUEST=GetLegendGraphic&LAYER=detailed_hydro_view&HEIGHT=10&WIDTH=10&FORMAT=image%2Fpng&TRANSPARENT=true&LEGEND_OPTIONS=forceLabels%3Afalse&EXCEPTIONS=application%2Fvnd.ogc.se_xml&RULE=Stream%20or%20River"]},{"name":"Pond","icons":["http://maps/geoserver/wms?&VERSION=1.1.0&REQUEST=GetLegendGraphic&LAYER=detailed_hydro_view&HEIGHT=10&WIDTH=10&FORMAT=image%2Fpng&TRANSPARENT=true&LEGEND_OPTIONS=forceLabels%3Afalse&EXCEPTIONS=application%2Fvnd.ogc.se_xml&RULE=Pond"]},{"name":"Lake Erie","icons":["http://maps/geoserver/wms?&VERSION=1.1.0&REQUEST=GetLegendGraphic&LAYER=detailed_hydro_view&HEIGHT=10&WIDTH=10&FORMAT=image%2Fpng&TRANSPARENT=true&LEGEND_OPTIONS=forceLabels%3Afalse&EXCEPTIONS=application%2Fvnd.ogc.se_xml&RULE=Lake"]},{"name":"Other Wet Areas","icons":["http://maps/geoserver/wms?&VERSION=1.1.0&REQUEST=GetLegendGraphic&LAYER=detailed_hydro_view&HEIGHT=10&WIDTH=10&FORMAT=image%2Fpng&TRANSPARENT=true&LEGEND_OPTIONS=forceLabels%3Afalse&EXCEPTIONS=application%2Fvnd.ogc.se_xml&RULE=Other%20Wet%20Areas"]}]},{"name":"Trails","classes":[{"name":"ADA, APT","icons":["http://maps/geoserver/wms?&VERSION=1.1.0&REQUEST=GetLegendGraphic&LAYER=cm_trails&HEIGHT=10&WIDTH=10&FORMAT=image%2Fpng&TRANSPARENT=true&LEGEND_OPTIONS=forceLabels%3Afalse&EXCEPTIONS=application%2Fvnd.ogc.se_xml&RULE=rule01"]},{"name":"Bridle","icons":["http://maps/geoserver/wms?&VERSION=1.1.0&REQUEST=GetLegendGraphic&LAYER=cm_trails&HEIGHT=10&WIDTH=10&FORMAT=image%2Fpng&TRANSPARENT=true&LEGEND_OPTIONS=forceLabels%3Afalse&EXCEPTIONS=application%2Fvnd.ogc.se_xml&RULE=rule02"]},{"name":"Hiking","icons":["http://maps/geoserver/wms?&VERSION=1.1.0&REQUEST=GetLegendGraphic&LAYER=cm_trails&HEIGHT=10&WIDTH=10&FORMAT=image%2Fpng&TRANSPARENT=true&LEGEND_OPTIONS=forceLabels%3Afalse&EXCEPTIONS=application%2Fvnd.ogc.se_xml&RULE=rule03"]},{"name":"Mountain Bike Trails","icons":["http://maps/geoserver/wms?&VERSION=1.1.0&REQUEST=GetLegendGraphic&LAYER=cm_trails&HEIGHT=10&WIDTH=10&FORMAT=image%2Fpng&TRANSPARENT=true&LEGEND_OPTIONS=forceLabels%3Afalse&EXCEPTIONS=application%2Fvnd.ogc.se_xml&RULE=rule04"]},{"name":"Connector Trail","icons":["http://maps/geoserver/wms?&VERSION=1.1.0&REQUEST=GetLegendGraphic&LAYER=cm_trails&HEIGHT=10&WIDTH=10&FORMAT=image%2Fpng&TRANSPARENT=true&LEGEND_OPTIONS=forceLabels%3Afalse&EXCEPTIONS=application%2Fvnd.ogc.se_xml&RULE=rule05"]}]}],"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}},{"baseURL":"http://maps/geoserver/wms?","opacity":1,"singleTile":false,"type":"WMS","layers":["reservation_bounds"],"format":"image/png","styles":[""],"customParams":{"TRANSPARENT":true,"TILED":false}},{"baseURL":"http://maps/geoserver/wms?","opacity":1,"singleTile":false,"type":"WMS","layers":["detailed_hydro_view"],"format":"image/png","styles":[""],"customParams":{"TRANSPARENT":true,"TILED":false}},{"baseURL":"http://maps/geoserver/wms?","opacity":1,"singleTile":false,"type":"WMS","layers":["cm_bridge_view"],"format":"image/png","styles":[""],"customParams":{"TRANSPARENT":true,"TILED":false}},{"baseURL":"http://maps/geoserver/wms?","opacity":1,"singleTile":false,"type":"WMS","layers":["cm_trails"],"format":"image/png","styles":[""],"customParams":{"TRANSPARENT":true,"TILED":false}},{"baseURL":"http://maps/geoserver/wms?","opacity":1,"singleTile":false,"type":"WMS","layers":["impervious_update","cm_buildings","cm_buildings_outline"],"format":"image/png","styles":[""],"customParams":{"TILED":"false","TRANSPARENT":true}},{"baseURL":"http://maps/geoserver/wms?","opacity":1,"singleTile":false,"type":"WMS","layers":["golf_view"],"format":"image/png","styles":[""],"customParams":{"TILED":"false","TRANSPARENT":true}},{"baseURL":"http://maps/geoserver/wms?","opacity":1,"singleTile":false,"type":"WMS","layers":["nhd_lake_erie"],"format":"image/png","styles":[""],"customParams":{"TILED":"false","TRANSPARENT":true}},{"baseURL":"http://maps/geoserver/wms?","opacity":1,"singleTile":false,"type":"WMS","layers":["reservation_boundaries_public_private_cm_dissolved_mask_gradien"],"format":"image/png","styles":[""],"customParams":{"TILED":"false","TRANSPARENT":true}},{"baseURL":"http://maps/geoserver/wms?","opacity":1,"singleTile":true,"type":"WMS","layers":["supplementary_shields","odot_interstate","odot_us_routes","odot_state_routes","planet_osm_line","cuyahoga_street_centerlines_labels","planet_osm_line_outside_cuy","detailed_hydro_labels","facilities_cm","facility_areas_cm"],"format":"image/png","styles":[""],"customParams":{"TILED":"false","TRANSPARENT":true}}],"pages":[{"center":[' || ST_X(ST_Centroid(the_geom)) || ',' || ST_Y(ST_Centroid(the_geom)) || '],"scale":2400,"rotation":0,"mapTitle":""}]}'::text
FROM loops
WHERE (ST_XMax(ST_Envelope(the_geom)) - ST_XMin(ST_Envelope(the_geom))) < (ST_YMax(ST_Envelope(the_geom)) - ST_YMin(ST_Envelope(the_geom)))
AND
(ST_YMax(ST_Envelope(the_geom)) - ST_YMin(ST_Envelope(the_geom))) < 7900


SELECT 'landscape, linear follow'::text, ST_Centroid(the_geom)
	FROM loops
		WHERE (ST_XMax(ST_Envelope(the_geom)) - ST_XMin(ST_Envelope(the_geom))) > (ST_YMax(ST_Envelope(the_geom)) - ST_YMin(ST_Envelope(the_geom)))
			AND
		(ST_XMax(ST_Envelope(the_geom)) - ST_XMin(ST_Envelope(the_geom))) > 15800

UNION ALL

SELECT 'landscape, quad page'::text, ST_Centroid(the_geom)
	FROM loops
		WHERE (ST_XMax(ST_Envelope(the_geom)) - ST_XMin(ST_Envelope(the_geom))) > (ST_YMax(ST_Envelope(the_geom)) - ST_YMin(ST_Envelope(the_geom)))
			AND
		(ST_XMax(ST_Envelope(the_geom)) - ST_XMin(ST_Envelope(the_geom))) > 7900
			AND
		(ST_XMax(ST_Envelope(the_geom)) - ST_XMin(ST_Envelope(the_geom))) <= 15800

UNION ALL

SELECT 'landscape, single page'::text, ST_Centroid(the_geom)
	FROM loops
		WHERE (ST_XMax(ST_Envelope(the_geom)) - ST_XMin(ST_Envelope(the_geom))) > (ST_YMax(ST_Envelope(the_geom)) - ST_YMin(ST_Envelope(the_geom)))
			AND
		(ST_XMax(ST_Envelope(the_geom)) - ST_XMin(ST_Envelope(the_geom))) < 7900

UNION ALL

SELECT 'portrait, quad page'::text, ST_Centroid(the_geom)
	FROM loops
		WHERE (ST_XMax(ST_Envelope(the_geom)) - ST_XMin(ST_Envelope(the_geom))) < (ST_YMax(ST_Envelope(the_geom)) - ST_YMin(ST_Envelope(the_geom)))
			AND
		(ST_YMax(ST_Envelope(the_geom)) - ST_YMin(ST_Envelope(the_geom))) > 7900
			AND
		(ST_YMax(ST_Envelope(the_geom)) - ST_YMin(ST_Envelope(the_geom))) <= 15800

UNION ALL

SELECT 'portrait, single page'::text, ST_Centroid(the_geom)
	FROM loops
		WHERE (ST_XMax(ST_Envelope(the_geom)) - ST_XMin(ST_Envelope(the_geom))) < (ST_YMax(ST_Envelope(the_geom)) - ST_YMin(ST_Envelope(the_geom)))
			AND
		(ST_YMax(ST_Envelope(the_geom)) - ST_YMin(ST_Envelope(the_geom))) < 7900

UNION ALL

SELECT 'portrait, linear follow'::text, ST_Centroid(the_geom)
	FROM loops
		WHERE (ST_XMax(ST_Envelope(the_geom)) - ST_XMin(ST_Envelope(the_geom))) < (ST_YMax(ST_Envelope(the_geom)) - ST_YMin(ST_Envelope(the_geom)))
			AND
		(ST_YMax(ST_Envelope(the_geom)) - ST_YMin(ST_Envelope(the_geom))) > 15800
;


				

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.

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
	    ;