Opened 6 months ago

Last modified 6 weeks ago

#4051 new defect

Strange result with raster ST_Intersects

Reported by: stev00 Owned by: Bborie Park
Priority: medium Milestone: PostGIS 2.2.8
Component: raster Version: 2.4.x
Keywords: Cc:

Description

I have odd behaviour when trying to use the ST_Intersects function with rasters, which I am using to filter raster tiles using an index for an ST_Value call.

On both Linux Red Hat (RDS) and Windows using Postgres 10.1 and Postgis 2.4 I get a null result for the filter call for a very small number of point geometries (in fact only one so far from around 1000 points). The only significant feature is that the x value of this point is a whole integer.

I am attaching a full sql file to reproduce a minimal version of the problem. If you look at the rasters in QGIS you should be able to see it covers the point, and adding a very small increment or decrement to the X value means a value is returned (but no difference with y value).

This combined with the whole number issue above means I wonder if it is something to do with the integer value.

I will continue looking into this, and it could obviously be a dependent library issue, but would be interested to see if others can reproduce it.

Attachments (1)

raster_intersects_potential_bug.sql (178.2 KB) - added by stev00 6 months ago.
minimally reproduce potential bug

Download all attachments as: .zip

Change History (5)

Changed 6 months ago by stev00

minimally reproduce potential bug

comment:1 Changed 6 months ago by robe

Component: postgisraster
Owner: changed from pramsey to Bborie Park

comment:2 Changed 6 months ago by stev00

I have done some more tests with the above data. The SQL commands run and results are shown below. I guess this is related to precision - it seems that ST_Worldtorastercoord gives different results to ST_Rastertoworldcoord. Will try and do a bit more digging.

with r as(
SELECT * from test_raster WHERE rast && ST_SetSRID(ST_GeomFromText('POINT(529004 182915.87)'),27700) AND fn_st_intersects_log(rast, ST_SetSRID(ST_GeomFromText('POINT(529004 182915.87)'),27700), 1) 
)

-- result columnx=51, rowy=26
select (ST_worldtorastercoord(rast, 529004, 182915.87)).* from r 

--result st_rastertoworldcoordx=529003.000000106 st_width=50 st_height=56 st_pixelwidth=1.00000000000698
select ST_RasterToWorldCoordX(rast, 50),st_width(rast), st_height(rast), ST_PixelWidth(rast) from r 

comment:3 Changed 6 months ago by robe

Milestone: PostGIS 2.4.4PostGIS 2.2.8

As much as I'd like to fix this for 2.4.4, I don't have enough time to investigate the issue before release or even if there is one.

I did do enough to determine that yes I got no ST_Value answer back for the ones you say don't work , but haven't confirmed if that's an issue or not or just there is no value for that part.

If it's an issue it's an issue in 2.2 as well so pushing to 2.2

Note: See TracTickets for help on using tickets.