<br><br><div class="gmail_quote">On Thu, Jun 30, 2011 at 2:16 PM, Pierre Racine <span dir="ltr"><<a href="mailto:Pierre.Racine@sbf.ulaval.ca">Pierre.Racine@sbf.ulaval.ca</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<br>
(I prefer to fight with you about your inefficient modular approach to seamless operators instead ;-) But this is another discussion...)<br><br></blockquote><div><br>Take a look at the patch on <a href="http://trac.osgeo.org/postgis/ticket/1058">http://trac.osgeo.org/postgis/ticket/1058</a> .  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).<br>
<br>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.<br>
<br>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.<br>
<br>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). <br>
<br>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. <br>
<br>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).<br>
<br>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.<br>
 </div></div>