[gdal-dev] GDAL v1.10.1 - Bad output from gdalwarp (works with 1.10.0)

[Jason Greenlaw - NOAA Affiliate]

When attempting to use a polygon shapefile (in custom web mercator
projection) with the gdalwarp cutline operation, operating on a GeoTIFF
(with same projection), the following error is output in both versions
1.10.0 and 1.10.1:

    ERROR 6: No translation for 'Mercator_Auxiliary_Sphere' to PROJ.4
format is known.

With GDAL 1.10.0, the operation still completes successfully with a zero
exit status, but with GDAL 1.10.1 the operation fails with a nonzero exit
status, and the resulting GeoTIFF has zeros for all pixels.

Here is the command I'm using:

    gdalwarp -co COMPRESS=LZW -cutline oceans_glakes_50m_3857.shp
-srcnodata 9999.0 -dstnodata 9999.0 -overwrite sst_global.tif
sst_global_cut.tif


The input shapefile was produced by reprojecting with ogr2ogr from
EPSG:4326 to my custom projection (WKT is included below), which is stored
in a *.prf file.  The source GeoTIFF is single-band with Float32 pixel type.

The custom WKT can be found here:  http://jgreenlaw.org/gdal/esri-3857.prf

The cutline shapefile can be downloaded at:
http://jgreenlaw.org/gdal/oceans_glakes_50m_3857.tar.gz (900kb)

The source GeoTIFF can be generated from this FLT/HDR raster:
http://jgreenlaw.org/gdal/sst_flt.tar.gz (11mb) by running these commands:

    gdalwarp -s_srs EPSG:4326 -t_srs esri-3857.prf -r bilinear -srcnodata
9999.0 -dstnodata 9999.0 -tr 4000 9000 -overwrite -co COMPRESS=LZW -co
TILED=yes sst.flt sst_projected.tif

    gdalwarp -co COMPRESS=LZW -te -20037507.0671618 -19971868.8804086
20037507.0671618 19971868.8804086 -srcnodata 9999.0 -dstnodata 9999.0 -tr
4000 9000 sst_projected.tif sst_global.tif

Does anyone know if this behavior is to be expected when the above error is
output?  Is there a better method for reprojecting to Web Mercator that
will align properly in ArcGIS?

Thanks
-Jason


Background:
The reason for using custom WKT for reprojecting to Web Mercator is because
I have experienced alignment problems displaying GDAL-produced GeoTIFFs in
ArcGIS when reprojecting using "EPSG:3857" as the target SRS (this seems to
be related to http://trac.osgeo.org/gdal/ticket/3962).

I generated the WKT below by copying the 102113 definition in
$GDAL_DATA/esri_extra.wkt and merging information from an ESRI Web Mercator
*.prj file.  Using this custom WKT as the target SRS with gdalwarp normally
results in rasters that align properly with ArcGIS Web Mercator data
(without on-the-fly datum transformation).


----------------------------------------------------------------------------------------------
Custom WKT for reprojecting to ESRI Web Mercator using GDAL (
http://jgreenlaw.org/gdal/esri-3857.prf):
----------------------------------------------------------------------------------------------

PROJCS['WGS_1984_Web_Mercator_Auxiliary_Sphere',
    GEOGCS['GCS_WGS_1984',
        DATUM['D_WGS_1984',
            SPHEROID['WGS_1984',6378137.0,298.257223563]
        ],
        PRIMEM['Greenwich',0.0],
        UNIT['Degree',0.0174532925199433]
    ],
    PROJECTION['Mercator_Auxiliary_Sphere'],
    PARAMETER['False_Easting',0.0],
    PARAMETER['False_Northing',0.0],
    PARAMETER['Central_Meridian',0.0],
    PARAMETER['latitude_of_origin',0.0],
    PARAMETER['Standard_Parallel_1',0.0],
    PARAMETER['Auxiliary_Sphere_Type',0.0],
    UNIT['Meter',1.0],
    EXTENSION["PROJ4",
        "+proj=merc +a=6378137 +b=6356752.314245179 +lat_ts=0.0 +lon_0=0.0
+x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
    ],
    AUTHORITY["EPSG","3857"]
]
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20140122/801f9788/attachment.html>