wiki:WKTRaster/SpecificationWorking03

Version 212 (modified by pracine, 12 years ago) ( diff )

—

PostGIS Raster Working Specifications for Future Versions

Quick Links

RoadMap (in order of importance)

  • Work on ST_MapAlgebra()
    • 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 MapAlgebra.
    • Optimized version of two rasters ST_MapAlgebra setting pixel areas when possible instead of just pixels by pixels
    • Multiband ST_MapAlgebra(raster, raster)
    • ST_MapAlgebraFctNbg working on a tiled coverage
  • Different variant of ST_SetValues()
  • ST_CreateOverview(), ST_IsRegularlyTiled(), ST_HasOverlaps(), ST_HasGaps(), ST_HasTileSameSize(), ST_HasTileAligned()
  • C version of the ST_Union(raster) aggregate
  • ST_Retile()
  • Multi-band ST_AddBand()
  • ST_UnionToRaster() and ST_BurnToRaster()
  • Set and update of the IsNodataValue flag indicating the whole band is nodata.
  • Subtiling of rasters (PostGIS 3.0 or create a new raster type)

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"

Objectives' Details

Note that accomplished objectives have been moved at the bottom of the page.

Objective FV.03 - Implement all the necessary versions of ST_MapAlgebra

Different versions of ST_MapAlgebra are still to do:

One raster versions:

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. Some versions were done:

Two rasters versions:

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.

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.

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.

One raster versions:

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…

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. 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.

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.

Two rasters versions. These versions must take into account different alignment, different extents, nodata and non-existent values.

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.

Parameters and details for each variant follow…

Details for 1) ST_MapAlgebraExpr(raster rast, integer band, text expression, text nodatavalueexpr, text pixeltype) → raster

This function is already implemented. See the specifications of Objective 2.0.02 and the ​documentation.

Details for 3) ST_MapAlgebraFctNgb()

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.

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:

SELECT ST_MapAlgebraFctNgb(rast, band, pixeltype, "ST_Mean", 2, 2, "ignore")

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."

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:

-"NULL": If any value is a nodata value, return NULL.
-"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.
-"value": Replace any nodata value with the value of the pixel being computed.
-a value: Replace any nodata value with this value and compute.

Open questions

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.

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.

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.

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).

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…

Users could write their own map algebra summarizing functions.

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:

SELECT ST_MapAlgebraNgb(rast, band, pixeltype, 2, 2, "ST_Mean", "ignore")

and the dimensions do not have to be passed to the summarizing functions since it could deduce them from ST_Width & ST_Height.

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.:

SELECT ST_MapAlgebraNgb("mycoveragetable", "myrastercolumn", band, pixeltype, 2, 2, "ST_Mean", "ignore")

Three difficulties must be solved to implement this function:

-The construction of the matrix must to be passed to the summarizing functions must be optimized when passing from one pixel to the other.
-We must see how it is possible to call a PL/pgSQL function from a C function
-We must see how to pass a variable number of parameter to a PL/pgSQL function

See also: Notes taken by David Zwarg during the Montreal Code Sprint 2011 and http://trac.osgeo.org/postgis/ticket/860

Details for 3) ST_MapAlgebraFctNgb()

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.

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'…

Details for 5) to 8) Two rasters versions

  • 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 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.
  • Both rasters must have the same SRID. If not, ST_MapAlgebra() should return an error: "ERROR: Operation on two geometries with different SRIDs".
  • 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.
  • 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."
  • The computation extent (extentexpr in the functions below) and hence the extent of the resulting raster can be:

-'FIRST' (the extent of the first raster) (default),
-'SECOND' (the extent of the first raster),
-'INTERSECTION' (the extent of the intersection of both rasters) (could also be default),
-'UNION' (the extent of the union of both rasters).

  • 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:

    SELECT ST_MapAlgebra(rast1, rast2, mathExpr) FROM mytable WHERE ST_Intersects(rast1, rast2)

    should be much faster and return many less empty rasters than:

    SELECT ST_MapAlgebra(rast1, rast2, mathExpr) FROM mytable
  • 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”)
  • Alternative expressions, evaluated in place of the main expression, can be provided as parameter to deal with nodata values:

-nodata1expr is evaluated when the first raster pixel value is nodata or when the first raster is absent from the pixel area being evaluated
-nodata2expr is evaluated when the second raster pixel value is nodata or when the second raster is absent from the pixel area being evaluated
-nodatanodataexpr is evaluated when both rasters pixel values are nodata or when the both rastera are absent from the pixel area being evaluated
-expression is evaluated when both raster pixel values are withdata.

All expressions default to NULL. When an expression is NULL every pixel filling the condition of this expression is set to NULL (nodata).

  • 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:

    'rast2', 'rast2', 'rast1'

    is simpler to understand than a single complex expression dealing with nodata like this:

    ‘CASE 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 three extra parameters may seams cumbersome.
  • 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)
    PR: I think this would greatly contribute to simplify the API.
  • List of functions implementable based on the two raster version of map algebra:

-ST_Intersection(raster, raster) See at the end of plpgsql/st_mapalgebra.sql
-ST_Union (raster, raster) not the aggregate one. See at the end of plpgsql/st_mapalgebra.sql
-ST_Difference() and ST_SymDifference() See at the end of plpgsql/st_mapalgebra.sql
-ST_Intersection(raster, raster) See at the end of plpgsql/st_mapalgebra.sql
-ST_Clip(raster, geometry) as a wrapper around ST_Intersection(raster, ST_AsRaster(geometry))
-ST_BurnToRaster(): specifications below.

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

See discussion in http://trac.osgeo.org/postgis/ticket/590

ST_SameAlignment(raster, raster) → boolean - done see below

Objective FV.04 - Being able to use "group by" to accumulate tiles to form a new raster.

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 base raster operations.

ST_Area(raster|geometry) → double

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

  • Not really useful since it is the same as ST_MapAlgebra(rast, 'CASE WHEN rast < x OR rast = y THEN rast ELSE NULL END') which is already implemented.

ST_Clip(raster|geometry,geometry) → same type as first argument
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.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)

Objective FV.10 - Being able to derive a raster layer from vector layer.

ST_Interpolate(points, pixelsize, method) → raster

Objective FV.11 - Being able to do on rasters most operations available on geometries

ST_Centroid(raster|geometry) → point geometry
ST_PointOnSurface(raster|geometry) → point geometry
ST_Buffer(raster|geometry, double) → same type as first arg.
ST_Difference(raster|geometry A, raster|geometry B) → same type as first argument
ST_SymDifference(raster|geometry,raster|geometry,'raster'|'geometry') → raster/geometry

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 textual name.

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 load 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") and to convert them witinot a CREATE TABLE with a ST_MakeBandInDB(rast, band)

Changes to the rt_band_get_data core function
ST_MakeRegisteredRaster(wildcardPath)
ST_SetPath(raster, band, string)
ST_MakeBandInDB(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);

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

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

Other functions

ST_AsBinary(raster, compression)
ST_RasterFromWKB(raster, [<srid>])
ST_RasterFromText(string, [<srid>])
ST_AsText(raster)

Accomplished objectives

Note: See TracWiki for help on using the wiki.