Changes between Version 65 and Version 66 of WKTRaster/SpecificationWorking02


Ignore:
Timestamp:
07/24/12 14:46:11 (12 years ago)
Author:
pracine
Comment:

--

Legend:

Unmodified
Added
Removed
Modified
  • WKTRaster/SpecificationWorking02

    v65 v66  
    12531253~~ Rasterize the geometry as a new raster (ST_AsRaster(geometry, pixeltype, val, nodataval, raster)) and then copy only pixels for which both raster bands have a value.
    12541254~~ Should be implemented as a wrapper around ST_MapAlgebra after rasterizing the geometry to a raster having the same alignment as the raster.
     1255
     1256----
     1257== '''Objective FV.03 - Implement all the necessary versions of ST_MapAlgebra''' ==
     1258
     1259Different versions of ST_MapAlgebra are planned:
     1260
     1261One raster versions:
     1262
     12631) '''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...
     1264
     12652) '''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.
     1266
     12673) '''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.
     1268
     12694) '''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.
     1270Some versions were done:
     1271
     1272Two rasters versions. These versions must take into account different alignment, different extents, nodata and non-existent values.
     1273
     12745) '''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.
     1275
     12766) '''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.
     1277
     12787) '''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.
     1279
     12808) '''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.
     1281
     1282Parameters and details for each variant follow...
     1283
     1284'''Details for 1) ST_MapAlgebraExpr(raster rast, integer band, text expression, text nodatavalueexpr, text pixeltype) -> raster'''
     1285
     1286 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].
     1287
     1288
     1289
     1290'''Details for 3) ST_MapAlgebraFctNgb()'''
     1291
     1292 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.
     1293
     1294 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:
     1295
     1296  SELECT ST_MapAlgebraFctNgb(rast, band, pixeltype, "ST_Mean", 2, 2, "ignore")
     1297
     1298 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."
     1299
     1300 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:
     1301
     1302  -"NULL": If any value is a nodata value, return NULL.[[BR]]
     1303  -"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]]
     1304  -"value": Replace any nodata value with the value of the pixel being computed.[[BR]]
     1305  -a value: Replace any nodata value with this value and compute.
     1306
     1307   '''Open questions'''
     1308
     1309   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.
     1310
     1311   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.
     1312
     1313   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.
     1314
     1315 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).
     1316
     1317 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...
     1318
     1319 Users could write their own map algebra summarizing functions.
     1320
     1321 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:
     1322
     1323  SELECT ST_MapAlgebraNgb(rast, band, pixeltype, 2, 2, "ST_Mean", "ignore")
     1324
     1325 and the dimensions do not have to be passed to the summarizing functions since it could deduce them from ST_Width & ST_Height.
     1326
     1327 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.:
     1328
     1329  SELECT ST_MapAlgebraNgb("mycoveragetable", "myrastercolumn", band, pixeltype, 2, 2, "ST_Mean", "ignore")
     1330
     1331 Three difficulties must be solved to implement this function:
     1332
     1333  -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]]
     1334  -We must see how it is possible to call a PL/pgSQL function from a C function[[BR]]
     1335  -We must see how to pass a variable number of parameter to a PL/pgSQL function
     1336
     1337 See also: [wiki:WKTRaster/MapAlgebra Notes taken by David Zwarg during the Montreal Code Sprint 2011] and http://trac.osgeo.org/postgis/ticket/860
     1338
     1339'''Details for 3) ST_MapAlgebraFctNgb()'''
     1340
     1341 This variant works on a tiled coverage. When computing values on the edge of a tile, it has the responsibility to provide out of range value contained in neighbour tiles to the user function.
     1342
     1343 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'...
     1344
     1345'''Details for 5) to 8) Two rasters versions'''
     1346
     1347 * A simple PL/pgSQL prototype of the two raster version of ST_MapAlgebra() exists in http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_mapalgebra.sql. This version iterates over every pixel of the unionized extent of two rasters even if there are large areas where both rasters are absent and hence are interpreted as nodata values.
     1348
     1349 * The prototype of an optimized version, trying to set those large areas of nodata values as well as areas where only one raster is present (this is the case when unioning two contiguous non-overlapping rasters) as a block (not pixel by pixel) is still in development. See http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_mapalgebra_optimized.sql The idea is to set blocks of the resulting raster using ST_SetValues() (described in [http://trac.osgeo.org/postgis/wiki/WKTRaster/SpecificationWorking02 Objective 2.0.05]) instead of processing one pixel at a time. This is somewhat similar to setting a large block of memory with memcpy() or memset() rather than setting a buffer one value at a time. The resulting raster is divided into rectangular block with the _MapAlgebraParts() function.
     1350
     1351 * Both rasters must have the same SRID. If not, ST_MapAlgebra() should return an error: "ERROR:  Operation on two geometries with different SRIDs".
     1352
     1353 * 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.
     1354
     1355 * 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."
     1356
     1357 * The computation extent (extentexpr in the functions below) and hence the extent of the resulting raster can be:
     1358
     1359  -'FIRST' (the extent of the first raster) (default),[[BR]]
     1360  -'SECOND' (the extent of the first raster),[[BR]]
     1361  -'INTERSECTION' (the extent of the intersection of both rasters) (could also be default),[[BR]]
     1362  -'UNION' (the extent of the union of both rasters).
     1363
     1364 * 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]]
     1365  SELECT ST_MapAlgebra(rast1, rast2, mathExpr) FROM mytable WHERE ST_Intersects(rast1, rast2)[[BR]][[BR]]
     1366 should be much faster and return many less empty rasters than:[[BR]][[BR]]
     1367  SELECT ST_MapAlgebra(rast1, rast2, mathExpr) FROM mytable
     1368
     1369 * 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”)
     1370
     1371 * Alternative expressions, evaluated in place of the main expression, can be provided as parameter to deal with nodata values:
     1372
     1373  -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]]
     1374  -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]]
     1375  -nodatanodataexpr is evaluated when both rasters pixel values are nodata or when the both rastera are absent from the pixel area being evaluated[[BR]]
     1376  -expression is evaluated when both raster pixel values are withdata.
     1377
     1378 All expressions default to NULL. When an expression is NULL every pixel filling the condition of this expression is set to NULL (nodata).
     1379
     1380 * 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.
     1381
     1382
     1383
     1384 * 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]]
     1385   PR: I think this would greatly contribute to simplify the API.
     1386
     1387 * List of functions implementable based on the two raster version of map algebra:
     1388
     1389 -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]]
     1390 -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]]
     1391 -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]]
     1392 -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]]
     1393 -ST_Clip(raster, geometry) as a wrapper around ST_Intersection(raster, ST_AsRaster(geometry))[[BR]]
     1394 -ST_BurnToRaster(): specifications below.
     1395
     1396'''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'''
     1397
     1398 See discussion in http://trac.osgeo.org/postgis/ticket/590
     1399
     1400'''ST_SameAlignment(raster, raster) -> boolean - done see below'''
    12551401
    12561402----