Saturday, 15 August 2015

geotiff tiff - Merging GDAL tiles and crop via bounding box?


I am using GDAL.


Assume I have four adjacent GDAL tiles in the formation of a square. I have a bounding box with a corner in each tile so that the rectangle formed spans all four tiles. My goal is to generate a DEM from those four adjacent GDAL files that is cropped to the bounding box. Also as an added bonus, assume these GDAL tiles are in folder of GDAL tiles. In this case, the DEM is a Geotile.


In essence I am asking how to perform three steps here:



  1. Find all GDAL tiles covered by the bounding box.

  2. Merge the tiles into one giant DEM

  3. Crop the DEM to the bounding box and output as GEOTIFF.



The methods I know to do this would be really efficient (ex: open each GDAL sequentially and find all the ones that contain a corner point, merge tiles, then crop, surely GDAL has a more efficient means of doing this).



Answer



This can be achieved with the help of GDAL's Virtual Raster Format. With this you can essentially skip the step of creating one giant DEM. The VRT will be handled by GDAL like a giant, merged DEM but is just a small XML file containing the file paths for each tile as well as some metadata. This can then be fed to gdalwarp together with a bounding box or a cutline.


This assumes that your tiles do not overlap. My guess is you are probably working with SRTM tiles for which this is the case.



  1. Merge the tiles into one giant, virtual DEM


gdalbuildvrt big_DEM.vrt tile_folder/*.tif




  1. Crop the DEM to your liking and write the result to a new file. In this case with a cropline vector dataset.


gdalwarp -of GTiff -cutline region_of_interest.kml -crop_to_cutline big_DEM.vrt small_DEM.tif


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...