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

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

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

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’”)

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

gdal_calc.py– Raster Calculator using Numby Functions

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 gdal_calc.py 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
gdal_calc.py -A input1.tif -B input2.tif --outfile=result.tif --calc="A+B"

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

# example 3 - set values of zero and below to null
gdal_calc.py -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:

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

GDAL, MrSid, and nearblack

Translating MrSid lossy compressed files into uncompressed imagery has its drawbacks, including licensing and artifacts.
Old versions of fwtools, which includes the GDAL utilities (and more), were compiled with a license that allowed for the translation of MrSid into e.g. Erdas Imagine images or GeoTiff. The licensing changed on that library, so FWTools and MS4W don’t do this anymore.

If you have a compiled version of the GDAL utilities with a legal license for converting, you still run into the issue of artifacts, e.g.

Nearly, but not quite black pixels at the edge of MrSid lossy compressed image

GDAL’s nearblack is a great utility for correcting this problem. If we wanted to batch this up on a windows system, it might look something like this:

FOR %f IN (*.sid) DO nearblack -o %~nf.img %f

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

One thing I’ve learned in the last few years is that there is no theoretical reason why (properly indexed and summarized) data cannot be displayed at all scales as quickly as at any scale. This is the principle at work behind the extraordinary efficiencies of delivering data and imagery through a slippy map interface like Google, Bing, and OpenLayers as well as efficient thick client interfaces like Google Earth. So, in principle, and largely in practice, serving spatial data for a whole county or the whole world shouldn’t be any more onerous than serving data for a particular site, so long as you have adequate storage for the pre-rendered and summarized data, and time to pre-render the data. As storage tends to be cheaper than processing and network speed, this is a no-brainer.

A number of great Open Source tools exist to help with serving large amounts of data efficiently, not the least of which is my favorite, GeoServer (paired with GeoWebCache). For serving imagery, in our case orthorectified 0.6-inch (0.1524 meter) aerial imagery, we have a few options. GeoServer does natively support GeoTiff, but for this large an area at this level of detail, we’d have to wade into the realm of BigTiff support through the GDAL extension, because we have 160GB imagery to serve. We could use wavelet compressed imagery, e.g. MrSid or ECW or Jpeg2000, but I don’t have a license to create a lossless version of these, and besides, storage is cheaper than processors– wavelet compressed imagery may be a good field solution, but for server side work, it doesn’t make a lot of sense unless it’s all you have available. Finally, there are two data source extensions to GeoServer meant for large imagery, the ImageMosaic Plugin, and the ImagePyramid Plugin. The ImageMosaic Plugin works well for serving large amounts of images, and has some great flexibility with respect to handling transparency and image overlap. The ImagePyramid extension is tuned for serving imagery at many scales. The latter is what we chose to deploy.

The ImagePyramid extension takes advantage of gdal_retile.py, a utility built as part of GDAL that takes an image or set of images and re-tiles them to a standardized size (e.g. 2048×2048) and creates overviews as separate images in a hierachy (here shown as outlines of the images):

But here’s the problem– for some reason, I can’t load all the images at once. If I do, only the low resolution pyramids (8-foot pixels and larger) load. If I break the area into smaller chunks, most of them fewer than 2000 images, they load fine.

1" = 50' scale snapshot

1:32,000 scale snapshot

Historical USGS Quads

This will be a very quick post today.  We wanted to use some historical USGS quads for analysis.  We had a ready source, used ArcGIS to georeference them to a grid, but then wanted something approaching a seamless product.

Enter gdal.

First step, create an image to put the mosaic into:


gdal_merge -o usgs_merge.tif -createonly input1.tif input2.tif ...


Then mosaic into that.  Since we have a shapefile with the shape of the bounds of each of the images in question, we can use a cutline to perform the mosaic, e.g.

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

In all honesty, this was scripted in Excel for expediency, but could easily be done in Windows Command with for-in-do, or BASH at a UNIX machine.

Final result:

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:

Mapping places unknown– free global datasets and FOSS GIS are a great combo

I wanted to put together a quick and dirty map of a biological reserve in Ecuador, sort of a laptop exploration of a place quite distant.  At first, I thought I’d use Shuttle Radar Topography Misssion data to get the elevation information.  Then I discovered the ASTER Global DEM which is 30m resolution for the whole world.  Wow.  Cool cool data. (I used the Japanese site to download, as I had trouble getting the US site to load).

So, first I downloaded and extracted the zipped geotiffs that come from the ASTER Global DEM.  I extracted the DEMs into the same directory, and mosaicked them into a single DEM:

gdalwarp AST*.tif ecuador_aster_dem.tif

Then I created a hillshade version of the DEM for mapping:

gdaldem hillshade ecuador_aster_dem.tif ecuador_aster_hillshade.tif

and then 100 and 50 meter contours for the site:

gdal_contour -a elev -3d -i 100 ecuador_aster_dem.tif ecuador_aster_contour_100.shp
gdal_contour -a elev -3d -i 50 ecuador_aster_dem.tif ecuador_aster_contour_50.shp

I loaded it all into Quantum GIS, symbolized it, and added some OpenStreetMap data.  To do this, I navigated to the location in OSM, chose the download tab, and chose OpenStreetMap XML Data as my output.  Using the OSM plug-in in QGIS, I imported the file, and walla– a nearly functional map.