[GIS] Gdal: How to reproject from equirectangular to mercator

coordinate systemgdalmercator

I work with raster ETOPO1, made a crop of my area of interest.

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'
# --- crop
gdal_translate -projwin -5 51 10 41 ETOPO1_Ice_g_geotiff.tif cropXL.tmp.tif
# --- resize for speed sake
gdalwarp -of GTiff -s_srs epsg:4326 -t_srs epsg:4326 -te -5 41 10 51 -ts 1200 0 cropXL.tmp.tif resized.tmp.tif
# --- reproject into mercator
gdalwarp -of GTiff -s_srs EPSG:4326 -t_srs EPSG:3857 resized.tmp.tif resized-mercator.tif

But this result into an incorrect file. What does I do wrong ? How to reproject correctly from EPSG:4326 to EPSG:3857 ?

Note: command inspired by Re-project raster image from Mercator to Equirectangular

Edit:

yug@yug-PC:~/projects/WikiAtlas_scripts/back$ gdalinfo ./01_reliefs/crop_resized.tmp.tif -hist
Driver: GTiff/GeoTIFF
Files: ./01_reliefs/crop_resized.tmp.tif
       ./01_reliefs/crop_resized.tmp.tif.aux.xml
Size is 1962, 1551
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    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"]],
    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"]],
    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","3857"]]
Origin = (11365720.009993232786655,1654650.287315490888432)
Pixel Size = (323.425157454019200,-323.425157454019200)
Metadata:
  AREA_OR_POINT=Area
  NC_GLOBAL#Conventions=COARDS/CF-1.0
  NC_GLOBAL#GMT_version=4.4.0
  NC_GLOBAL#history=grdreformat ETOPO1_Ice_g_gdal.grd ETOPO1_Ice_g_gmt4.grd=ni
  NC_GLOBAL#node_offset=0
  NC_GLOBAL#title=ETOPO1_Ice_g_gmt4.grd
  x#actual_range=-180, 180
  x#long_name=Longitude
  x#units=degrees
  y#actual_range=-90, 90
  y#long_name=Latitude
  y#units=degrees
  z#_FillValue=-2147483648
  z#actual_range=-10898, 8271
  z#long_name=z
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (11365720.010, 1654650.287) (102d 6' 0.00"E, 14d42' 0.00"N)
Lower Left  (11365720.010, 1153017.868) (102d 6' 0.00"E, 10d18' 6.40"N)
Upper Right (12000280.169, 1654650.287) (107d48' 1.26"E, 14d42' 0.00"N)
Lower Right (12000280.169, 1153017.868) (107d48' 1.26"E, 10d18' 6.40"N)
Center      (11683000.089, 1403834.078) (104d57' 0.63"E, 12d30'36.90"N)
Band 1 Block=1962x2 Type=Int16, ColorInterp=Gray
  Min=0.000 Max=0.000 
  Minimum=0.000, Maximum=0.000, Mean=0.000, StdDev=0.000
  1 buckets from 0 to 0:
  0 
  NoData Value=-2147483648
  Metadata:
    NETCDF_VARNAME=z
    STATISTICS_MAXIMUM=0
    STATISTICS_MEAN=0
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=0

Best Answer

This is odd, but by following exactly your workflow I get this as a result:

enter image description here

Gdalinfo with statistics for my version of resized-mercator.tif is here:

Corner Coordinates:
Upper Left  ( -556597.454, 6621293.723) (  5d 0' 0.00"W, 51d 0' 0.00"N)
Lower Left  ( -556597.454, 5011872.466) (  5d 0' 0.00"W, 40d59'48.55"N)
Upper Right ( 1113920.713, 6621293.723) ( 10d 0'23.47"E, 51d 0' 0.00"N)
Lower Right ( 1113920.713, 5011872.466) ( 10d 0'23.47"E, 40d59'48.55"N)
Center      (  278661.630, 5816583.095) (  2d30'11.74"E, 46d13'32.32"N)
Band 1 Block=1039x3 Type=Int16, ColorInterp=Gray

     Metadata:
       NETCDF_VARNAME=z
       STATISTICS_MAXIMUM=3569
       STATISTICS_MEAN=34.570234537257
       STATISTICS_MINIMUM=-4953
       STATISTICS_STDDEV=1037.5075489272
Related Question