# Raster Algebra¶

Django-raster has raster calculator functionality. The raster calculator allows rendering raster tiles based on algebraic formulas. The use is very similar to a standard z/x/y tile endpoint, but allows the evaluation of a broad range of algebraic expressions applied to existing pixel values. The z/x/y structure can be used directly in online mapping software such as OpenLayers or Leaflet.

Similar to the regular tiles endpoint, the django-raster url patterns need to be installed for the raster algebra endpoint to work. For the documentation we assume that the /raster/ base url is used as proposed in the Installation section.

## Raster algebra TMS endpoint¶

The raster algebra url base is used only to specify the z/x/y tile index. All the rest of the configuration is done through the query parameters. The input to the raster algebra is a named list of RasterLayer ids and a formula for evaluation. These values are passed to the backend through two required query parameters: layers and formula.

The layers query parameter identifies which raster layers to use for evaluation. It is a comma separated list of variable-name and RasterLayer id pairs. The variable names are matched with the names in the formula. An example is layers=a=2,b=4 which will match RasterLayer with id 2 to variable name a and the layer with id 4 with the variable name b.

The formula query parameter is a string specifying a formula for evaluation. The formula is an algebraic expression based on the names given to the layers in the layers query parameter. The formula has to be an expression that can be evaluated by the FormulaParser. It accepts a broad range of algebraic expressions. The endpoint supports most of the common mathematical operators (+, -, *, /, etc), functions (sin, cos, exp, etc.), and logical operators (&, !, >, =, etc.). It also has a set of predefined constants through reserved keywords such as pi PI or the Euler number E.

Putting it all together, an example request to the raster algebra endpoint could look like this:

/raster/algebra/{z}/{x}/{y}.png?layers=a=1,b=3,c=6&formula=log(a+b)*c&legend=5


In addition to the required query parameters: layers and formula, a Legend id can be specified using the legend query parameter. If specified, the legend will be used to interpret the result of the algebra expression. This is convenient to use predefined colormaps for the endpoint.

### Dynamic colormap¶

For a more dynamic rendering scheme, a dynamic colormap can be passed to the endpoint using the colormap query parameter. The following request would color all pixels that result in a value bigger than zero in red, and all other pixels in green.

/raster/algebra/{z}/{x}/{y}.png?layers=a=1,b=3,c=6&formula=log(a+b)*c&colormap={'x>0':'#FF0000','x<=0':'#00FF00'}


### Using specific bands¶

By default, the algebra and rgb endpoints use the first band in each layer specified. To use a specific band, use a 'variable:band' syntax, where variable is the name of the variable, and band is the band index. For example {'a:3': 23} would match band 3 of the RasterLayer with the id 23 to the variable name a.

### Encoding¶

Both the colormap and the formula should be properly url encoded. The examples here are not encoded and should be considered as instructive examples only.

### RGB endpoint¶

The algebra endpoint can also be used to render RGB images. For this, only three query parameters are expected: r, g, and b. If these three parameters are found in the list of query parameters, and no formula has been specified, the three input bands are interpreted as RGB channels of an RGB image. For example to use raster layer with id 1 as red, id 3 as green and id 6 as blue, the following url can be used:

/raster/algebra/{z}/{x}/{y}.png?layers=r=1,g=3,b=6


If the raw data in the tiles is not already scaled to the range [0, 255], an additional scaling factor can be specified, which will be used to rescale all three bands to the default RGB color range. For instance, the following query would assume that the input bands have values in the range of [5, 10000], and would rescale them to the RGB color space.

/raster/algebra/{z}/{x}/{y}.png?layers=r=1,g=3,b=6&scale=5,10000


An alpha channel can be activated by passing the alpha query parameter. The alpha parameter makes all the pixels transparent that have values equal to 0 in all three RGB channels.

For multi band rasters that have the rgb channels as bands and not in separate files, the band accessor syntax can be used. For instance, if the layer with id 23 is a 3-band RGB raster, the following would render the tiles as RGB using bands 0, 1, and 2:

/raster/algebra/{z}/{x}/{y}.png?layers=r:0=23,g:1=23,b:2=23


### Image Enhancement¶

The algebra and TMS endpoints support image enhancement using the ImageEnhance PIL module. The following query parameters arguments are passed to the corresponding image enhancers. The parameter value is passed to the enhancer as factor argument.

Enhancer query parameters.
Query Enhancer
enhance_color ImageEnhance.Color
enhance_contrast ImageEnhance.Contrast
enhance_brightness ImageEnhance.Brightness
enhance_sharpness ImageEnhance.Sharpness

The following example enhances the contrast of tiles from the RGB endpoint by a factor of 3:

/raster/algebra/{z}/{x}/{y}.png?layers=r=1,g=3,b=6&scale=5,10000&enhance_contrast=3


## Pixel Value Lookup¶

Single pixel values for raster algebra expressions can be looked up by coordinates. The endpoint works very similar to the raster algebra TMS endpoint, but instead of Z-X-Y tile indices, coordinates are passed through the url. The query parameters are analogue to the algebra TMS endpoint as described above.

The base url structure is

/raster/pixel/{xcoord}/{ycoord}/


For instance, the following request will return the pixel value of the requested raster algebra expressino for the coordinates xcoord = -9218229 and ycoord = 3229269. The coordinates must be provided in the web mercator projection (EPSG 3857).

/raster/pixel/-9218229/3229269/?layers=a=1,b=3,c=6&formula=log(a+b)*c


## Formula parser¶

At the heart of the raster calculator is the FormulaParser, which is based on the pyparsing package. The FormulaParser is a general purpose formula evaluation class. It is It does not know about rasters and operates with Numpy arrays directly. To use it, you need a dictionary with Numpy arrays of equal shape and a formula as string. The keys in the dictionary are the variable names and are used to match data to variables in the formula. Here are some examples of how to use the formula parser:

# Import parser and instantiate an instance.
>>> from raster.algebra.parser import FormulaParser
>>> parser = FormulaParser()
# Create a data dictionary and evaluate a simple sum.
>>> data = {'a': range(5), 'b': range(5)}
>>> formula = 'a + b'
>>> parser.evaluate(data, formula)
... array([0, 2, 4, 6, 8])
# Use the sin function and divide by b.
>>> formula = 'sin(a) / b'
>>> parser.evaluate(data, formula)
... array([ nan, 0.84147098, 0.45464871, 0.04704, -0.18920062])
# Use a logical array.
>>> data.update({'a_new_var': [True, False, False, True, False]})
>>> formula = '!a_new_var * a + 3'
>>> parser.evaluate(data, formula)
... array([ 3.,  4.,  5.,  3.,  7.])
# Use the PI keyword in a formula.
>>> formula = 'a * PI'
>>> parser.evaluate(data, formula)
>>> array([0. , 3.14159265, 6.28318531, 9.42477796, 12.56637061])


## Raster algebra parser¶

The RasterAlgebraParser class is a wrapper that can be used to apply the generic formula parser to raster objects directly. The use is identical to the generic case except that the objects in the data dictionary are expected to be :class:GDALRaster objects. The data arrays are extracted from the raster objects automatically and are passed to the formula parser. The result array is converted into a GDALRaster before returning.

By default, the first band is used for calculations, to specify a specific band to be used the syntax is 'variable:band', where variable is the name of the variable, and band is the band index. For example {'a:3': rst} would match band 3 of the GDALRaster rst to the variable name a.

Here is a complete example for how to use the RasterAlgebraParser.

>>> from raster.algebra.parser import RasterAlgebraParser
>>> parser = RasterAlgebraParser()
>>> base = {
>>>     'datatype': 1,
>>>     'driver': 'MEM',
>>>     'width': 2,
>>>     'height': 2,
>>>     'srid': 3086,
>>>     'origin': (500000, 400000),
>>>     'scale': (100, -100),
>>>     'skew': (0, 0),
>>>     'bands': [
>>>         {'nodata_value': 10},
>>>         {'nodata_value': 10},
>>>         {'nodata_value': 10},
>>>     ],
>>> }
>>> base['bands'][0]['data'] = range(20, 24)
>>> base['bands'][1]['data'] = range(10, 14)
>>> rast1 = GDALRaster(base)
>>> base['bands'][0]['data'] = [1, 1, 1, 1]
>>> rast2 = GDALRaster(base)
>>> base['bands'][0]['data'] = [30, 31, 32, 33]
>>> base['bands'][0]['nodata_value'] = 31
>>> rast3 = GDALRaster(base)
>>> data = dict(zip(['x:1', 'y:0', 'z'], [rast1, rast2, rast3]))
>>> rst = parser.evaluate_raster_algebra('x*(x>11) + 2*y + 3*z*(z==30)')
>>> rst.bands[0].data()
... array([[ 10.,  10.],
...        [ 14.,  15.]])


## Keywords, Operators and Functions¶

The following tables list the available operators, functions and reserved keywords from the FormulaParser and the corresponding raster calculator.

Keyword symbols
Keyword Symbol
Euler Number E
Pi PI
True Boolean TRUE
False Boolean FALSE
Null NULL
Infinite INF
Operator symbols
Operator Symbol
Add +
Substract -
Multiply *
Divide /
Power ^
Equal ==
Not Equal !=
Greater >
Greater or Equal >=
Less <
Less or Equal <=
Logical Or |
Logial And &
Logcal Not !
Fill Nodata Values ~
Unary And +
Unary Minus -
Unary Not !
Function symbols
Function Symbol
Sinus sin
Cosinus cos
Tangens tan
Natural Logarithm log
Exponential Function exp
Absolute Value abs
Integer int
Round round
Sign sign
Minimum min
Maximum max
Mean mean
Median median
Standard Deviation std
Sum sum