1

I'm trying to use gdal_warp to slice and convert some SRTM data. The data I'm interested in are in S30E152.hgt and S30E153.hgt. What I'd like to do is:

  1. Convert these DEM files from EPSG:4326 to EPSG:3857
  2. Select a subsection of the data at coordinates 17008000,-3590694 17051000,-3503623
  3. Write this subsection to a GeoTIFF file for further processing

Firstly, I add the two DEM files to a VRT:

gdalbuildvrt combined.vrt S30E152.hgt S30E153.hgt

Then I check the VRT to ensure the height values are sane:

gdalinfo -mm combined.vrt [snip] Band 1 Block=128x128 Type=Int16, ColorInterp=Undefined Computed Min/Max=-14.000,1501.000 [snip]

Then I use gdal_warp to convert and slice the VRT:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 -te 17008000 -3590694 17051000 -3503623 -r bilinear combined.vrt sliced.tiff

Then I check the GeoTIFF to ensure the height values are still sane:

gdalinfo -mm sliced.tiff [snip] Band 1 Block=449x9 Type=Int16, ColorInterp=Gray Computed Min/Max=0.000,0.000 [snip]

It appears that whenever I use gdalwarp with a VRT containing multiple DEMs, it always clobbers the height values. For comparison, I created a VRT with only a single DEM and tried the same operations (using a smaller subsection that fits within range).

gdalbuildvrt test.vrt S30E152.hgt gdalinfo -mm test.vrt (produces Min/Max=-14.000,1501.000) gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 -te 16931600 -3388300 16976200 -3439400 -r bilinear test.vrt test.tiff gdalinfo -mm test.tiff [snip] Band 1 Block=446x9 Type=Int16, ColorInterp=Gray Computed Min/Max=53.000,1471.000 [snip]

This seems to work OK. Maybe the problem was with my coordinates. So I thought I'd warp the entire VRT, and then select the subsection:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 combined.vrt converted.tiff gdalinfo -mm converted.tiff [snip] Band 1 Block=2328x1 Type=Int16, ColorInterp=Gray Computed Min/Max=-14.000,1501.000 [snip] gdaltranslate -projwin 17008000 -3503623 17051000 -3590694 converted.tiff sliced.tiff

Although I've double-checked the coordinates, I always get this warning:

Warning: Computed -srcwin 914 1337 449 910 falls partially outside raster extent. Going on however.

Now this gives me a file with reasonable height values:

gdalinfo -mm sliced.tiff [snip] Band 1 Block=449x9 Type=Int16, ColorInterp=Gray Computed Min/Max=0.000,369.000 [snip]

However when I check to see if the data is sane, I get an empty image:

gdal_translate -scale 0 369 0 255 -outsize 100 200 -of PNG sliced.tiff sliced.png

At this point, I'm totally lost. I think I have data, but am not sure.

  • Forgot to add, gdalwarp --version reports GDAL 1.10.1, released 2013/08/26. – simonhaines Sep 4 '14 at 5:20
  • The question is OK. It contains all information which is needed for reproducing the case (except links to source DEMs). Nothing less in not enough for solving this case. – user30184 Sep 4 '14 at 7:46
  • I can't reproduce the error with my own DEM images and GDAL 2.0-dev. Without access to S30E152.hgt and S30E153.hgt I can't do more. However, there may be something odd with gdalwarp and -te gis.stackexchange.com/questions/112395/…. It would be nice to catch it. – user30184 Sep 4 '14 at 8:47
  • @user30184 You should be able to get the data from the SRTM FTP site: dds.cr.usgs.gov/srtm/version2_1/SRTM3/Australia – WhiteboxDev Sep 4 '14 at 11:23
  • Convert the first .VRT into epsg:3857 and check the extents with gdalinfo. I get Lower Left (16920516.217,-3503628.939) Upper Right (17143287.264,-3375593.003) It means that your slice with S coordinates between -3590694 and -3503623 do not hit the data. – user30184 Sep 4 '14 at 11:48
0

Totally black output was a result from the extents which were used in the gdalwarp command. They were outside the image area and therefore the output was nothing else than nodata.

I recommend to use QGIS for debugging this kind of troubles. The VRT can be opened directly with QGIS as a raster layer and it can be reprojected on-the-fly to target projection which in this case was EPSG:3857. The target extents used as -te parameters can be placed on the QGIS map with the nice Quick WKT plugin. Convert -te 17008000 -3590694 17051000 -3503623 into two WKT points or build a corresponding polygon from the corners as below:

POINT (17008000 -3590694)
POINT (17051000 -3503623)

POLYGON ((
        17008000 -3590694, 
        17008000 -3503623, 
        17051000 -3503623, 
        17051000 -3590694, 
        17008000 -3590694
    ))

Having a look at the QGIS map windows reveals immediately the reason for the empty gdalwarp output.

DEM area and BBOX both in EPSG:3857

It is also possible to clip and extract a slice of the DEM with QGIS by using the Clipper tool that is in menu Raster - Extraction - Clipper. However, if the VRT file combines hundreds or thousands of original images it is more efficient to use direct GDAL commands because opening such VRT into QGIS can be slow.

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.