PostGIS Raster Working Specifications for Future Versions
Quick Links
RoadMap (in order of importance)(not every items have a description below)
1) Have a stable and fast GDAL driver (see driver's page, open tickets and working specifications page).
2) Break up RASTER_mapAlgebra2 so that the main engine is in rt_core instead of rt_pg. This is needed for C-based aggregate functions that run against ST_MapAlgebra.
3) Multiband input for ST_MapAlgebra(raster, raster).
4) C version of the ST_Union(raster) aggregate (Ticket #1364).
5) Optimized version of two rasters ST_MapAlgebra() setting pixel areas when possible instead of just pixels by pixels (Objective FV.25 below).
6) ST_MapAlgebraFctNbg() working on a tiled coverage (See Objective FV.24 below).
7) Different variant of ST_SetValues() (Objective FV.23 below and ticket #595).
8) Integrate the ST_CreateIndexRaster() plpgsql prototype (do not need to be implemented in C to be fast).
9) ST_Tile() (implemented in C).
10) Integrate the ST_AreaWeightedSummaryStats() plpgsql prototype (do not need to be implemented in C to be fast).
11) Integrate the ST_SummaryStatsAgg() plpgsql prototype (do not need to be implemented in C to be fast) or develop a C Aggregate versions (item 15).
12) Add some function interpolating a raster from a geometry layer (See Objective FV.28 below).
13) Integrate different custom functions, to be used with ST_MapAlgebraFct(), to extract of derive values from a coverage (raster or vector). (e.g. ST_Distance4ma resulting from the GSoC 2012) (See Objective FV.27 and FV.28 below)
14) ST_CreateOverview(), ST_IsRegularlyTiled(), ST_HasOverlaps(), ST_HasGaps(), ST_HasTileSameSize(), ST_HasTileAligned() (See Objective FV.20 below)
15) C Aggregate versions of stats functions (see the ST_SummaryStatsAgg() plpgsql prototype and ticket #1048).
16) Add raster support to the PostGIS GUI Loader.
17) Faster out-of-db data loading.
18) Set and keep up to date the IsNodataValue flag indicating that the whole band is nodata and add related optimizations. (See Objective FV.22 below and tickets #593, #594).
19) Multi-band ST_AddBand().
20) ST_UnionToRaster() and ST_BurnToRaster() (See Objective FV.19 below).
21) Have ST_Union() variants accepting a set of temporary expressions (necessary to add other aggregation options like RANGE…) (See Objective FV.04 below).
22) ST_Difference() and ST_SymDifference() (See Objective FV.26 below)
23) CUnit tests.
24) Subtiling of rasters (PostGIS 3.0 or create a new raster type) (See Objective FV.29 below).
Other important TODO
- Write a "Best practices for storing raster in PostGIS".
- Write a "Best practices for third party applications wishing to read raster stored in PostGIS".
- More tests with MapServer.
- More tests with GeoServer.
Objectives' Details
Note that accomplished objectives have been moved at the bottom of the page.
Objective FV.04 - Have ST_Union() variants accepting a set of temporary expressions
ST_Union(set of raster[, p_expression, p_nodata1expr, p_nodata2expr, p_nodatanodataexpr, t_expression, t_nodata1expr, t_nodata2expr, t_nodatanodataexpr, f_expression, f_nodata1expr, f_nodata2expr, f_nodatanodataexpr]) → raster
- ST_Union is an aggregate function which merges, one by one, all the selected rasters into a unique raster.
- It can be used to simply merge disjoint tiles into a bigger raster but also to blend overlapping areas together using different rules.
- As an aggregate function it is using a STATE and a FINAL function materialized by a state expression (p_expression) and a final expression (f_expression) and their respective nodata value alternatives (p_nodata1expr, p_nodata2expr, p_nodatanodataexpr, f_nodata1expr and f_nodata2expr, f_nodatanodataexpr). All expressions are SQL expressions that are evaluated by the internal two raster version of ST_MapAlgebra().
- An optional temporary expression (t_expression) and its nodata value alternatives (t_nodata1expr, t_nodata2expr, t_nodatanodataexpr) can be provided of used internally depending on which variant or which predefined expression is passed to ST_Union(). When it is used, MapAlgebra is called two times: one time to fill this temporary band using the t_expression and a second time to fill the first band using the p_expression. This temporary band is useful to implement some expressions needing two variables like the predefined ‘MEAN’ and ‘RANGE’ expressions.
- The different expressions resume like this:
- p_expression, p_nodata1expr, p_nodata2expr and p_nodatanodataexpr is the main set of state expressions accumulating the final values if no final expression are used.
- t_expression, t_nodata1expr, t_nodata2expr and t_nodatanodataexpr is the optional set of temporary expressions accumulating values in a temporary raster. t_expressions are evaluated before p_expressions.
- f_expression, f_nodata1expr, f_nodata2expr and p_nodatanodataexpr is the optional set of final expressions. Final expressions may refer to the raster resulting from the p_expression (rast1) and the t_expression (rast2) to determine the final pixel values.
- As explained in the ST_MapAlgebra() specifications, alternate expressions like p_nodata1expr and f_nodata2expr are used to simplify complex decision expression trying to deal with the presence of nodata value pixels. Having three short expressions like this:
'rast2', 'rast2', 'rast1', NULL
is simpler than a single complex expression dealing with nodata like this:
'CASE WHEN rast1 IS NULL AND rast2 IS NULL THEN NULL WHEN rast2 IS NULL THEN rast1 ELSE rast2 END'
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 so many parameters may seams cumbersome. One must also consider the fact that in many simple cases, the t_ set and the f_ set of expressions is optional as well as the alternate expressions of the p_ set. In these very simple cases only one p_ expression has to be passed to or used by ST_Union. Another factor reducing the apparent complexity is that in many cases, users will be happy with the predefined expressions, reducing a call to ST_Union to something like ST_Union(rast, 'SUM').
- Here are some commonly used predefined expressions and how they materialize as the four (4) "p_", "t_" and "f_" expressions passed to ST_Mapalgebra(rast1, rast2):
-LAST: 'rast2', 'rast2', 'rast1', NULL
-FIRST: 'rast1', 'rast2', 'rast1', NULL
-MIN: 'LEAST(rast1, rast2)', 'rast2', 'rast1', NULL
-MAX: 'GREATEST(rast1, rast2)', 'rast2', 'rast1', NULL
-COUNT: 'rast1 + 1', '1', 'rast1', '0'
-SUM: 'rast1 + rast2', 'rast2', 'rast1', NULL
-RANGE:
-p_expressions = 'LEAST(rast1, rast2)', 'rast2', 'rast1', NULL
-t_expressions = 'GREATEST(rast1, rast2)', 'rast2', 'rast1', NULL
-f_expressions = 'rast1 - rast2', NULL, NULL, NULL
-MEAN:
-p_expressions = 'rast1 + rast2', 'rast2', 'rast1', NULL sum
-t_expressions = 'rast1 + 1', '1', 'rast1', '0' count
-f_expressions = 'CASE WHEN rast2 > 0 THEN rast1 / rast2::float8 ELSE NULL END', NULL, NULL, NULL sum/count
-STANDARD_DEVIATION: NOT POSSIBLE, nead two passes. See how ST_SummaryStats() does it…
-MAX_LENGTH, MAX_AREA or MAX_ANYTHING (requires MapAlgebra being able to reference other bands then the first one):
-p_expressions = 'CASE WHEN rast1[2] > rast2[t] THEN rast1[1] ELSE rast2[1] END', 'rast2[1]', 'rast1[1]', NULL
-t_expressions = 'CASE WHEN rast1[2] > rast2[2] THEN rast1[2] ELSE rast2[2] END', 'rast2[2]', 'rast1[2]', NULL
-f_expressions = 'rast1', NULL, NULL, NULL
- A PL/pgSQL prototype of ST_Union exist in script/plpgsql/st_union.sql.
- Variants
1) ST_Union(raster, band)
2) ST_Union(raster, band, p_expression)
3) ST_Union(raster, band, p_expression, p_nodata1expr, p_nodata2expr, p_nodatanodataexpr)
4) ST_Union(raster, band, p_expression, t_expression)
5) ST_Union(raster, band, p_expression, p_nodata1expr, p_nodata2expr, p_nodatanodataexpr, t_expression, t_nodata1expr, t_nodata2expr, t_nodatanodataexpr)
6) ST_Union(raster, band, p_expression, t_expression, f_expression)
7) ST_Union(raster, band, p_expression, p_nodata1expr, p_nodata2expr, p_nodatanodataexpr, t_expression, t_nodata1expr, t_nodata2expr, t_nodatanodataexpr, f_expression, f_nodata1expr, f_nodata2expr, f_nodatanodataexpr)
ST_MultibandUnion(raster, band)
- Same as ST_Union but handle multi-band rasters.
ST_Collect(raster set|geometry set, 'raster'|'geometry') → raster/geometry
Objective FV.06 - Being able to do some basic raster operations.
ST_Flip(raster|geometry, 'vertical'|'horizontal') → same type as first argument
Objective FV.07 - Being able to convert a raster to standards formats.
ST_AsKML(raster|geometry) → string
ST_AsSVG(raster|geometry) → string
Objective FV.08 - Being able to control the validity of a raster.
ST_mem_size(raster|geometry) → integer
ST_isvalid(raster|geometry) → boolean
Objective FV.12 - Being able to use all the other topological operators
ST_Equals(raster|geometry, raster|geometry)
ST_Disjoint(raster|geometry, raster|geometry)
ST_Touches(raster|geometry, raster|geometry)
ST_Crosses(raster|geometry, raster|geometry)
ST_Covers(raster|geometry A, raster|geometry B)
ST_IsCoveredBy(raster|geometry A, raster|geometry B)
ST_Relate(raster|geometry, raster|geometry, intersectionPatternMatrix )
Objective FV.13 - Being able to edit a raster
ST_Affine(raster|geometry, …) → same type as input
ST_Translate(raster|geometry, …) → same type as input
ST_Scale(raster|geometry, …) → same type as input
ST_TransScale(raster|geometry, …) → same type as input
ST_RotateZ,Y,Z(raster|geometry, float8) → same type as input
Objective FV.15 - Being able to intersect two rasters to get a geometry.
ST_Intersection(raster, integer, raster, integer, 'geometry') → geometry - Returns a two bands raster with values only in the intersecting areas of both rasters. Integer parameters are the band number of the raster.
Variants
1) ST_Intersection(raster, integer, raster, integer, 'geometry') → geometry
2) ST_Intersection(raster, raster, integer, 'geometry') → geometry — default first raster to band # 1
3) ST_Intersection(raster, integer, raster, 'geometry') → geometry — default second raster to band # 1
4) ST_Intersection(raster, raster, 'geometry') → geometry — default both raster to band # 1
Objective FV.17 - Being able to refer to band by name. (PostGIS 3.0)
Add 8 digit string to each band in the base raster WKB format.
Adjust gdal2wktraster.py to be able to give names to each band when importing.
Adjust/overlaod every function to be able to refer to raster band by name.
Objective FV.18 - Being able to load rasters from SQL
The idea is to change the rt_band_get_data core function so it loads filesystem registered raster data using GDAL into the data base. This allow us to create a list of raster with a new ST_MakeRegisteredRaster("c:/temp/mytiff/*.tif") function and to convert them to internal rasters with a CREATE TABLE indbraster AS SELECT ST_MakeInDB(ST_MakeRegisteredRaster("c:/temp/mytiff/*.tif"), band) statement.
Changes to the rt_band_get_data core function
ST_MakeRegisteredRaster(wildcardPath)
ST_MakeInDB(rast, band)
Objective FV.19 - Being able burn geometries to existing raster
ST_BurnToRaster(rast, geometry[, value, expression, nodata1expr, nodata2expr, nodatanodataexpr]) → raster
- Burn a geometry in the band of an existing raster.
- No geoferencing metadata have to be provided since the ones of the provided raster are used.
- This function is mainly different from the ST_AsRaster(geometry) by the fact that it must provide some options to take into account the existing pixel values when burning the new value associated with the geometry.
- This function could also be called ST_SetValues(). However does not just "set" values. It set them in a more flexible way possibly taking existing values into account. This is why "Burn" might be a better name. Very open to other name propositions…
- It should be implemented as some wrappers around ST_MapAlgebra(raster, band, ST_AsRaster(geometry, val, raster, "CROP"), 'expression', 'nodata1expr', 'nodata2expr', 'nodatanodataexpr').
- 'expression', 'nodata1expr', 'nodata2expr', 'nodatanodataexpr' are SQL expressions passed to the underlying two raster map algebra function:
-'expression' is evaluated when both pixels have a value associated with them. The default for this expression is 'LAST'.
-'nodata1expr' is evaluated when rast1 is a nodata value and rast2 has a value associated with it. The default for this expression is NULL.
-'nodata2expr' is evaluated when rast2 is a nodata value and rast1 has a value associated with it. The default for this expression is NULL.
-'nodatanodataexpr' is evaluated when both pixels are nodata values. The default for this expression is NULL.
- Here are some commonly used predefined expressions and how they materialize as the ST_Mapalgebra(rast1, rast2) expression, nodata1expr and nodata2expr, and nodatanodataexpr.
-LAST: expression = 'rast2'. This is the default ST_BurnToRaster expression.
-COUNT: expression = '2', nodata1expr = '1', nodata2expr = '1', nodatanodataexpr = '0'
-MAX: expression = 'GREATEST(rast1, rast2)' , nodata1expr = 'rast2', nodata2expr = 'rast1', nodatanodataexpr = 'NULL'
-MIN: expression = 'LEAST(rast1, rast2)' , nodata1expr = 'rast2', nodata2expr = 'rast1', nodatanodataexpr = 'NULL'
-RANGE: expression = 'ABS(rast1 - rast2)' , nodata1expr = '0', nodata2expr = '0', nodatanodataexpr = 'NULL'
-SUM: expression = 'rast1 + rast2' , nodata1expr = 'rast2', nodata2expr = 'rast1', nodatanodataexpr = 'NULL'
-MEAN: expression = '(rast1 + rast2)/2' , nodata1expr = 'rast2', nodata2expr = 'rast1', nodatanodataexpr = 'NULL'
- When expression is 'LAST' the process can be implemented with GDALRasterizeLayers since, apparently GDALRasterizeLayers do not allow other burning rules.
ST_UnionToRaster(set of geometry, val[, pixeltype, nodataval, x float8, y float8, pixelsizex, pixelsizey, skewx, skewy, p_expression, p_nodata1expr, p_nodata2expr, p_nodatanodataexpr, t_expression, t_nodata1expr, t_nodata2expr, t_nodatanodataexpr, f_expression, f_nodata1expr, f_nodata2expr, f_nodatanodataexpr]) → raster
- ST_UnionToRaster is an aggregate function very similar to ST_Union(set of raster) which burn, one after the other, all the geometries of a table, or selected by a GROUP BY clause, into a unique raster.
- It is implemented as a wrapper around ST_Union(ST_AsRaster(geometry, value, pixeltype, nodataval, x float8, y float8, pixelsizex, pixelsizey, skewx, skewy), p_expression, p_nodata1expr, p_nodata2expr, p_nodatanodataexpr, t_expression, t_nodata1expr, t_nodata2expr, t_nodatanodataexpr, f_expression, f_nodata1expr, f_nodata2expr, f_nodatanodataexpr)) function.
- Variants
1) ST_UnionToRaster(set of geometry, val)
2) ST_UnionToRaster(set of geometry, val, pixeltype, nodataval, x float8, y float8, pixelsizex, pixelsizey, skewx, skewy, p_expression)
3) ST_UnionToRaster(set of geometry, val, pixeltype, nodataval, x float8, y float8, pixelsizex, pixelsizey, skewx, skewy, p_expression, p_nodata1expr, p_nodata2expr, p_nodatanodataexpr)
4) ST_UnionToRaster(set of geometry, val, pixeltype, nodataval, x float8, y float8, pixelsizex, pixelsizey, skewx, skewy, p_expression, t_expression)
5) ST_UnionToRaster(set of geometry, val, pixeltype, nodataval, x float8, y float8, pixelsizex, pixelsizey, skewx, skewy, p_expression, p_nodata1expr, p_nodata2expr, p_nodatanodataexpr, t_expression, t_nodata1expr, t_nodata2expr, t_nodatanodataexpr)
6) ST_UnionToRaster(set of geometry, val, pixeltype, nodataval, x float8, y float8, pixelsizex, pixelsizey, skewx, skewy, p_expression, t_expression, f_expression)
7) ST_UnionToRaster(set of geometry, val, pixeltype, nodataval, x float8, y float8, pixelsizex, pixelsizey, skewx, skewy, p_expression, p_nodata1expr, p_nodata2expr, p_nodatanodataexpr, t_expression, t_nodata1expr, t_nodata2expr, t_nodatanodataexpr, f_expression, f_nodata1expr, f_nodata2expr, f_nodatanodataexpr)
Objective FV.20 - Being able to determine topological characteristics of a coverage
ST_HasNoOverlaps(rasttable text, rastcolumn text)
ST_HasNoGaps(rasttable text, rastcolumn text)
ST_HasTileSameSize(rasttable text, rastcolumn text)
ST_HasTileAligned(rasttable text, rastcolumn text) — this is the actually the same as ST_SameAlignment()
ST_IsRegularlyTiled(rasttable text, rastcolumn text) if the four conditions above are fulfilled. The ST_HasTileSameSize() condition could be optional.
ST_HasOverlaps() and ST_HasGaps() might be interesting to apply on a geometry coverage as well.
ST_HasOverlaps(), ST_HasGaps(), ST_HasTileSameSize() and ST_HasTileAligned() could also be AGGREGATE functions so that they can be applied on the subset of a table. In this case they would accept a single 'rast' (or, even more useful, a 'geometry') parameter.
Objective FV.21 - Being able to create overviews in SQL
ST_CreateOverviews(rasttable text, rastcolumn text, overviewnumber int)
Creates the tiles corresponding to the lower resolution (or overviews) of a raster table. So that an overview table can be created like this: SELECT ST_CreateOverviews(rasttable, rast, 2);
This method is dependent on two helper functions: ST_SetValues (see FV.23) and ST_MakeEmptyCoverage (below)
Pseudo code:
IF ST_IsRegularlyTiled(rasttable text, rastcolumn text) THEN
— Check if the tiles can be divided into 2x2, 4x4, 8x8
IF 'SELECT ST_Width(rast) modulo 2 = (overviewnumber + 1) AND ST_Width(rast) modulo (overviewnumber + 1) = 0 FROM rasttable LIMIT 1'
ST_Resample(all the tiles)
Add a column with the number of the set of tile to union together
ST_Union(using this number)
ELSE
Add a column with the number of the set of tile to union together
ST_Union(using this number)
ST_Resample(all the tiles)
ELSE
Resample(all the tiles) — Without unioning them afterward
END IF
ST_MakeEmptyCoverage(tilewidth int, tileheight int, width int, height int, upperleftx float8, upperlefty float8, scalex float8, scaley float8, skewx float8, skewy float8)
See ticket #2249
Creates a set of rasters that regularly cover the area described by the coordinates. The net coverage of the returned set of rasters is the contiguous coverage over the area described by the coordinates. Each individual raster in the set will have all pixel values set to NODATA.
If the width of the empty coverage is not a multiple of the tilewidth, then the empty coverage will contain rasters on the right (highest upperleftx value) that are "thinner": they will be (width - (tilewidth * (width % tilewidth))) wide.
If the height of the empty coverage is not a multiple of the tileheight, then the empty coverage will contain rasters on the bottom (lowest upperlefty value) that are "shorter": they will be (height - (tileheight * (height % tileheight))) tall.
If the width and height are not multiples of tilewidth and tileheight, respectively, then the empty coverage will contain one raster at the bottom right corner (highest upperleftx value, lowest upperlefty value) that is both "thinner" and "shorter" by the calculations above for width and height, respectively.
- Variants
1) ST_MakeEmptyCoverage(tilewidth int, tileheight int, width int, height int, upperleftx float8, upperlefty float8, scalex float8, scaley float8, skewx float8, skewy float8, srid int4)
Variant 1 accepts an SRID parameter.
Objective FV.22 - Implement better support for NULL, Empty, HasNoBand(rast), HasNoBand(rast, band) and BandIsNoData rasters in all functions.
Each function should better handle NULL rasters, empty rasterm rasters with no band and bands only filled with nodata values.
1) Generally, when NULL rasters are provided, return NULL. If the function involves a second non NULL raster and something can be done with it, do it.
2) Generally, when empty rasters are provided (ST_IsEmpty, width=0 or height=0), return an empty raster. If the function involves a second non empty raster and something can be done with it, do it.
2) Generally, when a HasNoBand(rast) or a HasNoBand(rast, band) raster is provided return a raster with no band but with the right extent. If the function involves a second raster having a band or the band, treat the first missing band like a BandIsNoData.
4) A BandIsNoData raster is a normal raster but many functions can be optimized with this type of raster.
5) All functions having missing the requested information (about its extent when it is a NULL or an empty raster or about its band when it is a HasNoBand(rast) or a HasNoBand(rast, band) raster) should return NULL or a documented value.
6) Try as less as possible to return EXCEPTION (or ERROR).
See discussion in http://trac.osgeo.org/postgis/ticket/594
ST_IsEmpty(rast) → boolean
Returns TRUE if the raster width or the raster height is 0.
ST_HasNoBand
Variants
1) ST_HasNoBand(rast, band) → boolean
2) ST_HasNoBand(rast) → boolean
Returns TRUE if the the raster does not have this band.
Variant 2 returns TRUE if the the raster does not have any band.
ST_BandIsNoData
Variants
1) ST_BandIsNoData(rast, band) → boolean
2) ST_BandIsNoData(rast) → boolean
Returns TRUE if the specified raster band is only filled with no data value.
Variant 2 default to band 1.
Implementation details
This require a new flag to be set in the core at import and at edition. See discussion in http://trac.osgeo.org/postgis/ticket/593
Objective FV.23 - Being able to set and get values for part of a raster.
ST_SetValues
Set a set of pixel to a value.
Those functions are provided as a faster way to set pixels than the (future) two rasters version of ST_MapAlgebra which can also be used when more conditional flexibility is necessary. They should also be used in this version of ST_MapAlgebra to quicly set large areas to a single value.
Variants
The first series of variant set a defined area of the raster to a value.
1) ST_SetValues(rast raster, band int, x int, y int, width int, height int, val float8, keepdestnodata boolean)
2) ST_SetValues(rast raster, x int, y int, width int, height int, val float8, keepdestnodata boolean)
3) ST_SetValues(rast raster, x int, y int, width int, height int, val float8)
4) ST_SetValues(rast raster, band int, x int, y int, width int, height int, val float8)
The second series of variant set a defined area of the raster to corresponding values of another raster.
5) ST_SetValues(rast1 raster, band1 int, x int, y int, width int, height int, rast2 raster, band2 int, keepdestnodata boolean, keepsourcenodata boolean)
6) ST_SetValues(rast1 raster, x int, y int, width int, height int, rast2 raster, keepdestnodata boolean, keepsourcenodata boolean)
7) ST_SetValues(rast1 raster, x int, y int, width int, height int, rast2 raster)
x, y, width & height define the area of the raster to be edited. x & y are the raster coordinates of the upper left corner of this are. When x, y, width or height are out of the raster range, only the part of the range intersecting with the raster is set to the new value.
val is the new value to be set.
When keepdestnodata is TRUE, destination nodata value pixels are not set. Default is FALSE so that nodata value pixels are set.
When keepsourcenodata is TRUE, source nodata value pixels are not copied. Default is FALSE so that nodata value pixels are copied.
Variant 2 assume band to be 1.
Variant 3 assume band to be 1 and keepdestnodata to be FALSE.
Variant 4 assume only keepdestnodata to be FALSE.
Variant 6 assume both bands to be 1.
Variant 7 assume both bands to be 1, keepdestnodata to be FALSE and keepsourcenodata to be FALSE.
Implementation details
PL/pgSQL prototypes can be found in http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_setvalues.sql
Objective FV.24 - ST_MapAlgebraFctNbg() working on a tiled coverage
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. This version work transparently on a tiled coverage avoiding edges effects.
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).
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'…
See also: Notes taken by David Zwarg during the Montreal Code Sprint 2011 and ticket #860.
Objective FV.25 - Optimized version of two rasters ST_MapAlgebra()
- ST_MapAlgebra(raster, raster, 'UNION') can be optimized by not iterating on every pixels when it is not necessary. "Chunk" (or large areas) of pixels, outside the intersecting area of the two raster can be set in a "memset manner".
- 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 or in a "memset manner") 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 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.
Objective FV.26 - ST_Difference() and ST_SymDifference() based on ST_MapAlgebra(raster, raster)
-ST_Difference(raster1, band1, raster2, band2)
-ST_SymDifference(raster1, band1, raster2, band2)
See at the end of plpgsql/st_mapalgebra.sql
Objective FV.27 - Integrate different custom functions, to be used with ST_MapAlgebraFct(), to extract values from a coverage
- These are quick to implement ST_MapAlgebraFct() custom functions setting raster pixel values from any king of existing coverage. They all accept the name of the schema, the name of table and the name of raster (or geometry) column from which to extract values. They can also accept additional parameters.
- One way to organize/classify/make sense all those possible functions goes like this:
1) Functions searching for the pixel value in a raster coverage. This is the same principle as ST_Union() except that the rasters do not have to be aligned (nice advantage!).
2) Functions searching for the pixel value in a vector coverage. This is the same principle as ST_UnionToRaster() and ST_BurnToRaster().
3) Functions aggregating pixel values from a raster neighborhood. This is the same principle as ST_MapAlgebraFctNbg() except that the target raster do not have to be aligned with the source raster (nice advantage!).
4) Functions aggregating pixel values from a vector neighborhood. This is the only way to do that.
5) Functions deriving metric values from a raster or a geometry coverage. This is very similar to the two previous items but the result is not an aggregate value. e.g. Distance to something.
Category 1), not working on pixel neighbour area, hence only on aligned pixels from overlapping rasters, should eventually be replaced by the optimized version of ST_Union() (itself dependent on the optimized version of ST_MapAlgebra See Objective FV.25 above). Those alternative are less flexible than their future ST_Union alternative in that they only work on a whole table, not on a selection of a table (i.e. they are not aggregates like ST_Union is. One can always create a view on a table to work on a table subset though).
- Some examples of such custom functions exists:
- ST_FasterUnion('schemaname', 'tablename', 'rastercolumnname'), a temporary alternative to ST_Union(rast, 'LAST') - Merge all the tiles of a table together into a single raster. This function act as a wrapper around the ST_FirstRasterValue4ma() ST_MapAlgebraFct() custom function and is described here.
- ST_FirstGeomValue4ma() - Get the value from a specified column from the first geometry intersecting with the shape of the pixel (See this post to postgis-users)
- Other possible functions:
- ST_GeomToRaster('schemaname', 'tablename', 'geomolumnname', 'METHOD') (Category 2) - A generalization of ST_FirstGeomValue4ma() described above but accepting more methods of value extraction, depending on the type of coverage we want to extract from. This is very much like doing an intersection (actually more like an "identify" overlay operation) between a raster and a vector coverage.
This function can further be generalised by adding parameters to define a buffer area around the pixel centroid or the pixel boundary. This would provide a category 4 type of function. (e.g. ST_Density())
The first series of METHODs return statistics about the intersection of the footprint of each pixel (the intersecting area with polygons must be greater than 0.0000000001 to be considered):
- with any kind of coverage, returning value frequencies: COUNT_OF_VALUES, COUNT_OF_DISTINCT_VALUES, COUNT_OF_MOST_FREQUENT_VALUE, COUNT_OF_LEAST_FREQUENT_VALUE
- with any kind of coverage, returning brute (in opposition to weighted) stats on values: FIRST_VALUE, LAST_VALUE, MOST_FREQUENT_VALUE, LEAST_FREQUENT_VALUE, MIN_VALUE, MAX_VALUE, RANGE_OF_VALUES, SUM_OF_VALUES, MEAN_OF_VALUES, STDDEV_OF_VALUES
- with a line coverage, returning length stats: FIRST_LENGTH, LAST_LENGTH, MIN_LENGTH, MAX_LENGTH, RANGE_OF_LENGTHS, SUM_OF_LENGTHS, MEAN_OF_LENGTHS, STDDEV_OF_LENGTHS
- with a line coverage, returning stats on length of merged lines (having the same value): MIN_LENGTH_OF_MERGED_LINES, MAX_LENGTH_OF_MERGED_LINES, RANGE_OF_LENGTHS_OF_MERGED_LINES, MEAN_LENGTH_OF_MERGED_LINES, STDDEV_LENGTH_OF_MERGED_LINES, TOTAL_LENGTH_OF_LINES_HAVING_GREATEST_VALUE, TOTAL_LENGTH_OF_LINE_HAVING_SMALLEST_VALUE, TOTAL_LENGTH_OF_LINE_WITH_MOST_FREQUENT_VALUE, TOTAL_LENGTH_OF_LINE_WITH_LEAST_FREQUENT_VALUE
- with a line coverage, returning the value of lines having some special lengths: VALUE_OF_LONGEST, VALUE_OF_SHORTEST, VALUE_OF_MERGED_LONGEST, VALUE_OF_MERGED_SHORTEST
- with a line coverage, returning the value of lines weighted by their lengths: LENGTH_WEIGHTED_SUM_OF_VALUES, LENGTH_WEIGHTED_MEAN_OF_VALUES, LENGTH_WEIGHTED_STDDEV_OF_VALUES
- with a polygon coverage, returning area stats: FIRST_AREA, LAST_AREA, MIN_AREA, MAX_AREA, RANGE_OF_AREAS, SUM_OF_AREAS, MEAN_OF_AREAS, STDDEV_OF_AREAS
- with a polygon coverage, returning stats on area of merged polygons (having the same value): MIN_AREA_OF_MERGED_POLYGONS, MAX_AREA_OF_MERGED_POLYGONS, RANGE_OF_AREAS_OF_MERGED_POLYGONS, MEAN_AREA_OF_MERGED_POLYGONS, STDDEV_AREA_OF_MERGED_POLYGONS, TOTAL_AREA_OF_POLYGONS_HAVING_GREATEST_VALUE, TOTAL_AREA_OF_POLYGON_HAVING_SMALLEST_VALUE, TOTAL_AREA_OF_POLYGON_WITH_MOST_FREQUENT_VALUE, TOTAL_AREA_OF_POLYGON_WITH_LEAST_FREQUENT_VALUE
- with a polygon coverage, returning the value of polygons having some special areas: VALUE_OF_BIGGEST, VALUE_OF_SMALLEST, VALUE_OF_MERGED_BIGGEST, VALUE_OF_MERGED_SMALLEST
- with a polygon coverage, returning the value of polygons weighted by their areas: AREA_WEIGHTED_SUM_OF_VALUES, AREA_WEIGHTED_MEAN_OF_VALUES, AREA_WEIGHTED_STDDEV_OF_VALUES
The second series of METHODs return statistics about the intersection of the centroids of each pixel:
- with any kind of coverage, returning value frequencies: COUNT_OF_VALUES_AT_PIXEL_CENTROID, COUNT_OF_DISTINCT_VALUES_AT_PIXEL_CENTROID, COUNT_OF_MOST_FREQUENT_VALUE_AT_PIXEL_CENTROID, COUNT_OF_LEAST_FREQUENT_VALUE_AT_PIXEL_CENTROID
- with any kind of coverage, returning brute (in opposition to weighted) stats on values: FIRST_VALUE_AT_PIXEL_CENTROID, LAST_VALUE_AT_PIXEL_CENTROID, MOST_FREQUENT_VALUE_AT_PIXEL_CENTROID, LEAST_FREQUENT_VALUE_AT_PIXEL_CENTROID, MIN_VALUE_AT_PIXEL_CENTROID, MAX_VALUE_AT_PIXEL_CENTROID, RANGE_OF_VALUES_AT_PIXEL_CENTROID, SUM_OF_VALUES_AT_PIXEL_CENTROID, MEAN_OF_VALUES_AT_PIXEL_CENTROID, STDDEV_OF_VALUES_AT_PIXEL_CENTROID
- ST_GeomToRaster(rast, 'schemaname', 'tablename', 'geomolumnname', 'METHOD') - A variant of ST_GeomToRaster() accepting an already defined raster. This is the equivalent to a variant of ST_UnionToRaster(rast, set of geometry, …) described in Objective FV.19 but not an aggregate and faster while we wait for a faster ST_Union().
- ST_EuclidianDistance() (see the GSoC project) this function takes advantage of the KNN index of geometry to quickly find the geometry nearest to the pixel centroid and compute its distance.
Objective FV.28 - Add some function interpolating a raster from a geometry layer
- ST_Interpolate('schemaname', 'tablename', 'geomolumnname', 'METHOD') - Method could be 'NEAREST_NEIGHBOR', 'IDW', 'KRIGING', 'NATURAL_NEIGHBOR', 'SPLINE', 'REG_SPLINE_WITH_TENSION', etc…
Objective FV.29 - Subtiling of rasters (PostGIS 3.0 or create a new raster type)
Rewrite the RASTER type so that big rasters stored in one row (max 1GB) are internally tiled and quickly accessed.
Other functions
ST_AsBinary(raster, compression)
ST_RasterFromWKB(raster, [<srid>])
ST_RasterFromText(string, [<srid>])
ST_AsText(raster)
Cancelled Objectives
These now work by first converting the raster using ST_Polygon(raster)
ST_Area(raster|geometry) → double
ST_Centroid(raster|geometry) → point geometry
ST_PointOnSurface(raster|geometry) → point geometry
ST_Buffer(raster|geometry, double) → same type as first arg.
Accomplished Objectives
ST_ValueCount(raster, value) → integer - done see below
ST_ValuePercent(raster, value) → double precision - done see below
ST_Resample(raster, method, originx, originy, pixelsizex, pixelsizey) → raster - done see below
ST_SelectByValue(raster, 'expression') → same type as first argument
ST_Clip(raster|geometry,geometry) → same type as first argument
Objective FV.09 - Being able to use other major topological operators
ST_Within(raster|geometry A, raster|geometry B)
ST_Contains(raster|geometry A, raster|geometry B)
ST_Overlaps(raster|geometry, raster|geometry)
