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