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