libGrid, a terrain resampling library

The Grid Library (libGrid) is the resampling backend for the libMini terrain rendering engine.

The library aims to provide a generic resampling and manipulation tool for a wide variety of geospatial terrain data and imagery building on GDAL. Its main application is to produce terrain tilesets that can be displayed with the libMini terrain rendering engine. These tilesets are generated fully automatic.

The library also provides convenient methods to manipulate geospatial raster data. This includes the reprojection, transformation, interpolation and modification of the image channels or height fields. In particular, the explicit distinction between grid-centric and cell-centric data avoids one major cause of interpolation error. The library additionally utilizes image pyramids for all geospatial operations in a transparent way so that aliasing artifacts are intrinsically avoided.

Version 3.0 as of 25.May.2010
Copyright (c) 2010 by Stefan Roettger

Table of Contents

Terms of Usage

The libGrid software is licensed under the terms of the GNU GPL. No warranty WHATSOEVER is expressed; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE! See the GNU license for more details.

Back to top

General Information

The author's contact address is:

stefan:at:stereofx.org
www.stereofx.org

Back to top

Getting the Package

The latest version of libGrid is available via svn:

   svn co http://libgrid.googlecode.com/svn/libgrid/grid grid
Back to top

Compilation

The library compiles on a wide variety of platforms, but it requires the following dependencies to be installed on the system: libMini, ZLIB, Squish and GDAL. It is highly recommended to use GDAL versions 1.6 and higher with internal libtiff support.

If you checked out from svn, libMini and Squish are already contained in the source code and ZLIB is most likely part of the installation (e.g. on Linux/Ubuntu). So the only dependency that needs to be installed to compile the source code is GDAL.

The only build system supported so far is CMake. To compile on Windows, specify the libGrid directory as the source and binary directory in the CMake GUI and click the generator button. Then open the generated VC++ project and compile with Visual Studio. On Unix (and MacOS X), simply type the following on the command line:

   cmake . && make
If some of the required dependencies cannot be found, because they have not been installed at a standard location like /usr/local, you can specify the installation directory by overriding the variable LIBGRID_THIRDPARTY_DIR in the CMake configuration (via ccmake or the Windows CMake GUI). The path to libMini can be specified separately with the variable LIBMINI_PATH.

Optionally, libGrid allows the generation of S3TC compressed images via the squish library. To disable that feature, turn off the BUILD_GRID_WITH_SQUISH cmake option.

Back to top

Overview

The purpose of libGrid is to resample terrain data and create paged tilesets that can be displayed at real time with the libMini terrain rendering engine.

Terrain data usually comes as gridded elevation, bathymetry and satellite imagery at varying spatial resolution, location and geographic coordinate systems. A single piece of data is called a layer. In order to display such a heterogeneous data collection with libMini, the data needs to be transformed into a single consistent coordinate system. Data that does not fit into the graphics card in one piece needs to be broken up into smaller pieces. This process is called resampling and implemented in libGrid.

As opposed to commercial earth viewers like Google Earth, libGrid allows you to import, view and analyze whatever terrain data you like. You are not limited by policies, the limit is just the data itself.

Aside from massaging data for use with libMini, the library also serves as a generic tool for manipulating terrain data. For libMini related uses see the following Sections, else see the Programming API and the subsequent example use cases.

Back to top

Resampling

The whole resampling process can be described as follows:

After specifying a set of input layers and an area of interest, the terrain data is transformed into the respective output coordinate system. Then the area is broken up into evenly spaced tiles yielding the layout of the tileset. The data size of a generated tile depends on the resolution of the overlapping layers. The higher the resolution of the overlapping layers, the larger the size of the tile. The higher the number of tiles, the larger the extent of a tile and the smaller the size of it. libGrid automatically determines the number of tiles such that all the tiles fit into the graphics card in one piece. With those definitions, the resampling process is fully specified and will take some time, as it is a computationally expensive operation.

Once resampling is finished, the generated tileset is ready to be viewed with libMini. To do so, the terrain renderer will only page in all those tiles, which are visible for a specific point of view and viewing distance. This minimizes the amount of terrain data residing in graphics memory and the number of triangles being displayed. To further reduce the graphics memory foot print, libMini will also only page in the level of detail that is appropriate for a specific viewing distance. Thus, a tile in the distance will page in at a lower resolution than a tile that is close by. While the amount of available graphics memory limits the resolution of the visible tiles being rendered, the size and resolution of the entire tileset, however, is principally unlimited.

Back to top

Resampler Usage

In order to generate a tileset with libGrid, the layers need to be specified in a so called grid file. In that file, one layer per line is given. The first line denotes the name of the tileset. With respect to supported layer file formats, all the available formats of GDAL library are supported. In particular the GeoTIFF format is recommended.

To start the resampling process, simply invoke the "gridder" application with the name of the grid file as single argument. For more information on the additional capabilities and the syntax of the grid file format see here.

Back to top

Simple Resampling Example

In the simplest case, the set of input layers for the gridder just consists of two layers, a single elevation grid (height field) and a single geo-referenced image or satellite photo. For the Island of Oahu the libGrid package contains sample elevation data at 10m resolution and sample imagery at 25m resolution in the data directory. This directory needs to be checked out in addition to the source code we already checked out into our working directory:

   svn co http://libgrid.googlecode.com/svn/libgrid/data data
With that terrain data, a simple grid input file for the gridder could look as follows:

   "Sample-Tileset"     # tileset name
   path "data/tilesets" # output path

   "data/elev/Oahu-10m.tif" # input elevation layer 1
   "data/imag/Oahu-25m.tif" # input imagery layer 2
For convenience, the above file is located in the "data/tilesets" directory as sample grid file "Sample.grid". Now, if you start the gridder with that file, it will create a tileset in the "data/tilesets/Sample-Tileset" directory. On Linux, this corresponds to the following command line in the working directory: 
   grid/gridder data/tilesets/Sample.grid
Due to the size of the data this will take several minutes (5min on my Apple MacBook). Please be patient and have a Kona coffee.

To view the tileset, compile and start the libMini viewer using the following commands:

   svn co http://libmini.googlecode.com/svn/libmini/mini mini-viewer
   pushd mini-viewer; cmake -DBUILD_MINI_APPS=ON .; make; popd
   mini-viewer/viewer data/tilesets/Sample-Tileset
Before navigating around Oahu, press the 'f' key in the libMini viewer to increase the viewing distance (far plane) as needed. Also, press the 'h' key to bring up the HUD, where you see all available keyboard commands. For example, to jump above the terrain press 'j'. To increase the triangle mesh detail press 'T' multiple times (upper case, lower case will reduce the detail). Likewise, to increase texture detail press 'R'.

Back to top

Detail Textures

There is also the option to add a detail texture to a tileset. Currently the only detail texture file format supported is libMini's native db format. A converter utility is provided in the tools directory, which reads all common gdal image formats and exports to the db format. The exported textures are compressed with S3TC. For example, you can add hires compressed imagery at 1m resolution for Manana Island off the east coast of Oahu by passing the according generated db file as an additional argument to the libMini viewer:

   grid/tools/grid2db data/image/Oahu-Islands/MananaIsland.tif manana.db
   mini-viewer/viewer -m data/tilesets/Sample-Tileset manana.db
Back to top

Grid File Syntax

In the following, all the options of the grid file format will be described by means of a more sophisticated version of the sample tileset. We are basically going to add more layers to cover the bathymetry around Oahu. We are also adding high-resolution imagery of Manana Island off the eastern coast. And we are going to improve the appearance of the tileset by pre-shading the terrain. This is a precondition to render bathymetry, as there is no imagery naturally available for the sea floor. Before diving into the details, here is the resulting example grid file:

   "Oahu-Tileset"                # tileset name
   repo "data"                   # input repo
   path "data/tilesets"          # output path
   level 0                       # output level
   shade fill reproject compress # output attributes

   # elevation
   "elev/Oahu-10m.tif"
   # bathymetry
   "bthy/Oahu-50m-W.tif"
   "bthy/Oahu-50m-E.tif"
   "bthy/Hawaii-1km.tif" keep_bathy
   # islands
   "elev/Oahu-Islands/KapapaIsland.tif"
   "elev/Oahu-Islands/MananaIsland.tif"
   "elev/Oahu-Islands/MokuliiIsland.tif"
   "elev/Oahu-Islands/MokuluaIslands.tif"
   "elev/Oahu-Islands/MokumanuIsland.tif"
   "elev/Oahu-Islands/MokuoloeIsland.tif"
   "elev/Oahu-Islands/SandIsland.tif"
   # etc
   "elev/MakaiPier.tif" spacing 5

   # imagery
   "imag/Oahu-25m.tif" extents
   # islands
   "imag/Oahu-Islands/MananaIsland.tif" spacing 6 move 100 30
   "imag/Oahu-Islands/MokuluaIslands.jpgintif" spacing 1.5
   "imag/Oahu-Islands/FlatIsland.jpgintif" spacing 1.5
   "elev/Oahu-Islands/MokuoloeIsland.tif" apply_shading 4 0
After the tileset name in the top line of the grid file, a variety of attributes can be specified, which control the properties of the output tileset. The attributes are as follows:

repo dir: Specifies an output path prefix for the tileset.

path dir: Specifies an input path prefix for all layers.

level value: Specifies the output resolution of the tileset. Level 0 corresponds to the original resolution of the layers, level 1 corresponds to half of that and so forth. Can be used to quickly test the tileset output at a faster resampling lower resolution.

float: Output elevation as float, otherwise use 16bit integer. 16bit elevation may introduce quantization artifacts. In that case, bump up to float.

shade: Shade the imagery based on available elevation data to emphasize interesting structures like ridges. Used mainly to shade bathymetry.

fill: Fill undefined grid cells (elevation and imagery) by extrapolating from the neighbor cells. Automatically gets rid of small glitches in the input data.

reproject: Reproject all layers before resampling. May speed up resampling, but comes at the expense of an additional interpolation, which may degrade data crispness.

compress: Compress output imagery with the S3TC algorithm. Reduces size by a factor of 6. Note: Output imagery is always mip-mapped to ensure proper display at any level of detail at which the tiles are paged in. Side note: If the output imagery contains areas with undefined cells (no-data values), the output imagery will be amended with an alpha channel to mask those areas. Another side note: Elevation output is always compressed losslessly with the ZLIB. This has no performance impact.

Likewise, after each layer, a variety of layer attributes can be specified, which alter the appearance of each layer:

extents: Specifying that keyword after a layer means that the area of interest is the extents of the corresponding layer. If no extents are specified explicitly, the area of interest is set to the entire coverage of all layers.

spacing value: Overrides the cell spacing of a layer to be the specified value given in meters. Useful for oversampled data, where the native resolution is too high for display.

move off_x off_y: Moves the layer by the specified offset to correct certain displacements where the geo-referenced position does not correspond exactly with reality (or does not match other layers).

treat_black: Treats black as transparent, when black has been used to encode no-data values instead of naturally utilizing an alpha channel.

treat_white: Same as above, except for white being treated as transparent.

remove_bathy: Removes all elevation values below the sea level with no-data.

keep_bathy: Removes all elevation values above the sea level with no-data.

apply_shading: Transforms an elevation layer into an imagery layer by shading it. Useful to generate b/w imagery for areas where there is no native imagery available.

The above example grid file is supplied as "data/tilesets/Oahu.grid", so you can generate the corresponding tileset and check out the effects of the above options with the following command:

   grid/gridder data/tilesets/Oahu.grid
At the original output level the resampling process will take a while, so please be patient. For testing purposes, you can bump up the output level though:

   grid/gridder data/tilesets/Oahu.grid 1
The default output level is 0. For each higher level, the output resolution is halved so that processing time is reduced to a quarter. Have fun flying around Oahu,
Stefan

Back to top

Command Line Tools

The libGrid package also contains command line tools to query, copy and manipulate terrain data. For example, general geospatial information about tif or db files can be extracted analogue to the gdalinfo tool with the gridinfo tool. In order to print information about the Manana Island detail texture we use the command line

   grid/tools/gridinfo data/imag/Oahu-Islands/MananaIsland.tif
yielding the output

   grid_layer:
    valid=yes
    size_x=3082
    size_y=4917
    channels=3
    layout=cell-centered
    repo=""
    loadable=yes
    path="data/imag/Oahu-Islands/MananaIsland.tif"
   grid_extent:
    leftbottom=[ (637785,2.35537e+06,0) t=0 crs=UTM zone=4 datum=NAD83 ]
    rightbottom=[ (640404,2.35537e+06,0) t=0 crs=UTM zone=4 datum=NAD83 ]
    lefttop=[ (637785,2.35955e+06,0) t=0 crs=UTM zone=4 datum=NAD83 ]
    righttop=[ (640404,2.35955e+06,0) t=0 crs=UTM zone=4 datum=NAD83 ]
   grid_properties:
    type="byte"
    nodata=0
    scale=1
    bias=0
    spacing=0.850115(0.850173/0.850058)
    range=[0.0562092..0.956863]
Another handy tool is the gridcopy tool, which allows to copy and modify grid data. The most common usage is to read an uncompressed image and compress it either with LZW or the JPG-in-TIF option. The latter uses the YCbCr color space, which leads to a very reasonable compression ratio at a high image quality using just a regular JPEG compression embedded into a GeoTiff file. LZW is the recommended format for losslessly compressed data, whereas JPG-in-TIF is the recommended format for lossy compressions. Do not forget to keep the original uncompressed data around to prevent quality degradation from multiple lossy compressions. See the usage message of the gridcopy tool for further details on the available options.

Back to top

Programming API

Overview of the core modules and their purpose:

A grid_value encapsulates one cell value with multiple channels.

For single channel data like height fields, a cell value is created with the constructor grid_value(value,1) or as an abbreviation with data_value(value). Invalid or unspecified cell values are represented with nodata_value(). Multi-channel cell values can be created from single-channel cell values by appending additional channels via value.append(channel_value).

A cell value is undefined or said to be no-data if one ore more channel values are not-a-number (IEEE quiet NAN), which can be tested with is_nodata(value). Operators are defined for the standard arithmetic operations + - and == != as well as operators for linear combinations via * and /. As an example, the average of two rgb values can be computed as follows:

   #include <grid/grid_value.h>

   grid_value rgb1(0.0,3),rgb2(1.0,3);
   grid_value avg=(rgb1+rgb2)/2.0;
   std::cout << avg << std::endl;

   -> avg = (0.5,0.5,0.5)
A grid_data object contains multi-channel raster data of various types. Possible types are byte, unsigned short, signed short or float. The types come in two flavors with an absolute value range or a normalized range of 0..1. In both cases the values can be scaled and biased and for all types there is the possibility to identify a specific value with no-data. The raster data is initialized to no-data by default.

A typical rgb image with sx columns, sy rows and with black representing no-data values is setup with

   #include <grid/grid_data.h>

   grid_data image;
   image.init(sx,sy,3, // grid size and number of image channels
              GRID_TYPE_BYTE, // normalized unsigned byte
              TRUE, // cell-centric (image)
              1.0,0.0, // scale and bias
              0.0); // no-data value
A typical height field is setup with

   grid_data hfield;
   hfield.init(sx,sy,1, // single elevation channel
               GRID_TYPE_SIGNED_SHORT_INT, // absolute signed short
               FALSE, // grid-centric (height field)
               1.0,0.0, // scale and bias
               NAN); // no-data value
For all but the byte and unsigned short type, no-data values can also be represented with NAN. The float raster type stores the no-data values directly as NAN, whereas the value range of signed short raster data ranges from -32767 to 32767 with -32768 being a reserved value that is interpreted as NAN.

The data representation can be either cell-centric or grid-centric, indicating whether a cell value should be assumed to represent a sampling over the region of a cell (area orientated) or a lattice point (point orientated). The convention of libGrid is that imagery is cell-centric and height fields are grid-centric.

The following line sets an image pixel at raster position (x,y) to a constant brightness c. The raster positions are in the range 0≤x≤sx-1, 0≤y≤sy-1 with the origin at the bottom left corner.

   image.set(x,y,data_value(c));
The number of channels of the value parameter is expanded to fit the raster data. Likewise, the accessors to the raster data return a grid_value object that reflects the number of present channels:

   grid_value value=hfield.get(x,y);
The same is true for the bi-linear interpolator that returns a grid_value for a specific coordinate (u,v)

   grid_value value=hfield.interpolate(u,v);
except that it operates on a normalized coordinate range of [0..1,0..1]. For cell-centric data (that is imagery) the coordinate ((0.5+x)/sx,(0.5+y)/sy) maps to the center of the according cell at raster position (x,y). The coordinates range from (0.0,0.0) at the bottom left corner of the bottom left pixel to (1.0,1.0) at the top right corner of the top right pixel. For grid-centric data (that is elevation) the coordinate (x/(sx-1),y/(sy-1)) maps to the according lattice point (x,y). The coordinates range from (0.0,0.0) at the bottom left lattice point to (1.0,1.0) at the top right lattice point.

A grid_layer is a superset of a grid_data object. It additionally has a geospatial extent (grid_extent) which geo-references the raster data with 4 specific corner points on earth. In contrast to grid_data objects, the interpolators of a grid_layer do not work with normalized coordinates but with geo-referenced coordinates. Each geo-referenced coordinate is represented as a minicoord vector. The layer extent is defined to exactly enclose the

Please note that this is different from the GDAL convention that a geo-referenced extent is always area oriented as outlined in the GDAL data model. The distinction of libGrid between the two cell layouts enforces the correct interpretation of cell values.

The following common geospatial coordinate reference systems (CRS) are supported for defining the layer extents:

ECEF and LatLon coordinates assume to be referenced to the WGS84 datum. For UTM coordinates a variety of other datums (not to be mixed up with ellipsoids) are supported as listed below:

For other not so common coordinate systems or datums please refer to GDAL's utilities (in particular cs2cs and gdalwarp) to convert those to a supported CRS (UTM/WGS84 recommended).

The following code snippet illustrates how to setup a geo-referenced layer:

   #include <grid/grid_layer.h>

   minicoord::MINICOORD_TYPE type=minicoord::MINICOORD_UTM;
   int zone=4;
   minicoord::MINICOORD_DATUM datum=minicoord::MINICOORD_DATUM_NAD83;
   double base=0.0;

   minicoord leftbottom(miniv3d(601498.6948,2371702.575,base),type,zone,datum);
   minicoord rightbottom(miniv3d(603513.7721,2371702.575,base),type,zone,datum);
   minicoord lefttop(miniv3d(601498.6948,2373737.497,base),type,zone,datum);
   minicoord righttop(miniv3d(603513.7721,2373737.497,base),type,zone,datum);

   grid_extent ext(leftbottom,lefttop,rightbottom,righttop);

   grid_layer layer;
   layer.init(ext,
              size_x,size_y,channels,
              GRID_TYPE_SIGNED_SHORT,TRUE,
              1.0,0.0,
              NAN);
The layer API supports a variety of file formats through a generic interface. The most common used format is geotiff as supported by the GDAL implementation of the layer API, the grid_layer_gdal class. Besides geotiff the GDAL implementation also adds support for BT (Binary Terrain) and a slew of other common geospatial file formats. For a complete list of GDAL supported formats see here.

The libGrid package contains two other implementations of the layer API, that is grid_layer_db for the DB (DataBuf) format of libMini and grid_layer_pnm for the PNM (Portable Any Map) format with geospatial extensions as supported by libMini.

With the layer API it is almost a one-liner to load a layer from a file, for example with GDAL:

   #include <grid/grid_layer_gdal.h>

   grid_layer_gdal layer_gdal;
   bool success=layer.load("filename.tif");
The layer API will only load the header information when opening a file, unless the raster data of a layer is accessed explicitly.

Satellite raster data often comes as multiple files, that is a tileset, mosaic or multiple channels. In order to group those layers together logically and perform certain operations on them such as interpolation or resampling, the layers can be stored in a dynamic array, the grid_layers class. This is done by appending pointers to each layer to the array:

   #include <grid/grid_layer.h>
   #include <grid/grid_layers.h>

   grid_layer *a=new grid_layer;
   grid_layer *b=new grid_layer;

   grid_layers *layers=new grid_layers;

   layers->append(a);
   layers->append(b);

   layers->do_something(); // e.g. interpolation, resampling etc.

   delete layers;

   delete a;
   delete b;
A grid_layers object doesn't manage the stored layers, it just keeps references to them. If we need to store a collection of managed layers that are automatically released when going out of scope, we use the grid_list container class. Layers appended to the container, will be released when the container is destroyed.

   #include <grid/grid_list.h>

   grid_layer *a=new grid_layer;
   grid_layer *b=new grid_layer;

   grid_list *list=new grid_list;

   list->append(a); // list takes control over a
   list->append(b); // list takes control over b

   list->do_something();

   delete list; // implicitly releases managed layers
The grid_list class also offers a convenient loader wrapper for all known layer API implementations. As an example we load some height fields and interpolate the elevation at a well-known location given by a minicoord geo-referenced position:

   #include <grid/grid_list.h>

   grid_list list;

   // load elevation of Oahu
   list.load(path,"data/elev/MakaiPier.tif");
   list.load(path,"data/elev/Oahu-10m.tif");

   // define coordinate at T.C.'s Helicopter Platform on Makai Pier
   minicoord crd=minicoord(miniv3d(638128,2358101,0),minicoord::MINICOORD_UTM,4,minicoord::MINICOORD_DATUM_WGS84);

   // query elevation at coordinate
   double spacing=list.get_spacing();
   double elev=list.interpolate(crd,spacing);

   std::cout << "elevation=" << elev << "m" << std::endl;

   -> elevation=2.49258m
The interpolate method returns the interpolated value of the first layer that has valid elevation (no no-data) at the specified coordinate. The spacing parameter defines which level of detail of the underlying raster data pyramid (grid_pyramid) is accessed. A spacing of zero will access the highest available resolution. The data pyramid is generated automatically on demand.

The most common use case of the grid_list class is the resampling core of its superset the grid_tileset. With the resampling core one can reproject, resample and transform collections of raster data into tilesets that can be viewed with the libMini terrain renderer. For more information about resampling see the Section about the gridder application.

Another common use case is the merging core. For an axample usage, see the Section about Pansharpening.

Back to top

Programming Example

The following sample code shows how to load a layer, manipulate it, and save it back. It is pretty much self-explanatory, so the code is left to the user to explore on their own.

   #include <grid/grid.h>

   void main()
      {
      grid_resampler resampler;
      grid_layer *layer;

      // load layer
      layer=resampler.load("path","input_filename.tif");
      if (layer==NULL) exit(1);

      // check for elevation data
      if (!layer->grid_cell_centered) // images are cell-centered, height fields are grid-centered
         {
         // replace elevation values in the range [0..1000] with no-data
         layer->replace_range(data_value(0.0),data_value(100.0),nodata_value());

         // save to geotiff format (lzw compressed by default)
         if (!resampler.save("path","output_filename.tif",layer)) exit(1);
         }
      }
Also have a look at the applications in the tools directory, in particular the "gridinfo" tool, which illustrates how to retrieve basic properties such as size and spacing from a grid layer.

Back to top

Advanced Use Cases

As outlined in the Programming API, the libGrid core classes can be utilized to perform a variety of geospatially related tasks. In the following we have a look at some advanced use cases.

Pansharpening

Satellite imagery usually does not come as plain RGB image. Instead, it comes as multiple intensity bands, where each band represents a certain narrow frequency range. The bands typically range from thermal infrared over near infrared to visible red, green and blue. Most satellites do also carry a panchromatic (b/w) detector which delivers better spatial resolution at the expense of a wider frequency spectrum. For example the Landsat 7 (ETM+) satellite provides the visible bands with 30m resolution (spectrum spread ca. 100nm), but the panchromatic channel is provided with 15m resolution spanning the near infrared to blue spectrum (spectrum spread ca. 400nm).

Pansharpening is the task of merging color imagery with hi-res panchromatic imagery of the same region, so that the resolution of the color imagery is lifted to the level of the panchromatic imagery, hence it is panchromatically sharpened. The procedure is as follows:

Let r,g,b and n be the red green blue and nir intensities of the respective imagery bands. And let p be the panchromatic intensity. In the simplest case we just perform the so-called Brovey transformation, dividing the r,g and b intensities by the total intensity i=(r+g+b)/3 and multiplying with p.

In the case of Landsat imagery, we get unnatural colors though. This is due to the fact, that the panchromatic channel not only spans the visible spectrum but also covers the near infrared spectrum, thus we need to calculate the total intensity as i=(r+g+b+n)/4. Still we get unsatisfactory results, because the bands do not contribute evenly to the panchromatic channel and may exhibit additional scale and bias. To find out the exact contribution of each channel we assume a linear relationship given as

   p = alpha*n + beta*r + gamma*g + delta*b + kappa
The constants alpha, beta, gamma, delta and kappa can be determined by multiple regression. With those, the pansharpened intensities are computed as

   r|g|b' = r|g|b * p / ( alpha*n + beta*r + gamma*g + delta*b + kappa )
The merger application applies the above principle. The five input channels pan, nir, red, green and blue are specified on the command line and one gets the pansharpened imagery as output. For convenience, the data directory contains sample imagery in the data/landsat directory. The sample channels are small crops of the original huge Landsat 7 bands of Bavaria, Germany. The original imagery is available at the Global LandCover Facility (GLCF, path=193 row=26 year=1999, FTP). In Landsat notation, the channels panchro, nir, red, green and blue correspond to the bands 80, 40, 30, 20 and 10. For the particular sample data, the multiple regression analysis of the bands reveals that

p = 0.363*n + 0.181*r + 0.062*g + 0.124*b + -0.021 The command line to perform the according pansharpening operation with the merger tool is:

   grid/merger pan.tif nir.tif red.tif green.tif blue.tif pansharpened.tif
The resulting pansharpened output image still does not look natural, because the input channels are not normalized and pretty dark. To yield bright natural colors the output needs to be processed further with a standard image manipulation tool such as gimp or imagemagick. Typically, a brightness/gamma correction together with a tonal adjustment will usually suffice.

If a true-color background image is available, the pansharpened image can be matched with the background to yield the exact same naturally looking colors. This is described in the following section.

Back to top

Color-Matching

Color-matching is the task of balancing the colors of adjacent or overlapping imagery, so that no seams arise when composed together. The matcher application in the tools directory can be applied to accomplish this task. It balances the colors of a foreground image to match the colors of a set of background images.

The tool does so by calculating a coarse grid of linear regression coefficients per-channel. The regression coefficients constitute the one linear mapping between the image and the background that exhibits the lowest linear mapping error. The tool interpolates the regression coefficients per-pixel and applies the corresponding linear transformation to each pixel. This way the color-characteristics of the background are over-imposed onto the foreground imagery.

As a welcome side effect, the color-matching tool can also be used to modify the colors of a false-color image to yield true natural looking colors. This is achieved by matching the false-color image with a [low-resolution] true-color background image.

As an example, we can create true-color Landsat images by matching the Landsat false-color pan-sharpened images (15m resolution) with a Modis true-color image (500m resolution):

   grid/matcher landsat.tif modis.tif truecolor.tif
Back to top

NDVI

The NDVI, the Normalized Difference Vegetaton Index, is a measure for bioactivity. It builds on the fact that a healthy vegetation with active photosynthesis absorbs light in the red spectrum but almost no light in the near infrared spectrum. In contrast to that, barren land and water surfaces do not exhibit this characteristical absorption and absorb light more uniform across the spectrum. The non-uniformity of the absorption, that is the NDVI, is measured by the following formula:

   NDVI = (nir-red)/(nir+red)
Typical values range from -0.1 for no bioactivity over 0.2-0.5 for average to 0.5-0.8 and higher for an increased photosynthetic activity.

The merger application (as introduced in this section) supports the calculation of the NDVI as an option. To compute the NDVI from the Landsat nir (channel 40) and red (channel 30) bands, we use the following command line:

   grid/merger nir.tif red.tif ndvi.tif
Back to top eof