1.. _rfc-62: 2 3======================================================================================= 4RFC 62 : Raster algebra 5======================================================================================= 6 7Author: Ari Jolma 8 9Contact: ari.jolma at gmail.com 10 11Status: Development 12 13Implementation in version: 14 15Summary 16------- 17 18It is proposed that a set of functions or methods are written for raster 19band objects to support "raster algebra", i.e., a set of operations, 20which modify bands or compute values from bands. An example of a 21modification is adding a value to all the cells of the band. An example 22of a computation is the maximum cell value in the band. Operations may 23or may not take arguments, in addition to the band itself, and if they 24take, the argument may be a numeric value, a data structure, or another 25band. Similarly, the computed value may be a simple numeric value, a 26data structure, or another band. 27 28Rationale 29--------- 30 31Raster algebra is a well known branch of geospatial science and 32technology and an often needed tool. Currently GDAL does not have 33comprehensive support for raster algebra in core. 34 35Related work 36------------ 37 38Python bindings: raster bands or parts of raster bands can be read into 39/ written from numpy array objects. Huge rasters can be iterated block 40by block. numpy methods allow efficient implementation of many raster 41algebra methods using python. There is some support for parallel 42processing using numpy. 43 44Perl bindings: raster bands or parts of raster bands can be read into / 45written from PDL objects. Huge rasters can be iterated block by block. 46PDL methods allow efficient implementation of many raster algebra 47methods using perl. There is some support for parallel processing using 48PDL. 49 50QGIS raster calculator: raster calculation string is parsed into an 51expression tree and output band is computed row by row from the inputs. 52All computations are done with double precision floating point numbers. 53The calculation does not support parallel processing. 54 55PostGIS raster: raster algebra is supported by callback functions. 56 57There is existing research, which may be exploited: 58 59Parallel Priority-Flood depression filling for trillion cell digital 60elevation models on desktops or clusters 61`http://www.sciencedirect.com/science/article/pii/S0098300416301704 <http://www.sciencedirect.com/science/article/pii/S0098300416301704>`__ 62 63Parallel Non-divergent Flow Accumulation For Trillion Cell Digital 64Elevation Models On Desktops Or Clusters 65`https://arxiv.org/abs/1608.04431 <https://arxiv.org/abs/1608.04431>`__ 66 67Requirements (Goals) 68-------------------- 69 70The implementation should be data type aware. This may mean code written 71with templates. 72 73The implementation should be parallel processing friendly. 74 75The implementation should allow a relatively easy to use C++ / C API. 76This may mean interface, which does not use templates. 77 78The implementation should allow arbitrary functions on cell values. 79I.e., be extensible by the user. 80 81The implementation should allow focal methods. I.e., methods, where the 82value of a cell depends on its neighborhood. 83 84Approach 85-------- 86 87The implementation does not need to be tightly integrated with the core. 88This means an "add-on" type solution is ok. 89 90GDAL design sets some constraints/requirements to raster algebra 91implementation: 1) the access to data is based on blocks, 2) GDAL 92supports several datatypes, even complex values, 3) there is no 93immediate support for the not-simple data structures needed by some 94methods (by "method" I refer to functions of raster algebra in this 95text), 4) data can be read from a band in parallel but writing needs to 96be exclusive. 97 98Changes 99------- 100 101Drivers 102------- 103 104Drivers are not affected. 105 106Bindings 107-------- 108 109The functionality will be added to the bindings. 110 111Utilities 112--------- 113 114Existing utilities are not affected but new utilities may be written 115taking advantage of the new functionality. 116 117Documentation 118------------- 119 120Must be written. 121 122Test Suite 123---------- 124 125Must be written. 126 127Compatibility Issues 128-------------------- 129 130Related tickets 131--------------- 132 133Implementation 134-------------- 135 136A proposed implementation is developed at 137`https://github.com/ajolma/raster_algebra <https://github.com/ajolma/raster_algebra>`__ 138 139This code attempts to solve the problem as follows. (The source is in 140transition from an old approach, which was based on operators as 141methods, while the new approach is based on operator classes) 142 143Classes 'operand' and 'operator' are defined. An operand is an object, 144which holds data and an operator is an object, which computes a result 145(essentially an operand) from operands. 146 147Raster algebra computation is a tree of operand and operator objects, 148which is executed in a recursive fashion. 149 150There are interface classes and templated concrete classes. The concrete 151classes inherit from the interface classes. 152 153Two operand classes are defined: a number and a band. There is a need 154for other types of operands. For example a classifier would map integer 155values or real number ranges into numbers. Code for such exists in the 156source but it is not organized to reflect the new approach. 157 158A central method is 'compute' in band class, which is basically the 159effective block loop code presented in the documentation for 160GDALRasterBand::ReadBlock. 161 162Multiple data types are supported by template concrete class for bands 163and by overloaded get_value method, which returns the value in required 164data type. 165 166Voting history 167-------------- 168 169