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

# Source Code for Module osgeo.utils.gdal_sieve

      1  #!/usr/bin/env python3
      2  # ******************************************************************************
      3  #  $Id$
      4  #
      5  #  Project:  GDAL Python Interface
      6  #  Purpose:  Application for applying sieve filter to raster data.
      7  #  Author:   Frank Warmerdam, warmerdam@pobox.com
      8  #
      9  # ******************************************************************************
     10  #  Copyright (c) 2008, Frank Warmerdam
     11  #  Copyright (c) 2009-2010, 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 os.path
     33  import sys
     34
     35  from osgeo import gdal
     36
     37
     38 -def Usage():
     39      print("""
     40  gdal_sieve [-q] [-st threshold] [-4] [-8] [-o name=value]
     41             srcfile [-nomask] [-mask filename] [-of format] [dstfile]
     42  """)
     43      sys.exit(1)
     44
     45
     46 -def DoesDriverHandleExtension(drv, ext):
     47      exts = drv.GetMetadataItem(gdal.DMD_EXTENSIONS)
     48      return exts is not None and exts.lower().find(ext.lower()) >= 0
     49
     50
     51 -def GetExtension(filename):
     52      ext = os.path.splitext(filename)[1]
     53      if ext.startswith('.'):
     54          ext = ext[1:]
     55      return ext
     56
     57
     58 -def GetOutputDriversFor(filename):
     59      drv_list = []
     60      ext = GetExtension(filename)
     61      for i in range(gdal.GetDriverCount()):
     62          drv = gdal.GetDriver(i)
     63          if (drv.GetMetadataItem(gdal.DCAP_CREATE) is not None or
     64              drv.GetMetadataItem(gdal.DCAP_CREATECOPY) is not None) and \
     65             drv.GetMetadataItem(gdal.DCAP_RASTER) is not None:
     66              if ext and DoesDriverHandleExtension(drv, ext):
     67                  drv_list.append(drv.ShortName)
     68              else:
     69                  prefix = drv.GetMetadataItem(gdal.DMD_CONNECTION_PREFIX)
     70                  if prefix is not None and filename.lower().startswith(prefix.lower()):
     71                      drv_list.append(drv.ShortName)
     72
     73      # GMT is registered before netCDF for opening reasons, but we want
     74      # netCDF to be used by default for output.
     75      if ext.lower() == 'nc' and not drv_list and \
     76         drv_list[0].upper() == 'GMT' and drv_list[1].upper() == 'NETCDF':
     77          drv_list = ['NETCDF', 'GMT']
     78
     79      return drv_list
     80
     81
     82 -def GetOutputDriverFor(filename):
     83      drv_list = GetOutputDriversFor(filename)
     84      ext = GetExtension(filename)
     85      if not drv_list:
     86          if not ext:
     87              return 'GTiff'
     88          else:
     89              raise Exception("Cannot guess driver for %s" % filename)
     90      elif len(drv_list) > 1:
     91          print("Several drivers matching %s extension. Using %s" % (ext if ext else '', drv_list[0]))
     92      return drv_list[0]
     93
     94
     95 -def main(argv):
     96      threshold = 2
     97      connectedness = 4
     98      quiet_flag = 0
     99      src_filename = None
    100
    101      dst_filename = None
    102      frmt = None
    103
    104      mask = 'default'
    105
    106      gdal.AllRegister()
    107      argv = gdal.GeneralCmdLineProcessor(argv)
    108      if argv is None:
    109          sys.exit(0)
    110
    111      # Parse command line arguments.
    112      i = 1
    113      while i < len(argv):
    114          arg = argv[i]
    115
    116          if arg == '-of' or arg == '-f':
    117              i = i + 1
    118              frmt = argv[i]
    119
    120          elif arg == '-4':
    121              connectedness = 4
    122
    123          elif arg == '-8':
    124              connectedness = 8
    125
    126          elif arg == '-q' or arg == '-quiet':
    127              quiet_flag = 1
    128
    129          elif arg == '-st':
    130              i = i + 1
    131              threshold = int(argv[i])
    132
    133          elif arg == '-nomask':
    134              mask = 'none'
    135
    136          elif arg == '-mask':
    137              i = i + 1
    138              mask = argv[i]
    139
    140          elif arg == '-mask':
    141              i = i + 1
    142              mask = argv[i]
    143
    144          elif arg[:2] == '-h':
    145              Usage()
    146
    147          elif src_filename is None:
    148              src_filename = argv[i]
    149
    150          elif dst_filename is None:
    151              dst_filename = argv[i]
    152
    153          else:
    154              Usage()
    155
    156          i = i + 1
    157
    158      if src_filename is None:
    159          Usage()
    160
    161      # =============================================================================
    162      #   Verify we have next gen bindings with the sievefilter method.
    163      # =============================================================================
    164      try:
    165          gdal.SieveFilter
    166      except AttributeError:
    167          print('')
    168          print('gdal.SieveFilter() not available.  You are likely using "old gen"')
    169          print('bindings or an older version of the next gen bindings.')
    170          print('')
    171          sys.exit(1)
    172
    173      # =============================================================================
    174      # Open source file
    175      # =============================================================================
    176
    177      if dst_filename is None:
    178          src_ds = gdal.Open(src_filename, gdal.GA_Update)
    179      else:
    180          src_ds = gdal.Open(src_filename, gdal.GA_ReadOnly)
    181
    182      if src_ds is None:
    183          print('Unable to open %s ' % src_filename)
    184          sys.exit(1)
    185
    186      srcband = src_ds.GetRasterBand(1)
    187
    188      if mask == 'default':
    189          maskband = srcband.GetMaskBand()
    190      elif mask == 'none':
    191          maskband = None
    192      else:
    193          mask_ds = gdal.Open(mask)
    194          maskband = mask_ds.GetRasterBand(1)
    195
    196      # =============================================================================
    197      #       Create output file if one is specified.
    198      # =============================================================================
    199
    200      if dst_filename is not None:
    201          if frmt is None:
    202              frmt = GetOutputDriverFor(dst_filename)
    203
    204          drv = gdal.GetDriverByName(frmt)
    205          dst_ds = drv.Create(dst_filename, src_ds.RasterXSize, src_ds.RasterYSize, 1,
    206                              srcband.DataType)
    207          wkt = src_ds.GetProjection()
    208          if wkt != '':
    209              dst_ds.SetProjection(wkt)
    210          gt = src_ds.GetGeoTransform(can_return_null=True)
    211          if gt is not None:
    212              dst_ds.SetGeoTransform(gt)
    213
    214          dstband = dst_ds.GetRasterBand(1)
    215      else:
    216          dstband = srcband
    217
    218      # =============================================================================
    219      # Invoke algorithm.
    220      # =============================================================================
    221
    222      if quiet_flag:
    223          prog_func = None
    224      else:
    225          prog_func = gdal.TermProgress_nocb
    226
    227      result = gdal.SieveFilter(srcband, maskband, dstband,
    228                                threshold, connectedness,
    229                                callback=prog_func)
    230
    231      src_ds = None
    232      dst_ds = None
    233      mask_ds = None
    234
    235      return result
    236
    237  if __name__ == '__main__':
    238      sys.exit(main(sys.argv))
    239