Package osgeo :: Package utils :: Module gdal_merge |
[hide private]
[frames] | no frames]

# Source Code for Module osgeo.utils.gdal_merge

      1  #!/usr/bin/env python3
      2  ###############################################################################
      3  # $Id$
      4  #
      5  # Project:  InSAR Peppers
      6  # Purpose:  Module to extract data from many rasters into one output.
      7  # Author:   Frank Warmerdam, warmerdam@pobox.com
      8  #
      9  ###############################################################################
     10  # Copyright (c) 2000, Atlantis Scientific Inc. (www.atlsci.com)
     11  # Copyright (c) 2009-2011, Even Rouault <even dot rouault at spatialys.com>
     12  #
     13  # This library is free software; you can redistribute it and/or
     14  # modify it under the terms of the GNU Library General Public
     15  # License as published by the Free Software Foundation; either
     16  # version 2 of the License, or (at your option) any later version.
     17  #
     18  # This library is distributed in the hope that it will be useful,
     19  # but WITHOUT ANY WARRANTY; without even the implied warranty of
     20  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     21  # Library General Public License for more details.
     22  #
     23  # You should have received a copy of the GNU Library General Public
     24  # License along with this library; if not, write to the
     25  # Free Software Foundation, Inc., 59 Temple Place - Suite 330,
     26  # Boston, MA 02111-1307, USA.
     27  ###############################################################################
     28  # changes 29Apr2011
     29  # If the input image is a multi-band one, use all the channels in
     30  # building the stack.
     31  # anssi.pekkarinen@fao.org
     32
     33  import math
     34  import os.path
     35  import sys
     36  import time
     37
     38  from osgeo import gdal
     39
     40  progress = gdal.TermProgress_nocb
     41
     42  __version__ = '$id$'[5:-1]
     43
     44
     45 -def DoesDriverHandleExtension(drv, ext):
     46      exts = drv.GetMetadataItem(gdal.DMD_EXTENSIONS)
     47      return exts is not None and exts.lower().find(ext.lower()) >= 0
     48
     49
     50 -def GetExtension(filename):
     51      ext = os.path.splitext(filename)[1]
     52      if ext.startswith('.'):
     53          ext = ext[1:]
     54      return ext
     55
     56
     57 -def GetOutputDriversFor(filename):
     58      drv_list = []
     59      ext = GetExtension(filename)
     60      for i in range(gdal.GetDriverCount()):
     61          drv = gdal.GetDriver(i)
     62          if (drv.GetMetadataItem(gdal.DCAP_CREATE) is not None or
     63              drv.GetMetadataItem(gdal.DCAP_CREATECOPY) is not None) and \
     64             drv.GetMetadataItem(gdal.DCAP_RASTER) is not None:
     65              if ext and DoesDriverHandleExtension(drv, ext):
     66                  drv_list.append(drv.ShortName)
     67              else:
     68                  prefix = drv.GetMetadataItem(gdal.DMD_CONNECTION_PREFIX)
     69                  if prefix is not None and filename.lower().startswith(prefix.lower()):
     70                      drv_list.append(drv.ShortName)
     71
     72      # GMT is registered before netCDF for opening reasons, but we want
     73      # netCDF to be used by default for output.
     74      if ext.lower() == 'nc' and not drv_list and \
     75         drv_list[0].upper() == 'GMT' and drv_list[1].upper() == 'NETCDF':
     76          drv_list = ['NETCDF', 'GMT']
     77
     78      return drv_list
     79
     80
     81 -def GetOutputDriverFor(filename):
     82      drv_list = GetOutputDriversFor(filename)
     83      ext = GetExtension(filename)
     84      if not drv_list:
     85          if not ext:
     86              return 'GTiff'
     87          else:
     88              raise Exception("Cannot guess driver for %s" % filename)
     89      elif len(drv_list) > 1:
     90          print("Several drivers matching %s extension. Using %s" % (ext if ext else '', drv_list[0]))
     91      return drv_list[0]
     92
     93
     94  # =============================================================================
     95 -def raster_copy(s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
     96                  t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
     97                  nodata=None, verbose=0):
     98
     99      if verbose != 0:
    100          print('Copy %d,%d,%d,%d to %d,%d,%d,%d.'
    101                % (s_xoff, s_yoff, s_xsize, s_ysize,
    102                   t_xoff, t_yoff, t_xsize, t_ysize))
    103
    104      if nodata is not None:
    105          return raster_copy_with_nodata(
    106              s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
    107              t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
    108              nodata)
    109
    110      s_band = s_fh.GetRasterBand(s_band_n)
    111      m_band = None
    112      # Works only in binary mode and doesn't take into account
    113      # intermediate transparency values for compositing.
    114      if s_band.GetMaskFlags() != gdal.GMF_ALL_VALID:
    115          m_band = s_band.GetMaskBand()
    116      elif s_band.GetColorInterpretation() == gdal.GCI_AlphaBand:
    117          m_band = s_band
    118      if m_band is not None:
    119          return raster_copy_with_mask(
    120              s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
    121              t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
    122              m_band)
    123
    124      s_band = s_fh.GetRasterBand(s_band_n)
    125      t_band = t_fh.GetRasterBand(t_band_n)
    126
    127      data = s_band.ReadRaster(s_xoff, s_yoff, s_xsize, s_ysize,
    128                               t_xsize, t_ysize, t_band.DataType)
    129      t_band.WriteRaster(t_xoff, t_yoff, t_xsize, t_ysize,
    130                         data, t_xsize, t_ysize, t_band.DataType)
    131
    132      return 0
    133
    134  # =============================================================================
    135
    136
    137 -def raster_copy_with_nodata(s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
    138                              t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
    139                              nodata):
    140      try:
    141          import numpy as Numeric
    142      except ImportError:
    143          import Numeric
    144
    145      s_band = s_fh.GetRasterBand(s_band_n)
    146      t_band = t_fh.GetRasterBand(t_band_n)
    147
    148      data_src = s_band.ReadAsArray(s_xoff, s_yoff, s_xsize, s_ysize,
    149                                    t_xsize, t_ysize)
    150      data_dst = t_band.ReadAsArray(t_xoff, t_yoff, t_xsize, t_ysize)
    151
    152      if not Numeric.isnan(nodata):
    153          nodata_test = Numeric.equal(data_src, nodata)
    154      else:
    155          nodata_test = Numeric.isnan(data_src)
    156
    157      to_write = Numeric.choose(nodata_test, (data_src, data_dst))
    158
    159      t_band.WriteArray(to_write, t_xoff, t_yoff)
    160
    161      return 0
    162
    163  # =============================================================================
    164
    165
    166 -def raster_copy_with_mask(s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
    167                            t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
    168                            m_band):
    169      try:
    170          import numpy as Numeric
    171      except ImportError:
    172          import Numeric
    173
    174      s_band = s_fh.GetRasterBand(s_band_n)
    175      t_band = t_fh.GetRasterBand(t_band_n)
    176
    177      data_src = s_band.ReadAsArray(s_xoff, s_yoff, s_xsize, s_ysize,
    178                                    t_xsize, t_ysize)
    179      data_mask = m_band.ReadAsArray(s_xoff, s_yoff, s_xsize, s_ysize,
    180                                     t_xsize, t_ysize)
    181      data_dst = t_band.ReadAsArray(t_xoff, t_yoff, t_xsize, t_ysize)
    182
    183      mask_test = Numeric.equal(data_mask, 0)
    184      to_write = Numeric.choose(mask_test, (data_src, data_dst))
    185
    186      t_band.WriteArray(to_write, t_xoff, t_yoff)
    187
    188      return 0
    189
    190  # =============================================================================
    191
    192
    193 -def names_to_fileinfos(names):
    194      """
    195      Translate a list of GDAL filenames, into file_info objects.
    196
    197      names -- list of valid GDAL dataset names.
    198
    199      Returns a list of file_info objects.  There may be less file_info objects
    200      than names if some of the names could not be opened as GDAL files.
    201      """
    202
    203      file_infos = []
    204      for name in names:
    205          fi = file_info()
    206          if fi.init_from_name(name) == 1:
    207              file_infos.append(fi)
    208
    209      return file_infos
    210
    211  # *****************************************************************************
    212
    213
    214 -class file_info(object):
    215      """A class holding information about a GDAL file."""
    216
    217 -    def __init__(self):
    218          self.band_type = None
    219          self.bands = None
    220          self.ct = None
    221          self.filename = None
    222          self.geotransform = None
    223          self.lrx = None
    224          self.lry = None
    225          self.projection = None
    226          self.ulx = None
    227          self.uly = None
    228          self.xsize = None
    229          self.ysize = None
    230
    231 -    def init_from_name(self, filename):
    232          """
    233          Initialize file_info from filename
    234
    235          filename -- Name of file to read.
    236
    237          Returns 1 on success or 0 if the file can't be opened.
    238          """
    239          fh = gdal.Open(filename)
    240          if fh is None:
    241              return 0
    242
    243          self.filename = filename
    244          self.bands = fh.RasterCount
    245          self.xsize = fh.RasterXSize
    246          self.ysize = fh.RasterYSize
    247          self.band_type = fh.GetRasterBand(1).DataType
    248          self.projection = fh.GetProjection()
    249          self.geotransform = fh.GetGeoTransform()
    250          self.ulx = self.geotransform[0]
    251          self.uly = self.geotransform[3]
    252          self.lrx = self.ulx + self.geotransform[1] * self.xsize
    253          self.lry = self.uly + self.geotransform[5] * self.ysize
    254
    255          ct = fh.GetRasterBand(1).GetRasterColorTable()
    256          if ct is not None:
    257              self.ct = ct.Clone()
    258          else:
    259              self.ct = None
    260
    261          return 1
    262
    263 -    def report(self):
    264          print('Filename: ' + self.filename)
    265          print('File Size: %dx%dx%d'
    266                % (self.xsize, self.ysize, self.bands))
    267          print('Pixel Size: %f x %f'
    268                % (self.geotransform[1], self.geotransform[5]))
    269          print('UL:(%f,%f)   LR:(%f,%f)'
    270                % (self.ulx, self.uly, self.lrx, self.lry))
    271
    272 -    def copy_into(self, t_fh, s_band=1, t_band=1, nodata_arg=None, verbose=0):
    273          """
    274          Copy this files image into target file.
    275
    276          This method will compute the overlap area of the file_info objects
    277          file, and the target gdal.Dataset object, and copy the image data
    278          for the common window area.  It is assumed that the files are in
    279          a compatible projection ... no checking or warping is done.  However,
    280          if the destination file is a different resolution, or different
    281          image pixel type, the appropriate resampling and conversions will
    282          be done (using normal GDAL promotion/demotion rules).
    283
    284          t_fh -- gdal.Dataset object for the file into which some or all
    285          of this file may be copied.
    286
    287          Returns 1 on success (or if nothing needs to be copied), and zero one
    288          failure.
    289          """
    290          t_geotransform = t_fh.GetGeoTransform()
    291          t_ulx = t_geotransform[0]
    292          t_uly = t_geotransform[3]
    293          t_lrx = t_geotransform[0] + t_fh.RasterXSize * t_geotransform[1]
    294          t_lry = t_geotransform[3] + t_fh.RasterYSize * t_geotransform[5]
    295
    296          # figure out intersection region
    297          tgw_ulx = max(t_ulx, self.ulx)
    298          tgw_lrx = min(t_lrx, self.lrx)
    299          if t_geotransform[5] < 0:
    300              tgw_uly = min(t_uly, self.uly)
    301              tgw_lry = max(t_lry, self.lry)
    302          else:
    303              tgw_uly = max(t_uly, self.uly)
    304              tgw_lry = min(t_lry, self.lry)
    305
    306          # do they even intersect?
    307          if tgw_ulx >= tgw_lrx:
    308              return 1
    309          if t_geotransform[5] < 0 and tgw_uly <= tgw_lry:
    310              return 1
    311          if t_geotransform[5] > 0 and tgw_uly >= tgw_lry:
    312              return 1
    313
    314          # compute target window in pixel coordinates.
    315          tw_xoff = int((tgw_ulx - t_geotransform[0]) / t_geotransform[1] + 0.1)
    316          tw_yoff = int((tgw_uly - t_geotransform[3]) / t_geotransform[5] + 0.1)
    317          tw_xsize = int((tgw_lrx - t_geotransform[0]) / t_geotransform[1] + 0.5) \
    318              - tw_xoff
    319          tw_ysize = int((tgw_lry - t_geotransform[3]) / t_geotransform[5] + 0.5) \
    320              - tw_yoff
    321
    322          if tw_xsize < 1 or tw_ysize < 1:
    323              return 1
    324
    325          # Compute source window in pixel coordinates.
    326          sw_xoff = int((tgw_ulx - self.geotransform[0]) / self.geotransform[1] + 0.1)
    327          sw_yoff = int((tgw_uly - self.geotransform[3]) / self.geotransform[5] + 0.1)
    328          sw_xsize = int((tgw_lrx - self.geotransform[0]) /
    329                         self.geotransform[1] + 0.5) - sw_xoff
    330          sw_ysize = int((tgw_lry - self.geotransform[3]) /
    331                         self.geotransform[5] + 0.5) - sw_yoff
    332
    333          if sw_xsize < 1 or sw_ysize < 1:
    334              return 1
    335
    336          # Open the source file, and copy the selected region.
    337          s_fh = gdal.Open(self.filename)
    338
    339          return raster_copy(s_fh, sw_xoff, sw_yoff, sw_xsize, sw_ysize, s_band,
    340                             t_fh, tw_xoff, tw_yoff, tw_xsize, tw_ysize, t_band,
    341                             nodata_arg, verbose)
    342
    343
    344  # =============================================================================
    345 -def Usage():
    346      print('Usage: gdal_merge.py [-o out_filename] [-of out_format] [-co NAME=VALUE]*')
    347      print('                     [-ps pixelsize_x pixelsize_y] [-tap] [-separate] [-q] [-v] [-pct]')
    348      print('                     [-ul_lr ulx uly lrx lry] [-init "value [value...]"]')
    349      print('                     [-n nodata_value] [-a_nodata output_nodata_value]')
    350      print('                     [-ot datatype] [-createonly] input_files')
    351      print('                     [--help-general]')
    352      print('')
    353
    354
    355 -def main(argv=None):
    356      verbose = 0
    357      quiet = 0
    358      names = []
    359      frmt = None
    360      out_file = 'out.tif'
    361
    362      ulx = None
    363      psize_x = None
    364      separate = 0
    365      copy_pct = 0
    366      nodata = None
    367      a_nodata = None
    368      create_options = []
    369      pre_init = []
    370      band_type = None
    371      createonly = 0
    372      bTargetAlignedPixels = False
    373      start_time = time.time()
    374
    375      gdal.AllRegister()
    376      if argv is None:
    377          argv = argv
    378      argv = gdal.GeneralCmdLineProcessor(argv)
    379      if argv is None:
    380          sys.exit(0)
    381
    382      # Parse command line arguments.
    383      i = 1
    384      while i < len(argv):
    385          arg = argv[i]
    386
    387          if arg == '-o':
    388              i = i + 1
    389              out_file = argv[i]
    390
    391          elif arg == '-v':
    392              verbose = 1
    393
    394          elif arg == '-q' or arg == '-quiet':
    395              quiet = 1
    396
    397          elif arg == '-createonly':
    398              createonly = 1
    399
    400          elif arg == '-separate':
    401              separate = 1
    402
    403          elif arg == '-seperate':
    404              separate = 1
    405
    406          elif arg == '-pct':
    407              copy_pct = 1
    408
    409          elif arg == '-ot':
    410              i = i + 1
    411              band_type = gdal.GetDataTypeByName(argv[i])
    412              if band_type == gdal.GDT_Unknown:
    413                  print('Unknown GDAL data type: %s' % argv[i])
    414                  sys.exit(1)
    415
    416          elif arg == '-init':
    417              i = i + 1
    418              str_pre_init = argv[i].split()
    419              for x in str_pre_init:
    420                  pre_init.append(float(x))
    421
    422          elif arg == '-n':
    423              i = i + 1
    424              nodata = float(argv[i])
    425
    426          elif arg == '-a_nodata':
    427              i = i + 1
    428              a_nodata = float(argv[i])
    429
    430          elif arg == '-f' or arg == '-of':
    431              i = i + 1
    432              frmt = argv[i]
    433
    434          elif arg == '-co':
    435              i = i + 1
    436              create_options.append(argv[i])
    437
    438          elif arg == '-ps':
    439              psize_x = float(argv[i + 1])
    440              psize_y = -1 * abs(float(argv[i + 2]))
    441              i = i + 2
    442
    443          elif arg == '-tap':
    444              bTargetAlignedPixels = True
    445
    446          elif arg == '-ul_lr':
    447              ulx = float(argv[i + 1])
    448              uly = float(argv[i + 2])
    449              lrx = float(argv[i + 3])
    450              lry = float(argv[i + 4])
    451              i = i + 4
    452
    453          elif arg[:1] == '-':
    454              print('Unrecognized command option: %s' % arg)
    455              Usage()
    456              sys.exit(1)
    457
    458          else:
    459              names.append(arg)
    460
    461          i = i + 1
    462
    463      if not names:
    464          print('No input files selected.')
    465          Usage()
    466          sys.exit(1)
    467
    468      if frmt is None:
    469          frmt = GetOutputDriverFor(out_file)
    470
    471      Driver = gdal.GetDriverByName(frmt)
    472      if Driver is None:
    473          print('Format driver %s not found, pick a supported driver.' % frmt)
    474          sys.exit(1)
    475
    476      DriverMD = Driver.GetMetadata()
    477      if 'DCAP_CREATE' not in DriverMD:
    478          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).' % frmt)
    479          sys.exit(1)
    480
    481      # Collect information on all the source files.
    482      file_infos = names_to_fileinfos(names)
    483
    484      if ulx is None:
    485          ulx = file_infos[0].ulx
    486          uly = file_infos[0].uly
    487          lrx = file_infos[0].lrx
    488          lry = file_infos[0].lry
    489
    490          for fi in file_infos:
    491              ulx = min(ulx, fi.ulx)
    492              uly = max(uly, fi.uly)
    493              lrx = max(lrx, fi.lrx)
    494              lry = min(lry, fi.lry)
    495
    496      if psize_x is None:
    497          psize_x = file_infos[0].geotransform[1]
    498          psize_y = file_infos[0].geotransform[5]
    499
    500      if band_type is None:
    501          band_type = file_infos[0].band_type
    502
    503      # Try opening as an existing file.
    504      gdal.PushErrorHandler('CPLQuietErrorHandler')
    505      t_fh = gdal.Open(out_file, gdal.GA_Update)
    506      gdal.PopErrorHandler()
    507
    508      # Create output file if it does not already exist.
    509      if t_fh is None:
    510
    511          if bTargetAlignedPixels:
    512              ulx = math.floor(ulx / psize_x) * psize_x
    513              lrx = math.ceil(lrx / psize_x) * psize_x
    514              lry = math.floor(lry / -psize_y) * -psize_y
    515              uly = math.ceil(uly / -psize_y) * -psize_y
    516
    517          geotransform = [ulx, psize_x, 0, uly, 0, psize_y]
    518
    519          xsize = int((lrx - ulx) / geotransform[1] + 0.5)
    520          ysize = int((lry - uly) / geotransform[5] + 0.5)
    521
    522          if separate != 0:
    523              bands = 0
    524
    525              for fi in file_infos:
    526                  bands = bands + fi.bands
    527          else:
    528              bands = file_infos[0].bands
    529
    530          t_fh = Driver.Create(out_file, xsize, ysize, bands,
    531                               band_type, create_options)
    532          if t_fh is None:
    533              print('Creation failed, terminating gdal_merge.')
    534              sys.exit(1)
    535
    536          t_fh.SetGeoTransform(geotransform)
    537          t_fh.SetProjection(file_infos[0].projection)
    538
    539          if copy_pct:
    540              t_fh.GetRasterBand(1).SetRasterColorTable(file_infos[0].ct)
    541      else:
    542          if separate != 0:
    543              bands = 0
    544              for fi in file_infos:
    545                  bands = bands + fi.bands
    546              if t_fh.RasterCount < bands:
    547                  print('Existing output file has less bands than the input files. You should delete it before. Terminating gdal_merge.')
    548                  sys.exit(1)
    549          else:
    550              bands = min(file_infos[0].bands, t_fh.RasterCount)
    551
    552      # Do we need to set nodata value ?
    553      if a_nodata is not None:
    554          for i in range(t_fh.RasterCount):
    555              t_fh.GetRasterBand(i + 1).SetNoDataValue(a_nodata)
    556
    557      # Do we need to pre-initialize the whole mosaic file to some value?
    558      if pre_init is not None:
    559          if t_fh.RasterCount <= len(pre_init):
    560              for i in range(t_fh.RasterCount):
    561                  t_fh.GetRasterBand(i + 1).Fill(pre_init[i])
    562          elif len(pre_init) == 1:
    563              for i in range(t_fh.RasterCount):
    564                  t_fh.GetRasterBand(i + 1).Fill(pre_init[0])
    565
    566      # Copy data from source files into output file.
    567      t_band = 1
    568
    569      if quiet == 0 and verbose == 0:
    570          progress(0.0)
    571      fi_processed = 0
    572
    573      for fi in file_infos:
    574          if createonly != 0:
    575              continue
    576
    577          if verbose != 0:
    578              print("")
    579              print("Processing file %5d of %5d, %6.3f%% completed in %d minutes."
    580                    % (fi_processed + 1, len(file_infos),
    581                       fi_processed * 100.0 / len(file_infos),
    582                       int(round((time.time() - start_time) / 60.0))))
    583              fi.report()
    584
    585          if separate == 0:
    586              for band in range(1, bands + 1):
    587                  fi.copy_into(t_fh, band, band, nodata, verbose)
    588          else:
    589              for band in range(1, fi.bands + 1):
    590                  fi.copy_into(t_fh, band, t_band, nodata, verbose)
    591                  t_band = t_band + 1
    592
    593          fi_processed = fi_processed + 1
    594          if quiet == 0 and verbose == 0:
    595              progress(fi_processed / float(len(file_infos)))
    596
    597      # Force file to be closed.
    598      t_fh = None
    599
    600
    601  if __name__ == '__main__':
    602      sys.exit(main(sys.argv))
    603