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;

 

 

GDAL Slopes– Local Hydrologic Slope vs. the Standard Approach

Open Source software is not, of course just about warm and fuzzies, great support, rapid development cycles, shared costs, etc., it’s also about getting your hands dirty with someone else’s code and implementing stuff more quickly and more intelligently because of it, and hopefully learning something in the process.  You don’t have to poke under the hood to drive the car, but sometimes it’s nice to anyway.

The two normal approaches to calculating slope in GIS raster tools (used in applications from ESRI’s Spatial Analyst and GDAL’s gdaldem) are Horn’s formula and Zevenbergen & Thorne’s formula, both integrated calculations of slope that do a bit of generalization in the process of calculating slope.  This is particularly useful if you are using slope for more than a localized calculations, for example in doing estimates of local light conditions.  That said, if you want to match local hydrologic slope (something which will more closely approximate a field measurement of slope), you are better off using a different formula altogether.  An article a colleague sent me: “When GIS Slope is Not What You Think” in this month’s The Forestry Source tipped me off to this, and I thought I’d do some quick and dirty digging in the GDAL code to see if I could implement something a little different.

Here is the standard Horn approach as currently implemented in GDAL:

float GDALSlopeHornAlg (float* afWin, float fDstNoDataValue, void* pData)
{
    const double radiansToDegrees = 180.0 / M_PI;
    GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData;
    double dx, dy, key;

    using std::max;
    dx = (max((afWin[0], afWin[3], afWin[3], afWin[6])) -
          max((afWin[2], afWin[5], afWin[5], afWin[8])))/psData->ewres;

    dy = ((afWin[6], afWin[7], afWin[7], afWin[8]) -
          max((afWin[0], afWin[1], afWin[1], afWin[2])))/psData->nsres;

    key = (dx * dx + dy * dy);

    if (psData->slopeFormat == 1)
        return (float) (atan(sqrt(key) / (8*psData->scale)) * radiansToDegrees);
    else
        return (float) (100*(sqrt(key) / (8*psData->scale)));
}

Borrowing from gdaldem roughness we implement this:

float GDALSlopeHydroAlg (float* afWin, float fDstNoDataValue, void* pData)

{
    // Hydrologic Slope is the max
    //  local slope btw center cell and adjacent cells

    const double radiansToDegrees = 180.0 / M_PI;
    GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData;
                double key;

    float pafLocalMin = afWin[0];
     float pafLocalMax = afWin[0];

    for ( int k = 1; k < 9; k++)
    {
        if (afWin[k] < pafLocalMin)
        {
            pafLocalMin=afWin[k];
        }

        if (afWin[k] > pafLocalMax)
        {
            pafLocalMax=afWin[k];
        }
    }

    key = afWin[4] - pafLocalMin;

    if (psData->slopeFormat == 1)
        return (float) (atan(sqrt(key) / (2*psData->scale)) * radiansToDegrees);
    else
        return (float) (100*(sqrt(key) / (2*psData->scale)));
}

Output is rougher than the original, as expected:

Local Hydrologic Slope

Zevenbergen & Thorne’s Slope

And looking at the stats we see a larger standard deviation for local hydrologic slope:

Z:\>gdalinfo Hydrologic_slope.tif
Driver: GTiff/GeoTIFF
Files: Hydrologic_slope.tif
...

Metadata:
STATISTICS_MINIMUM=0
STATISTICS_MAXIMUM=46.168109893799
STATISTICS_MEAN=10.206419514979
STATISTICS_MEDIAN=2.8041767686283e-266
STATISTICS_MODE=1.4463930749435e-307
STATISTICS_STDDEV=5.3794018830967
LAYER_TYPE=athematic

Z:\>gdalinfo Zevenbergen_Thorne.tif
Driver: GTiff/GeoTIFF
Files: Zevenbergen_Thorne.tif
...
Metadata:
STATISTICS_MINIMUM=0
STATISTICS_MAXIMUM=52.424499511719
STATISTICS_MEAN=3.0666069784504
STATISTICS_MEDIAN=4.331070843283e-184
STATISTICS_MODE=-4
STATISTICS_STDDEV=3.567882216098
LAYER_TYPE=athematic

A couple of additional notes.  First, this is 2.5ft resolution data.  At this resolution, the differences in algorithms may be somewhat academic.  Also, of note– the range for the Zevenbergen-Thorn slope algorithm is greater than that for the local hydrologic slope.  This is an outcome I did not anticipate.

 

Update:

Hmmm, looks like I have more revisions and questions from the GDAL dev board. May have to revisit this… .  Forgot about the corners of my 3×3 window… .

Landscape Position: Conclusion? (part 2)

From earlier post:

“I’ve managed to pilot most of a fast high resolution landscape position workflow with PovRay as my magic tool. The final steps I hope to pipe through PostGIS Raster. In the meantime a screenshot and description: blues are riparian, raw ocre, etc upland categories, grey is mostly flat lake plain and mid slopes, all derived from just a high res DEM input (no hydro lines as input, so it works on areas where you don’t know where the streams may be). There will be more categories in the final, in the meantime, welcome to December.”

I didn’t use PostGIS Raster, but the incredibly useful gdal_calc.py to finish out my landscape position calculations.  I will endeavor to post my code to github soon, but the parts and pieces include using gdal and PovRay.  PovRay helps with the sampling (really!) of nearby neighbors in the raster DEM, and gdal does the averaging and differencing of those to get relative landscape position.  I spent some time yesterday scratching my head over how to show all the new landscape position information on a readable and useable map, and after discussion with collegues, decided to use it to divide the world into two categories– riparian & lowland + headwaters & upland (not reflected yet in the labels).  To find out more about landscape position, follow this link, or better yet this one.  (BTW, the green is park land, blue is riparian/lowland/stream networks, purple is the basin boundary).

Riparian map based on landscape position, calculated using PovRay and GDAL.

Landscape Position: Conclusion?

I’ve managed to pilot most of a fast high resolution landscape position workflow with PovRay as my magic tool. The final steps I hope to pipe through PostGIS Raster. In the meantime a screenshot and description: blues are riparian, raw ocre, etc upland categories, grey is mostly flat lake plain and mid slopes, all derived from just a high res DEM input (no hydro lines as input, so it works on areas where you don’t know where the streams may be). There will be more categories in the final, in the meantime, welcome to December.

EcoHackNYC– Cool projects, fun new ideas, human waste, CartoDB and other flotsam.

I took a bus to New York City this weekend to enjoy the company of fellow hackers at EcoHackNYC, organized by Javier Torre and Andrew Hill of Vizzuality and Robin Kraft of REDD Metrics.  Due to delays in Pittsburg, I missed the ignite talks on Friday, arriving on NYU’s campus on Saturday morning. Several groups formed around topics and we started hacking. I worked on Robin Kraft’s crew– helping to put together a draft visualization of deforestation data for Indonesia. I suspect that something interesting will continue to evolve out of that project. Robin envisions visualizations for tropical deforestation on a monthly basis globally, and he’s not far away from that. It should be really great to see. We also got some wild-eyed data visualization/HTML5/data transfer protocol ideas out of that project, thinking about how to stuff any sort of data into a PNG tile.  The cool thing is that while our wild-eyed plans would require special data prep, they would not necessarily require changes to existing middlewhere/tile-rendering software initially, just client level magic. More on that later, unless we discover we’re re-inventing aluminium rims and someone has already built a hot rod.

It turns out, when I stayed in a hotel in Midtown Manhattan, had it been raining hard enough to activate the CSOs, my waste would have come out just upstream from the UN. True story.

For now though, I’ll highlight one of the more polished products presented at the end of the hacking event the Don’t Flush Me project (warning– potty humor and mild curse words involved). I did not work on this project. With so many pun opportunities I probably would have derailed the project with linguistic abuses of monumental puniness, so maybe it’s for the best. I’m holding back right now. I just want you all to appreciate that. I like potty puns so much, they are my first and second most favorite pun type. That said, the best pun of the night was by a gentleman named Francois who was with the big-carbon-footprint group (meaning they flew in from Brazil for the conference– it was nice to have the international contingent). I’m always particularly impressed with puns done in an individual’s second language. It shows true mastery. Of course, now I can’t remember the pun. Anyway, I digress substantially more than usual.

Quick 3rd party description of the project: Don’t Flush Me took the combined sewer overflow (CSO) pipe diagrams and outfall data for Manhattan and created an interface that geocodes from address or IP and returns a polygon that shows the shared sewershed and overflow location for the input location.

Additionally, it has the Wunderground API wrapped in, in order to report whether there has been rain in the last day for your sewershed. Further work on the project would model capacity vs. rainfall and report whether you should wait to flush, and other such features. There are some small browser-dependent geocoding api bugs, but that’s forgivable considering this was put together in fewer than 8 hours. Also, the code is reportedly not particularly well organized yet. But who cares– the functionality is awesome. As an aside, if memory serves me, the back end PostGIS services are hosted in a CartoDB instance.

I think Don’t Flush Me could serve as a great model for reporting mechanisms for Sewer Districts with CSOs. Brilliant, well-scoped, and well-executed work, also fun to use. If you are not in Manhattan, here are a couple of addresses you can feed into the geocoding engine:

11 West 53 Street  New York, NY 10019

1 Wall Street, New York, NY

and an homage to the kindness of strangers:

Chelsea Park, New York, NY.

Postgis for breakfast: ST_Donut — Revision

A commenter on my last ST_Donut post pointed out that we were essentially using not one but two buffers and and ST_Intersection to test where a point lay.  Bad Spatial SQL. Very bad Spatial SQL. Actually, at Paul Ramsey’s PostGIS for Power Users presentation at FOSS4G this year, I think he mentioned not doing that very thing, and I chuckled to myself, “Well, even I know better than that.” So much for false pride.

So here’s the revision, adapted as before from Sorokine:

Revised Oct 28th

Postgis for breakfast: ST_Donut

This post typed into my iPod as an homage. Assisted today by my collegue, J. Stein.

Moderate obfuscation of locations is an important technique for the protection of data, say something sensitive like the nesting locations of the very rare and strange fuzzy-bellied gnat catcher. We still want to display the data, but want to make it slightly wrong.

A naive approach would place it a random distance away, i.e. Location +- random()*500ft. Unfortunately, this could be a distance as small as 0 ft away, not an outcome we want.

So, instead we constrain our placement of random points to a minimum and maximum distance away. Enter ST_Donut. We need to construct a donut, and then place a random, ahem, sprinkle on it, and return a point:

CREATE OR REPLACE FUNCTION ST_Donut (
  geom Geometry, dough numeric, nut numeric, geom2 Geometry 
 )
 RETURNS Geometry
 AS $$
DECLARE

  donut Geometry;
  
BEGIN  
  donut = ST_SetSRID(ST_Difference(ST_Difference(ST_Buffer(geom, dough), ST_Buffer(geom, nut)), ST_SRID(geom)),geom2);

RETURN donut;
END;

$$ LANGUAGE plpgsql;

 

Thanks to sorokine:

 CREATE OR REPLACE FUNCTION RandomPoint ( geom Geometry ) RETURNS Geometry AS $$ DECLARE maxiter INTEGER := 1000; i INTEGER := 0; x0 DOUBLE PRECISION; dx DOUBLE PRECISION; y0 DOUBLE PRECISION; dy DOUBLE PRECISION; xp DOUBLE PRECISION; yp DOUBLE PRECISION; rpoint Geometry; BEGIN -- find envelope x0 = ST_XMin(geom); dx = (ST_XMax(geom) - x0); y0 = ST_YMin(geom); dy = (ST_YMax(geom) - y0); WHILE i maxiter THEN RAISE NOTICE 'number of interations exceeded max'; END IF; RETURN rpoint; END; $$ LANGUAGE plpgsql; 

Mmmm, donut…

Edit: October 8, 11:06

The real function is:

CREATE OR REPLACE FUNCTION ST_Donut (
  geom Geometry, dough numeric, nut numeric, geom2 Geometry
 )
 RETURNS Geometry
 AS $$
DECLARE

  donut Geometry;

BEGIN
  donut = ST_SetSRID(ST_Difference(ST_Buffer(geom, dough), ST_Buffer(geom, nut)), ST_SRID(geom));

The code I had in the block at the top is the overloaded version. In this case it can take in two geometries; the second geometry is a limiting geometry– say a boundary outside of which you don’t want to place the randomized point. PostgreSQL let’s you overload functions, so if two functions have the same name, but different variable input, whether by type or number or both, it will run the function that matches the name and inputs. This allows us to create two functions with the same name, but slightly different purposes. BTW, this overdetermined function should probably check to see if the second geometry is a polygon, but it does not. Buyer beware… .

Landscape Position Continued– absolutely relative position calculation <– Pics!

Input:

Output:

Landscape Position Continued– Median and ImageMagick

Highlighting ridges with 250ft buffer (on 2.5ft DEM) with just ImageMagick:


convert lscape_posit.png -median 100 median100.png
composite -compose difference lscape_posit.png median100.png difference_median100.png


Input:

Output:

BTW, median calculations of this size are slow, even in ImageMagick.

Landscape Position Continued– absolutely relative position calculation

I apologize in advance– this first post will be heavy on code, short on explanation. Landscape position, e.g. previous posts, can be trivial to calculate, but to make the calculations scalable to a large area, some batching is necessary. In this case, instead of a McNab index, we’re calculating the traditional GIS landscape position. Enter my favorite non-geographic tool, PovRay… . To the difference between, we’ll use PovRay in combination with imagemagick:

#include "transforms.inc"
#version 3.6;

#declare widthx=3425;
#declare heighty=1707.5;

#declare Aperture=300;
#declare Laps=20;

#declare Ind=1-pow(1-clock,1.2);
#declare PosX=cos(2*pi*Laps*Ind)*Aperture*Ind;
#declare PosY=sin(2*pi*Laps*Ind)*Aperture*Ind;

#declare Camera_Location = ;
#declare Camera_Lookat = ;

camera {
	orthographic
	location Camera_Location
	look_at Camera_Lookat
	right widthx*x
	up heighty*y
}

background {color  }

union {

	height_field {
		png "dem.png"
		scale ;
		translate ;
	}

	pigment {
		gradient x color_map {
			[0 color rgb 1]
			[1 color rgb 0]
			}

		scale 
		Reorient_Trans(x, Camera_Lookat-Camera_Location)
		translate Camera_Location
	}

	finish {ambient 1}
}

povray +Ilscape_posit.pov +Olscape_posit.png +FN16 +W1370 +H683 +KFI1 +KFF99 -D
convert lscape_posit??.png -average average.png
composite -compose difference lscape_posit.png average.png difference.png

In truth, I think I could do all of this in imagemagick, but it might not be fast enough. More testing to follow... .