====== Raster Operations ======

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

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


===== Arithmetic Operations =====

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.

==== Example of use ====

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.

<code c++>
// 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;
</code>

===== Resampling =====

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 (4x4) to interpolate, creating a smooth result. Is the slowest method.
  * Bilinear: this method uses the 4 nearest pixels (2x2) 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.

==== Example of use ====

<code c++>
// 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;
</code>

===== Cropping =====

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

For example, input Raster and a rectangle to crop:

{{:wiki:designimplementation:raster:input-crop.png?200|}}

Output Raster:

{{:wiki:designimplementation:raster:output-crop.png?100|}}

==== Example of use ====

The following code shows how to crop one raster.

<code c++>
// 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;
</code>

===== Reprojection =====

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

==== Example of use ====

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

<code c++>
// 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;
</code>