prl900 · November 7, 2018 23:34
diff --git a/gdal_merge.py b/gdal_merge.py
 #!/g/data/xc0/software/python/miniconda2/bin/python
 ###############################################################################
 # $Id: gdal_merge.py 33790 2016-03-26 12:42:12Z goatbar $
 #
 # Project:  InSAR Peppers
 # Purpose:  Module to extract data from many rasters into one output.
 # Author:   Frank Warmerdam, [email protected]
 #
 ###############################################################################
 # 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.
 # [email protected]

 import math
 import sys
 import time

 from osgeo import gdal

 try:
    progress = gdal.TermProgress_nocb
 except:
    progress = gdal.TermProgress

 __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 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 ))

    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 )

    s_band = s_fh.GetRasterBand( s_band_n )
    m_band = None
    # Works only in binary mode and doesn't take into account
    # intermediate transparency values for compositing.
    if s_band.GetMaskFlags() != gdal.GMF_ALL_VALID:
        m_band = s_band.GetMaskBand()
    elif s_band.GetColorInterpretation() == gdal.GCI_AlphaBand:
        m_band = s_band
    if m_band is not None:
        return raster_copy_with_mask(
            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,
            m_band )

    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

    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 raster_copy_with_mask( 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,
                           m_band ):
    try:
        import numpy as Numeric
    except ImportError:
        import Numeric

    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_mask = m_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 )

    mask_test = Numeric.equal(data_mask, 0)
    to_write = Numeric.choose( mask_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
    start_time = time.time()

    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('Unrecognized 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 in %d minutes."
                  % (fi_processed+1,len(file_infos),
                     fi_processed * 100.0 / len(file_infos),
                     int(round((time.time() - start_time)/60.0)) ))
            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())
	#!/g/data/xc0/software/python/miniconda2/bin/python
	###############################################################################
	# $Id: gdal_merge.py 33790 2016-03-26 12:42:12Z goatbar $
	#
	# Project: InSAR Peppers
	# Purpose: Module to extract data from many rasters into one output.
	# Author: Frank Warmerdam, [email protected]
	#
	###############################################################################
	# 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.
	# [email protected]

	import math
	import sys
	import time

	from osgeo import gdal

	try:
	progress = gdal.TermProgress_nocb
	except:
	progress = gdal.TermProgress

	__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 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 ))

	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 )

	s_band = s_fh.GetRasterBand( s_band_n )
	m_band = None
	# Works only in binary mode and doesn't take into account
	# intermediate transparency values for compositing.
	if s_band.GetMaskFlags() != gdal.GMF_ALL_VALID:
	m_band = s_band.GetMaskBand()
	elif s_band.GetColorInterpretation() == gdal.GCI_AlphaBand:
	m_band = s_band
	if m_band is not None:
	return raster_copy_with_mask(
	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,
	m_band )

	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

	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 raster_copy_with_mask( 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,
	m_band ):
	try:
	import numpy as Numeric
	except ImportError:
	import Numeric

	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_mask = m_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 )

	mask_test = Numeric.equal(data_mask, 0)
	to_write = Numeric.choose( mask_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
	start_time = time.time()

	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('Unrecognized 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 in %d minutes."
	% (fi_processed+1,len(file_infos),
	fi_processed * 100.0 / len(file_infos),
	int(round((time.time() - start_time)/60.0)) ))
	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())