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;

 

 

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.

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

Landscape Position and McNab Indices (cont.2)

In one and two previous posts, I talked about McNab indices and what they mean and how to compute them.  This is a short post just showing another screenshot of a McNab mesoscale.  The previous image was from a stream valley running through the glaciated Allegheny Plateau.  This image is a stream cut through the soft shale of the Lake Plain of the Eastern Basin of Lake Erie:

 

Landscape Position and McNab Indices (cont.)

I typed that last one too quickly– too many typos, but my wife says I’m not supposed to revise blogs, but move on… .

So, for clarity, let’s talk a little more about McNab indices.  Field-derived McNab indices are a measure of average angle from the observer to the horizon (mesoscale landform index), or from the observer to another field person a set distance away, e.g. 10m or 18m, or whatever the plot size or settled upon (minor landform index).  Typically these are done in the field at a discreet number of directions, e.g. 4 cardinal directions, or 4 cardinal plus 4 ordinal (NE,NW SE, SW, SE) directions. The landform image in this post and the last is a calculation of mesoscale landform, which is harder to do in a classic GIS than minor landform (I’ll have a follow-up post on minor landform, probably using ArcGIS).

To calculate the values for mesoscale landform computationally, we require calculating the angle to the horizon for a certain number of directions for each point in the landscape.  This has the potential to be computationally intensive and complicated.  If, however, we conceive of the problem differently, as a 3D calculation we can perform in PovRay, arriving at the answer is simplified by the coding already done by the PovRay programmers.

Essentially, McNab mesoscale indices are a field proxy for the steradian of a site.

What does that mean?  Well, essentially, a map of steradians is a proxy for the shading of a white uneven surface on a cloudy day– or how much diffuse light is available to a given site.  Used in combination with site aspect, this is enough information to determine most of of the light conditions of a site, which is why McNab indices in combination with other factors are a good predictor of site productivity, and correlated with different plant communities across the landscape.

How does an uneven white surface shade on a cloudy day? Steeper areas with more of the sky occluded are darker while wide open spaces, like the bottom of river valleys and the tops of ridges and plateaus are brighter.  If you want to witness this effect, look to snow on the ground on a cloudy day (and what a great winter to do it).  (The only difference with snow is subsurface scattering which evens out the light quite a bit.)  You’ll notice the divots in the snow to be darker than the peaks, and the edges of the divots to be darkest of all.

The question then, is how to we compute this within PovRay?  We could use radiance as a global illumination model, but the calculation of inter-reflectivity that is at the core of radiance, while an important real factor in the landscape, would fail to replicate the original mesoscale landform index.  Instead, we set up a series of lights in a dome to illuminate the whole sky sphere, blur them a bit, and call it a day, a technique developed for HDRI rendering by Ben Weston, whose code I borrowed heavily.  The more lights the element of the landscape is exposed to the brighter, and vice versa, essentially making the brightness of an element proportional to the exposure to the sky sphere.  Unlike Ben, I used a simple white sphere to get even lighting.

Rendered at 16-bit resolution, we have a possible range of values from 0-65535.  Let’s assume that a linear transformation of these values will result in values representing the proportion of the sky sphere.  From there, transformation to steradians and then to solid angles in degrees is trivial.  Once it’s solid angles in degrees, it represents the same kind of value as a mesoscale landform index would give us (more later…).

Landscape Position and McNab Indices

Just a quick teaser post for our forestry/ecology readers out there.  I have a methodology developed for calculating McNab indices that directly corresponds with the field technique (unlike, as far as I know, any previous GIS-based techniques– which are probably adequate proxies).

What is a McNab index?  Well there are two kinds, the minor landforms and mesoscale landforms that are field-measured topographic position or terrain shape indices inform the location of ecological processes across the landscape.  So, for example, some plant forest types like ridges, and some like ravines.  The question is, quantitatively, how raviney or ridgey should it be for a given species, association, or alliance, and how is it measured?  Basically either average angle to horizon, or average angle to the local landscape, e.g. 10m away are the two McNab indices.  See i.e. http://www.treesearch.fs.fed.us/pubs/1150 and http://www.treesearch.fs.fed.us/pubs/24472.

So here’s the output (more to come), including code.  The darker the shade, the lower the relative position, the lighter the shade, the higher the valley, e.g. ridges and planes.  I know, it just looks like a hillshade, but there’s deeper stuff happening here:

Topographic Position Index and Ecological Land Type (warning completely unrefined not quite Geologic dribble– with bad maps :) …)

Warning.  What follows is somewhat informed, but I’m no geologist.  I just play one on wordpress.

Understanding the basic underlying geology and associated topography plus site history helps us achieve a basic understanding of a sites ecological potential.  At the most basic level, we expect different wildlife and vegetation dynamics in a floodplain vs. a mountain ridge.  Classification of digital elevation models can be done with a concept called topographic position index (of which Jensen Enterprises has a pretty good explanation).

I’ve been playing with TPI for two areas in the Allegheny Plateau in Ohio– one in an unglaciated portion, and one in a glaciated portion.  The Allegheny Plateau developed as the once very tall Appalachians wore down to the nubs they are today.  That plateau has since eroded into little hill-like remains, with stream valleys redefining a hill and hollow landscape, and most of the ridge-tops being near the same elevation.  In the glaciated portion, this is overlayed with glacial till, a re-flattening of these “hilltops”, and the definition of gorges in places where sub-glacial rivers subdivided the terrain.

The TPI attempt I did for the unglaciated portion roughly defines the ridge tops in red, stream valleys and slope bases in blue, and slopes as yellow/yellow green.  The map that follows is the overlay of TPI performed in gdaldem at multiple scales, but we are somewhat limited in what we can do with gdaldem, because we can’t control our window size.

Then I performed a proper TPI analysis at 80 feet and 300 feet for a glaciated portion of the plateau in Northern Ohio, following Jensen’s guidlines.  What follows is a photo of a printed map, cause I forgot to bring an jpeg of pdf home… .  Here, upland flat portions are dark green.  They are adjacent to light blue areas that are also upland, but exceed 5% slopes.  We can see upland areas adjacent to cliff areas as well as peaks in red, lower areas of high relief in orange, gorge valleys in dark blue, U-shaped valleys in a lighter blue (stream polygons and lines are also overlayed, adding to color confusion– hence the warning in the title).

What intrigues me about this map are the areas around the streams just northwest of the center of the map.  Capturing the gorges and upland areas successfully was not a surprise.  What is surprising is that we capture and map the areas that seem to me likely contribute directly to surface flow into the streams as distinct areas from the upland plains.  Next, I think I’ll compare this analysis with the soil map for the area, as well as the detailed vegetation cover we have for the area.