Raster Operations

Some basic operations involving two rasters (or two bands) are available in TerraLib.

They include arithmetic operations, resampling, cropping, and reprojection.

The 4 basic arithmetic operations (+,-,/,*) are possible to be performed between 2 rasters (or 2 bands). When using 2 rasters, they must fit in terms of number of bands, and number of rows x columns.

A basic use of arithmetic operations between rasters is the calculation of the NDVI (Normalized Difference Vegetation Index), that follows this equation:

$$B_{NDVI} = \frac{B_{IRed} - B_{Red}}{B_{IRed} + B_{Red}}$$

where $B_{NDVI}$ is a single band derived from Red and Infra-Red bands.

// load input data
std::map<std::string, std::string> rinfo;
rinfo["URI"] = "input.tif";
te::rst::Raster* rin = te::rst::RasterFactory::open(rinfo);
 
// in this raster, band 0 is the red band, and band 1 is the infra-red band
te::rst::Band* bred = rin->getBand(0);
te::rst::Band* bired = rin->getBand(1);
 
// create output raster with 2 bands (rout[0] - ndvi, rout[1] - tmp)
std::map<std::string, std::string> orinfo;
orinfo["URI"] = "output-ndvi.tif";
te::rst::Grid* grid = new te::rst::Grid(*rin->getGrid());
std::vector<te::rst::BandProperty*> bprops;
bprops.push_back(new te::rst::BandProperty(0, te::dt::DOUBLE_TYPE, "ndvi band"));
bprops.push_back(new te::rst::BandProperty(0, te::dt::DOUBLE_TYPE, "tmp band"));
te::rst::Raster* rout = te::rst::RasterFactory::make(grid, bprops, orinfo);
 
// use rout[0] as numerator
te::rst::Band* numerator = rout->getBand(0);
te::rst::Copy(*bired, *numerator);
*numerator -= *bred;
 
// create denominator band
te::rst::Band* denominator = rout->getBand(1);
te::rst::Copy(*bired, *denominator);
*denominator += *bred;
 
// divide numerator by denominator, than rout[0] contains the ndvi
*numerator /= *denominator;
 
// clean up
delete rin;
delete rout;

TerraLib contains 3 methods for raster resampling, which creates a copy of the raster with a different spatial resolution.

The available methods in TerraLib are:

• Bicubic: it considers 16 pixels (4×4) to interpolate, creating a smooth result. Is the slowest method.
• Bilinear: this method uses the 4 nearest pixels (2×2) and interpolates using a weighted average. The weights are defined by the distance to the neighbor.
• Nearest Neighbor: this method interpolates a new pixel value by copying the nearest neighbor pixel. Is the fastest method, but can create artifacts.

// load input raster
std::map<std::string, std::string> rinfo;
rinfo["URI"] = "input.tif";
te::rst::Raster* inraster = te::rst::RasterFactory::open(rinfo);
 
// define output raster name
std::map<std::string, std::string> orinfo;
orinfo["URI"] = "output-resampling-nn.tif";
 
// applies the interpolation using neares neighbor technique
te::rst::Raster* nnraster = inraster->resample(te::rst::Interpolator::NearestNeighbor, 2, orinfo);
 
// clean up
delete inraster;
delete nnraster;

This operation allows to crop one raster, creating a sub-image, through the operation Trim.

For example, input Raster and a rectangle to crop:

Output Raster:

The following code shows how to crop one raster.

// load input raster
std::map<std::string, std::string> rinfo;
rinfo["URI"] = "input.tif";
te::rst::Raster* inraster = te::rst::RasterFactory::open(rinfo);
 
// set envelope based on line/column 
te::gm::Coord2D cur = inraster->getGrid()->gridToGeo(299, 150);
te::gm::Coord2D cll = inraster->getGrid()->gridToGeo(100, 299);
te::gm::Envelope* blenv = new te::gm::Envelope(cll.x, cll.y, cur.x, cur.y);
 
// define output raster parameters
std::map<std::string, std::string> orinfo;
orinfo["URI"] = "output-crop.tif";
 
// perform trim operation
te::rst::Raster* outraster = inraster->trim(blenv, orinfo);
 
// clean up
delete inraster;
delete outraster;

This operation reprojects a raster to another SRS (Spatial Reference System).

The following code shows how to reproject one raster with a given SRID.

// load the input raster
std::map<std::string, std::string> rinfo;
rinfo["URI"] = "input.tif";
te::rst::Raster* inraster = te::rst::RasterFactory::open(rinfo);
 
// define raster output parameters
std::map<std::string, std::string> orinfo;
orinfo["URI"] = "output-reprojection.tif";
 
// call reprojection with changing SRID to 4326 (WGS84)
te::rst::Raster* outraster = te::rst::Reproject(inraster, 4326, orinfo);
 
// clean up
delete inraster;
delete outraster;