2

I'm using PostGIS to calculate statistics on vector and raster data. The raster data is stored tables with one 128x128 pixel tile per row. I'm seeing an error when comparing clipping certain polygons to certain raster tiles. In a dataset of roughly 10 million polygons, the error occurs a few dozen times. Why is this error occurring? What's the fix?

I've created a minimal example that reproduces the error:

with a (rast) as (values('0100000100756F1820CBD6E23E1D701820CBD6E2BEA14D8090536055C09D958109C1BA3E4000000000000000000000000000000000AD10000001002200020001010101010101010101010101010101010101010101010101010101010101010101'::raster)),
     b (geom) as (select st_setsrid(st_geomfromtext('POLYGON((-85.5046457277349 30.7311595527075,-85.5046492027152 30.7294850816075,-85.5051002810699 30.7294907086655,-85.5050985678614 30.7303174886974,-85.5050968186801 30.7311624161489,-85.5046457277349 30.7311595527075))'), 4269))
select st_clip(rast, geom, true)
from a, b;

This yields the following error:

ERROR:  rt_raster_iterator: The set of rasters provided (custom extent included, if appropriate) do not have the same alignment
CONTEXT:  PL/pgSQL function st_clip(raster,integer[],geometry,double precision[],boolean) line 8 at RETURN

The example raster data was generated with ST_MakeEmptyRaster, and ST_AddBand, so I don't think the problem is related to corrupted data. Running other operations on this example raster (and the original data) works fine. For example, here's the output from ST_Summary:

Raster of 1x34 pixels has 1 band and extent of BOX(-85.5051003697358 30.7292025857055,-85.505091386583 30.7295080129021)
band 1 of pixtype 4BUI is in-db with no NODATA value

At least some of the time, the error goes away if a raster tile is unioned with neighboring tiles.

This occurs on both OS X and Windows, and in CARTO. PostGIS version info:

POSTGIS="2.4.3 r16312" PGSQL="100"
GEOS="3.6.2-CAPI-1.10.2 4d2925d" PROJ="Rel. 4.9.3, 15 August 2016"
GDAL="GDAL 2.2.3, released 2017/11/20"
LIBXML="2.7.8" LIBJSON="0.12" LIBPROTOBUF="1.2.1" RASTER
  • I am guessing it is an alignment issue. As your raster is only one pixel wide and the scale is 8.98315284119493e-06 and your polygon is 0.000454553335003993 wide, even though they do intersect, there is nothing to actually clip. I think the raster would have to be wider than the geometry to be clippable by it. – John Powell May 31 '18 at 19:55
  • @JohnPowellakaBarça The issue occurs with a bigger raster, too. I just made the smallest raster that would reproduce the issue. Clipping a small raster to a big polygon and vis versa generally works fine. – Neil F May 31 '18 at 19:57
  • Perhaps you could provide an example of a bigger raster demonstrating this problem too. Do you really require such levels of precision? Often, weird errors happen in Postgis due to arbitrary precision and rounding errors in intersections (with polygons and GEOS, anyway), which can often be solved by using ST_SnapToGrid or its raster sister, RT_ST_SnapToGrid – John Powell May 31 '18 at 20:01
  • Using ST_SnapToGrid on either the raster or the geometry appears to fix things! Thanks for the suggestion. If you want to make it an answer I'll accept it. – Neil F May 31 '18 at 20:31
  • Sure. I have seen similar questions before and it is only recently I made the (educated) guess that this was the reason :D. – John Powell May 31 '18 at 21:03
2

It is quite common to see errors like this due to rounding errors with coordinates with arbitrary precision. It happens with vector/vector operations (which come from GEOS), with errors like:

Self-intersection at or near point

and it happens with raster/vector operations, albeit with a different error, such as:

The set of rasters provided (custom extent included, if appropriate) do not have the same alignment.

In both cases, reducing the precision by using ST_SnapToGrid or the raster equivalent, RT_ST_SnapToGrid, often solves the problem.

I guessed this was your issue as your polygon is 0.000454553335003993 wide and your raster 34x1, with the 1 representing a scale of 8.98315284119493e-06, even though they do intersect and neither contains the other, as you can check with:

WITH 
  a (rast) AS (VALUES('0100000100756F1820CBD6E23E1D701820CBD6E2BEA14D8090536055C09D958109C1BA3E4000000000000000000000000000000000AD10000001002200020001010101010101010101010101010101010101010101010101010101010101010101'::raster)),
  b (geom) AS (
     SELECT ST_Setsrid(ST_Geomfromtext('POLYGON((-85.5046457277349 30.7311595527075,-85.5046492027152 30.7294850816075,-85.5051002810699 30.7294907086655,-85.5050985678614 30.7303174886974,-85.5050968186801 30.7311624161489,-85.5046457277349 30.7311595527075))'), 4269))
SELECT ST_Contains (geom, ST_Convexhull(rast))
  FROM a, b;

which returns false.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.