PostGIS WKT Raster Beta 0.1 Final Specifications

Quick Links

• WKT Raster Home Page
• Planning & Funding

• Beta 0.1.6 Working Specifications
• Beta 0.1.6 Documentation

• Beta 0.2.4 Working Specifications
• Working Specifications for Future Versions

---

Objective 0.1.6c - Being able to load rasters in the database

gdal2wktraster.py - A prototype of the translation utility implemented in Python and with use of ​GDAL and its bindings to Python. Create an SQL commands output to create a table of raster. As input raster file, all GDAL formats are accepted. The script is available as gdal2wktraster.py script.

USAGE:
gdal2wktraster.py -r <RASTERFILE> -t [<SCHEMA>].<TABLE> [<OPTIONS>]

-r <RASTERFILE> Specifies input raster file. Multiple -r options can be specified for a number of input files or wildcards can be used (? and *).
-t [<SCHEMA>.]<TABLE> Name of destination table in with or without target schema specified.

OPTIONS:

Raster processing: Optional parameters used to manipulate input raster dataset

-s <SRID> Set the SRID field of output raster. Default is -1.
-b <NBBAND> Specify the band number of the band to be extracted from input raster file.
-k <BLOCK_SIZE> Specify the size of the block of the input raster, assuming regular blocking mode. Must be specified as WIDTHxHEIGHT.
-R Register the raster as a filesystem (out-db) raster. Only the absolute path to the raster and its georeferencing informations are stored. Raster pixel values stay in the filesystem raster file and are not stored in the database.
-l <OVERVIEW_LEVEL> Create overview tables named as o_<LEVEL>_<RASTER_TABLE> and populate it with GDAL-provided overviews. Regular blocking only.

Database processing: Optional parameters used to manipulate database objects

-c Creates a new table and populates it with raster input file(s), this is the default if you do not specify any options.
-d Drops the table, then recreates it and populates it with current raster file(s) data.
-f <COLUMN> Name of target column for raster data. Default column name is rast.
-F Add a "filename" column containing the original name of the loaded raster file.
-I Create a GiST index on the raster column.
-M Issue VACUUM command against all generated tables.
-V Create RASTER_OVERVIEWS table used to store overviews metadata.

Miscellanous: Other optional parameters

-e <ENDIAN> Control endianness of generated binary output of raster; specify 0 for XDR (big-endian) and 1 for NDR (little-endian). Only NDR output is supported right now.
-w <VERSION> Specify version of WKT Raster protocol. Default is 0; only this value is supported right now.
-o <FILE> Output file for generated SQL commands. If not specified, stdout is assumed.
-v Switch on excessively verbose mode, useful for debugging.
-h Display this help screen.

---

Objective 0.1.6d - Being able to get all the properties of a raster (all the “Get” functions).

ST_SRID(raster|geometry) -> integer
Return the SRID associated with the raster.

ST_Width(raster) -> integer
Return the width of the raster.

ST_Height(raster) -> integer
Return the height of the raster.

ST_PixelSizeX(raster) -> float64
Return the georeference's X pixel size of the raster. ​See.

ST_PixelSizeY(raster) -> float64
Return the georeference's Y pixel size of the raster. ​See.

ST_RotationX(raster) -> float64
Return the georeference's X rotation.

ST_RotationY(raster) -> float64
Return the georeference's Y rotation. ​See.

ST_UpperLeftX(raster) -> float64
Return the georeference's X-coordinate of the upper left corner of the upper left pixel. ​See.

ST_UpperLeftY(raster) -> float64
Return the georeference's Y-coordinate of the upper left corner of the upper left pixel. ​See.

ST_Georeference(raster, text) -> string
Returns the georeference meta data in GDAL or ESRI format as commonly seen in a world file. For ESRI, ​see for GDAL, ​see.

ST_NumBands(raster) -> integer
Return the number of band included in the raster.

ST_BandPixelType(raster, integer) -> string
Return the pixel type of the specified 1-based Nth band of raster. Band index is 1-based. The function returns one of the following values:

• 1BB - 1-bit boolean
• 2BUI - 2-bit unsigned integer
• 4BUI - 4-bit unsigned integer
• 8BSI - 8-bit signed integer
• 8BUI - 8-bit unsigned integer
• 16BSI - 16-bit signed integer
• 16BUI - 16-bit unsigned integer
• 32BSI - 32-bit signed integer
• 32BUI - 32-bit unsigned integer
• 32BF - 32-bit float
• 64BF - 64-bit float

ST_BandNoDataValue(raster, integer) -> float32
Return the NoDataValue of the specified 1-based Nth band of raster. Band index is 1-based. The value is always returned as a float32 even if the pixel type is integer.

---

Objective 0.1.6e - Being able to intersect geometries and rasters to produce geometries .

ST_Raster_to_box2d(raster) should be renamed ST_box2d(raster).

ST_Box2D(raster) -> BOX2D - Returns a BOX2D representing the extent of the raster.

If the raster is rotated, the result is a box2d enclosing the rotated raster.

The only difference with ST_Envelope(raster) is that ST_Box2D(raster) returns a BOX2D, not a geometry.

Implementation details

This function is not anymore available in PostGIS 1.5 and is therefore probably not worth implementing.

This function replaces ST_Raster_to_box2d() and should be used for the "raster AS box2d" cast. It should also be used instead of ST_raster_envelope() when creating indexes in gdal2wktraster.py (to be tested).

ST_ConvexHull(raster) -> polygon geometry - Returns the minimum convex geometry that encloses every pixel from the raster.

The resulting polygon DOES NOT TAKE the NODATA values into account. The resulting polygon enclose every pixel, even those containing NODATA values. The resulting polygon can be rotated or not depending on the rotation of the raster. If the raster is not rotated ST_ConvexHull(geometry) is (almost) equivalent to ST_Envelope(raster) (ST_Envelope floor the coordinates and hence add a little buffer around the raster).

Implementation details

This function is actually already implemented by the ST_raster_envelope() function which is wrongly named. There is also a ticket (#348) related to this function: http://trac.osgeo.org/postgis/ticket/348

ST_Envelope(raster) -> polygon geometry - Returns the minimum bounding box for the supplied raster, as a geometry.

If the raster is rotated, the envelope is a non-rotated box enclosing the rotated raster.

The only difference with ST_Box2D(raster) is that ST_Envelope(raster) returns a geometry, not a BOX2D.

Implementation details

This function should be implemented as a SQL or a PL/pgSQL function looking like: 'SELECT ST_envelope(ST_ConvexHull(raster))'

There is already a ST_raster_envelope() function but this has actually the behavior of the ST_ConvexHull(raster) described below. There is also a ticket (#348) related to this function: http://trac.osgeo.org/postgis/ticket/348

ST_Polygon(raster, integer) -> polygon geometry - Returns a geometry encomparsing every pixel having a significant value (different than the NODATA value) for the provided band.

This polygon geometry might contain holes if some internal areas of the raster contain pixels with NODATA values.

Variant 1: ST_Polygon(raster) -> polygon geometry -- default to band # 1

Implementation details

Could also be named ST_Outline().

This function could be roughly implemented as a SQL function looking like 'SELECT ST_Collect((ST_DumpAsPolygons(raster)).geom)'. For sure a more specialised version could be faster than ST_AsPolygon.

ST_DumpAsPolygons(raster, integer) -> geomval set - Returns a set of "geomval" value, one for each contiguous group of pixel sharing the same value for the provided band.

This is a set-returning function (SRF). A "geomval" value is a complex type composed of a polygon geometry (one for each contiguous group of pixel sharing the same value) and the value associated with this geometry. These values are always returned as a value of type double precision. The shape of each polygon follow pixels edges.

Variant 1: ST_DumpAsPolygons(raster) -> geomval set -- default to band # 1

This function should be used with precaution on raster with float pixel type as it could return one "geomval" (or polygon) per pixel. This kind of raster should be reclassified (using the planned ST_Reclassify(raster,...) function) before using ST_DumpAsPolygons().

Implementation details

This function is at the base of the first version of ST_Intersection which compute the intersection between a raster and a geometry.

To avoid linking directly with PostGIS (see "Why avoid to link with PostGIS?" below) It should be implemented as a SQL wrapper around DumpAsWKTPolygons() looking something like this:

CREATE TYPE geomval AS (geom geometry, val float8);
CREATE TYPE wktgeomval AS (wktgeom text, val float8);

CREATE OR REPLACE FUNCTION ST_DumpAsPolygons(rast raster) RETURNS SETOF geomval AS $$

SELECT ST_GeomFromText(wktgeomval.wktgeom), wktgeomval.val FROM DumpAsWKTPolygons(%1) AS wktgeomval;

$$ LANGUAGE SQL;

It should then be used like this:

SELECT (gv).val, (gv).geom FROM (SELECT ST_DumpAsPolygons(rast) gv FROM sometable) foo

DumpAsWKTPolygons(raster, integer) -> wktgeomval set - Returns a set of "geomval" value, one for each group of pixel sharing the same value for the provided band.

Variant 1: DumpAsWKTPolygons(raster) -> wktgeomval set -- default to band # 1

This is a set-returning function (SRF). A "wktgeomval " value is a complex type composed of a the wkt representation of a geometry (one for each group of pixel sharing the same value) and the value associated with this geometry. These values are always returned as a value of type double precision. The shape of each polygon follow pixels edges. Areas with NODATA values are not included in the resulting set.

This function should be used with precaution on raster with float pixel type as it could return one "geomval" (or polygon) per pixel. This kind of raster should be reclassified (using the planned ST_Reclassify(raster,...) function) before using DumpAsWKTPolygons().

DumpAsWKTPolygons() should not be used directly. Use ST_DumpAsPolygons() instead.

Implementation details

This function construct the WKT strings representing the geometries grouping pixels of a raster having the same value. It does not directly construct geometries and therefore prevent WKT Raster from having to link with liblwgeom.a (see "Why avoid to link with PostGIS?" below) It should also return the value associated with this WKT polygon.

It should be implemented as a C function passing the raster to GDAL and converting each polygon produced by the ​GDALPolygonize function to a WKT string accompagned by its corresponding value in a wktgeomval type. See ​http://www.postgresql.org/docs/8.4/interactive/xfunc-c.html#AEN44970 on how to return rows and sets in a C function.

ST_Intersects(raster, integer, geometry) -> boolean - Returns TRUE if the geometry and the raster "spatially intersect" - (share any portion of space) and FALSE if they don't (they are disjoint).

Variant 1: ST_Intersects(geometry, raster, integer) -> boolean -- the third parameter is the raster band number

Variant 2: ST_Intersects(geometry, raster) -> boolean -- default to band # 1

Variant 3: ST_Intersects(raster, geometry) -> boolean -- default to band # 1

Variant 4: ST_Intersects(raster, raster) -> boolean

This function permform a full intersection test and proceed in three steps to determine if the raster intersects with the geometry:

1) First, it checks if the bounding box of the raster intersects with the bounding box of the geometry.

2) If the first test returns TRUE, it checks if the geometry returned by ST_ConvexHull(raster) intersects with the geometry.

3) If the second test returns TRUE, it checks if the geometry returned by ST_Polygon(raster, integer) intersects with the geometry. This test is slower since it involve the computation of the raster shape and it might involve the geometry shape. This test takes NODATA values into account. i.e. A geometry intersecting only with a NODATA value area of a raster is NOT actually not intersecting with this raster. This behavior may be controled by limiting the test to certain conditions as stated below.

Variant 4 proceeds in a very similar way except that convex hulls of both rasters are computed and compared in step 2) and both shape of raster are computed and compared in step 3).

If you want to limit the intersection test to the first condition, simply use the && operator. The raster and the geometry will be casted to their respective bounding box (box2d objects):

raster && geometry

If you want to limit the intersection test to the second condition you can write:

ST_Intersects(ST_ConvexHull(raster), geometry)

Implementation details

This function should be implemented as a pl/pgSQL function performing the three test described one after and conditionally to the other.

It might be faster to skip test 2) if this test is not signicantly faster than test 3).

ST_Intersection(raster, integer, geometry) -> geometry - Returns a set of "geomval" value representing the shared portion of the geometry and the raster band areas sharing a common meaningfull values. The integer parameter is the band number of the raster.

Variant 1:

ST_Intersection(raster, geometry) -> geometry -- default to band # 1
ST_Intersection(geometry, raster, integer) -> geometry -- the integer parameter is the band number of the raster
ST_Intersection(geometry, raster) -> geometry -- default to band # 1

Variant 2:

ST_Intersection(raster, integer, geometry, 'raster') -> raster -- the integer parameter is the band number of the raster
ST_Intersection(raster, geometry, 'raster') -> raster -- default to band # 1

Variant 3:

ST_Intersection(geometry, raster, integer, 'raster') -> raster -- the integer parameter is the band number of the raster
ST_Intersection(geometry, raster, 'raster') -> raster -- default to band # 1

Variant 4:

ST_Intersection(raster, integer, raster, integer) -> raster -- the integer parameters are the band number of the rasters
ST_Intersection(raster, raster, integer) -> raster -- default first raster to band # 1
ST_Intersection(raster, integer, raster) -> raster -- default second raster to band # 1
ST_Intersection(raster, raster) -> raster -- default both rasters to band # 1

Variant 5:

ST_Intersection(raster, integer, raster, integer, 'geometry') -> geometry
ST_Intersection(raster, raster, integer, 'geometry') -> geometry -- default first raster to band # 1
ST_Intersection(raster, integer, raster, 'geometry') -> geometry -- default second raster to band # 1
ST_Intersection(raster, raster, 'geometry') -> geometry -- default both raster to band # 1

This is a set-returning function (SRF). It returns a set of "geomval" representing the point set intersection of the geometry and the areas of the raster sharing a common meaningfull value (NODATA values ARE taken into account). The raster is first polygonised using ST_DumpAsPolygons(raster) and ST_Intersection(geometry, geometry) is then performed between the provided geometry and all the geometries resulting from the polygonisation of the raster.

If the geometries do not share any space (are disjoint), then an empty geometry collection is returned.

ST_Intersection in conjunction with ST_Intersects is very useful for clipping geometries such as in bounding box, buffer, region queries where you only want to return that portion of a geometry that sits in a country or region of interest.

If you only want to compute the intersection between the convex hull of the raster without polygonising it to group of pixels sharing a common value, do:

ST_Intersection(ST_ConvexHull(raster), geometry)

If you only want to compute the intersection between the shape of the raster without polygonising it to group of pixels sharing a common value, do:

ST_Intersection(ST_Polygon(raster), geometry)

Implementation details

This function should be implemented as a pl/pgSQL function (possibly only a SQL function) performing the intersection between the provided geometry and the table generated by ST_DumpAsPolygons(raster).

Priority 1 should be given to the raster/geometry variants (original and variant 1) returning a geometry.

Priority 2 should be given to the raster/raster variants (5) returning a geometry.

Priority 3 should be given to the other variants (2, 3 & 4) returning a raster. Variant 2 & 3 should be implemented after the planned ST_AsRaster(geometry) function. What should be the raster result of the intersection of two rasters is still to be specified.

---

Objective 0.1.6i - Being able to quickly get metadata for raster and band.

ST_Metadata(raster) -> record
Returns all the metadata associated with a raster (upper left x, upper left y, width, height, pixelsize x, pixelsize y, skew x, skew y, srid, numbands in this order). Does not include any band metadata. If you want band metadata, add ST_BandMetadata(raster, integer) to your query.

It should be used like this: "SELECT (md).* FROM (SELECT ST_Metadata(rast) AS md FROM myrastertable) foo"

ST_BandMetadata(raster, integer) -> record
Returns all the metadata associated with a raster band (pixeltype, has nodata value, nodata value, is stored out-db, path, in this order). Does not include embedding raster metadata. If you want embedding raster metadata, add ST_Metadata(raster) to your query.

It should be used like this: "SELECT (bmd).* FROM (SELECT ST_BandMetadata(rast) AS bmd FROM myrastertable) foo"

---

RASTER_COLUMNS Metadata Table

The following metadata table provides external applications and possibly internal tools the ability to discover tables with raster data, and to learn various information about those datasets.

Column Name	Type	Constraints	Comments
r_table_catalog	character varying(256)	NOT NULL
r_table_schema	character varying(256)	NOT NULL
r_table_name	character varying(256)	NOT NULL
r_column	character varying(256)	NOT NULL
srid	integer	NOT NULL	Use 0 for unknown SRID
pixel_types	ARRAY=VARCHAR=	NOT NULL	An array of pixeltypes, one per band. The band count is implicit in the size of this array.
out_db	ARRAY=BOOLEAN=	NOT NULL	An array of boolean values, one per band. "false" if internal tiles, "true" if tiles are references to raster files outside the database.
regular_blocking	boolean	NOT NULL	False by default, true if all blocks are equal sized, abutted and non-overlapping, started at top left origin (see below).
nodata_values	ARRAY=DOUBLE=		An array of nodata values, one per band. The entry may be NULL to indicate no nodata values.
pixelsize_x	double		Width of a pixel in geounits (per SRID).
pixelsize_y	double		Height of a pixel in geounits (per SRID).
blocksize_x	integer		The width of a block in pixels (NULL if irregular).
blocksize_y	integer		The height of a block in pixels (NULL if irregular).
extent	GEOMETRY		A polygon geometry containing all raster tiles, or NULL if predefined bounds are not known. For "regular_blocking" cases this geometry will be a simple rectangle. In other cases it might be an irregular polygon.

If the regular_blocking field is true a number of restrictions are placed on the raster column that is defined:

1. All tiles must have the same size (blocksize_x and blocksize_y).
2. All tiles must be non-overlapping, and appear on regular block grid.
3. The top left block must start at the top left corner of the extent.
4. The right most column, and bottom row of blocks may have portions that extend beyond the raster extent. These areas will be assumed to be nodata and not part of the described raster.
5. The extent field must be a simple rectangle (non-rotated).

It is permissable to for regular_blocking rasters to omit some blocks/tiles (sparse rasters) in which case the missing tiles should be assumed to be all nodata or zero valued.

---