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

Jason Greenlaw - NOAA Affiliate jason.greenlaw at noaa.gov
Wed Jan 22 13:40:38 PST 2014 

Yes, ESRI does seem to have their own way of doing things, so it's not too
surprising.

What's strange is that an ESRI-generated .prj file for EPSG:3857 looks
something like this (no mention of the semi-minor axis):

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["Standard_Parallel_1",0.0],
PARAMETER["Auxiliary_Sphere_Type",0.0],
UNIT["Meter",1.0],
AUTHORITY["EPSG",3857]
]



But in the ArcGIS GUI, the details of the projection are reported as the
following, with the modified semi-minor axis (this is what I used to modify
the proj4 string in my WKT):

WGS_1984_Web_Mercator_Auxiliary_Sphere
WKID: 3857 Authority: EPSG

Projection: Mercator_Auxiliary_Sphere
False_Easting: 0.0
False_Northing: 0.0
Central_Meridian: 0.0
Standard_Parallel_1: 0.0
Auxiliary_Sphere_Type: 0.0
Linear Unit: Meter (1.0)

Geographic Coordinate System: GCS_WGS_1984
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_WGS_1984
  Spheroid: WGS_1984
    Semimajor Axis: 6378137.0
    Semiminor Axis: 6356752.314245179
    Inverse Flattening: 298.257223563


Jason


On Wed, Jan 22, 2014 at 4:18 PM, Even Rouault
<even.rouault at mines-paris.org>wrote:

> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20140122/86117b87/attachment.html>

More information about the gdal-dev mailing list