gdal - gdalwarp `-te -5.8 41 10 51.5 -ts 1980 0` not working with mercator?

I was working with native lat/log projection on etopo1:

#donwload:
curl -L -C - 'http://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/grid_registered/georeferenced_tiff/ETOPO1_Ice_g_geotiff.zip' -o ./ETOPO1.zip
unzip -n ./ETOPO1.zip '*.tif'

#commands:
gdal_translate -projwin -5.8 51.5 10 41 ../data/noaa/ETOPO1_Ice_g_geotiff.tif cropXL.tmp.tif
gdalwarp -of GTiff -s_srs EPSG:4326 -t_srs EPSG:4326 cropXL.tmp.tif reproj.tmp.tif
gdalwarp  -te -5.8 41 10 51.5 -ts 1980 0   reproj.tmp.tif resized.tmp.tif

all working fine.

I now wish to work in mercator. Also, I made a tiny change to reproject to mercator -t_srs EPSG:3857, but it doesn't work:

gdal_translate -projwin -5.8 51.5 10 41 ../data/noaa/ETOPO1_Ice_g_geotiff.tif cropXL.tmp.tif
gdalwarp -of GTiff -s_srs EPSG:4326 -t_srs EPSG:3857 cropXL.tmp.tif reproj.tmp.tif
gdalwarp -te -5.8 41 10 51.5  -ts 1980 0   reproj.tmp.tif resized.tmp.tif

The reprojection (2) seems fine. The last command (3) runs but resized.tmp.tif is now flat, all black. No error message. But the values x,y,z seems non-sense.

What am I doing wrong ?

Answer

Solution

The part that is going wrong is where you assumed that the coordinates in 3857 are degrees (like -5.8, 41). The units of 3857 are metres. So you're asking for something tiny, off the map (near the origin).

Lets look at a conversion, using pyproj:

from pyproj import Proj, transform

inProj = Proj(init='epsg:4326')

outProj = Proj(init='epsg:3857')
x1,y1 = -5.8, 41
x2,y2 = transform(inProj,outProj,x1,y1)
print x2,y2

That is going to convert the input to EPSG:3857, with an output of:

-645653.046601 5012341.66385

When you project your data, you need to project the bounding box too.

gdalsrsinfo

You can check the projection properties on epsg.io (epsg.io/4326 and epsg.io/3857) or via gdalsrsinfo :

$ gdalsrsinfo EPSG:4326

PROJ.4 : '+proj=longlat +datum=WGS84 +no_defs '
OGC WKT :
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],

    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,                    # <========= units !
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]

and

$gdalsrsinfo EPSG:3395

PROJ.4 : '+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs '

OGC WKT :
PROJCS["WGS 84 / World Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,

            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,                   # <========= units !
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],

    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","3395"]]