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

# Source Code for Module osgeo.utils.gdal_pansharpen

      1  #!/usr/bin/env python3
      2  # -*- coding: utf-8 -*-
      3  ###############################################################################
      4  # $Id$
      5  #
      6  #  Project:  GDAL scripts
      7  #  Purpose:  Perform a pansharpening operation
      8  #  Author:   Even Rouault <even.rouault at spatialys.com>
      9  #
     10  ###############################################################################
     11  #  Copyright (c) 2015, Even Rouault <even.rouault at spatialys.com>
     12  #
     13  #  Permission is hereby granted, free of charge, to any person obtaining a
     14  #  copy of this software and associated documentation files (the "Software"),
     15  #  to deal in the Software without restriction, including without limitation
     16  #  the rights to use, copy, modify, merge, publish, distribute, sublicense,
     17  #  and/or sell copies of the Software, and to permit persons to whom the
     18  #  Software is furnished to do so, subject to the following conditions:
     19  #
     20  #  The above copyright notice and this permission notice shall be included
     21  #  in all copies or substantial portions of the Software.
     22  #
     23  #  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
     24  #  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
     25  #  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
     26  #  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
     27  #  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
     28  #  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
     29  #  DEALINGS IN THE SOFTWARE.
     30  ###############################################################################
     31
     32  import os
     33  import os.path
     34  import sys
     35  from osgeo import gdal
     36
     37
     38 -def DoesDriverHandleExtension(drv, ext):
     39      exts = drv.GetMetadataItem(gdal.DMD_EXTENSIONS)
     40      return exts is not None and exts.lower().find(ext.lower()) >= 0
     41
     42
     43 -def GetExtension(filename):
     44      ext = os.path.splitext(filename)[1]
     45      if ext.startswith('.'):
     46          ext = ext[1:]
     47      return ext
     48
     49
     50 -def GetOutputDriversFor(filename):
     51      drv_list = []
     52      ext = GetExtension(filename)
     53      for i in range(gdal.GetDriverCount()):
     54          drv = gdal.GetDriver(i)
     55          if (drv.GetMetadataItem(gdal.DCAP_CREATE) is not None or
     56              drv.GetMetadataItem(gdal.DCAP_CREATECOPY) is not None) and \
     57             drv.GetMetadataItem(gdal.DCAP_RASTER) is not None:
     58              if ext and DoesDriverHandleExtension(drv, ext):
     59                  drv_list.append(drv.ShortName)
     60              else:
     61                  prefix = drv.GetMetadataItem(gdal.DMD_CONNECTION_PREFIX)
     62                  if prefix is not None and filename.lower().startswith(prefix.lower()):
     63                      drv_list.append(drv.ShortName)
     64
     65      # GMT is registered before netCDF for opening reasons, but we want
     66      # netCDF to be used by default for output.
     67      if ext.lower() == 'nc' and not drv_list and \
     68         drv_list[0].upper() == 'GMT' and drv_list[1].upper() == 'NETCDF':
     69          drv_list = ['NETCDF', 'GMT']
     70
     71      return drv_list
     72
     73
     74 -def GetOutputDriverFor(filename):
     75      drv_list = GetOutputDriversFor(filename)
     76      ext = GetExtension(filename)
     77      if not drv_list:
     78          if not ext:
     79              return 'GTiff'
     80          else:
     81              raise Exception("Cannot guess driver for %s" % filename)
     82      elif len(drv_list) > 1:
     83          print("Several drivers matching %s extension. Using %s" % (ext if ext else '', drv_list[0]))
     84      return drv_list[0]
     85
     86
     87 -def Usage():
     88      print('Usage: gdal_pansharpen [--help-general] pan_dataset {spectral_dataset[,band=num]}+ out_dataset')
     89      print('                       [-of format] [-b band]* [-w weight]*')
     90      print('                       [-r {nearest,bilinear,cubic,cubicspline,lanczos,average}]')
     91      print('                       [-threads {ALL_CPUS|number}] [-bitdepth val] [-nodata val]')
     92      print('                       [-spat_adjust {union,intersection,none,nonewithoutwarning}]')
     93      print('                       [-verbose_vrt] [-co NAME=VALUE]* [-q]')
     94      print('')
     95      print('Create a dataset resulting from a pansharpening operation.')
     96      return -1
     97
     98
     99 -def gdal_pansharpen(argv):
    100
    101      argv = gdal.GeneralCmdLineProcessor(argv)
    102      if argv is None:
    103          return -1
    104
    105      pan_name = None
    106      last_name = None
    107      spectral_ds = []
    108      spectral_bands = []
    109      out_name = None
    110      bands = []
    111      weights = []
    112      frmt = None
    113      creation_options = []
    114      callback = gdal.TermProgress_nocb
    115      resampling = None
    116      spat_adjust = None
    117      verbose_vrt = False
    118      num_threads = None
    119      bitdepth = None
    120      nodata = None
    121
    122      i = 1
    123      argc = len(argv)
    124      while i < argc:
    125          if (argv[i] == '-of' or argv[i] == '-f') and i < len(argv) - 1:
    126              frmt = argv[i + 1]
    127              i = i + 1
    128          elif argv[i] == '-r' and i < len(argv) - 1:
    129              resampling = argv[i + 1]
    130              i = i + 1
    131          elif argv[i] == '-spat_adjust' and i < len(argv) - 1:
    132              spat_adjust = argv[i + 1]
    133              i = i + 1
    134          elif argv[i] == '-b' and i < len(argv) - 1:
    135              bands.append(int(argv[i + 1]))
    136              i = i + 1
    137          elif argv[i] == '-w' and i < len(argv) - 1:
    138              weights.append(float(argv[i + 1]))
    139              i = i + 1
    140          elif argv[i] == '-co' and i < len(argv) - 1:
    141              creation_options.append(argv[i + 1])
    142              i = i + 1
    143          elif argv[i] == '-threads' and i < len(argv) - 1:
    144              num_threads = argv[i + 1]
    145              i = i + 1
    146          elif argv[i] == '-bitdepth' and i < len(argv) - 1:
    147              bitdepth = argv[i + 1]
    148              i = i + 1
    149          elif argv[i] == '-nodata' and i < len(argv) - 1:
    150              nodata = argv[i + 1]
    151              i = i + 1
    152          elif argv[i] == '-q':
    153              callback = None
    154          elif argv[i] == '-verbose_vrt':
    155              verbose_vrt = True
    156          elif argv[i][0] == '-':
    157              sys.stderr.write('Unrecognized option : %s\n' % argv[i])
    158              return Usage()
    159          elif pan_name is None:
    160              pan_name = argv[i]
    161              pan_ds = gdal.Open(pan_name)
    162              if pan_ds is None:
    163                  return 1
    164          else:
    165              if last_name is not None:
    166                  pos = last_name.find(',band=')
    167                  if pos > 0:
    168                      spectral_name = last_name[0:pos]
    169                      ds = gdal.Open(spectral_name)
    170                      if ds is None:
    171                          return 1
    172                      band_num = int(last_name[pos + len(',band='):])
    173                      band = ds.GetRasterBand(band_num)
    174                      spectral_ds.append(ds)
    175                      spectral_bands.append(band)
    176                  else:
    177                      spectral_name = last_name
    178                      ds = gdal.Open(spectral_name)
    179                      if ds is None:
    180                          return 1
    181                      for j in range(ds.RasterCount):
    182                          spectral_ds.append(ds)
    183                          spectral_bands.append(ds.GetRasterBand(j + 1))
    184
    185              last_name = argv[i]
    186
    187          i = i + 1
    188
    189      if pan_name is None or not spectral_bands:
    190          return Usage()
    191      out_name = last_name
    192
    193      if frmt is None:
    194          frmt = GetOutputDriverFor(out_name)
    195
    196      if not bands:
    197          bands = [j + 1 for j in range(len(spectral_bands))]
    198      else:
    199          for band in bands:
    200              if band < 0 or band > len(spectral_bands):
    201                  print('Invalid band number in -b: %d' % band)
    202                  return 1
    203
    204      if weights and len(weights) != len(spectral_bands):
    205          print('There must be as many -w values specified as input spectral bands')
    206          return 1
    207
    208      vrt_xml = """<VRTDataset subClass="VRTPansharpenedDataset">\n"""
    209      if bands != [j + 1 for j in range(len(spectral_bands))]:
    210          for i, band in enumerate(bands):
    211              sband = spectral_bands[band - 1]
    212              datatype = gdal.GetDataTypeName(sband.DataType)
    213              colorname = gdal.GetColorInterpretationName(sband.GetColorInterpretation())
    214              vrt_xml += """  <VRTRasterBand dataType="%s" band="%d" subClass="VRTPansharpenedRasterBand">
    215        <ColorInterp>%s</ColorInterp>
    216    </VRTRasterBand>\n""" % (datatype, i + 1, colorname)
    217
    218      vrt_xml += """  <PansharpeningOptions>\n"""
    219
    220      if weights:
    221          vrt_xml += """      <AlgorithmOptions>\n"""
    222          vrt_xml += """        <Weights>"""
    223          for i, weight in enumerate(weights):
    224              if i > 0:
    225                  vrt_xml += ","
    226              vrt_xml += "%.16g" % weight
    227          vrt_xml += "</Weights>\n"
    228          vrt_xml += """      </AlgorithmOptions>\n"""
    229
    230      if resampling is not None:
    231          vrt_xml += '      <Resampling>%s</Resampling>\n' % resampling
    232
    233      if num_threads is not None:
    234          vrt_xml += '      <NumThreads>%s</NumThreads>\n' % num_threads
    235
    236      if bitdepth is not None:
    237          vrt_xml += '      <BitDepth>%s</BitDepth>\n' % bitdepth
    238
    239      if nodata is not None:
    240          vrt_xml += '      <NoData>%s</NoData>\n' % nodata
    241
    242      if spat_adjust is not None:
    243          vrt_xml += '      <SpatialExtentAdjustment>%s</SpatialExtentAdjustment>\n' % spat_adjust
    244
    245      pan_relative = '0'
    246      if frmt.upper() == 'VRT':
    247          if not os.path.isabs(pan_name):
    248              pan_relative = '1'
    249              pan_name = os.path.relpath(pan_name, os.path.dirname(out_name))
    250
    251      vrt_xml += """    <PanchroBand>
    252        <SourceFilename relativeToVRT="%s">%s</SourceFilename>
    253        <SourceBand>1</SourceBand>
    254      </PanchroBand>\n""" % (pan_relative, pan_name)
    255
    256      for i, sband in enumerate(spectral_bands):
    257          dstband = ''
    258          for j, band in enumerate(bands):
    259              if i + 1 == band:
    260                  dstband = ' dstBand="%d"' % (j + 1)
    261                  break
    262
    263          ms_relative = '0'
    264          ms_name = spectral_ds[i].GetDescription()
    265          if frmt.upper() == 'VRT':
    266              if not os.path.isabs(ms_name):
    267                  ms_relative = '1'
    268                  ms_name = os.path.relpath(ms_name, os.path.dirname(out_name))
    269
    270          vrt_xml += """    <SpectralBand%s>
    271        <SourceFilename relativeToVRT="%s">%s</SourceFilename>
    272        <SourceBand>%d</SourceBand>
    273      </SpectralBand>\n""" % (dstband, ms_relative, ms_name, sband.GetBand())
    274
    275      vrt_xml += """  </PansharpeningOptions>\n"""
    276      vrt_xml += """</VRTDataset>\n"""
    277
    278      if frmt.upper() == 'VRT':
    279          f = gdal.VSIFOpenL(out_name, 'wb')
    280          if f is None:
    281              print('Cannot create %s' % out_name)
    282              return 1
    283          gdal.VSIFWriteL(vrt_xml, 1, len(vrt_xml), f)
    284          gdal.VSIFCloseL(f)
    285          if verbose_vrt:
    286              vrt_ds = gdal.Open(out_name, gdal.GA_Update)
    287              vrt_ds.SetMetadata(vrt_ds.GetMetadata())
    288          else:
    289              vrt_ds = gdal.Open(out_name)
    290          if vrt_ds is None:
    291              return 1
    292
    293          return 0
    294
    295      vrt_ds = gdal.Open(vrt_xml)
    296      out_ds = gdal.GetDriverByName(frmt).CreateCopy(out_name, vrt_ds, 0, creation_options, callback=callback)
    297      if out_ds is None:
    298          return 1
    299      return 0
    300
    301
    302 -def main(argv):
    303      return gdal_pansharpen(argv)
    304
    305
    306  if __name__ == '__main__':
    307      sys.exit(main(sys.argv))
    308