| 77 | | == '''Objective FV.03 - Implement all the necessary versions of ST_MapAlgebra''' == |
| 78 | | |
| 79 | | Different versions of ST_MapAlgebra are still to do: |
| 80 | | |
| 81 | | One raster versions: |
| 82 | | |
| 83 | | 4) '''ST_MapAlgebraFctNgb('''rasttable text, rastcolumn text, band int, pixeltype text, radius int, funcname text[, funcargs text]''') -''' A one raster version taking a '''table name and a raster column name''' (in order to work on a '''tiled coverage''') and a '''user defined function''' (with optional parameters) of the set of first, second, etc... '''neighbours''' of a pixel. The passed matrix should include values for out of bound pixels taken from neighbour tiles. |
| 84 | | Some versions were done: |
| 85 | | |
| 86 | | Two rasters versions: |
| 87 | | |
| 88 | | 6) '''ST_MapAlgebraExpr('''rast1table text, rast1column text, band1 integer, rast2table text, rast2column text, band2 integer, pixeltype text, extentexpr text, expression text, nodata1expr text, nodata2expr text, nodatanodatadaexpr text''') -''' A two rasters version taking '''two table names and two raster column name''' (in order to work on a '''tiled coverage''') and an expression of two pixels at a time, one from each rasters. |
| 89 | | |
| 90 | | 7) '''ST_MapAlgebraFct('''rast1 raster, band1 integer, rast2 raster, band2 integer, pixeltype text, extentexpr text, funcname text[, funcargs text]''') -''' A two rasters version taking a '''user defined function''' (with optional parameters) of two pixels at a time, one from each rasters. The function is a user defined PL/pgSQL function taking two float8 values and returning one value. |
| 91 | | |
| 92 | | 8) '''ST_MapAlgebraFct('''rast1table text, rast1column text, band1 integer, rast2table text, rast2column text, band2 integer, pixeltype text, extentexpr text, funcname text[, funcargs text]''') -''' A two rasters version taking a '''table name and a raster column name''' (in order to work on a '''tiled coverage''') and a '''user defined function''' (with optional parameters) of two pixels at a time, one from each rasters. The function is a user defined PL/pgSQL function taking two float8 values and returning one value. Non-existent values are found in the neighbour tiles when possible. |
| 93 | | |
| 94 | | One raster versions: |
| 95 | | |
| 96 | | 1) '''ST_MapAlgebraExpr('''rast raster, band int, pixeltype text, expression text, nodataexpr text''') -''' A one raster version taking an '''expression''' of one pixel at a time. Already implemented... |
| 97 | | |
| 98 | | 2) '''ST_MapAlgebraFct('''rast raster, band int, pixeltype text, funcname text[, funcargs text]''') -''' A one raster version taking a '''user defined function''' (with optional parameters) of one pixel at a time. The function is a user defined PL/pgSQL function taking a float8 value and returning a value. [http://trac.osgeo.org/postgis/ticket/860 Code was developped by David Zwarg], needs to be integrated in the source tree. This version is much faster than 1) but requires the user to write a PL/pgSQL function. |
| 99 | | |
| 100 | | 3) '''ST_MapAlgebraFctNgb('''rast raster, band int, pixeltype text, radius int, funcname text[, funcargs text]''') -''' A one raster version taking a '''user defined function''' (with optional parameters) of the set of first, second, etc... '''neighbours''' of a pixel. The function is a user defined PL/pgSQL function taking a matrix containing the neighbour values and returning one value. Code do not exist yet but will be very much similar to 2). Out of bound pixels values are set to NULL. This version requires the user to write a PL/pgSQL function. Many predefined function should be delivered. |
| 101 | | |
| 102 | | Two rasters versions. These versions must take into account different alignment, different extents, nodata and non-existent values. |
| 103 | | |
| 104 | | 5) '''ST_MapAlgebraExpr('''rast1 raster, band1 integer, rast2 raster, band2 integer, pixeltype text, extentexpr text, expression text, nodata1expr text, nodata2expr text, nodatanodatadaexpr text''') -''' A two rasters version taking an expression of two pixels at a time, one from each rasters. |
| 105 | | |
| 106 | | Parameters and details for each variant follow... |
| 107 | | |
| 108 | | '''Details for 1) ST_MapAlgebraExpr(raster rast, integer band, text expression, text nodatavalueexpr, text pixeltype) -> raster''' |
| 109 | | |
| 110 | | This function is already implemented. See the [http://trac.osgeo.org/postgis/wiki/WKTRaster/SpecificationWorking02 specifications of Objective 2.0.02] and the [http://postgis.refractions.net/documentation/manual-svn/RT_ST_MapAlgebra.html documentation]. |
| 111 | | |
| 112 | | |
| 113 | | |
| 114 | | '''Details for 3) ST_MapAlgebraFctNgb()''' |
| 115 | | |
| 116 | | For now ST_MapAlgebra expressions refer only to the pixel being computed. e.g. "rast * 100". The original plan was to allow refering to neighbour pixels using two coordinated relative to the pixel being computed. e.g. "rast[-1,0] * 100" where rast[-1,0] refer to the value of the pixel one pixel to the left of the pixel being computed. However this syntax might prove to be hard to use when many neighbours are to be used. |
| 117 | | |
| 118 | | An alternative syntax would involve another function name (e.g. ST_MapAlgebraNgb or ST_MovingWindow) and a way to define a neighbour rectangular region around the computed pixel (e.g.: "2,2" meaning a rectangle encompassing the two neighbour pixels in each direction) and a function to call with this matrix of pixel values. A complete example might look like: |
| 119 | | |
| 120 | | SELECT ST_MapAlgebraFctNgb(rast, band, pixeltype, "ST_Mean", 2, 2, "ignore") |
| 121 | | |
| 122 | | So this would mean "for each pixel, compute the average of all the 1 + 8 + 16 = 25 pixels surrounding the current pixel and "ignore" pixels with nodata values." |
| 123 | | |
| 124 | | The "ST_Mean" summarizing function should accept three parameters: an array of float8 values, a X and a Y dimension, and optionnally a "what to do with nodata values". The possible value for this last parameter could be: |
| 125 | | |
| 126 | | -"NULL": If any value is a nodata value, return NULL.[[BR]] |
| 127 | | -"ignore": Ignore any nodata value so that if 4 pixels on 25 are nodata values, do the computation with the 21 remaining. This is the default.[[BR]] |
| 128 | | -"value": Replace any nodata value with the value of the pixel being computed.[[BR]] |
| 129 | | -a value: Replace any nodata value with this value and compute. |
| 130 | | |
| 131 | | '''Open questions''' |
| 132 | | |
| 133 | | DZ: Since this accepts a user function, the user function should determine how to handle NODATA values inside of the neighborhood. These different modes of operating should be arguments to the userfunctions, and not to ST_MapAlgebraFctNgb. |
| 134 | | |
| 135 | | Pierre: They ARE arguments to the user function as are 2 and 2 in this example. But all user function MUST accept at least these three parameters. |
| 136 | | |
| 137 | | DZ: The user function can determine the dimensions of the neighborhood from the incoming 2-dimensional array, so the user function must accept: neighborhood float[][], nodataflag text, variadic args text[]. It is not necessary to pass the neighborhood dimensions to the user function. They are still required in the main ST_MapAlgebraFctNgb function, as it must construct the neighborhood array. |
| 138 | | |
| 139 | | Any remaining parameters to ST_MapAlgebraNgb could be passed to the summarizing functions for its own need (e.g. "round" to specify that only the pixel forming a circle should be used in the computing). |
| 140 | | |
| 141 | | A number of predefined summarizing function could be delivered: ST_Max, ST_Min, ST_Sum, ST_Distinct, ST_Mean, ST_STD, ST_Range, ST_Quantile, ST_Median, ST_Majority, ST_Minority, ST_Slope, ST_Aspect, and more... |
| 142 | | |
| 143 | | Users could write their own map algebra summarizing functions. |
| 144 | | |
| 145 | | A more sophisticated version would pass a georeferenced raster instead of just a value matrix so that summarizing function could use this geoereference (e.g. to determine a value from the whole coverage (with ST_Value) when the neighbours are out of the bound of the raster). Passing a raster would allow existing raster functions (like the summarizing function which are to come). Only the optional "what to do with nodata values" could be needed and some additional parameters. In this case the example become: |
| 146 | | |
| 147 | | SELECT ST_MapAlgebraNgb(rast, band, pixeltype, 2, 2, "ST_Mean", "ignore") |
| 148 | | |
| 149 | | and the dimensions do not have to be passed to the summarizing functions since it could deduce them from ST_Width & ST_Height. |
| 150 | | |
| 151 | | An even more sophisticated version should get a raster table and a raster column as parameters and try to search for neighbour in the whole raster coverage when out of bound pixels are part of the neighbourhood. e.g.: |
| 152 | | |
| 153 | | SELECT ST_MapAlgebraNgb("mycoveragetable", "myrastercolumn", band, pixeltype, 2, 2, "ST_Mean", "ignore") |
| 154 | | |
| 155 | | Three difficulties must be solved to implement this function: |
| 156 | | |
| 157 | | -The construction of the matrix must to be passed to the summarizing functions must be optimized when passing from one pixel to the other.[[BR]] |
| 158 | | -We must see how it is possible to call a PL/pgSQL function from a C function[[BR]] |
| 159 | | -We must see how to pass a variable number of parameter to a PL/pgSQL function |
| | 79 | == '''Objective FV.24 - ST_MapAlgebraFctNbg() working on a tiled coverage''' == |
| | 80 | |
| | 81 | |
| | 82 | '''ST_MapAlgebraFctNgb('''rasttable text, rastcolumn text, band int, pixeltype text, radius int, funcname text[, funcargs text]''') -''' |
| | 83 | |
| | 84 | A one raster version taking a '''table name and a raster column name''' (in order to work on a '''tiled coverage''') and a '''user defined function''' (with optional parameters) of the set of first, second, etc... '''neighbours''' of a pixel. This version work transparently on a tiled coverage avoiding edges effects. |
| | 85 | |
| | 86 | The function could 1) search itself for out of bound pixel values or 2) pass the name of the table to search for those values and let this job to the custom function. (I think 1) is a better option since the actual predefined custom functions would work tranparently). |
| | 87 | |
| | 88 | Open question: What to do when there are more than one pixel (because they are overlapping) value neighbouring the current pixel? Could be another parameter saying take the 'MIN', the 'MAX', the 'FIRST', the 'LAST', the 'MEAN'... |
| 175 | | * Both rasters must have the same SRID. If not, ST_MapAlgebra() should return an error: "ERROR: Operation on two geometries with different SRIDs". |
| 176 | | |
| 177 | | * An optional option could be added to allow specifying which raster, between the first and the second, is the MASTER raster. The MASTER raster is used to determine the pixeltype, the alignment and the nodatavalue of the resulting raster. If the MASTER raster do not have a nodata value defined and it is necessary to set one in the resulting raster then the minimal possible value for the pixeltype should be used as nodata value. In case this option would not be added or in case it is not passed to the function, the first raster should be used to determine those informations. |
| 178 | | |
| 179 | | * Both raster must be aligned. This is determined using the ST_SameAlignment() function described below. If a resampling is necessary they use the planned ST_Resample() function to resample the second raster to the first one before processing (or the first to the second if a MASTER parameter is provided). If ST_Resample() is not yet implemented when these functions are implemented, just return an error message: "ERROR: MapAlgebra on rasters with different alignment not yet implemented." |
| 180 | | |
| 181 | | * The computation extent (extentexpr in the functions below) and hence the extent of the resulting raster can be: |
| 182 | | |
| 183 | | -'FIRST' (the extent of the first raster) (default),[[BR]] |
| 184 | | -'SECOND' (the extent of the first raster),[[BR]] |
| 185 | | -'INTERSECTION' (the extent of the intersection of both rasters) (could also be default),[[BR]] |
| 186 | | -'UNION' (the extent of the union of both rasters). |
| 187 | | |
| 188 | | * In certain cases, for example when the computation extent is set to 'INTERSECTION', it is better to reduce the execution of ST_Mapalgebra(rast1, rast2) to overlapping rasters. This reduction is performed by adding a 'WHERE ST_Intersects(raster, raster)' clause to the query (to be implemented). If such a clause is not used and the rasters don't overlap an empty raster is returned. This behaviour is similar to the one of ST_Intersection(). Hence:[[BR]][[BR]] |
| 189 | | SELECT ST_MapAlgebra(rast1, rast2, mathExpr) FROM mytable WHERE ST_Intersects(rast1, rast2)[[BR]][[BR]] |
| 190 | | should be much faster and return many less empty rasters than:[[BR]][[BR]] |
| 191 | | SELECT ST_MapAlgebra(rast1, rast2, mathExpr) FROM mytable |
| 192 | | |
| 193 | | * The expression are evaluated by calling the equivalent or an EXECUTE SQL statement. “EXECUTE” in pl/pgsql or SPI_execute in C. This mechanism ensures that users can use any PostgreSQL function and operators in the expression as well as their own custom pl/pgsql functions. Raster values in the expressions are refered by “rast1” and “rast2”. Ex.: ST_MapAlgebra(raster1, band1, raster2, band2, “rast1 + cos(rast2) + 4”) |
| 194 | | |
| 195 | | * Alternative expressions, evaluated in place of the main expression, can be provided as parameter to deal with nodata values: |
| 196 | | |
| 197 | | -nodata1expr is evaluated when the first raster pixel value is nodata or when the first raster is absent from the pixel area being evaluated[[BR]] |
| 198 | | -nodata2expr is evaluated when the second raster pixel value is nodata or when the second raster is absent from the pixel area being evaluated[[BR]] |
| 199 | | -nodatanodataexpr is evaluated when both rasters pixel values are nodata or when the both rastera are absent from the pixel area being evaluated[[BR]] |
| 200 | | -expression is evaluated when both raster pixel values are withdata. |
| 201 | | |
| 202 | | All expressions default to NULL. When an expression is NULL every pixel filling the condition of this expression is set to NULL (nodata). |
| 203 | | |
| 204 | | * Alternate expressions like nodata1expr and nodata2expr are used to simplify complex decision expression trying to deal with the presence of nodata value pixels. Having three short expressions like this:[[BR]][[BR]] 'rast2', 'rast2', 'rast1'[[BR]][[BR]]is simpler to understand than a single complex expression dealing with nodata like this:[[BR]][[BR]] ‘CASE WHEN rast2 IS NULL THEN rast1 ELSE rast2 END’[[BR]][[BR]]This is a simple case. In more complex cases, expressions can quickly get incomprehensible and alternate expressions greatly simplify the task of writing expressions, even if having three extra parameters may seams cumbersome. |
| 205 | | |
| 206 | | |
| 207 | | |
| 208 | | * Open Question 1 (DZ): Should ST_MapAlgebra resample the resulting rasters internally, or should we let the users resample it afterward with ST_Resample: ST_Resample(ST_MapAlgebra(), originx, originy, pixelsizex, pixelsizey)[[BR]] |
| 209 | | PR: I think this would greatly contribute to simplify the API. |
| 210 | | |
| 211 | | * List of functions implementable based on the two raster version of map algebra: |
| 212 | | |
| 213 | | -ST_Intersection(raster, raster) See at the end of [http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_mapalgebra.sql plpgsql/st_mapalgebra.sql][[BR]] |
| 214 | | -ST_Union (raster, raster) not the aggregate one. See at the end of [http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_mapalgebra.sql plpgsql/st_mapalgebra.sql][[BR]] |
| 215 | | -ST_Difference() and ST_SymDifference() See at the end of [http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_mapalgebra.sql plpgsql/st_mapalgebra.sql][[BR]] |
| 216 | | -ST_Intersection(raster, raster) See at the end of [http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_mapalgebra.sql plpgsql/st_mapalgebra.sql][[BR]] |
| 217 | | -ST_Clip(raster, geometry) as a wrapper around ST_Intersection(raster, ST_AsRaster(geometry))[[BR]] |
| 218 | | -ST_BurnToRaster(): specifications below. |
| 219 | | |
| 220 | | '''Details for 5) ST_MapAlgebraExpr(rast1 raster, band1 integer, rast2 raster, band2 integer, expression text, pixeltype text, extentexpr text, nodata1expr text, nodata2expr text, nodatanodataexpr text) -> raster''' |
| 221 | | |
| 222 | | See discussion in http://trac.osgeo.org/postgis/ticket/590 |
| 223 | | |
| 224 | | '''ST_SameAlignment(raster, raster) -> boolean - done see below''' |
| 225 | | |
| 226 | | ---- |
| 227 | | == '''Objective FV.04 - Being able to use "group by" to accumulate tiles to form a new raster.''' == |
| | 97 | == '''Objective FV.26 - ST_Difference() and ST_SymDifference() based on ST_MapAlgebra(raster, raster)''' == |
| | 98 | |
| | 99 | -ST_Difference(raster1, band1, raster2, band2) |
| | 100 | |
| | 101 | -ST_SymDifference(raster1, band1, raster2, band2) |
| | 102 | |
| | 103 | See at the end of [http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_mapalgebra.sql plpgsql/st_mapalgebra.sql][[BR]] |
| | 104 | |
| | 105 | ---- |
| | 106 | == '''Objective FV.04 - Have ST_Union() variants accepting a set of temporary expressions''' == |
