| 3231 | | -- Get the intersection between with the geometry. |
| 3232 | | -- We will search for withvalue pixel only in this area. |
| 3233 | | geomintersect := st_intersection(geom, convexhull); |
| 3234 | | |
| 3235 | | --RAISE NOTICE 'geomintersect=%', st_astext(geomintersect); |
| 3236 | | |
| 3237 | | -- If the intersection is empty, return false |
| 3238 | | IF st_isempty(geomintersect) THEN |
| 3239 | | RETURN FALSE; |
| 3240 | | END IF; |
| 3241 | | |
| 3242 | | -- We create a minimalistic buffer around the intersection in order to scan every pixels |
| 3243 | | -- that would touch the edge or intersect with the geometry |
| 3244 | | SELECT sqrt(scalex * scalex + skewy * skewy), width, height INTO scale, w, h FROM ST_Metadata(rast); |
| 3245 | | IF scale != 0 THEN |
| 3246 | | geomintersect := st_buffer(geomintersect, scale / 1000000); |
| 3247 | | END IF; |
| 3248 | | |
| 3249 | | --RAISE NOTICE 'geomintersect2=%', st_astext(geomintersect); |
| 3250 | | |
| 3251 | | -- Find the world coordinates of the bounding box of the intersecting area |
| 3252 | | x1w := st_xmin(geomintersect); |
| 3253 | | y1w := st_ymin(geomintersect); |
| 3254 | | x2w := st_xmax(geomintersect); |
| 3255 | | y2w := st_ymax(geomintersect); |
| 3256 | | nodata := st_bandnodatavalue(rast, nband); |
| 3257 | | |
| 3258 | | --RAISE NOTICE 'x1w=%, y1w=%, x2w=%, y2w=%', x1w, y1w, x2w, y2w; |
| 3259 | | |
| 3260 | | -- Convert world coordinates to raster coordinates |
| 3261 | | x1 := st_world2rastercoordx(rast, x1w, y1w); |
| 3262 | | y1 := st_world2rastercoordy(rast, x1w, y1w); |
| 3263 | | x2 := st_world2rastercoordx(rast, x2w, y1w); |
| 3264 | | y2 := st_world2rastercoordy(rast, x2w, y1w); |
| 3265 | | x3 := st_world2rastercoordx(rast, x1w, y2w); |
| 3266 | | y3 := st_world2rastercoordy(rast, x1w, y2w); |
| 3267 | | x4 := st_world2rastercoordx(rast, x2w, y2w); |
| 3268 | | y4 := st_world2rastercoordy(rast, x2w, y2w); |
| 3269 | | |
| 3270 | | --RAISE NOTICE 'x1=%, y1=%, x2=%, y2=%, x3=%, y3=%, x4=%, y4=%', x1, y1, x2, y2, x3, y3, x4, y4; |
| 3271 | | |
| 3272 | | -- Order the raster coordinates for the upcoming FOR loop. |
| 3273 | | x1 := int4smaller(int4smaller(int4smaller(x1, x2), x3), x4); |
| 3274 | | y1 := int4smaller(int4smaller(int4smaller(y1, y2), y3), y4); |
| 3275 | | x2 := int4larger(int4larger(int4larger(x1, x2), x3), x4); |
| 3276 | | y2 := int4larger(int4larger(int4larger(y1, y2), y3), y4); |
| 3277 | | |
| 3278 | | -- Make sure the range is not lower than 1. |
| 3279 | | -- This can happen when world coordinate are exactly on the left border |
| 3280 | | -- of the raster and that they do not span on more than one pixel. |
| 3281 | | x1 := int4smaller(int4larger(x1, 1), w); |
| 3282 | | y1 := int4smaller(int4larger(y1, 1), h); |
| 3283 | | |
| 3284 | | -- Also make sure the range does not exceed the width and height of the raster. |
| 3285 | | -- This can happen when world coordinate are exactly on the lower right border |
| 3286 | | -- of the raster. |
| 3287 | | x2 := int4smaller(x2, w); |
| 3288 | | y2 := int4smaller(y2, h); |
| 3289 | | |
| 3290 | | --RAISE NOTICE 'x1=%, y1=%, x2=%, y2=%', x1, y1, x2, y2; |
| 3291 | | |
| 3292 | | -- Search exhaustively for withvalue pixel on a moving 3x3 grid |
| 3293 | | -- (very often more efficient than searching on a mere 1x1 grid) |
| 3294 | | FOR xinc in 0..2 LOOP |
| 3295 | | FOR yinc in 0..2 LOOP |
| 3296 | | FOR x IN x1+xinc..x2 BY 3 LOOP |
| 3297 | | FOR y IN y1+yinc..y2 BY 3 LOOP |
| 3298 | | -- Check first if the pixel intersects with the geometry. Often many won't. |
| 3299 | | bintersect := NOT st_isempty(st_intersection(st_pixelaspolygon(rast, x, y), geom)); |
| 3300 | | |
| 3301 | | IF bintersect THEN |
| 3302 | | -- If the pixel really intersects, check its value. Return TRUE if with value. |
| 3303 | | pixelval := st_value(rast, nband, x, y); |
| 3304 | | IF pixelval != nodata THEN |
| 3305 | | RETURN TRUE; |
| 3306 | | END IF; |
| 3307 | | END IF; |
| 3308 | | END LOOP; |
| 3309 | | END LOOP; |
| 3310 | | END LOOP; |
| 3311 | | END LOOP; |
| | 3208 | -- get band surface |
| | 3209 | surface := ST_BandSurface(rast, nband); |
| | 3210 | |
| | 3211 | IF surface IS NOT NULL THEN |
| | 3212 | RETURN ST_Intersects(geom, surface); |
| | 3213 | END IF; |
