text: #!/usr/bin/env python
  ###############################################################################
  # $Id$
  #
  # Project:  InSAR Peppers
  # Purpose:  Module to extract data from many rasters into one output.
  # Author:   Frank Warmerdam, warmerdam@pobox.com
  #
  ###############################################################################
  # Copyright (c) 2000, Atlantis Scientific Inc. (www.atlsci.com)
  # Copyright (c) 2009-2011, Even Rouault <even dot rouault at mines-paris dot org>
  # 
  # This library is free software; you can redistribute it and/or
  # modify it under the terms of the GNU Library General Public
  # License as published by the Free Software Foundation; either
  # version 2 of the License, or (at your option) any later version.
  # 
  # This library is distributed in the hope that it will be useful,
  # but WITHOUT ANY WARRANTY; without even the implied warranty of
  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  # Library General Public License for more details.
  # 
  # You should have received a copy of the GNU Library General Public
  # License along with this library; if not, write to the
  # Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  # Boston, MA 02111-1307, USA.
  ###############################################################################
  # changes 29Apr2011
  # if the input image is a multi-band one, use all 
  # the channels in building the stack
  # anssi.pekkarinen@fao.org
  
  
  try:
      from osgeo import gdal
  except ImportError:
      import gdal
  
  try:
      progress = gdal.TermProgress_nocb
  except:
      progress = gdal.TermProgress
  
  
  import sys
  import math
  
  __version__ = '$id$'[5:-1]
  verbose = 0
  quiet = 0
  
  
  # =============================================================================
  def raster_copy( s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
                   t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
                   nodata=None ):
  
      if nodata is not None:
          return raster_copy_with_nodata(
              s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
              t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
              nodata )
  
      if verbose != 0:
          print('Copy %d,%d,%d,%d to %d,%d,%d,%d.' \
                % (s_xoff, s_yoff, s_xsize, s_ysize,
               t_xoff, t_yoff, t_xsize, t_ysize ))
  
      s_band = s_fh.GetRasterBand( s_band_n )
      t_band = t_fh.GetRasterBand( t_band_n )
  
      data = s_band.ReadRaster( s_xoff, s_yoff, s_xsize, s_ysize,
                               t_xsize, t_ysize, t_band.DataType )
      t_band.WriteRaster( t_xoff, t_yoff, t_xsize, t_ysize,
                          data, t_xsize, t_ysize, t_band.DataType )
          
  
      return 0
      
  # =============================================================================
  def raster_copy_with_nodata( s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
                               t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
                               nodata ):
      try:
          import numpy as Numeric
      except ImportError:
          import Numeric
      
      if verbose != 0:
          print('Copy %d,%d,%d,%d to %d,%d,%d,%d.' \
                % (s_xoff, s_yoff, s_xsize, s_ysize,
               t_xoff, t_yoff, t_xsize, t_ysize ))
  
      s_band = s_fh.GetRasterBand( s_band_n )
      t_band = t_fh.GetRasterBand( t_band_n )
  
      data_src = s_band.ReadAsArray( s_xoff, s_yoff, s_xsize, s_ysize,
                                     t_xsize, t_ysize )
      data_dst = t_band.ReadAsArray( t_xoff, t_yoff, t_xsize, t_ysize )
  
      nodata_test = Numeric.equal(data_src,nodata)
      to_write = Numeric.choose( nodata_test, (data_src, data_dst) )
  
      t_band.WriteArray( to_write, t_xoff, t_yoff )
  
      return 0
      
  # =============================================================================
  def names_to_fileinfos( names ):
      """
      Translate a list of GDAL filenames, into file_info objects.
  
      names -- list of valid GDAL dataset names.
  
      Returns a list of file_info objects.  There may be less file_info objects
      than names if some of the names could not be opened as GDAL files.
      """
      
      file_infos = []
      for name in names:
          fi = file_info()
          if fi.init_from_name( name ) == 1:
              file_infos.append( fi )
  
      return file_infos
  
  # *****************************************************************************
  class file_info:
      """A class holding information about a GDAL file."""
  
      def init_from_name(self, filename):
          """
          Initialize file_info from filename
  
          filename -- Name of file to read.
  
          Returns 1 on success or 0 if the file can't be opened.
          """
          fh = gdal.Open( filename )
          if fh is None:
              return 0
  
          self.filename = filename
          self.bands = fh.RasterCount
          self.xsize = fh.RasterXSize
          self.ysize = fh.RasterYSize
          self.band_type = fh.GetRasterBand(1).DataType
          self.projection = fh.GetProjection()
          self.geotransform = fh.GetGeoTransform()
          self.ulx = self.geotransform[0]
          self.uly = self.geotransform[3]
          self.lrx = self.ulx + self.geotransform[1] * self.xsize
          self.lry = self.uly + self.geotransform[5] * self.ysize
  
          ct = fh.GetRasterBand(1).GetRasterColorTable()
          if ct is not None:
              self.ct = ct.Clone()
          else:
              self.ct = None
  
          return 1
  
      def report( self ):
          print('Filename: '+ self.filename)
          print('File Size: %dx%dx%d' \
                % (self.xsize, self.ysize, self.bands))
          print('Pixel Size: %f x %f' \
                % (self.geotransform[1],self.geotransform[5]))
          print('UL:(%f,%f)   LR:(%f,%f)' \
                % (self.ulx,self.uly,self.lrx,self.lry))
  
      def copy_into( self, t_fh, s_band = 1, t_band = 1, nodata_arg=None ):
          """
          Copy this files image into target file.
  
          This method will compute the overlap area of the file_info objects
          file, and the target gdal.Dataset object, and copy the image data
          for the common window area.  It is assumed that the files are in
          a compatible projection ... no checking or warping is done.  However,
          if the destination file is a different resolution, or different
          image pixel type, the appropriate resampling and conversions will
          be done (using normal GDAL promotion/demotion rules).
  
          t_fh -- gdal.Dataset object for the file into which some or all
          of this file may be copied.
  
          Returns 1 on success (or if nothing needs to be copied), and zero one
          failure.
          """
          t_geotransform = t_fh.GetGeoTransform()
          t_ulx = t_geotransform[0]
          t_uly = t_geotransform[3]
          t_lrx = t_geotransform[0] + t_fh.RasterXSize * t_geotransform[1]
          t_lry = t_geotransform[3] + t_fh.RasterYSize * t_geotransform[5]
  
          # figure out intersection region
          tgw_ulx = max(t_ulx,self.ulx)
          tgw_lrx = min(t_lrx,self.lrx)
          if t_geotransform[5] < 0:
              tgw_uly = min(t_uly,self.uly)
              tgw_lry = max(t_lry,self.lry)
          else:
              tgw_uly = max(t_uly,self.uly)
              tgw_lry = min(t_lry,self.lry)
          
          # do they even intersect?
          if tgw_ulx >= tgw_lrx:
              return 1
          if t_geotransform[5] < 0 and tgw_uly <= tgw_lry:
              return 1
          if t_geotransform[5] > 0 and tgw_uly >= tgw_lry:
              return 1
              
          # compute target window in pixel coordinates.
          tw_xoff = int((tgw_ulx - t_geotransform[0]) / t_geotransform[1] + 0.1)
          tw_yoff = int((tgw_uly - t_geotransform[3]) / t_geotransform[5] + 0.1)
          tw_xsize = int((tgw_lrx - t_geotransform[0])/t_geotransform[1] + 0.5) \
                     - tw_xoff
          tw_ysize = int((tgw_lry - t_geotransform[3])/t_geotransform[5] + 0.5) \
                     - tw_yoff
  
          if tw_xsize < 1 or tw_ysize < 1:
              return 1
  
          # Compute source window in pixel coordinates.
          sw_xoff = int((tgw_ulx - self.geotransform[0]) / self.geotransform[1])
          sw_yoff = int((tgw_uly - self.geotransform[3]) / self.geotransform[5])
          sw_xsize = int((tgw_lrx - self.geotransform[0]) \
                         / self.geotransform[1] + 0.5) - sw_xoff
          sw_ysize = int((tgw_lry - self.geotransform[3]) \
                         / self.geotransform[5] + 0.5) - sw_yoff
  
          if sw_xsize < 1 or sw_ysize < 1:
              return 1
  
          # Open the source file, and copy the selected region.
          s_fh = gdal.Open( self.filename )
  
          return \
              raster_copy( s_fh, sw_xoff, sw_yoff, sw_xsize, sw_ysize, s_band,
                           t_fh, tw_xoff, tw_yoff, tw_xsize, tw_ysize, t_band,
                           nodata_arg )
  
  
  # =============================================================================
  def Usage():
      print('Usage: gdal_merge.py [-o out_filename] [-of out_format] [-co NAME=VALUE]*')
      print('                     [-ps pixelsize_x pixelsize_y] [-tap] [-separate] [-q] [-v] [-pct]')
      print('                     [-ul_lr ulx uly lrx lry] [-init "value [value...]"]')
      print('                     [-n nodata_value] [-a_nodata output_nodata_value]')
      print('                     [-ot datatype] [-createonly] input_files')
      print('                     [--help-general]')
      print('')
  
  # =============================================================================
  #
  # Program mainline.
  #
  
  def main( argv=None ):
  
      global verbose, quiet
      verbose = 0
      quiet = 0
      names = []
      format = 'GTiff'
      out_file = 'out.tif'
  
      ulx = None
      psize_x = None
      separate = 0
      copy_pct = 0
      nodata = None
      a_nodata = None
      create_options = []
      pre_init = []
      band_type = None
      createonly = 0
      bTargetAlignedPixels = False
      
      gdal.AllRegister()
      if argv is None:
          argv = sys.argv
      argv = gdal.GeneralCmdLineProcessor( argv )
      if argv is None:
          sys.exit( 0 )
  
      # Parse command line arguments.
      i = 1
      while i < len(argv):
          arg = argv[i]
  
          if arg == '-o':
              i = i + 1
              out_file = argv[i]
  
          elif arg == '-v':
              verbose = 1
  
          elif arg == '-q' or arg == '-quiet':
              quiet = 1
  
          elif arg == '-createonly':
              createonly = 1
  
          elif arg == '-separate':
              separate = 1
  
          elif arg == '-seperate':
              separate = 1
  
          elif arg == '-pct':
              copy_pct = 1
  
          elif arg == '-ot':
              i = i + 1
              band_type = gdal.GetDataTypeByName( argv[i] )
              if band_type == gdal.GDT_Unknown:
                  print('Unknown GDAL data type: %s' % argv[i])
                  sys.exit( 1 )
  
          elif arg == '-init':
              i = i + 1
              str_pre_init = argv[i].split()
              for x in str_pre_init:
                  pre_init.append(float(x))
  
          elif arg == '-n':
              i = i + 1
              nodata = float(argv[i])
  
          elif arg == '-a_nodata':
              i = i + 1
              a_nodata = float(argv[i])
  
          elif arg == '-f':
              # for backward compatibility.
              i = i + 1
              format = argv[i]
  
          elif arg == '-of':
              i = i + 1
              format = argv[i]
  
          elif arg == '-co':
              i = i + 1
              create_options.append( argv[i] )
  
          elif arg == '-ps':
              psize_x = float(argv[i+1])
              psize_y = -1 * abs(float(argv[i+2]))
              i = i + 2
  
          elif arg == '-tap':
              bTargetAlignedPixels = True
  
          elif arg == '-ul_lr':
              ulx = float(argv[i+1])
              uly = float(argv[i+2])
              lrx = float(argv[i+3])
              lry = float(argv[i+4])
              i = i + 4
  
          elif arg[:1] == '-':
              print('Unrecognised command option: %s' % arg)
              Usage()
              sys.exit( 1 )
  
          else:
              names.append(arg)
  
          i = i + 1
  
      if len(names) == 0:
          print('No input files selected.')
          Usage()
          sys.exit( 1 )
  
      Driver = gdal.GetDriverByName(format)
      if Driver is None:
          print('Format driver %s not found, pick a supported driver.' % format)
          sys.exit( 1 )
  
      DriverMD = Driver.GetMetadata()
      if 'DCAP_CREATE' not in DriverMD:
          print('Format driver %s does not support creation and piecewise writing.\nPlease select a format that does, such as GTiff (the default) or HFA (Erdas Imagine).' % format)
          sys.exit( 1 )
  
      # Collect information on all the source files.
      file_infos = names_to_fileinfos( names )
  
      if ulx is None:
          ulx = file_infos[0].ulx
          uly = file_infos[0].uly
          lrx = file_infos[0].lrx
          lry = file_infos[0].lry
          
          for fi in file_infos:
              ulx = min(ulx, fi.ulx)
              uly = max(uly, fi.uly)
              lrx = max(lrx, fi.lrx)
              lry = min(lry, fi.lry)
  
      if psize_x is None:
          psize_x = file_infos[0].geotransform[1]
          psize_y = file_infos[0].geotransform[5]
  
      if band_type is None:
          band_type = file_infos[0].band_type
  
      # Try opening as an existing file.
      gdal.PushErrorHandler( 'CPLQuietErrorHandler' )
      t_fh = gdal.Open( out_file, gdal.GA_Update )
      gdal.PopErrorHandler()
      
      # Create output file if it does not already exist.
      if t_fh is None:
      
          if bTargetAlignedPixels:
              ulx = math.floor(ulx / psize_x) * psize_x
              lrx = math.ceil(lrx / psize_x) * psize_x
              lry = math.floor(lry / -psize_y) * -psize_y
              uly = math.ceil(uly / -psize_y) * -psize_y
      
          geotransform = [ulx, psize_x, 0, uly, 0, psize_y]
  
          xsize = int((lrx - ulx) / geotransform[1] + 0.5)
          ysize = int((lry - uly) / geotransform[5] + 0.5)
  
  
          if separate != 0:
              bands=0
  
              for fi in file_infos:
                  bands=bands + fi.bands
          else:
              bands = file_infos[0].bands
  
  
          t_fh = Driver.Create( out_file, xsize, ysize, bands,
                                band_type, create_options )
          if t_fh is None:
              print('Creation failed, terminating gdal_merge.')
              sys.exit( 1 )
              
          t_fh.SetGeoTransform( geotransform )
          t_fh.SetProjection( file_infos[0].projection )
  
          if copy_pct:
              t_fh.GetRasterBand(1).SetRasterColorTable(file_infos[0].ct)
      else:
          if separate != 0:
              bands=0
              for fi in file_infos:
                  bands=bands + fi.bands            
              if t_fh.RasterCount < bands :
                  print('Existing output file has less bands than the input files. You should delete it before. Terminating gdal_merge.')
                  sys.exit( 1 )
          else:
              bands = min(file_infos[0].bands,t_fh.RasterCount)
  
      # Do we need to set nodata value ?
      if a_nodata is not None:
          for i in range(t_fh.RasterCount):
              t_fh.GetRasterBand(i+1).SetNoDataValue(a_nodata)
  
      # Do we need to pre-initialize the whole mosaic file to some value?
      if pre_init is not None:
          if t_fh.RasterCount <= len(pre_init):
              for i in range(t_fh.RasterCount):
                  t_fh.GetRasterBand(i+1).Fill( pre_init[i] )
          elif len(pre_init) == 1:
              for i in range(t_fh.RasterCount):
                  t_fh.GetRasterBand(i+1).Fill( pre_init[0] )
  
      # Copy data from source files into output file.
      t_band = 1
  
      if quiet == 0 and verbose == 0:
          progress( 0.0 )
      fi_processed = 0
      
      for fi in file_infos:
          if createonly != 0:
              continue
          
          if verbose != 0:
              print("")
              print("Processing file %5d of %5d, %6.3f%% completed." \
                    % (fi_processed+1,len(file_infos),
                       fi_processed * 100.0 / len(file_infos)) )
              fi.report()
  
          if separate == 0 :
              for band in range(1, bands+1):
                  fi.copy_into( t_fh, band, band, nodata )
          else:
              for band in range(1, fi.bands+1):
                  fi.copy_into( t_fh, band, t_band, nodata )
                  t_band = t_band+1
              
          fi_processed = fi_processed+1
          if quiet == 0 and verbose == 0:
              progress( fi_processed / float(len(file_infos))  )
      
      # Force file to be closed.
      t_fh = None
  
  if __name__ == '__main__':
      sys.exit(main())
  