Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Archive for December, 2011

2011 in review — Laziest post possible

Posted by smathermather on December 31, 2011

Happy New Year everyone!  In a beautifully global fashion, the first New Year wishes I received this year were out of Africa from two folks I follow on twitter.  In that same vein, my blog is comfortably global, with readers on all the continents, Antarctica notwithstanding.  As I used to do polar research, I suppose I should change that 6 out of 7 to 7 out of 7, but it could be that it might be hard to get decent numbers down there.  WordPress may also not tabulate them.  Regardless, hello to summer residents of Antarctica, all 5-thousand or so of you.  I love polar stereographic.  UTM wouldn’t be the same without its partner UPS.

Anyway, I am humbled by my readership, people who keep me in line with smart comments, and guide my growth in the geospatial world.  Report excerpt and link as follows:

“The stats helper monkeys prepared a 2011 annual report for this blog.

Here’s an excerpt:

The concert hall at the Syndey Opera House holds 2,700 people. This blog was viewed about 14,000 times in 2011. If it were a concert at Sydney Opera House, it would take about 5 sold-out performances for that many people to see it.

Click here to see the complete report.

Posted in Other | 2 Comments »

gdal_warp, cutlines, and cwhere– simple tip for use on Linux

Posted by smathermather on December 17, 2011

Mini GDAL tip of the day:

gdalwarp, especially in combination with gdal_merge, is a powerful tool for doing all sorts on nice aggregation (read: mosaic’ing) of spatial raster data.  Unfortunately, at least with a google search, there’s very little to be found on demonstrating the use of queries in conjunction with cutlines, probably because in general these queries are not difficult to figure out.

In a previous post, I demonstrated using this to stitch together USGS quads.  I found out when I tried to run this code on a Linux machine, the command like code was a little different.  Instead of:

gdal_merge -o usgs_merge.tif -createonly input1.tif input2.tif ...
gdalwarp -cutline quadrangle_index_24k.shp -cwhere "quadname_ = West_View" West_View.tif usgs_merge.tif

I found I needed some extra single quotes to get my string literal to not be interpreted as a field:

gdalwarp -cutline quadrangle_index_24k.shp -cwhere "quadname_ = 'West_View'" West_View.tif usgs_merge.tif

Ah, that’s better.  Now I have a nice big mosaic of DEMs for my region with perfect cuts where the overlaps used to be.  (The change is hard to see– look for single quotes around West_View in “quadname_ = ‘West_View'”)

Posted in GDAL | Tagged: , , , | Leave a Comment »

GDAL Slopes– Local Hydrologic Slope vs. the Standard Approach

Posted by smathermather on December 15, 2011

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);
        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)

        if (afWin[k] > pafLocalMax)

    key = afWin[4] - pafLocalMin;

    if (psData->slopeFormat == 1)
        return (float) (atan(sqrt(key) / (2*psData->scale)) * radiansToDegrees);
        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


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

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.



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

Posted in Analysis, Ecology, GDAL | 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 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 »

Leaflet, GeoServer, and Open Source Software Ramblings

Posted by smathermather on December 6, 2011


I posted this post about an apparent problem in rendering of GeoJSON in Leaflet, and now it’s fixed a week later.  Why gush now, when for a few years developers on the GeoServer board have been fixing bugs I’ve found (often by the next day)?  Well for one, I very rarely find bugs in GeoServer– I just think I do, document everything, send it on, and then find out I’ve made a mistake.  It’s hard to gush about making a mistake and then having the GeoServer community kind enough to show me the errors (and they are kind about it). I should gush anyway, but ego sometimes prevents it. For example, the last time I found a “bug” in GeoServer, one of the developers fixed it so no one else could make the same mistake (the mistake included forgetting to install fonts) and did it while laid up with an injury.

With Leaflet, I guess what’s really nice is the kindness of strangers, so it comes as a surprise.  The GeoServer listserve has come to feel like family.  They take great care of me and I take them for granted.  I don’t know the Leaflet developers from Adam, so it feels quite warm and fuzzy to have my bugs bugs I’ve discovered squashed anyway, all for pride in a codebase that should function right.

The short and long of it: warm and fuzzies all around. Yay, Open Source.


Posted in Database, GeoServer | Tagged: , , , | Leave a Comment »– Raster Calculator using Numby Functions

Posted by smathermather on December 1, 2011

gdal_calc is a python script that makes it easy to do band math and logical operations with gdal/numby. This ranges from basic arithemtic operations to logical functions. And while has been around since 2008, it is but is a recent and revelatory discovery for me. I had just today resigned myself to properly learning python so as to use the gdal bindings. But, my laziness may continue unabated. Not learning Python. Now that is lazy.

Anyway, want to do band math? It’s as simple as follows (straight from the python script itself):

# example 1 - add two files together -A input1.tif -B input2.tif --outfile=result.tif --calc="A+B"

# example 2 - average of two layers -A input.tif -B input2.tif --outfile=result.tif --calc="(A+B)/2"

# example 3 - set values of zero and below to null -A input.tif --outfile=result.tif --calc="A*(A>0)" --NoDataValue=0

So, I had little modification to average my set of images from PovRay: -A input.tif -B input2.tif -C input3.tif -D input4.tif -E input5.tif --outfile=result.tif --calc="(A+B+C+D+E)/5"

I might write some bash wrappers to create built in functions for basic things I’ll do repeatedly, like average a set of images based on a directory listing. Of course, a little voice inside says “You should just do that in Python.” We’ll see.

Posted in Analysis, GDAL, Image Processing | Tagged: , , , , | 2 Comments »