[postgis-devel] pgraster: st_resample error

[BJ Park]

On Thu, Dec 22, 2011 at 7:08 AM, Tom van Tilburg
<tom.van.tilburg at gmail.com> wrote:
> Hi all,
>
> Can anyone tell me if it's me or postgis making the mistake.
> Bottom line is that st_resample from epsg 28992 to 900913 gives me offsets
> of around 20 kilometers in Y direction.
>
> ===SQL=======
> WITH
> raster As (
>    SELECT  ST_AddBand(ST_MakeEmptyRaster(100, 100, 195000, 450000, 1, 1, 0,
> 0, 28992),'32BUI')  As rast
>    ),
> sub As (
>    SELECT ST_Transform(ST_SetSrid(ST_Extent(ST_Resample(rast,
> 900913)),900913),28992) geom FROM raster
>    ),
> sub2 As (
>    SELECT ST_Transform(ST_Transform(ST_SetSrid(ST_Extent(rast), 28992),
> 900913), 28992) geom FROM raster
> )
> SELECT 'st_resample' As class, ST_Extent(geom)::Text FROM sub
> UNION
> SELECT 'st_transform' As class, ST_Extent(geom)::Text From sub2
> UNION
> SELECT 'original' As class, St_Extent(rast)::Text FROM raster
>
> ===RESULT======
> "st_resample";"BOX(195165.982320333 429216.850820453,195268.212183354
> 429318.822637769)"
> "st_transform";"BOX(195000.000157893 450000.000438874,195100.000157985
> 450100.000438936)"
> "original";"BOX(195000 450000,195100 450100)"
>
> ===spatial_ref_sys====
> 28992: "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889
> +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel
> +towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812 +units=m
> +no_defs "
>
> 900913: "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0
> +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"

Tom,

I simplified your query to transform only from 28992 to 900913 instead
of back and forth and used ST_ConvexHull instead of ST_Extent.

{{{
WITH
raster As (
   SELECT  ST_AddBand(ST_MakeEmptyRaster(100, 100, 195000, 450000, 1,
1, 0, 0, 28992),'32BUI')  As rast
   ),
sub As (
   SELECT ST_ConvexHull(ST_Resample(rast, 900913)) geom FROM raster
   ),
sub2 As (
   SELECT ST_Transform(ST_ConvexHull(rast), 900913) geom FROM raster
)
SELECT 'st_resample' As class, ST_AsText(geom) FROM sub
UNION
SELECT 'st_transform' As class, ST_AsText(geom) From sub2
UNION
SELECT 'original' As class, ST_AsText(ST_ConvexHull(rast)) FROM raster
}}}

I would expect that sub and sub2 have the same polygon but they're
not. After checking the upperleft corner with the gdaltransform
utility, there definitely is something with ST_Resample.

I'll take a look at it later today.

-bborie