[raster] n-raster ST_MapAlgebra

[Bborie Park]

This is to replace ALL ST_MapAlgebra functions, 1 and 2 raster variants. The question I need to address is whether or not passing a single raster multiple times to the C side is efficient (multiple points to the same memory address). If the answer is yes, the supporting multiple bands from one raster is easy...

[pracine]

Hopefully this will not affect the performance of the two raster ST_MapAlgebra()...

[Bborie Park]

The performance should be about the same. I don't plan on replacing the two raster variant immediately as the n-raster variant will require a different signature.

[Bborie Park]

Another question that needs answering is:

Should Map Algebra support function-based, expression-based or both?

[Bborie Park]

Core signature for the n-raster ST_MapAlgebra...

With 1 or more rasters being passed to the function, a new composite type is needed. This type should be usable for many future applications as it is just denoting a very fundamental idea: Band b of Raster r.

CREATE TYPE rastbandarg AS (
	rast raster,
	nband integer
);

The main function calling the C back-end is a non-user function.

CREATE OR REPLACE FUNCTION _st_mapalgebra(
	rastbandargset rastbandarg[],
	callbackfunc regprocedure,
	pixeltype text DEFAULT NULL,
	distancex integer DEFAULT 0, distancey integer DEFAULT 0,
	extenttype text DEFAULT 'INTERSECTION',
	customextent raster DEFAULT NULL,
	VARIADIC userargs text[] DEFAULT NULL
)
	RETURNS boolean
	AS 'MODULE_PATHNAME', 'RASTER_nMapAlgebra'
	LANGUAGE 'c' STABLE;

Still haven't thought of actual user function signatures but some actual use cases that come to mind are below.

Do neighborhood operation on a tile taking into account neighboring tiles. The output raster would have the extent of the input custom extent, which is the tile of interest.

SELECT ST_MapAlgebra(
	ARRAY[ROW(ST_Union(f2.rast), 1)]::rastbandarg[],
	'ngbfunc(float[][][], boolean[][][], int[][][], text[])'::regprocedure,
	NULL,
	1, 1,
	'CUSTOM',
	f1.rast
)
FROM foo f1
JOIN foo f2
	ON ST_Intersects(f1.rast, f2.rast)
WHERE f1.rid = 1

A multi-band operation, such as for computing EVI or NDVI.

SELECT ST_MapAlgebra(
	ARRAY[ROW(rast, 1), ROW(rast, 4)]::rastbandarg[],
	'evi(float[][][], boolean[][][], int[][][], text[])'::regprocedure
)
FROM foo

If you noticed above, the callback function is different than in PostGIS 2.0. A correct callback should look like...

callbackfunc(value float[][][], nodata boolean[][][], pos int[][], variadic userargs text[])

The parameters "value" and "nodata" are now 3D arrays:

value[RASTER_NUM][ROW_Y][COLUMN_X] = pixel value

nodata[RASTER_NUM][ROW_Y][COLUMN_X] = pixel is NODATA (TRUE or FALSE)

where RASTER_NUM, ROW_Y, COLUMN_X are 1-based, starting from 1.

pos[][] is a 2D array with 2 elements in the second dimension of the X, Y position of the source pixel of interest.

pos[RASTER_NUM] = [COLUMN_X, ROW_Y]

[pracine]

Some questions:

What are distancex, distancey and customextent?

Also: You say "This is to replace ALL ST_MapAlgebra functions". You mean only ST_MapAlgebraFct ones or ST_MapAlgebraExpr ones as well?

[Bborie Park]

distancex and distancey are for passing neighborhoods to the callback, not just the pixel value of interest.

customextent is so that the user can specify an extent that they want. It is in essence a filter and as a side-affect permits correct coverage handling.

I'm not a fan of the expression function variants as they have way too many issues with parameter substitution and the most obvious place for user error. So, my statement "This is to replace ALL ST_MapAlgebra functions" should be amended to "This is to replace ALL ST_MapAlgebra functions using callback functions". If anything, I'm for leaving the existing expression functions as is with no new additions to them to provide multi-raster/multi-band support.

[pracine]

Replying to dustymugs:

distancex and distancey are for passing neighborhoods to the callback, not just the pixel value of interest.

So you will replace ST_MapAlgebraFct (1 band version), ST_MapAlgebraFct (2 band version) and ST_MapAlgebraFctNgb (1 band version). Right?

Are the existing signatures going to be mapped to the new function or we will have to rewrite all our custom functions?

customextent is so that the user can specify an extent that they want. It is in essence a filter and as a side-affect permits correct coverage handling.

Can we call it a "mask"? In essence your multiband MapAlgebra is really what it is, a multiband map algebra. But I would add "intra-raster" multiband meaning the main goal is to do a computation involving many bands generally from the same raster. I can not think, and you did not show an application, where, say, three rasters would come from three different raster tables involving two joins and a multitude of possible resulting extents. My point here is to ask: Why not making a simpler multiband mapalgebra working only on the bands of a same raster? No need to pass many rasters (just an array of band index) and hence no need to create a new type. A much simpler function signature. No hassle with the possible extent since they are always identical. No braking of the existing mapalgebra (which work very well right now)... I generally tend not to complexify things if I don't have a useful user case. And I don't think it would be efficient and wise to join three raster coverages.

[robe]

hmm -- what about the case where you want to overlay a geometry or perhaps more than one on a raster --- e.g. in doc example: ​http://www.postgis.org/documentation/manual-svn/RT_ST_MapAlgebraExpr2.html

You could end up with a multitude of rasters not coming from the same raster.

[pracine]

Replying to robe:

hmm -- what about the case where you want to overlay a geometry or perhaps more than one on a raster --- e.g. in doc example: ​http://www.postgis.org/documentation/manual-svn/RT_ST_MapAlgebraExpr2.html

Which example are you refering to? Are three tables involved?

[Bborie Park]

Replying to pracine:

Replying to dustymugs:

distancex and distancey are for passing neighborhoods to the callback, not just the pixel value of interest.

So you will replace ST_MapAlgebraFct (1 band version), ST_MapAlgebraFct (2 band version) and ST_MapAlgebraFctNgb (1 band version). Right?

Correct

Are the existing signatures going to be mapped to the new function or we will have to rewrite all our custom functions?

Existing signatures will still stay the same. The custom functions will need to be rewritten though.

customextent is so that the user can specify an extent that they want. It is in essence a filter and as a side-affect permits correct coverage handling.

Can we call it a "mask"? In essence your multiband MapAlgebra is really what it is, a multiband map algebra. But I would add "intra-raster" multiband meaning the main goal is to do a computation involving many bands generally from the same raster. I can not think, and you did not show an application, where, say, three rasters would come from three different raster tables involving two joins and a multitude of possible resulting extents. My point here is to ask: Why not making a simpler multiband mapalgebra working only on the bands of a same raster? No need to pass many rasters (just an array of band index) and hence no need to create a new type. A much simpler function signature. No hassle with the possible extent since they are always identical. No braking of the existing mapalgebra (which work very well right now)... I generally tend not to complexify things if I don't have a useful user case. And I don't think it would be efficient and wise to join three raster coverages.

Multiband (intra-raster) does work with my current signature. I've tested the memory usage of passing the same raster twice in the same function call. Both arguments point to the same memory space. Adding another function signature to do multiple bands of one raster would just be a wrapper around the one signature.

The reality is that none of us know all possible use-cases. The core iterator function (rt_raster_iterator()) is written to be generic and be as flexible as the can be. Hence why I'm implementing the SQL functions to maintain that flexibility.

[Bborie Park]

A change to the callback signature. Completely forgot about the magical NULL in SQL...

double precision allbackfunc(value double precision[][][], pos int[][], variadic userargs text[])

The nodata matrix is not needed due to the NULL. In C, a separate matrix is needed.

[Bborie Park]

I should elaborate on the custom functions. My statement "The custom functions will need to be rewritten though.". What I mean is that the signature can stay the same. PostgreSQL treats matrix[][] the same as matrix[] or matrix[][][][]. Bracket sets beyond the first set are more for user understanding but have no real meaning beyond that. A matrix defined as matrix[] can have one or more dimensions.

What will change with the custom functions are checks to require three dimensions and how to handle those three dimensions.

[robe]

Replying to pracine:

Replying to robe:

hmm -- what about the case where you want to overlay a geometry or perhaps more than one on a raster --- e.g. in doc example: ​http://www.postgis.org/documentation/manual-svn/RT_ST_MapAlgebraExpr2.html

Which example are you refering to? Are three tables involved?

The last one with the parcel boundaries. I may very well have a temperature raster and elevation too I might want to do things with that I got from different sources.

[Bborie Park]

You should be able to whatever you want with the in-development signature.

[Bborie Park]

Added in r10392