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

[oliver.klages at ksti.de]

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?

i could just modify the coordinates to match the world (longitudinal) using gdal_translate but i fear i am overlooking something else here which might just come back to bite me.


Bet regards,


Oliver

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20160817/4539c65d/attachment.html>