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

Even Rouault even.rouault at mines-paris.org
Wed Jan 22 13:18:07 PST 2014 

Le mercredi 22 janvier 2014 22:09:33, Jason Greenlaw - NOAA Affiliate a écrit :
> Even,
> 
> Thank you for catching the quote problem.  I've fixed my WKT accordingly.
> 
> Unfortunately altering the SRS to GDAL's default EPSG:3857 definition
> results in the same alignment problem when viewed in ArcGIS Desktop/Server,
> so I altered it back to my custom WKT (using double quotes this time) - now
> it's aligning properly again, but I'm still seeing the same error:
> 
>     ERROR 6: No translation for Mercator_Auxiliary_Sphere to PROJ.4 format
> is known.

Yes, this is an issue. Nobody agrees on how to represent WebMercator in WKT. 
There's a ticket pending in Trac about that.

While looking at the .prj, I see a difference in the EXTENSION["PROJ4" node 
w.r.t what GDAL does.

The .prj has :
"+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"

Whereas gdalsrsinfo EPSG:3857 shows :
"+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"

So the ESRI stuff uses a non-spheric ellipsoid, which I believe is wrong...

You can probably get the same result as ESRI by using the "+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" as -a_srs value.


> 
> It still produces valid output with 1.10.0, but with 1.10.1 it still exits
> prematurely with status 1 and produces an all-zero image.
> 
> It seems that something was changed between 1.10.0 and 1.10.1 that made
> this a fatal error.

I don't see anything obvious in the log that explain the behavioural change, 
but erroring out make sense when there's an error... So 1.10.1 behaviour seems 
OK.

> 
> Thanks,
> Jason
> 
> 
> On Wed, Jan 22, 2014 at 3:07 PM, Even Rouault
> 
> <even.rouault at mines-paris.org>wrote:
> > 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

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

More information about the gdal-dev mailing list