RFC 62 : Raster algebra

Author: Ari Jolma

Contact: ari.jolma at gmail.com

Status: Development

Implementation in version:

Summary

It is proposed that a set of functions or methods are written for raster band objects to support "raster algebra", i.e., a set of operations, which modify bands or compute values from bands. An example of a modification is adding a value to all the cells of the band. An example of a computation is the maximum cell value in the band. Operations may or may not take arguments, in addition to the band itself, and if they take, the argument may be a numeric value, a data structure, or another band. Similarly, the computed value may be a simple numeric value, a data structure, or another band.

Rationale

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.

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.

1) 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.

2) 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.

3) 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);

Changes

Drivers

Drivers are not affected.

Bindings

The functionality will be added to the bindings.

Utilities

A new utility will be written to take advantage of the functionality.

Documentation

Will be written.

Test Suite

Will be written.

Compatibility Issues

Related tickets

Implementation

The implementation will be done by Ari Jolma.

The proposed implementation is developed on ​https://github.com/ajolma/raster_algebra

Voting history