105 | | 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. |
106 | | |
107 | | 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: |
108 | | |
109 | | SELECT ST_MapAlgebraNgb(rast, band, pixeltype, "ST_Mean", 2, 2, "ignore") |
110 | | |
111 | | 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." |
112 | | |
113 | | 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: |
| 105 | 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. |
| 106 | |
| 107 | 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: |
| 108 | |
| 109 | SELECT ST_MapAlgebraNgb(rast, band, pixeltype, "ST_Mean", 2, 2, "ignore") |
| 110 | |
| 111 | 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." |
| 112 | |
| 113 | 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: |
120 | | 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). |
121 | | |
122 | | 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... |
123 | | |
124 | | Users could write their own map algebra summarizing functions. |
125 | | |
126 | | 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: |
127 | | |
128 | | SELECT ST_MapAlgebraNgb(rast, band, pixeltype, 2, 2, "ST_Mean", "ignore") |
129 | | |
130 | | and the dimensions do not have to be passed to the summarizing functions since it could deduce them from ST_Width & ST_Height. |
131 | | |
132 | | 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.: |
133 | | |
134 | | SELECT ST_MapAlgebraNgb("mycoveragetable", "myrastercolumn", band, pixeltype, 2, 2, "ST_Mean", "ignore") |
135 | | |
136 | | Three difficulties must be solved to implement this function: |
| 120 | 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). |
| 121 | |
| 122 | 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... |
| 123 | |
| 124 | Users could write their own map algebra summarizing functions. |
| 125 | |
| 126 | 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: |
| 127 | |
| 128 | SELECT ST_MapAlgebraNgb(rast, band, pixeltype, 2, 2, "ST_Mean", "ignore") |
| 129 | |
| 130 | and the dimensions do not have to be passed to the summarizing functions since it could deduce them from ST_Width & ST_Height. |
| 131 | |
| 132 | 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.: |
| 133 | |
| 134 | SELECT ST_MapAlgebraNgb("mycoveragetable", "myrastercolumn", band, pixeltype, 2, 2, "ST_Mean", "ignore") |
| 135 | |
| 136 | Three difficulties must be solved to implement this function: |
148 | | These function must first determine if a resampling of one raster is necessary: |
149 | | |
150 | | '''ST_SameAlignment(raster, raster)''' |
151 | | |
152 | | This function returns true if both rasters' grids are aligned. |
153 | | |
154 | | Two rasters grid are aligned if: |
155 | | |
156 | | -They share the same pixel scales and skews |
157 | | |
158 | | -At least one of any of the four corner of any pixel of one raster fall on any corner of the grid of the other raster. |
159 | | |
160 | | Alignment is not the same concept as georeference. Two rasters can be aligned but not have the same georeference. |
161 | | |
162 | | Rotation is important here since two rasters grid might look perfectly aligned but are not if their rotated are different. |
163 | | |
164 | | Knowing if two rasters share the same alignment is useful when it is time to decide if one of them must be resampled before doing other operations (like ST_MapAlgebra on two rasters). |
165 | | |
166 | | '''Variants''' |
167 | | |
168 | | 1) ST_SameAlignment(ulx1, uly1, scalex1, scaley1, skewx1, skewy1, ulx2, uly2, scalex2, scaley2, skewx2, skewy2) |
169 | | |
170 | | 2) ST_SameAlignment(rast1 raster, rast2 raster) |
171 | | |
172 | | The first variant is useful to PL/pgSQL function which have already get the values of the parameters. |
173 | | |
174 | | '''Implementation Details''' |
175 | | |
176 | | Only the first variant should be implemented in C. The second one is a PL/pgSQL variants. The C implementation should follow the PL/pgSQL version implemented in http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_mapalgebra.sql |
177 | | |
178 | | It is not clear if this PL/pgSQL implementation works when raster are rotated. To verify. |
179 | | |
180 | | See discussion in http://trac.osgeo.org/postgis/ticket/589 |
181 | | |
182 | | |
183 | | '''ST_MapAlgebra on two rasters''' |
184 | | |
185 | | '''ST_MapAlgebraExpr(raster rast, integer band, text expression, text nodatavalueexpr, text pixeltype) -> raster''' |
186 | | |
187 | | ST_MapAlgebra performs the specified mathematical expression on the input raster, returning a raster. Both rasters must have the same SRID. If both rasters do not have the same SRID, ST_MapAlgebra should return an error: "ERROR: Operation on two geometries with different SRIDs" |
188 | | |
189 | | The first raster passed to ST_MapAlgebra is the 'master' raster, unless either: |
190 | | - one (1) additional parameter specifies the index (in the parameter list) of the 'master' raster.[[BR]] |
191 | | - one (1) additional parameter specifies a raster whose origin and cell size should be used to compute the output raster.[[BR]] |
192 | | - four (4) additional parameters are passed specifying the origin, cell size, and raster rotation. |
193 | | |
194 | | This function implicitly calls ST_Intersects(raster, raster]) to detect if the rasters are overlapping before performing any computation. Additionally, the resulting raster will have the same extent as the result of ST_Intersection(raster, integer, geometry). |
195 | | |
196 | | One of the implications of the ST_Intersects inclusion is that: |
| 148 | Both rasters must have the same SRID. If not, ST_MapAlgebra should return an error: "ERROR: Operation on two geometries with different SRIDs" |
| 149 | |
| 150 | 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 an optional 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." |
| 151 | |
| 152 | Both raster must overlap. This is determined using ST_Intersects(raster, raster) (to be implemented). If they don't overlap an empty raster is returned. |
| 153 | |
| 154 | One of the implications of this is that users should add a WHERE clause to avoid useless computation with non-overlapping rasters. This is the same principle as with ST_Intersection(). Hence: |
| 178 | |
| 179 | '''ST_SameAlignment(raster, raster)''' |
| 180 | |
| 181 | This function returns true if both rasters' grids are aligned. |
| 182 | |
| 183 | Two rasters grid are aligned if: |
| 184 | |
| 185 | -They share the same pixel scales and skews |
| 186 | |
| 187 | -At least one of any of the four corner of any pixel of one raster fall on any corner of the grid of the other raster. |
| 188 | |
| 189 | Alignment is not the same concept as georeference. Two rasters can be aligned but not have the same georeference. |
| 190 | |
| 191 | Rotation is important here since two rasters grid might look perfectly aligned but are not if their rotated are different. |
| 192 | |
| 193 | Knowing if two rasters share the same alignment is useful when it is time to decide if one of them must be resampled before doing other operations (like ST_MapAlgebra on two rasters). |
| 194 | |
| 195 | '''Variants''' |
| 196 | |
| 197 | 1) ST_SameAlignment(ulx1, uly1, scalex1, scaley1, skewx1, skewy1, ulx2, uly2, scalex2, scaley2, skewx2, skewy2) |
| 198 | |
| 199 | 2) ST_SameAlignment(rast1 raster, rast2 raster) |
| 200 | |
| 201 | The first variant is useful to PL/pgSQL function which have already get the values of the parameters. |
| 202 | |
| 203 | '''Implementation Details''' |
| 204 | |
| 205 | Only the first variant should be implemented in C. The second one is a PL/pgSQL variants. The C implementation should follow the PL/pgSQL version implemented in http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_mapalgebra.sql |
| 206 | |
| 207 | It is not clear if this PL/pgSQL implementation works when raster are rotated. To verify. |
| 208 | |
| 209 | See discussion in http://trac.osgeo.org/postgis/ticket/589 |