[postgis-devel] pgraster: st_resample error

[BJ Park]

On Thu, Dec 22, 2011 at 7:55 AM, BJ Park <dustymugs at gmail.com> wrote:
> 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

I've filed a ticket #1402 for this problem.