[raster] ST_Resample producing incorrect transformation

[dustymugs]

As per email ​http://postgis.refractions.net/pipermail/postgis-devel/2011-December/016960.html, ST_Resample is returning raster with incorrect coordinates.

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

Double checking the upper-left corner coordinates with gdaltrasform proves that something is amiss.

{{{gdaltransform -s_srs EPSG:28992 -t_srs "+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" 195000 450000 664595.851266449 6806797.3239794 43.5106153006745 }}}

[dustymugs]

Fascinating. Depending on the format passed to GDAL's gdaltransform, the resulting answer varies

gdaltransform -s_srs 'PROJCS["Amersfoort / RD New",GEOGCS["Amersfoort",DATUM["Amersfoort",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],AUTHORITY["EPSG","6289"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4289"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Oblique_Stereographic"],PARAMETER["latitude_of_origin",52.15616055555555],PARAMETER["central_meridian",5.38763888888889],PARAMETER["scale_factor",0.9999079],PARAMETER["false_easting",155000],PARAMETER["false_northing",463000],AUTHORITY["EPSG","28992"],AXIS["X",EAST],AXIS["Y",NORTH]]' -t_srs 'PROJCS["Popular Visualisation CRS / Mercator (deprecated)",GEOGCS["Popular Visualisation CRS",DATUM["Popular_Visualisation_Datum",SPHEROID["Popular Visualisation Sphere",6378137,0,AUTHORITY["EPSG","7059"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6055"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4055"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3785"],AXIS["X",EAST],AXIS["Y",NORTH]]'
195000 450000
664595.85121482 6773064.24012838 -13206.6219004256

gdaltransform -s_srs "+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" -t_srs "+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"
195000 450000
664595.859810406 6806797.32452805 43.5179952728836

gdaltransform -s_srs EPSG:28992 -t_srs 'PROJCS["Popular Visualisation CRS / Mercator (deprecated)",GEOGCS["Popular Visualisation CRS",DATUM["Popular_Visualisation_Datum",SPHEROID["Popular Visualisation Sphere",6378137,0,AUTHORITY["EPSG","7059"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6055"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4055"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3785"],AXIS["X",EAST],AXIS["Y",NORTH]]'
195000 450000
664595.85126645 6773064.23984736 -13206.6222764212

[dustymugs]

Sample raster file of same spatial attributes as the raster in PostGIS

[dustymugs]

After testing a raster file (attached to ticket) of the same spatial attributes as that in PostGIS with gdalwarp, the resulting raster matches that returning from ST_Resample. So, the algorithm used by PostGIS Raster (rt_raster_gdal_warp) behaves the same way as gdalwarp.

Because of the warping output matching and the differing outputs from gdaltransform for the same srs expressed in different formats, I can't decide if this is an issue or not. If it is an issue, it isn't one resolvable with the current implementation. The alternative would be to write our own warp API, but that is a task that shouldn't be taken lightly.

[rouault]

dustymugs,

all the results you observe with gdaltransform can be perfectly explained. If you look closely, none of the formulations of the SRS translates exactly to same proj.4 string.

You can play with the testepsg / gdalsrsinfo (in GDAL 1.9.0) utilities to examine the subtle differences.

Concerning Amersfoort, the expension of EPSG:28992 using GDAL data file of GDAL 1.8 or 1.9 doesn't have the same TOWGS84 values as the explicit proj4 definition you are using in second transform.

But I believe the main error comes from the google mercator definition. As it is an unusual SRS, using its WKT form as of the first line cannot work, because its expansion to proj4 will lack the "+nadgrids=@null" which is crucial. So use the proj4 form, or EPSG:900913 (or better EPSG:3857), or the full WKT version with the PROJ4 EXTENSION node :

PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+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 +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]

I'm loosely following postgis activity, so asking on gdal-dev ML can be a better place in case you have questions on GDAL working.

Best regards,

Even

[dustymugs]

Thanks Even for the informative details. I'll ask in the GDAL ML regarding methods to improve the use of the GDAL Warp API.

[dustymugs]

warped version of test.png into EPSG:3857 from EPSG:28992

[dustymugs]

Attached test.tif is output from gdalwarp using command

gdalwarp -t_srs "EPSG:3857" test.png test.tif

The warped tiff has the same drift seen in PostGIS Raster. Time to ask the GDAL ML about whether or not this is an issue…

[dustymugs]

After testing, the transformed rasters are correct spatially upon checking in qgis. Expected features (boundary edges and terrain features) are correct.