[postgis-devel] [PostGIS] #1402: [raster] ST_Resample producing incorrect transformation

[PostGIS]

#1402: [raster] ST_Resample producing incorrect transformation
-----------------------------+----------------------------------------------
  Reporter:  dustymugs       |       Owner:  dustymugs    
      Type:  defect          |      Status:  closed       
  Priority:  blocker         |   Milestone:  PostGIS 2.0.0
 Component:  postgis raster  |     Version:  trunk        
Resolution:  wontfix         |    Keywords:               
-----------------------------+----------------------------------------------

Comment(by 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

-- 
Ticket URL: <http://trac.osgeo.org/postgis/ticket/1402#comment:3>
PostGIS <http://trac.osgeo.org/postgis/>
The PostGIS Trac is used for bug, enhancement & task tracking, a user and developer wiki, and a view into the subversion code repository of PostGIS project.