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

# Source Code for Module osgeo.utils.gdal_fillnodata

      1  #!/usr/bin/env python3
      2  # ******************************************************************************
      3  #  $Id$
      4  #
      5  #  Project:  GDAL Python Interface
      6  #  Purpose:  Application for filling nodata areas in a raster by interpolation
      7  #  Author:   Frank Warmerdam, warmerdam@pobox.com
      8  #
      9  # ******************************************************************************
     10  #  Copyright (c) 2008, Frank Warmerdam
     11  #  Copyright (c) 2009-2011, Even Rouault <even dot 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 sys
     33
     34  from osgeo import gdal
     35
     36
     37 -def CopyBand(srcband, dstband):
     38      for line in range(srcband.YSize):
     39          line_data = srcband.ReadRaster(0, line, srcband.XSize, 1)
     40          dstband.WriteRaster(0, line, srcband.XSize, 1, line_data,
     41                              buf_type=srcband.DataType)
     42
     43
     44 -def Usage():
     45      print("""
     46  gdal_fillnodata [-q] [-md max_distance] [-si smooth_iterations]
     47                  [-o name=value] [-b band]
     48                  srcfile [-nomask] [-mask filename] [-of format] [-co name=value]* [dstfile]
     49  """)
     50      sys.exit(1)
     51
     52
     53 -def main(argv):
     54      max_distance = 100
     55      smoothing_iterations = 0
     56      options = []
     57      quiet_flag = 0
     58      src_filename = None
     59      src_band = 1
     60
     61      dst_filename = None
     62      frmt = 'GTiff'
     63      creation_options = []
     64
     65      mask = 'default'
     66
     67      gdal.AllRegister()
     68      argv = gdal.GeneralCmdLineProcessor(argv)
     69      if argv is None:
     70          sys.exit(0)
     71
     72      # Parse command line arguments.
     73      i = 1
     74      while i < len(argv):
     75          arg = argv[i]
     76
     77          if arg == '-of' or arg == '-f':
     78              i = i + 1
     79              frmt = argv[i]
     80
     81          elif arg == '-co':
     82              i = i + 1
     83              creation_options.append(argv[i])
     84
     85          elif arg == '-q' or arg == '-quiet':
     86              quiet_flag = 1
     87
     88          elif arg == '-si':
     89              i = i + 1
     90              smoothing_iterations = int(argv[i])
     91
     92          elif arg == '-b':
     93              i = i + 1
     94              src_band = int(argv[i])
     95
     96          elif arg == '-md':
     97              i = i + 1
     98              max_distance = float(argv[i])
     99
    100          elif arg == '-nomask':
    101              mask = 'none'
    102
    103          elif arg == '-mask':
    104              i = i + 1
    105              mask = argv[i]
    106
    107          elif arg == '-mask':
    108              i = i + 1
    109              mask = argv[i]
    110
    111          elif arg[:2] == '-h':
    112              Usage()
    113
    114          elif src_filename is None:
    115              src_filename = argv[i]
    116
    117          elif dst_filename is None:
    118              dst_filename = argv[i]
    119
    120          else:
    121              Usage()
    122
    123          i = i + 1
    124
    125      if src_filename is None:
    126          Usage()
    127
    128      # =============================================================================
    129      #   Verify we have next gen bindings with the sievefilter method.
    130      # =============================================================================
    131      try:
    132          gdal.FillNodata
    133      except AttributeError:
    134          print('')
    135          print('gdal.FillNodata() not available.  You are likely using "old gen"')
    136          print('bindings or an older version of the next gen bindings.')
    137          print('')
    138          sys.exit(1)
    139
    140      # =============================================================================
    141      # Open source file
    142      # =============================================================================
    143
    144      if dst_filename is None:
    145          src_ds = gdal.Open(src_filename, gdal.GA_Update)
    146      else:
    147          src_ds = gdal.Open(src_filename, gdal.GA_ReadOnly)
    148
    149      if src_ds is None:
    150          print('Unable to open %s' % src_filename)
    151          sys.exit(1)
    152
    153      srcband = src_ds.GetRasterBand(src_band)
    154
    155      if mask == 'default':
    156          maskband = srcband.GetMaskBand()
    157      elif mask == 'none':
    158          maskband = None
    159      else:
    160          mask_ds = gdal.Open(mask)
    161          maskband = mask_ds.GetRasterBand(1)
    162
    163      # =============================================================================
    164      #       Create output file if one is specified.
    165      # =============================================================================
    166
    167      if dst_filename is not None:
    168
    169          drv = gdal.GetDriverByName(frmt)
    170          dst_ds = drv.Create(dst_filename, src_ds.RasterXSize, src_ds.RasterYSize, 1,
    171                              srcband.DataType, creation_options)
    172          wkt = src_ds.GetProjection()
    173          if wkt != '':
    174              dst_ds.SetProjection(wkt)
    175          gt = src_ds.GetGeoTransform(can_return_null=True)
    176          if gt:
    177              dst_ds.SetGeoTransform(gt)
    178
    179          dstband = dst_ds.GetRasterBand(1)
    180          CopyBand(srcband, dstband)
    181          ndv = srcband.GetNoDataValue()
    182          if ndv is not None:
    183              dstband.SetNoDataValue(ndv)
    184
    185          color_interp = srcband.GetColorInterpretation()
    186          dstband.SetColorInterpretation(color_interp)
    187          if color_interp == gdal.GCI_PaletteIndex:
    188              color_table = srcband.GetColorTable()
    189              dstband.SetColorTable(color_table)
    190
    191      else:
    192          dstband = srcband
    193
    194      # =============================================================================
    195      # Invoke algorithm.
    196      # =============================================================================
    197
    198      if quiet_flag:
    199          prog_func = None
    200      else:
    201          prog_func = gdal.TermProgress_nocb
    202
    203      result = gdal.FillNodata(dstband, maskband,
    204                               max_distance, smoothing_iterations, options,
    205                               callback=prog_func)
    206
    207
    208      src_ds = None
    209      dst_ds = None
    210      mask_ds = None
    211
    212      return result
    213
    214
    215  if __name__ == '__main__':
    216      sys.exit(main(sys.argv))
    217