[postgis-devel] [PostGIS] #1618: [raster] ST_Transform to SRID 3978 produce wrong results

[PostGIS]

#1618: [raster] ST_Transform to SRID 3978 produce wrong results
---------------------+------------------------------------------------------
 Reporter:  pracine  |       Owner:  dustymugs    
     Type:  defect   |      Status:  assigned     
 Priority:  blocker  |   Milestone:  PostGIS 2.0.0
Component:  raster   |     Version:  trunk        
 Keywords:           |  
---------------------+------------------------------------------------------

Comment(by dustymugs):

 There is definitely something going on with GDAL's handling of the proj4
 text passed to it.  It looks like it is using a spatial reference totally
 different than 3978.

 By default, we pass GDAL proj4 text (based upon advice from GDAL folks).
 If you have GDAL 1.9, you can use the gdalsrsinfo utility to see what I'm
 talking about.

 When gdalsrsinfo is called for EPSG:3978

 {{{
 root at slack64:/home/software/geo# gdalsrsinfo EPSG:3978

 PROJ.4 : '+proj=lcc +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +x_0=0 +y_0=0
 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs '

 OGC WKT :
 PROJCS["NAD83 / Canada Atlas Lambert",
     GEOGCS["NAD83",
         DATUM["North_American_Datum_1983",
             SPHEROID["GRS 1980",6378137,298.257222101,
                 AUTHORITY["EPSG","7019"]],
             TOWGS84[0,0,0,0,0,0,0],
             AUTHORITY["EPSG","6269"]],
         PRIMEM["Greenwich",0,
             AUTHORITY["EPSG","8901"]],
         UNIT["degree",0.0174532925199433,
             AUTHORITY["EPSG","9122"]],
         AUTHORITY["EPSG","4269"]],
     PROJECTION["Lambert_Conformal_Conic_2SP"],
     PARAMETER["standard_parallel_1",49],
     PARAMETER["standard_parallel_2",77],
     PARAMETER["latitude_of_origin",49],
     PARAMETER["central_meridian",-95],
     PARAMETER["false_easting",0],
     PARAMETER["false_northing",0],
     UNIT["metre",1,
         AUTHORITY["EPSG","9001"]],
     AXIS["Easting",EAST],
     AXIS["Northing",NORTH],
     AUTHORITY["EPSG","3978"]]
 }}}

 Now, if we took that PROJ.4 text outputted above and passed that into
 gdalsrsinfo

 {{{root at slack64:/home/software/geo# gdalsrsinfo '+proj=lcc +lat_1=49
 +lat_2=77 +lat_0=49 +lon_0=-95 +x_0=0 +y_0=0 +ellps=GRS80
 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs '

 PROJ.4 : '+proj=lcc +lat_1=49 +lat_0=49 +lon_0=-95 +k_0=1 +x_0=0 +y_0=0
 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs '

 OGC WKT :
 PROJCS["unnamed",
     GEOGCS["GRS 1980(IUGG, 1980)",
         DATUM["unknown",
             SPHEROID["GRS80",6378137,298.257222101],
             TOWGS84[0,0,0,0,0,0,0]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433]],
     PROJECTION["Lambert_Conformal_Conic_1SP"],
     PARAMETER["latitude_of_origin",49],
     PARAMETER["central_meridian",-95],
     PARAMETER["scale_factor",1],
     PARAMETER["false_easting",0],
     PARAMETER["false_northing",0],
     UNIT["Meter",1]]
 }}}

 As you can see, the OGC WKT is different and is causing the projection
 issues.  I'll need to think about how to fix this...

-- 
Ticket URL: <http://trac.osgeo.org/postgis/ticket/1618#comment:5>
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.