Context Navigation

← Previous Change
Wiki History
Next Change →

Changes between Version 7 and Version 8 of rfc62_raster_algebra

Timestamp:: Sep 24, 2016, 12:37:32 AM (8 years ago)
Author:: Ari Jolma
Comment:: Rewrite. Add related work and requirements.

Legend:

: Unmodified
: Added
: Removed
: Modified

rfc62_raster_algebra

-              v7
+              v8
 Raster algebra is a well known branch of geospatial science and technology and an often needed tool. Currently GDAL does not have comprehensive support for raster algebra in core.
+== Related work ==
+Python bindings: raster bands or parts of raster bands can be read into / written from numpy array objects. Huge rasters can be iterated block by block. numpy methods allow efficient implementation of many raster algebra methods using python. There is some support for parallel processing using numpy.
+Perl bindings: raster bands or parts of raster bands can be read into / written from PDL objects. Huge rasters can be iterated block by block. PDL methods allow efficient implementation of many raster algebra methods using perl. There is some support for parallel processing using PDL.
+QGIS raster calculator: raster calculation string is parsed into an expression tree and output band is computed row by row from the inputs. All computations are done with double precision floating point numbers. The calculation does not support parallel processing.
+PostGIS raster: raster algebra is supported by callback functions.
+There is existing research, which may be exploited:
+Parallel Priority-Flood depression filling for trillion cell digital elevation models on desktops or clusters http://www.sciencedirect.com/science/article/pii/S0098300416301704
+Parallel Non-divergent Flow Accumulation For Trillion Cell Digital Elevation Models On Desktops Or Clusters https://arxiv.org/abs/1608.04431
+== Requirements (Goals) ==
+The implementation should be data type aware. This may mean code written with templates.
+The implementation should be parallel processing friendly.
+The implementation should allow a relatively easy to use C++ / C API. This may mean interface, which does not use templates.
+The implementation should allow arbitrary functions on cell values. I.e., be extensible by the user.
+The implementation should allow focal methods. I.e., methods, where the value of a cell depends on its neighborhood.
 == Approach ==
 GDAL presents a few problems to raster algebra implementation: 1) the access to data is based on blocks, 2) GDAL supports several datatypes, even complex values, and 3) there is no immediate support for the not-simple data structures needed by some methods (by "method" I refer to functions of raster algebra in this text). It is proposed to solve these problems as follows.
+The implementation does not need to be tightly integrated with the core. This means an "add-on" type solution is ok.
+) The code is organized into two levels: the first one, which loops through all blocks in band(s) and maintains a cache of needed blocks, and the second one, which loops through all cells in a block. The first level is generic one and the method and its arguments just pass through it. The method is passed as a callback and arguments as pointers to objects of base class. The second one is method specific, it is where the method is really implemented. The method callbacks must have similar prototype and their return value is used to, e.g., determine whether the block needs to be written to disk etc.
+) C++ templates will be used to generically support multiple cell datatypes. However, an interface class will be written, which is not templated, so the user is not distracted with the datatype too much.
+) An interface class structure is used for arguments and return values (non-band). The interface classes are not templated, they shadow actual classes, which are templated and use STL as much as possible.
+== The API ==
+A proposal is to have only a C++ interface, and also base the bindings on this C++ interface.
+{{{
+// the API:
+// The C++ API (these are all interface classes, so everything is virtual as the actual code is in hidden template classes
+class gma_object_t {
+public:
+    virtual ~gma_object_t();
+    virtual gma_class_t get_class();
+};
+class gma_number_t : public gma_object_t {
+public:
+    virtual GDALDataType datatype();
+    virtual void set_value(double value);
+    virtual void set_value(int value);
+    virtual int value_as_int();
+    virtual double value_as_double();
+    virtual bool is_defined();
+    virtual void set_inf(int inf); // -1 to minus inf, 0 to not inf, 1 to plus inf
+    virtual bool is_inf();
+    virtual bool is_integer();
+    virtual bool is_float();
+};
+class gma_pair_t : public gma_object_t {
+public:
+    virtual void set_first(gma_object_t *first);
+    virtual void set_second(gma_object_t *second);
+    virtual gma_object_t *first();
+    virtual gma_object_t *second();
+};
+class gma_bins_t : public gma_object_t {
+public:
+    virtual GDALDataType datatype();
+    virtual unsigned int size();
+    virtual void push(int value);
+    virtual void push(double value);
+};
+class gma_histogram_t : public gma_object_t {
+public:
+    virtual GDALDataType datatype();
+    virtual unsigned int size();
+    virtual gma_object_t *at(unsigned int i);
+    virtual void print();
+};
+class gma_classifier_t : public gma_object_t {
+public:
+    virtual GDALDataType datatype();
+    virtual unsigned int size();
+    virtual void add_class(gma_number_t *interval_max, gma_number_t *value);
+    virtual void add_value(gma_number_t *old_value, gma_number_t *new_value);
+    virtual void add_default(gma_number_t *default_value);
+};
+class gma_cell_t  : public gma_object_t {
+public:
+    virtual GDALDataType datatype();
+    virtual int& x();
+    virtual int& y();
+    virtual void set_value(double value);
+    virtual void set_value(int value);
+    virtual int value_as_int();
+    virtual double value_as_double();
+};
+/*
+  Return value 0 interrupts, 1 denotes ok, and 2 denotes ok and a need
+  to save the cell value back to band.
+*/
+typedef int (*gma_cell_callback_f)(gma_cell_t*, gma_object_t*);
+class gma_cell_callback_t : public gma_object_t {
+public:
+    virtual void set_callback(gma_cell_callback_f callback);
+    virtual void set_user_data(gma_object_t*);
+};
+class gma_logical_operation_t  : public gma_object_t {
+public:
+    virtual GDALDataType datatype();
+    virtual void set_operation(gma_operator_t);
+    virtual void set_value(int value);
+    virtual void set_value(double value);
+};
+class gma_hash_t : public gma_object_t {
+public:
+    virtual GDALDataType datatype();
+    virtual int size();
+    virtual std::vector<gma_number_t*> *keys_sorted();
+    virtual gma_object_t *get(gma_number_t *key);
+};
+class gma_band_t : public gma_object_t {
+public:
+    virtual void update();
+    virtual GDALRasterBand *band();
+    virtual GDALDataset *dataset();
+    virtual GDALDriver *driver();
+    virtual GDALDataType datatype();
+    virtual bool datatype_is_integer();
+    virtual bool datatype_is_float();
+    virtual int w();
+    virtual int h();
+    virtual void set_progress_fct(GDALProgressFunc progress, void * progress_arg);
+    virtual gma_band_t *new_band(const char *name, GDALDataType datatype);
+    virtual gma_number_t *new_number();
+    virtual gma_number_t *new_int(int value);
+    virtual gma_pair_t *new_pair();
+    virtual gma_pair_t *new_range();
+    virtual gma_bins_t *new_bins();
+    virtual gma_cell_t *new_cell();
+    virtual gma_classifier_t *new_classifier();
+    virtual gma_cell_callback_t *new_cell_callback();
+    virtual gma_logical_operation_t *new_logical_operation();
+    virtual void print();
+    virtual void rand();
+    virtual void abs();
+    virtual void exp();
+    virtual void log();
+    virtual void log10();
+    virtual void sqrt();
+    virtual void sin();
+    virtual void cos();
+    virtual void tan();
+    virtual void ceil();
+    virtual void floor();
+    virtual void assign(int value);
+    virtual void assign_all(int value);
+    virtual void add(int summand);
+    virtual void subtract(int);
+    virtual void multiply(int);
+    virtual void divide(int);
+    virtual void modulus(int divisor);
+    virtual void assign(double value);
+    virtual void assign_all(double value);
+    virtual void add(double summand);
+    virtual void subtract(double);
+    virtual void multiply(double);
+    virtual void divide(double);
+    virtual void classify(gma_classifier_t*);
+    virtual void cell_callback(gma_cell_callback_t*);
+    // arg = NULL, pair:(n,pair:(min,max)), or bins; returns histogram
+    virtual gma_histogram_t *histogram(gma_object_t *arg = NULL);
+    // returns hash of a hashes, keys are zone numbers
+    virtual gma_hash_t *zonal_neighbors();
+    virtual gma_number_t *get_min();
+    virtual gma_number_t *get_max();
+    // returns a pair of numbers
+    virtual gma_pair_t *get_range();
+    virtual std::vector<gma_cell_t*> *cells();
+    // op can be used to make the operation conditional,
+    // the test is made against the value of the parameter band
+    virtual void assign(gma_band_t *, gma_logical_operation_t *op = NULL);
+    virtual void add(gma_band_t *, gma_logical_operation_t *op = NULL);
+    virtual void subtract(gma_band_t *, gma_logical_operation_t *op = NULL);
+    virtual void multiply(gma_band_t *, gma_logical_operation_t *op = NULL);
+    virtual void divide(gma_band_t *, gma_logical_operation_t *op = NULL);
+    virtual void modulus(gma_band_t *, gma_logical_operation_t *op = NULL);
+    // this = value where decision is true
+    // the decision band must be of type uint8_t
+    virtual void decision(gma_band_t *value, gma_band_t *decision);
+    virtual gma_hash_t *zonal_min(gma_band_t *zones);
+    virtual gma_hash_t *zonal_max(gma_band_t *zones);
+    virtual void rim_by8(gma_band_t *areas);
+    virtual void fill_depressions(gma_band_t *dem);
+    virtual void D8(gma_band_t *dem);
+    virtual void route_flats(gma_band_t *dem);
+    virtual void upstream_area(gma_band_t *);
+    virtual void catchment(gma_band_t *, gma_cell_t *);
+};
+// a band object factory
+gma_band_t *gma_new_band(GDALRasterBand *b);
+}}}
+GDAL design sets some constraints/requirements to raster algebra implementation: 1) the access to data is based on blocks, 2) GDAL supports several datatypes, even complex values, 3) there is no immediate support for the not-simple data structures needed by some methods (by "method" I refer to functions of raster algebra in this text), 4) data can be read from a band in parallel but writing needs to be exclusive.
 == Changes ==
 …
 == Utilities ==
 A new utility will be written to take advantage of the functionality.
+Existing utilities are not affected but new utilities may be written taking advantage of the new functionality.
 == Documentation ==
 Will be written.
+Must be written.
 == Test Suite ==
 Will be written.
+Must be written.
 == Compatibility Issues ==
 …
 == Implementation ==
+The implementation will be done by Ari Jolma.
+A proposed implementation is developed at https://github.com/ajolma/raster_algebra
+The proposed implementation is developed on https://github.com/ajolma/raster_algebra
+This code attempts to solve the problem as follows. (The source is in transition from an old approach, which was based on operators as methods, while the new approach is based on operator classes)
+Classes 'operand' and 'operator' are defined. An operand is an object, which holds data and an operator is an object, which computes a result (essentially an operand) from operands.
+Raster algebra computation is a tree of operand and operator objects, which is executed in a recursive fashion.
+There are interface classes and templated concrete classes. The concrete classes inherit from the interface classes.
+Two operand classes are defined: a number and a band. There is a need for other types of operands. For example a classifier would map integer values or real number ranges into numbers. Code for such exists in the source but it is not organized to reflect the new approach.
+A central method is 'compute' in band class, which is basically the effective block loop code presented in the documentation for GDALRasterBand::ReadBlock.
+Multiple data types are supported by template concrete class for bands and by overloaded get_value method, which returns the value in required data type.
 == Voting history