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

[Even Rouault]

Le mercredi 22 janvier 2014 20:27:19, Jason Greenlaw - NOAA Affiliate a écrit :
> 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.

The issue is that esri-3857.prf contains single quote characters instead of 
double quotes, and oceans_glakes_50m_3857.prj contain single quotes inside 
double quotes. Both are not valid AFAIK

You should be able to overcome this with :

ogr2ogr -a_srs EPSG:3857 oceans_glakes_50m_3857.shp 
oceans_glakes_50m_3857_fixed.shp

and

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

> 
> 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"]
> ]

-- 
Geospatial professional services
http://even.rouault.free.fr/services.html