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 > > <[email protected]>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 _______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
