Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Archive for the ‘Landscape Position’ Category

Landscape Position using GDAL — PT 3

Posted by smathermather on November 25, 2014

More landscape position pictures — just showing riparianess. See also

https://smathermather.wordpress.com/2014/11/22/landscape-position-using-gdal/

and

https://smathermather.wordpress.com/2014/11/24/landscape-position-using-gdal-pt-2/

valleyz3 valleyz2 valleyz1 valleyz

Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL

Posted in Analysis, Ecology, GDAL, Landscape Position, Other, POV-Ray | Tagged: , , , , , , , , , , | Leave a Comment »

Landscape Position using GDAL — PT 2

Posted by smathermather on November 24, 2014

Just pretty pics today of estimated riparianess. If you prefer a bit of code, see previous post https://smathermather.wordpress.com/2014/11/22/landscape-position-using-gdal/

valleys_improved valleys_improved_1 valleys_improved_2 valleys_improved_3 valleys_improved_3.1 valleys_improved_4 valleys_improved_5

Posted in Analysis, Ecology, GDAL, Landscape Position, Other, POV-Ray | Tagged: , , , , , , , , , , | 4 Comments »

Landscape Position using GDAL

Posted by smathermather on November 22, 2014

Hat tip again to Seth Fitzsimmons. I’ve been looking for a good, easy to use smoothing algorithm for rasters. Preferably something so easy, I don’t even need to write a little python, and so efficient I can run it on 30GB+ datasets and have it complete before I get distracted again by the next shiny project (a few hours).

Seth’s solution? Downsample to a low resolution using GDAL, then sample back up to a higher resolution in order to smooth the raster. My innovation to his approach? Use Lanczos resampling to keep location static, and get a great smooth model:

Unsmoothed DEM

Unsmoothed DEM

Smoothed DEM

Smoothed DEM

Code to do this in gdal follows. “-tr” sets our resamping resolution, “-r lanczos” sets our resampling algorithm, and the “-co” flags are not strictly necessary, but I’ve got a 30GB dataset, so it helps to chop up the inside of the TIFF in little squares to optimize subsequent processing.

gdalwarp -tr 50 50 -srcnodata "0 -32767" -r lanczos  -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" oh_leap_dem.tif oh_leap_dem_50.tif
gdalwarp -tr 10 50 -srcnodata "0 -32767" -r lanczos  -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" oh_leap_dem_50.tif oh_leap_dem_10-50.tif

At first this excited me for cartographic reasons. We can use this to simplify contours, and then use simplified contours at different zoom levels for maps:

But, we can also use this for analyses. For example, if we difference these smoothed images with our original digital elevation model, we get a measurement of local elevation difference, the first step in establishing where valleys, ridges, and other land forms are.

# Resample to lower resolution
gdalwarp -tr 328.0523587211646 328.0523587211646 -srcnodata "0 -32767" -r lanczos  -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" oh_leap_dem.tif oh_leap_dem_328.tif
# Upsample again to get nicely smoothed data
gdalwarp -tr 3.048293887897243 3.048293887897243 -srcnodata "0 -32767" -r lanczos  -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" oh_leap_dem_328.tif oh_leap_dem_3-328.tif
# Merge two datasets together into single image as separate bands to ensure they are the same dimensions
# (gdal_calc, as a wrapper for numpy requires this)
gdal_merge -separate -o oh_leap_dem_3-328_m.tif oh_leap_dem.tif oh_leap_dem_3-328.tif
# And now we'll use gdal_calc to difference our elevation model with the smoothed one to get relative elevation 
gdal_calc -A oh_leap_dem_3-328_m.tif -B oh_leap_dem_3-328_m.tif --A_band=1 --B_band=2 --outfile=oh_leap_dem_lp_328.tif --calc="A-B"

So, if we want a good proxy for riparian zones, we can use a technique like this, instead of buffering our streams and rivers a fixed distance (in this case, I’ve used 4 different distances:

Map of landscape position estimated valleys in Cuyahoga County, Ohio

Map of landscape position estimated valleys in Cuyahoga County, Ohio

Pretty snazzy riparian finder. It seems to also find upland headwater wetlands (most of them historic and drained for Cuyahoga County). I am now running on 4 million acres of Ohio at a 10ft (~3 meter) resolution. It’s that efficient.

Addendum: It also finds escarpment edges, like the Portage Escarpment in the above, so it is a mix of a few landforms. Darn handy nonetheless.

Posted in Analysis, Ecology, GDAL, Landscape Position, Other, POV-Ray | Tagged: , , , , , , , , , , | 2 Comments »

Kicking the tires of PostGIS 2.0 — Testing ST_MakeValid

Posted by smathermather on April 20, 2012

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;

 

 

Posted in Database, Ecology, Landscape Position, PostGIS, PostgreSQL | Tagged: , , , , , , , , | Leave a Comment »

Landscape Position: Conclusion? (part 2)

Posted by smathermather on December 7, 2011

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.

Posted in Analysis, Ecology, Landscape Position, Other, POV-Ray | Tagged: , , , , , , , , , , , | Leave a Comment »

Landscape Position: Conclusion?

Posted by smathermather on November 30, 2011

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.

20111130-221953.jpg

Posted in Analysis, Ecology, Landscape Position, Other, POV-Ray | Tagged: , , , , , , , , , , , | 1 Comment »

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

Posted by smathermather on July 14, 2011

Input:

Output:

Posted in Ecology, Image Processing, ImageMagick, Landscape Position, POV-Ray | Tagged: , , , , | Leave a Comment »

Landscape Position Continued– Median and ImageMagick

Posted by smathermather on July 14, 2011

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.

Posted in Ecology, Image Processing, ImageMagick, Landscape Position | Tagged: , , | Leave a Comment »

Landscape Position Continued– absolutely relative position calculation

Posted by smathermather on July 14, 2011

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

Posted in Ecology, Image Processing, ImageMagick, Landscape Position, POV-Ray | Tagged: , , , , | 1 Comment »

Landscape Position and McNab Indices (cont.2)

Posted by smathermather on February 1, 2011

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:

 

Posted in Landscape Position | Leave a Comment »