Merging tiles using gdal_merge.py

The geospatial data abstraction library is awesome!

With gdal_merge.py you can merge tiles of geospatial data together so that instead of multiple images you have one image. See the two images we want to merge below.

Find a link to two example data sets here. Extract the files into one folder.

In QGIS

Step 1: Go to Raster -> miscellaneous -> Merge

Step 2: Choose “input directory instead of files” and select the folder with the files you want to merge.

Step 3: Fill out the name of the output file.

Step 4: We add the nodata value (9999.990234375) by hand using the edit button.

Step 5: Press the Ok button.

To clean the data up we add a gdalwarp command.

Step 6: Go to raster -> projections -> warp (reproject)

Step 7: Fill out the inputfile, output file, source SRS, Target SRS and nodata values

Command line

You can of course also directly run the code constructed in QGIS in your terminal.

gdal_merge.py -n 9999.990234375 -of GTiff -o /path/to/myoutputfile.tif /path/to/inputfile1.tif /path/to/inputfile2.tif /path/to/inputfileX.tif

gdalwarp -s_srs EPSG:4326 -t_srs epsg:4326 -dstnodata 0 -of GTiff /path/to/myoutputfile.tif /path/to/newoutputfile.tif

When you have a large number of files it is better to construct a string of the command in python and execute this automatically:

import os

# three horizontal tiles
horizontaltiles = [27,28,29]

# horizontal tiles have different set of vertical tiles
verticaltiles =  [[6,7,8],[6,7,8,9],[6,7,8,9]]

# location of the files
pathname = "/path/to/"
# the basename
basename ="ETensemble250m"


# for all years
for year in range(2003,2015,1):
   # for all months
   for month in range(1,13,1):
      
      # the gdal_merge command
      command = "gdal_merge.py -n 9999.990234375 -of GTiff -o " 
      
      # the output name
      outputname = pathname +  "ET_ensemble_" + str(year) + "-" + str(month) + "-01_temp.tif " 
      
      # set the string here. We want the command to run for each month in all years
      command = command + outputname
      
      # for all horizontal tiles horizontal tiles
      for h in horizontaltiles:
         
         # get index of hortizontal tile
         myindex = horizontaltiles.index(h)
         
         # for all vertica tiles
         for v in verticaltiles[myindex]:
            tile = "h" + str(h) + "v" + str(v)
            thedate = str(year) + "-" + str(month) + "-01"
            command = command + pathname +   basename + "-" + tile + "-"+ thedate + ".tif "
      print " "  
      print command
      
      # run gdal_merge.py tool with os call
      os.system(command)
      
      # warp to clean up data
      newoutputname = pathname +  "ET_ensemble_" + str(year) + "-" + str(month) + "-01.tif "
      warpcommand  = "gdalwarp -s_srs EPSG:4326 -t_srs epsg:4326 -dstnodata 0 -of GTiff " + outputname  + " " + newoutputname
      
      print " "
      print warpcommand
      
      # execute gdalwarp command with os call
      os.system(warpcommand)