[postgis-devel] pgraster: st_resample error

[dustymugs]

On 12/22/2011 07:59 AM, BJ Park wrote:
> 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've done some digging around and the answers provided by ST_Resample 
match that provided by GDAL's gdalwarp utility.  I added some notes and 
testing results to the ticket #1402 regarding this.  I did test a PNG 
raster matching the spatial attributes of the raster created in PostGIS 
Raster and the answers are identical.

I suppose one test to check the outputs of GDAL and PostGIS Raster is to 
see how the warped output differs from a completely different program 
(ArcGIS, ERDAS, ENVI Imagine).

You could say that GDAL is incorrect but I'm not going to say that as 
that project has had years of inspections by users to make sure the 
outputs are correct.

-bborie