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