| 112 | |
| 113 | ---- |
| 114 | == '''Objective 0.1.6e - Being able to intersect vector and raster to produce vector.''' == |
| 115 | |
| 116 | List of PostGIS functions similar or related to ST_GetBBox(), ST_Envelope() and ST_Polygon(): |
| 117 | |
| 118 | [http://postgis.refractions.net/documentation/manual-1.4/ST_Boundary.html ST_Boundary(geometry)] (not really - always return a geometry a dimension lower - i.e. the boundary of a polygon is a polyline.) |
| 119 | |
| 120 | ST_box(geometry) (return a PostgreSQL box object) |
| 121 | |
| 122 | [http://postgis.refractions.net/documentation/manual-1.4/ST_Box2D.html ST_box2d(geometry)] (return a box2d object) |
| 123 | |
| 124 | [http://postgis.refractions.net/documentation/manual-1.4/ST_Box3D.html ST_box3d(geometry)] (return a box3d object) |
| 125 | |
| 126 | getBBOX(geometry) (Deprecation in 1.2.3) |
| 127 | |
| 128 | [http://postgis.refractions.net/documentation/manual-1.4/ST_Envelope.html ST_Envelope(geometry)] Returns the minimum bounding box for the supplied geometry, as a geometry. The polygon is defined by the corner points of the bounding box ((MINX, MINY), (MINX, MAXY), (MAXX, MAXY), (MAXX, MINY), (MINX, MINY)). (PostGIS will add a ZMIN/ZMAX coordinate as well). In PostGIS, the bounding box of a geometry is represented internally using float4s instead of float8s that are used to store geometries. The bounding box coordinates are floored, guarenteeing that the geometry is contained entirely within its bounds. This has the advantage that a geometry's bounding box is half the size as the minimum bounding rectangle, which means significantly faster indexes and general performance. But it also means that the bounding box is NOT the same as the minimum bounding rectangle that bounds the geometry. |
| 129 | |
| 130 | envelope(geometry) (Deprecation in 1.2.3) |
| 131 | |
| 132 | [http://postgis.refractions.net/documentation/manual-1.4/ST_Extent.html ST_extent(geometry set)] |
| 133 | |
| 134 | [http://postgis.refractions.net/documentation/manual-1.4/ST_ExteriorRing.html ST_ExteriorRing(geometry)] |
| 135 | |
| 136 | [http://postgis.refractions.net/documentation/manual-1.4/ST_ConvexHull.html ST_ConvexHull(geometry)] |
| 137 | |
| 138 | '''Function definitions for WKT Raster 0.1.6''' |
| 139 | |
| 140 | ST_GetBBox(raster) (not yet implemented but planned) is replaced with ST_Envelope() since there is no more an equivalent ST_GetBBox function in PostGIS. |
| 141 | |
| 142 | ST_Raster_to_box2d(raster) should be renamed ST_box2d(raster). |
| 143 | |
| 144 | '''ST_Box2D(raster) -> BOX2D''' - Returns a BOX2D representing the extent of the raster. |
| 145 | |
| 146 | If the raster is rotated, the result is a box2d enclosing the rotated raster. |
| 147 | |
| 148 | The only difference with ST_Envelope(raster) is that ST_Box2D(raster) returns a BOX2D, not a geometry. |
| 149 | |
| 150 | '''Implementation details''' |
| 151 | |
| 152 | This function is not anymore available in PostGIS 1.5 and is therefore probably not worth implementing. |
| 153 | |
| 154 | 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). |
| 155 | |
| 156 | '''ST_Envelope(raster) -> polygon geometry''' - Returns the minimum bounding box for the supplied raster, as a geometry. |
| 157 | |
| 158 | If the raster is rotated, the envelope is a non-rotated box enclosing the rotated raster. |
| 159 | |
| 160 | The only difference with ST_Box2D(raster) is that ST_Envelope(raster) returns a geometry, not a BOX2D. |
| 161 | |
| 162 | |
| 163 | '''Implementation details''' |
| 164 | |
| 165 | This function should be implemented as a SQL or a PL/pgSQL function looking like: 'SELECT ST_geometry(ST_Box2D(rast))' |
| 166 | |
| 167 | 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 |
| 168 | |
| 169 | [[Image(WKTRasterEnvelopeConvexHullAndPolygon.gif)]] |
| 170 | |
| 171 | '''ST_ConvexHull(raster) -> polygon geometry''' - Returns the minimum convex geometry that encloses every pixel from the raster. |
| 172 | |
| 173 | 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). |
| 174 | |
| 175 | '''Implementation details''' |
| 176 | |
| 177 | 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 |
| 178 | |
| 179 | '''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. |
| 180 | |
| 181 | This polygon geometry might contain holes if some internal areas of the raster contain pixels with NODATA values. |
| 182 | |
| 183 | Variant 1: ST_Polygon(raster) -> polygon geometry -- default to band # 1 |
| 184 | |
| 185 | '''Implementation details''' |
| 186 | |
| 187 | Could also be named ST_Outline(). |
| 188 | |
| 189 | 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. |
| 190 | |
| 191 | '''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. |
| 192 | |
| 193 | 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. |
| 194 | |
| 195 | Variant 1: ST_DumpAsPolygons(raster) -> geomval set -- default to band # 1 |
| 196 | |
| 197 | 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(). |
| 198 | |
| 199 | '''Implementation details''' |
| 200 | |
| 201 | This function is at the base of the first version of ST_Intersection which compute the intersection between a raster and a geometry. |
| 202 | |
| 203 | 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: |
| 204 | |
| 205 | CREATE TYPE geomval AS (geom geometry, val float8);[[BR]] |
| 206 | CREATE TYPE wktgeomval AS (wktgeom text, val float8); |
| 207 | |
| 208 | CREATE OR REPLACE FUNCTION ST_DumpAsPolygons(rast raster) RETURNS SETOF geomval AS $$[[BR]] |
| 209 | SELECT ST_GeomFromText(wktgeomval.wktgeom), wktgeomval.val FROM DumpAsWKTPolygons(%1) AS wktgeomval;[[BR]] |
| 210 | $$ LANGUAGE SQL; |
| 211 | |
| 212 | It should then be used like this: |
| 213 | |
| 214 | SELECT (gv).val, (gv).geom FROM (SELECT ST_DumpAsPolygons(rast) gv FROM sometable) foo |
| 215 | |
| 216 | [[Image(WKTRasterAsPolygon.gif)]] |
| 217 | |
| 218 | '''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. |
| 219 | |
| 220 | Variant 1: DumpAsWKTPolygons(raster) -> wktgeomval set -- default to band # 1 |
| 221 | |
| 222 | 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. |
| 223 | |
| 224 | 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(). |
| 225 | |
| 226 | DumpAsWKTPolygons() should not be used directly. Use ST_DumpAsPolygons() instead. |
| 227 | |
| 228 | '''Implementation details''' |
| 229 | |
| 230 | 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. |
| 231 | |
| 232 | It should be implemented as a C function passing the raster to GDAL and converting each polygon produced by the [http://www.gdal.org/gdal__alg_8h.html 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. |
| 233 | |
| 234 | |
| 235 | '''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). |
| 236 | |
| 237 | Variant 1: ST_Intersects(geometry, raster, integer) -> boolean -- the third parameter is the raster band number |
| 238 | |
| 239 | Variant 2: ST_Intersects(geometry, raster) -> boolean -- default to band # 1 |
| 240 | |
| 241 | Variant 3: ST_Intersects(raster, geometry) -> boolean -- default to band # 1 |
| 242 | |
| 243 | Variant 4: ST_Intersects(raster, raster) -> boolean |
| 244 | |
| 245 | This function permform a full intersection test and proceed in three steps to determine if the raster intersects with the geometry: |
| 246 | |
| 247 | 1) First, it checks if the bounding box of the raster intersects with the bounding box of the geometry. |
| 248 | |
| 249 | 2) If the first test returns TRUE, it checks if the geometry returned by ST_ConvexHull(raster) intersects with the geometry. |
| 250 | |
| 251 | 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. |
| 252 | |
| 253 | 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). |
| 254 | |
| 255 | 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): |
| 256 | |
| 257 | raster && geometry |
| 258 | |
| 259 | If you want to limit the intersection test to the second condition you can write: |
| 260 | |
| 261 | ST_Intersects(ST_ConvexHull(raster), geometry) |
| 262 | |
| 263 | '''Implementation details''' |
| 264 | |
| 265 | This function should be implemented as a pl/pgSQL function performing the three test described one after and conditionally to the other. |
| 266 | |
| 267 | It might be faster to skip test 2) if this test is not signicantly faster than test 3). |
| 268 | |
| 269 | '''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. |
| 270 | |
| 271 | Variant 1: |
| 272 | |
| 273 | ST_Intersection(raster, geometry) -> geometry -- default to band # 1[[BR]] |
| 274 | ST_Intersection(geometry, raster, integer) -> geometry -- the integer parameter is the band number of the raster[[BR]] |
| 275 | ST_Intersection(geometry, raster) -> geometry -- default to band # 1 |
| 276 | |
| 277 | Variant 2: |
| 278 | |
| 279 | ST_Intersection(raster, integer, geometry, 'raster') -> raster -- the integer parameter is the band number of the raster[[BR]] |
| 280 | ST_Intersection(raster, geometry, 'raster') -> raster -- default to band # 1 |
| 281 | |
| 282 | Variant 3: |
| 283 | |
| 284 | ST_Intersection(geometry, raster, integer, 'raster') -> raster -- the integer parameter is the band number of the raster[[BR]] |
| 285 | ST_Intersection(geometry, raster, 'raster') -> raster -- default to band # 1 |
| 286 | |
| 287 | Variant 4: |
| 288 | |
| 289 | ST_Intersection(raster, integer, raster, integer) -> raster -- the integer parameters are the band number of the rasters[[BR]] |
| 290 | ST_Intersection(raster, raster, integer) -> raster -- default first raster to band # 1[[BR]] |
| 291 | ST_Intersection(raster, integer, raster) -> raster -- default second raster to band # 1[[BR]] |
| 292 | ST_Intersection(raster, raster) -> raster -- default both rasters to band # 1 |
| 293 | |
| 294 | Variant 5: |
| 295 | |
| 296 | ST_Intersection(raster, integer, raster, integer, 'geometry') -> geometry[[BR]] |
| 297 | ST_Intersection(raster, raster, integer, 'geometry') -> geometry -- default first raster to band # 1[[BR]] |
| 298 | ST_Intersection(raster, integer, raster, 'geometry') -> geometry -- default second raster to band # 1[[BR]] |
| 299 | ST_Intersection(raster, raster, 'geometry') -> geometry -- default both raster to band # 1 |
| 300 | |
| 301 | 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. |
| 302 | |
| 303 | If the geometries do not share any space (are disjoint), then an empty geometry collection is returned. |
| 304 | |
| 305 | 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. |
| 306 | |
| 307 | 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: |
| 308 | |
| 309 | ST_Intersection(ST_ConvexHull(raster), geometry) |
| 310 | |
| 311 | 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: |
| 312 | |
| 313 | ST_Intersection(ST_Polygon(raster), geometry) |
| 314 | |
| 315 | '''Implementation details''' |
| 316 | |
| 317 | 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). |
| 318 | |
| 319 | Priority 1 should be given to the raster/geometry variants (original and variant 1) returning a geometry. |
| 320 | |
| 321 | Priority 2 should be given to the raster/raster variants (5) returning a geometry. |
| 322 | |
| 323 | 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. |
| 324 | |