[gdal-dev] Problem with Gdal transforming "world" wgs84 to EPSG:3857

[Even Rouault]

Le mercredi 17 août 2016 09:07:28, oliver.klages at ksti.de a écrit :
> down
> votefavorite<http://stackoverflow.com/questions/38953211/using-gdalwarp-to
> -transform-from-wgs84-to-epsg3857-does-not-cover-world#>
> 
> 
> i am using gdal 1.10 and 2.1.1.
> 
> i have a VRT datasource in WGS84 where i forced the corner coordinates to
> the min/max values of EPSG:3857 (-180,85.5,180,-85.5).
> 
> gdalinfo output for this VRT looks like this:
> 
> Size is 1296001, 601200
> Coordinate System is:
> GEOGCS["WGS 84",
>     DATUM["WGS_1984",
>         SPHEROID["WGS 84",6378137,298.257223563,
>             AUTHORITY["EPSG","7030"]],
>         AUTHORITY["EPSG","6326"]],
>     PRIMEM["Greenwich",0],
>     UNIT["degree",0.0174532925199433],
>     AUTHORITY["EPSG","4326"]]
> Origin = (-180.000000000000000,85.500000000000000)
> Pixel Size = (0.000277777563443,-0.000284431137725)
> Corner Coordinates:
> Upper Left  (-180.0000000,  85.5000000) (180d 0' 0.00"W, 85d30' 0.00"N)
> Lower Left  (-180.0000000, -85.5000000) (180d 0' 0.00"W, 85d30' 0.00"S)
> Upper Right ( 180.0000000,  85.5000000) (180d 0' 0.00"E, 85d30' 0.00"N)
> Lower Right ( 180.0000000, -85.5000000) (180d 0' 0.00"E, 85d30' 0.00"S)
> Center      (   0.0000000,  -0.0000000) (  0d 0' 0.01"E,  0d 0' 0.00"S)
> Band 1 Block=128x128 Type=Int16, ColorInterp=Gray
> 
> 
> Basically, i have the world minus the poles.
> 
> Now i want to convert this to EPSG:3857.
> 
> I use gdalwarp for this, using bilinear interpolation:
> 
> ./gdalwarp -of VRT -co TILED=YES -srcnodata 9999 -t_srs 'EPSG:3785'  -multi
>  wgs84.vrt  wmerc.vrt  -overwrite -r bilinear
> 
> 
> running gdalinfo on wmerc then gives this:
> 
> Size is 995026, 1025175
> Coordinate System is:
> PROJCS["Popular Visualisation CRS / Mercator (deprecated)",
>     GEOGCS["Popular Visualisation CRS",
>         DATUM["Popular_Visualisation_Datum",
>             SPHEROID["Popular Visualisation Sphere",6378137,0,
>                 AUTHORITY["EPSG","7059"]],
>             TOWGS84[0,0,0,0,0,0,0],
>             AUTHORITY["EPSG","6055"]],
>         PRIMEM["Greenwich",0,
>             AUTHORITY["EPSG","8901"]],
>         UNIT["degree",0.0174532925199433,
>             AUTHORITY["EPSG","9122"]],
>         AUTHORITY["EPSG","4055"]],
>     PROJECTION["Mercator_1SP"],
>     PARAMETER["central_meridian",0],
>     PARAMETER["scale_factor",1],
>     PARAMETER["false_easting",0],
>     PARAMETER["false_northing",0],
>     UNIT["metre",1,
>         AUTHORITY["EPSG","9001"]],
>     AXIS["X",EAST],
>     AXIS["Y",NORTH],
>     EXTENSION["PROJ4","+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"], AUTHORITY["EPSG","3785"]]
> Origin = (-20037508.342789247632027,20644642.363762538880110)
> Pixel Size = (40.275313937642778,-40.275354414328227)
> Corner Coordinates:
> Upper Left  (-20037508.343,20644642.364) (180d 0' 0.00"E, 85d28'11.27"N)
> Lower Left  (-20037508.343,-20644644.098) (180d 0' 0.00"E, 85d28'11.28"S)
> Upper Right (20037476.183,20644642.364) (179d59'58.96"E, 85d28'11.27"N)
> Lower Right (20037476.183,-20644644.098) (179d59'58.96"E, 85d28'11.28"S)
> Center      ( -16.0797308,  -0.8670919) (  0d 0' 0.52"W,  0d 0' 0.03"S)
> Band 1 Block=512x128 Type=Int16, ColorInterp=Gray
>   NoData Value=9999
> 
> 
> Note that the corner coordinates for upper/lower left look correct, but the
> corner coordinates for upper/lower right (the longitude) is missing is
> missing 32 units. Broadly put, i am missing a sliver on the right side,
> but only in regards to the coordinates. The data is there, but the
> coordinates to the right seem wrong.
> 
> Why is that?

When gdalwarp must guess the output extent and resolution, this involves a 
number of things. One of those is that gdalwarp will compute the output 
resolution to be square pixels so it is not generally possible to have both 
exact extent and square pixels at the same time.

But you can force the target extent manually with the -te parameters

As we have:

$ echo "-180 -85.5" | gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857
-20037508.3427892 -20644642.3637625 0

$ echo "180 85.5" | gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857
20037508.3427892 20644642.3637625 0

You can use :
-te -20037508.3427892 -20644642.3637625 20037508.3427892 20644642.3637625

The pixel size will be slightly altered from perfect square pixels to almost 
square pixels to enforce the specified extent.

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com