[raster] ST_Transform to SRID 3978 produce wrong results

[pracine]

I have a raster in srid=4269:

-- Raster unprojected
SELECT ST_AsBinary(ST_AddBand(ST_MakeEmptyRaster(10, 10, -168.0, 85.0, 0.083, -0.083, 0, 0, 4269), 1, '8BUI', 1, 0)::geometry);

I project the geometric extent of this raster to srid=3978:

-- Raster extent projected
SELECT ST_AsBinary(ST_Transform(ST_AddBand(ST_MakeEmptyRaster(10, 10, -168.0, 85.0, 0.083, -0.083, 0, 0, 4269), 1, '8BUI', 1, 0)::geometry, 3978));

When I project the raster (not the vectorization of its extent) to the same srid=3978, it is wrongly projected:

--Raster projected
SELECT ST_AsBinary(ST_Transform(ST_AddBand(ST_MakeEmptyRaster(10, 10, -168.0, 85.0, 0.083, -0.083, 0, 0, 4269), 1, '8BUI', 1, 0), 3978)::geometry);

[dustymugs]

So, I can't seem to get QGIS to behave right with SRID 3978 (It doesn't have it built in and gives me crazy results when I do). I'm testing with the SRID 3347 instead and everything looks about right.

The test I use is...

ROP TABLE IF EXISTS test_1618;
WITH foo AS (
	SELECT ST_AddBand(ST_MakeEmptyRaster(10, 10, -168.0, 85.0, 0.083, -0.083, 0, 0, 4269), 1, '8BUI', 1, 0) AS rast
)
SELECT
	gid,
	geom
INTO test_1618
FROM ((
	SELECT
		0 AS gid,
		rast::geometry AS geom
	FROM foo
	) UNION ALL (
	SELECT
		1 AS gid,
		ST_Transform(rast::geometry, 3347) AS geom
	FROM foo
	) UNION ALL (
	SELECT
		2 AS gid,
		ST_Transform(rast, 3347)::geometry AS geom
	FROM foo
	)
) bar

I've checked the debug output for ST_Transform(raster) and it matches with what you'd get from GDAL's gdalwarp utility. This is especially true considering that the ST_Transform(raster, 3347) is absolutely basic and uses everything that GDAL suggests.

[pracine]

How does GDAL behave with srid=3978? Can you correctly gdalwarp a raster to 3978?

[pracine]

Could we make a test trying to identify if there is any other srid to which we can not transform properly? Like looping over every srid and comparing the resulting extent with the equivalent geometry transformation? I know it would be hard to create realistic extents for every srid but maybe there is a way to do it without too much hassle.

I can't believe I choose the only rebel srid to build my big Canadian database...

[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@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@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]],

PROJECTIONLambert_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...

[dustymugs]

Shoot... The second part formatted...

root@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]]

[dustymugs]

Try r9386 as I've tweaked what spatial reference information is used by GDAL.

[pracine]

Did you try the first example I provided? Still gives the same thing for me...

[dustymugs]

Yes. The answer coming out is exactly the same answer as that from gdalwarp. I've attached a TIFF matching the raster in your example. If you warp that image with gdalwarp

gdalwarp -t_srs EPSG:3978 test4269.tif test3978.tif

The output matches that of ST_Transform.

[dustymugs]

Sample source raster

[dustymugs]

Output warped raster

[pracine]

So we agree that this is a GDAL bug?

[dustymugs]

I don't see a bug. The transformed raster's extent is correct as it does cover the transformed area. If you feel that it is a bug, you are more than welcome to discuss it with Even or Frank.