Opened 14 years ago

Closed 12 years ago

Last modified 11 years ago

#3569 closed enhancement (fixed)

PostGIS Raster driver: Check query that fetchs data from db

Reported by: jorgearevalo Owned by: jorgearevalo
Priority: highest Milestone:
Component: GDAL_Raster Version: svn-trunk
Severity: normal Keywords: postgis_raster
Cc:

Description

Currently, in IReadBlock, the query to fetch a block of data from database is:

SELECT <rast_column> FROM <table> WHERE <rast_column> ~ ST_SetSRID(ST_MakeBox2D(ST_Point(x1, y1), ST_Point(x2, y2)), srid)

We should use "&&" (=intersects) instead of "~" (=contains), because "~" will only work with regularly blocked rasters, but "&&" doesn't seem to work with this type of rasters just now (r19624)

Change History (9)

comment:1 by jorgearevalo, 13 years ago

Status: newassigned
Summary: WKT Raster driver: Check query that fetchs data from dbPostGIS Raster driver: Check query that fetchs data from db

Update: Andrew Simpson sent a bug report to PostGIS trac on this topic:

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

comment:3 by jorgearevalo, 13 years ago

Keywords: postgis added; wkt removed

comment:4 by jorgearevalo, 13 years ago

Priority: normalhigh
Type: defectenhancement

The query has changed. Now I'm implementing IRasterIO. The new query looks something like this:

SELECT rid, <rast_column>, ST_Width(<rast_column>), ST_Height(<rast_column>) FROM <table> WHERE ST_Intersects(<rast_column>, ST_PolygonFromText('<envelope>', <srid>))

The <envelope> is built in this way (IRasterIO function, at Dataset level):

    GetGeoTransform(adfTransform);
    ulx = nXOff;
    uly = nYOff;
    lrx = nXOff + nXSize;
    lry = nYOff + nYSize;

    
    /* Calculate the intersection polygon */
    adfProjWin[0] = adfTransform[0] +
        ulx * adfTransform[1] +
        uly * adfTransform[2];
    
    adfProjWin[1] = adfTransform[3] +
        ulx * adfTransform[4] +
        uly * adfTransform[5];

    adfProjWin[2] = adfTransform[0] +
        lrx * adfTransform[1] +
        uly * adfTransform[2];
 
    adfProjWin[3] = adfTransform[3] + 
        lrx * adfTransform[4] +
        uly * adfTransform[5];
    
    adfProjWin[4] = adfTransform[0] +
        lrx * adfTransform[1] +
        lry * adfTransform[2];

    adfProjWin[5] = adfTransform[3] +
        lrx * adfTransform[4] +
        lry * adfTransform[5];

    adfProjWin[6] = adfTransform[0] +
        ulx * adfTransform[1] +
        lry * adfTransform[2];

    adfProjWin[7] = adfTransform[3] +
        ulx * adfTransform[4] +
        lry * adfTransform[5];

    osCommand.Printf("SELECT rid, %s, ST_Width(%s), ST_Height(%s) FROM %s.%s WHERE
 ST_Intersects(%s, ST_PolygonFromText('POLYGON((%f %f, %f %f, %f %f, %f %f, %f 
%f))', %d))", pszColumn, pszColumn, pszColumn, pszSchema, pszTable, pszColumn, 
adfProjWin[0], adfProjWin[1], adfProjWin[2], adfProjWin[3], adfProjWin[4], 
adfProjWin[5], adfProjWin[6], adfProjWin[7], adfProjWin[0], adfProjWin[1], nSrid);

comment:5 by jorgearevalo, 13 years ago

Milestone: 1.8.11.8.2

comment:6 by jorgearevalo, 13 years ago

Priority: highhighest

comment:7 by jorgearevalo, 12 years ago

Resolution: fixed
Status: assignedclosed

The query is working fine. I resolve this case.

comment:8 by pracine, 12 years ago

Keywords: postgis_raster added; postgis raster removed

comment:9 by Even Rouault, 11 years ago

Milestone: 1.8.2

Milestone 1.8.2 deleted

Note: See TracTickets for help on using tickets.