[postgis-devel] Raster todo list

[Pierre Racine]

> Well since you asked, here's how I would approach the problem:
> 
> Make ST_Intersection(raster,raster) return a mask (e.g., a one-band raster
> having only true/false values). The mask represents the common area.
> MapAlgebra uses this mask to determine the size and properties of the resultant
> raster (corner location, # pixels in x and y direction, size of pixels, etc.).
> MapAlgebra adds the ability to compute a value using a user-specified arbitrary
> expression, but only does so over the grid cells having a "true" value in the mask.
> 
> Better, MapAlgebra is not concerned with how the mask is generated, and so it
> can use ST_Difference, ST_SymDifference, ST_Union, or ST_Intersection to
> determine the spatial extent of the result with little or no change of code. Better
> yet, since all MapAlgebra needs is a mask, it can accept a user-specified mask
> which was generated by an arbitrarily complex process. Even better yet, if users
> are allowed to specify the mask, the mask can be re-used for multiple calls to
> MapAlgebra and does not need to be re-computed.
> 
> Clearly, a one-raster MapAlgebra variant is all that is needed if ST_Difference or
> ST_SymDifference is used to define the area for which results are desired, as the
> "valid area" is guaranteed to only have data from one raster (not necessarily the
> same raster).
> 
>  Using this philosophy, there is a clear and logical separation between the spatial
> operation (ST_{Difference,SymDifference,Union,Intersection}) and the value
> computation part (MapAlgebra). More modular, easier to understand, easier to
> re-use, easier to maintain.

Clearly we are speaking about almost the same thing here but you see it more simple and practical as a modular process and I see it more practical and efficient as a one complex process.

In my two rasters prototype of MapAlgebra(raster, raster) there is also this separation between the spatial operation and the value computation part. The spatial operation determine the area of the FIRST, SECOND, INTERSECTION or UNION of both raster. This is a first step toward your mask with the difference that this one is strictly geometric, represent only the extent to be iterated over and hence do not require iterating over all the pixels.

The second part of MapAlgebra(raster, raster) determine values for each pixel using the right expression. This is your mapalgebra part but more flexible because it deal with pixels values not just with a mask.

You said the one raster version of mapalgebra should accept a mask. This is fundamentally a two raster mapalgebra. The first raster is the one you want to perform calculation on and the second one is the mask. What if both raster are not aligned? There is also no guarantee that the mask will contain only true/false values. You are clearly building a more restrictive MapAlgebra(raster, raster). For me this is a waste of flexibility. 

First concern: you CAN NOT do 'rast1 * rast2' (nor 'rast1 + rast2') with a one raster mapalgebra accepting only a raster and a mask... You could do this only if all the pixels of rast2 would have the same value and then act as a constant.

Second concern: let say that you create a MapAlpgebra(raster, raster) taking into account that the second raster is not just a mask but a raster containing values that can be used in the expression (so you can resolve expression like 'rast1 + rast2') and that you do not mind about the spatial operation which are performed by a first function ST_{Difference,SymDifference,Union,Intersection}

In this case:

First, if you really want to implement 'rast1 + rast2' then what you need is actually a three raster version of MapAlgebra(rast1, rast2, mask)... What if none of them is aligned? Did you speak about simplicity?

Second: Performance: How many time do you have to iterate over all the pixels of one raster or the other to do 'rast1 + rast2 + 120' on the intersection of two rasters?

-one pass to compute the intersection mask with ST_Intersection (you might have to resample one time if the two rasters are not aligned)
-one pass to create the resulting raster with ST_MapAlgebra (you might have to resample again)

Worst: you have to iterate over ALL the pixels since the mask must have the same extent as the two rasters. 

In my prototype of mapalgebra the mask extent computation is geometric thus VERY fast (no iteration) and the value computation require only one pixel iteration over only the intersecting extent and this only once for the whole 'rast1 + rast2 + 120' expression. You can event optimize in many case when the expression resume to 'rast1' or 'rast2' by copying chunk of pixel values (this is the case for ST_Collect). This is what the optimized version is trying to do.

Clearly in the end with the modular approach we will have to write a two raster mapalgebra (or worst a three raster one) and it will perform bad because it will have to iterate over too many pixels too many times.

There is implicitly a mask in the combination of both rasters. You can determine this mask at the same time as you process pixel values. This is much more efficient than precomputing the mask and then computing values.

I join to this mail, to be executed in this order:

1) A fixed version the the two rasters version of ST_Mapalgebra (mapalgebrafixed.sql)

2) A mapalgebra implementation of all your set function (ST_{Difference,SymDifference,Union,Intersection}) and more (setOperatorsWithMapAlgebra.sql). I wrote them between 8 and 10 this morning... Actually it took me more time to fix ST_MapAlgebra following so many changes since many months than writing the necessary expressions...

3) The same set functions with everything you need to display the results in OpenJump (setOperatorsWithMapAlgebraInOpenJump.sql).

Beeing able to do MapAlpgebra(raster, raster) efficiently is fundamental to implement any raster modelisation process. It is the most fundamental operator of any respectful raster analysis package. You can easily implement any vector like operator with a flexible enough ST_MapAlpgebra(raster, ST_AsRaster(geometry)).and much more.

I won't have time to comment much more on this as I have to get ready for Paris.

Pierre
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mapalgebrafixed.sql
Type: application/octet-stream
Size: 33978 bytes
Desc: mapalgebrafixed.sql
URL: <http://lists.osgeo.org/pipermail/postgis-devel/attachments/20110617/b995a344/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: setOperatorsWithMapAlgebra.sql
Type: application/octet-stream
Size: 2122 bytes
Desc: setOperatorsWithMapAlgebra.sql
URL: <http://lists.osgeo.org/pipermail/postgis-devel/attachments/20110617/b995a344/attachment-0001.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: setOperatorsWithMapAlgebraInOpenJump.sql
Type: application/octet-stream
Size: 2975 bytes
Desc: setOperatorsWithMapAlgebraInOpenJump.sql
URL: <http://lists.osgeo.org/pipermail/postgis-devel/attachments/20110617/b995a344/attachment-0002.obj>