Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Archive for May, 2012

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

Posted by smathermather on May 18, 2012

Just pretty pictures today.

Posted in Cartography, Database, Optics, Other, PostGIS, PostgreSQL, POV-Ray | 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 »

GeoServer and efficient delivery of raster data (image pyramid layer) (update)

Posted by smathermather on May 11, 2012

A perennial favorite on this blog is “GeoServer and efficient delivery of raster data (image pyramid layer)“. I am neither the last nor the first authority on this topic (check the GeoSolutions blog for authoritative work on GeoServer and raster, also look to the GeoServer documentation), but I’ve had some good experiences with serving rasters in GeoServer, especially using image pyramid layers

Read the original, as this will just augment, but here are some targets to hit with the retiling necessary for larger datasets. This is the command I currently use for the retiling:


gdal_retile.py -v -r bilinear -levels 4 -ps 6144 6144 -co "TILED=YES" -co "BLOCKXSIZE=256" -co "BLOCKYSIZE=256" -s_srs  EPSG:3734 -targetDir aerial_2011 --optfile list.txt

Don’t be afraid of big block sizes. Bump your memory up in your application container, stop worrying and learn to love the larger tif. I keep my total number of output tifs to no more than 2000, where I start to see performance issues in my implementation.
Also, give the image pyramid code a break. After retiling, do this:

mkdir 0
mv *.tif 0/.

Posted in GeoServer, GeoWebCache | Tagged: , , , , | 1 Comment »

GeoExt/Mapfish MultiPage Prints

Posted by smathermather on May 11, 2012

I can’t remember if I posted on this, but there is a cool demo of a GeoExt multipage print here:

http://dev.geoext.org/sandbox/pgiraud/playground/print-multipage/geoext.ux/ux/MultiPagePrint/examples/MultiPagePrint.html

Definitely an interface for advanced users, but nice nonetheless.

Posted in Other | 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: , , , , , | Leave a Comment »

Advanced Cartography in #GeoServer– #SLDs can make for very pretty maps

Posted by smathermather on May 10, 2012

I said it in the title, I’ll say it again, SLDs are not just for ugly maps.

I’ve heard from several professionals in this geospatial sector that CartoCSS is so cool– “look at the great cartography you can get out of it”. And, truth be told, creating a standard for cartography that aligns so well with existing web development standards is a brilliant way to get designers in on the map development process (as demonstrated aptly by e.g. Stamen), something sorely needed since Geography decided that cartographers weren’t needed anymore, some time in the late 90s/ early double aughts. That said, and at the risk of sounding parochial, beautiful maps can be produced using SLDs in GeoServer too. The below map is a snapshot from our GeoServer instance, and evolved from a set of cartographic conventions we’ve been developing in our shop, adapted in 4 hours by a darn smart part time employee (Tom Kraft) who had never built an SLD before.

Without further addo:

Posted in GeoServer | Tagged: , , , , | 6 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 »

OpenLayers, GeoExt, GeoServer, and GetFeatureInfo

Posted by smathermather on May 6, 2012

I wrote an earlier post on using GetFeatureInfo through OpenLayers to bring back a formatted html document with pictures, and formated tables, etc.  It wasn’t sophisticated, but got the job done.  Since around that time, as I’ve been building out our services, the speed with which a GetFeatureInfo request returns has g o t t  e  n  p    r     o     g       r       e        s         s            i             v              e             l            y     slower, to the point where the feature is essentially unusable.  Why?  Poor client coding.  Here’s my getfeatureinfo codeblock for OpenLayers:


////// WMS GetFeatureInfo
var info = new OpenLayers.Control.WMSGetFeatureInfo({
             drillDown : false,
queryVisible : true,
panMapIfOutOfView : false,
             url : GeoserverWMS,
layerUrls : [GeowebcacheURL],
eventListeners : {
getfeatureinfo : function (event) {
popup = new OpenLayers.Popup.FramedCloud(
                             "popinfo",
map.getLonLatFromPixel(event.xy),
null,
event.text,
null,
                             true);
map.addPopup(popup, true);
}
}
});

map.addControl(info);
info.activate();

//  end of popup code

The problem?  This code will query all visible layers.   Lots of visible layers results in a GetFeatureInfo request that looks like this:
localhost:8080/geoserver/wms?&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&LAYERS=summer_aerial_2,summer_aerial_1,reservation_boundaries_public_private_cm_dissolved_mask_gradien,odot_interstate,odot_us_routes,odot_state_routes,reservation_bounds,detailed_hydro_view,cm_bridge_view,cm_trails,impervious_update,cm_buildings,cm_buildings_outline,golf_view,nhd_lake_erie,supplementary_shields,planet_osm_line,cuyahoga_street_centerlines_labels,planet_osm_line_outside_cuy,detailed_hydro_labels,facilities_cm,facility_areas_cm&QUERY_LAYERS=summer_aerial_2,summer_aerial_1,reservation_boundaries_public_private_cm_dissolved_mask_gradien,odot_interstate,odot_us_routes,odot_state_routes,reservation_bounds,detailed_hydro_view,cm_bridge_view,cm_trails,impervious_update,cm_buildings,cm_buildings_outline,golf_view,nhd_lake_erie,supplementary_shields,planet_osm_line,cuyahoga_street_centerlines_labels,planet_osm_line_outside_cuy,detailed_hydro_labels,facilities_cm,facility_areas_cm&STYLES=,,,,,,,,,,,,,,,,,,,,,&BBOX=2232424.250515%2C625357.428203%2C2233363.463113%2C625592.231353&FEATURE_COUNT=10&HEIGHT=213&WIDTH=852&FORMAT=image%2Fpng&INFO_FORMAT=text%2Fhtml&SRS=EPSG%3A3734&X=320&Y=91
My naive attempts to fix this were simply adding a layers key/value pair, e.g.:

////// WMS GetFeatureInfo
     var info = new OpenLayers.Control.WMSGetFeatureInfo({
drillDown : false,
queryVisible : true,
panMapIfOutOfView : false,
             url : GeoserverWMS,
layers : [reservation_bounds],
layerUrls : [GeowebcacheURL],
eventListeners : {
getfeatureinfo : function (event) {
                     popup = new OpenLayers.Popup.FramedCloud(
                             "popinfo",
map.getLonLatFromPixel(event.xy),
null,
event.text,
null,
                             true);
map.addPopup(popup, true);
}
}
});

map.addControl(info);
info.activate();

//  end of popup code
… but for whatever reason, maybe the way I instantiated the OpenLayers.Layer.WMS instances using Ext.each (ahem, the way I modified a vendor’s cleverly written code…), the layers aren’t recognized as objects by the name I might expect.  So, I hack and create a duplicate layer by a different name.  I’ll have a non-hack version here soon (I hope) which will use the existing array of WMS instances, loop through them, and only add the appropriate ones to the layers array.  In the mean time, this should work.
var reservation_bounds = new OpenLayers.Layer.WMS("Reservation Boundaries", GeoserverWMS,
{'layers': 'base:reservation_bounds', transparent: true, format: 'image/png'},
             {isBaseLayer: false}
);

////// WMS GetFeatureInfo
var info = new OpenLayers.Control.WMSGetFeatureInfo({
drillDown : false,
queryVisible : true,
panMapIfOutOfView : false,
             url : GeoserverWMS,
layers : [reservation_bounds],
layerUrls : [GeowebcacheURL],
eventListeners : {
getfeatureinfo : function (event) {
popup = new OpenLayers.Popup.FramedCloud(
                             "popinfo",
map.getLonLatFromPixel(event.xy),
null,
event.text,
null,
                             true);
map.addPopup(popup, true);
}
}
});

map.addControl(info);
info.activate();

//  end of popup code

Posted in GeoExt, GeoServer, Javascript, OpenLayers | Tagged: , , , , | Leave a Comment »