[raster] Raster ST_Contains

[andrewsimpson]

Greetings,

I have loaded 1 band 8bit data to a raster table with tile size 500x500, and am trying to render through MapServer?/GDAL. Some of the tiles render, but most do not. After turning logging verbose on PostgreSQL I see that GDAL wktrasterraterband.cpp is calling for data a tile at a time, and many of the calls are not returning any data. Below is the metadata for one such tile, and the SQL that GDAL executes. It appears the ST_Contains thinks some part of the Box2D is outside the rast. I modified the wktrasterraterband.cpp IReadBlock function to buffer(Box2d,-1) and the all the tiles render.

I know this is not a good solution, any information on what I may do to remedy this would be great.

Thanks, Drew

PostgreSQL 8.4.4 on i686-redhat-linux-gnu, compiled by GCC gcc (GCC) 4.1.2 20080704 (Red Hat 4.1.2-48), 32-bit Postgis-1.5.1 WKTRaster-0.1.6d GDAL-1.7.2 Python-2.7 Proj-4.7.0 Geos-3.2.2

SELECT rid, (md).*, (bmd).*

FROM (SELECT rid, ST_Metadata(rast) AS md,

ST_BandMetadata(rast) AS bmd

FROM psi_tiled where rid in ( 87 )

) foo;

rid | upperleftx | upperlefty | width | height | pixelsizex | pixelsizey | skewx | skewy | srid | numbands | pixeltype | hasnodatavalue |

nodatavalue | isoutdb | path

---

+-------------+---------+------

87 | 504979.725998 | 2929894.254002 | 500 | 500 | 0.1016 | -0.1016 | 0 | 0 | 26917 | 1 | 8BUI | f |

0 | f | (1 row)

SELECT rid, rast FROM public.psi_tiled WHERE _ST_Contains(rast, ST_SetSRID( ST_MakeBox2D(ST_Point(504979.725998, 2929843.454002), ST_Point(505030.525998, 2929894.254002)), 26917))

rid | rast

---

(0 rows)

[jorgearevalo]

Hi Andrew,

ST_Contains considers a polygon doesn't contain its boundaries, and that's what's happening here, I think.

The BOX2D constructed to make the query is based on the pixel corners of each block (=row in raster table in tiled rasters) transformed in map coordinates using geotransform array. Making a buffer around this box could solve the problem, yes, but I think getting a tile per server round it's an improvable strategy.

I have a new version of GDAL driver, with some improvements, but this part hasn't be improved yet. Suggestions are welcome before committing it.

Many thanks for the report, Andrew.

[pracine]

This is a GDAL driver bug, no? Can we change the status?

[jorgearevalo]

Yes, we can. I've added this to the GDAL trac:  http://trac.osgeo.org/gdal/ticket/3569