[postgis-devel] [raster] C Implementation map algebra

[Bryce L Nordgren]

On Thu, Jun 30, 2011 at 2:16 PM, Pierre Racine
<Pierre.Racine at sbf.ulaval.ca>wrote:

>
> (I prefer to fight with you about your inefficient modular approach to
> seamless operators instead ;-) But this is another discussion...)
>
>
Take a look at the patch on http://trac.osgeo.org/postgis/ticket/1058 .
This code is posted mostly for review by humans (you and bborie, among
others), to spark just such a discussion. :) It's main function is to give
you alternatives. It won't be usable in compiled form until we resolve the
whole "raster now depends on postgis" issue (a.k.a. boo hoo I can't link my
way out of a wet paper bag).

I essentially refactored your SQL code into two functions: one to calculate
the extent of the result; and the other to calculate values for each pixel.
Both have "extension points", and this is where modularity comes in. Note
that there is no extension point for "map algebra" at this point, but this
is described below. The iteration engine itself should never change.

Consider this a "rev 0.1" edition. The notion of how data types differ from
information types is handled clumsily. For example, a raster can represent
either a mask (spatial info only) or something with meaningful values
(spatial+value). I didn't really work that out in my head until after
writing this first draft. With the conceptual toolkit I developed in the
architecture document, I believe it is possible to generalize the extension
points such that this engine can be reused for any operation returning a
raster and taking two inputs: the inputs can be either geometries, geomvals
or rasters.

End users don't interact with this function. I'd expect that various end
user functions would use this as a back end. For instance, to implement a
mapalgebra function taking two spatial+value rasters and not taking a mask,
you'd calculate the extent with rt_raster_raster_prep(), providing the two
rasters and the spatial operation you want (union, difference, intersection,
etc.) You'd then make an empty raster to contain the result. Populate it
with values using rt_raster_raster_op_engine(), giving it extension points
relevant to the call: you'd re-use one of the spatial operations (to
determine whether a cell has data), and write an "expression evaluator"
conforming to two_raster_value_op. The only new bit is the expression
evaluator. Everything else is re-used. The program flow is the same as your
SQL variant (only one iteration over all the pixels).

To support using "masks", the only change you'd make is to NOT call
rt_raster_raster_prep(). Create a result raster using the extent, projection
and geometry of the mask. Additionally, you'd provide a spatial operation
called "isMaskTrue" (not yet written, but trivial) to the
rt_raster_raster_op_engine() call.

And note: if I really can generalize the inputs to this function such that
they can be either geometry-values or rasters, you get more than just a "two
raster map algebra". You get seamless raster and vector input handling (for
all cases which return a raster).

The important code is at the end of rt_api.c, with extension point
signatures modeled as typedefs in rt_api.h. I defined a token two-raster
intersection method as RASTER_intersection_rr2r() in rt_pg.c as an example
of how all this comes together.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-devel/attachments/20110630/a6df3909/attachment.html>