0

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
  • In what way the result is incorrect? – user30184 Sep 1 '14 at 16:16
  • I expanded the post to the test is easily reproducible. I'am digging in to explain why it's incorrect. Basically, altitude is lost, i have a reprojected image which is all black in Qgis. Without data, I can't process it further. – Hugolpz Sep 1 '14 at 16:20
  • Looks good on my QGIS. Check what options you use for rendering from Layer properties - Style. Play with contrast enhancement options - Stretch to MinMax is a good starting point. You can also check the histogram with gdalinfo -hist – user30184 Sep 1 '14 at 16:43
  • I added my gdal info -hist for my file an other area of interest and other output name, but same process. – Hugolpz Sep 1 '14 at 17:27
  • What is your GDAL version? I also suggest to run gdalinfo for each interim tiff you create. Statistics with zeros as min and max prove that your image is really all black. I sent my version as an answer because I want to encourage you with a screen capture. – user30184 Sep 1 '14 at 17:57
2

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
  • Ok, I refactored my code to write down this question more clearly, I don't exactly know were was my error / trouble, I suspect the sequence of command (reproject then resize = bug, resize then reproject = ok), but it now works as you told me. Good ! – Hugolpz Sep 2 '14 at 9:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.