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

Source Code for Module osgeo.utils.gdal_pansharpen

  1  #!/usr/bin/env python3 
  2  # -*- coding: utf-8 -*- 
  3  ############################################################################### 
  4  # $Id$ 
  5  # 
  6  #  Project:  GDAL scripts 
  7  #  Purpose:  Perform a pansharpening operation 
  8  #  Author:   Even Rouault <even.rouault at spatialys.com> 
  9  # 
 10  ############################################################################### 
 11  #  Copyright (c) 2015, Even Rouault <even.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 
 33  import os.path 
 34  import sys 
 35  from osgeo import gdal 
 36   
 37   
38 -def DoesDriverHandleExtension(drv, ext):
39 exts = drv.GetMetadataItem(gdal.DMD_EXTENSIONS) 40 return exts is not None and exts.lower().find(ext.lower()) >= 0
41 42
43 -def GetExtension(filename):
44 ext = os.path.splitext(filename)[1] 45 if ext.startswith('.'): 46 ext = ext[1:] 47 return ext
48 49
50 -def GetOutputDriversFor(filename):
51 drv_list = [] 52 ext = GetExtension(filename) 53 for i in range(gdal.GetDriverCount()): 54 drv = gdal.GetDriver(i) 55 if (drv.GetMetadataItem(gdal.DCAP_CREATE) is not None or 56 drv.GetMetadataItem(gdal.DCAP_CREATECOPY) is not None) and \ 57 drv.GetMetadataItem(gdal.DCAP_RASTER) is not None: 58 if ext and DoesDriverHandleExtension(drv, ext): 59 drv_list.append(drv.ShortName) 60 else: 61 prefix = drv.GetMetadataItem(gdal.DMD_CONNECTION_PREFIX) 62 if prefix is not None and filename.lower().startswith(prefix.lower()): 63 drv_list.append(drv.ShortName) 64 65 # GMT is registered before netCDF for opening reasons, but we want 66 # netCDF to be used by default for output. 67 if ext.lower() == 'nc' and not drv_list and \ 68 drv_list[0].upper() == 'GMT' and drv_list[1].upper() == 'NETCDF': 69 drv_list = ['NETCDF', 'GMT'] 70 71 return drv_list
72 73
74 -def GetOutputDriverFor(filename):
75 drv_list = GetOutputDriversFor(filename) 76 ext = GetExtension(filename) 77 if not drv_list: 78 if not ext: 79 return 'GTiff' 80 else: 81 raise Exception("Cannot guess driver for %s" % filename) 82 elif len(drv_list) > 1: 83 print("Several drivers matching %s extension. Using %s" % (ext if ext else '', drv_list[0])) 84 return drv_list[0]
85 86
87 -def Usage():
88 print('Usage: gdal_pansharpen [--help-general] pan_dataset {spectral_dataset[,band=num]}+ out_dataset') 89 print(' [-of format] [-b band]* [-w weight]*') 90 print(' [-r {nearest,bilinear,cubic,cubicspline,lanczos,average}]') 91 print(' [-threads {ALL_CPUS|number}] [-bitdepth val] [-nodata val]') 92 print(' [-spat_adjust {union,intersection,none,nonewithoutwarning}]') 93 print(' [-verbose_vrt] [-co NAME=VALUE]* [-q]') 94 print('') 95 print('Create a dataset resulting from a pansharpening operation.') 96 return -1
97 98
99 -def gdal_pansharpen(argv):
100 101 argv = gdal.GeneralCmdLineProcessor(argv) 102 if argv is None: 103 return -1 104 105 pan_name = None 106 last_name = None 107 spectral_ds = [] 108 spectral_bands = [] 109 out_name = None 110 bands = [] 111 weights = [] 112 frmt = None 113 creation_options = [] 114 callback = gdal.TermProgress_nocb 115 resampling = None 116 spat_adjust = None 117 verbose_vrt = False 118 num_threads = None 119 bitdepth = None 120 nodata = None 121 122 i = 1 123 argc = len(argv) 124 while i < argc: 125 if (argv[i] == '-of' or argv[i] == '-f') and i < len(argv) - 1: 126 frmt = argv[i + 1] 127 i = i + 1 128 elif argv[i] == '-r' and i < len(argv) - 1: 129 resampling = argv[i + 1] 130 i = i + 1 131 elif argv[i] == '-spat_adjust' and i < len(argv) - 1: 132 spat_adjust = argv[i + 1] 133 i = i + 1 134 elif argv[i] == '-b' and i < len(argv) - 1: 135 bands.append(int(argv[i + 1])) 136 i = i + 1 137 elif argv[i] == '-w' and i < len(argv) - 1: 138 weights.append(float(argv[i + 1])) 139 i = i + 1 140 elif argv[i] == '-co' and i < len(argv) - 1: 141 creation_options.append(argv[i + 1]) 142 i = i + 1 143 elif argv[i] == '-threads' and i < len(argv) - 1: 144 num_threads = argv[i + 1] 145 i = i + 1 146 elif argv[i] == '-bitdepth' and i < len(argv) - 1: 147 bitdepth = argv[i + 1] 148 i = i + 1 149 elif argv[i] == '-nodata' and i < len(argv) - 1: 150 nodata = argv[i + 1] 151 i = i + 1 152 elif argv[i] == '-q': 153 callback = None 154 elif argv[i] == '-verbose_vrt': 155 verbose_vrt = True 156 elif argv[i][0] == '-': 157 sys.stderr.write('Unrecognized option : %s\n' % argv[i]) 158 return Usage() 159 elif pan_name is None: 160 pan_name = argv[i] 161 pan_ds = gdal.Open(pan_name) 162 if pan_ds is None: 163 return 1 164 else: 165 if last_name is not None: 166 pos = last_name.find(',band=') 167 if pos > 0: 168 spectral_name = last_name[0:pos] 169 ds = gdal.Open(spectral_name) 170 if ds is None: 171 return 1 172 band_num = int(last_name[pos + len(',band='):]) 173 band = ds.GetRasterBand(band_num) 174 spectral_ds.append(ds) 175 spectral_bands.append(band) 176 else: 177 spectral_name = last_name 178 ds = gdal.Open(spectral_name) 179 if ds is None: 180 return 1 181 for j in range(ds.RasterCount): 182 spectral_ds.append(ds) 183 spectral_bands.append(ds.GetRasterBand(j + 1)) 184 185 last_name = argv[i] 186 187 i = i + 1 188 189 if pan_name is None or not spectral_bands: 190 return Usage() 191 out_name = last_name 192 193 if frmt is None: 194 frmt = GetOutputDriverFor(out_name) 195 196 if not bands: 197 bands = [j + 1 for j in range(len(spectral_bands))] 198 else: 199 for band in bands: 200 if band < 0 or band > len(spectral_bands): 201 print('Invalid band number in -b: %d' % band) 202 return 1 203 204 if weights and len(weights) != len(spectral_bands): 205 print('There must be as many -w values specified as input spectral bands') 206 return 1 207 208 vrt_xml = """<VRTDataset subClass="VRTPansharpenedDataset">\n""" 209 if bands != [j + 1 for j in range(len(spectral_bands))]: 210 for i, band in enumerate(bands): 211 sband = spectral_bands[band - 1] 212 datatype = gdal.GetDataTypeName(sband.DataType) 213 colorname = gdal.GetColorInterpretationName(sband.GetColorInterpretation()) 214 vrt_xml += """ <VRTRasterBand dataType="%s" band="%d" subClass="VRTPansharpenedRasterBand"> 215 <ColorInterp>%s</ColorInterp> 216 </VRTRasterBand>\n""" % (datatype, i + 1, colorname) 217 218 vrt_xml += """ <PansharpeningOptions>\n""" 219 220 if weights: 221 vrt_xml += """ <AlgorithmOptions>\n""" 222 vrt_xml += """ <Weights>""" 223 for i, weight in enumerate(weights): 224 if i > 0: 225 vrt_xml += "," 226 vrt_xml += "%.16g" % weight 227 vrt_xml += "</Weights>\n" 228 vrt_xml += """ </AlgorithmOptions>\n""" 229 230 if resampling is not None: 231 vrt_xml += ' <Resampling>%s</Resampling>\n' % resampling 232 233 if num_threads is not None: 234 vrt_xml += ' <NumThreads>%s</NumThreads>\n' % num_threads 235 236 if bitdepth is not None: 237 vrt_xml += ' <BitDepth>%s</BitDepth>\n' % bitdepth 238 239 if nodata is not None: 240 vrt_xml += ' <NoData>%s</NoData>\n' % nodata 241 242 if spat_adjust is not None: 243 vrt_xml += ' <SpatialExtentAdjustment>%s</SpatialExtentAdjustment>\n' % spat_adjust 244 245 pan_relative = '0' 246 if frmt.upper() == 'VRT': 247 if not os.path.isabs(pan_name): 248 pan_relative = '1' 249 pan_name = os.path.relpath(pan_name, os.path.dirname(out_name)) 250 251 vrt_xml += """ <PanchroBand> 252 <SourceFilename relativeToVRT="%s">%s</SourceFilename> 253 <SourceBand>1</SourceBand> 254 </PanchroBand>\n""" % (pan_relative, pan_name) 255 256 for i, sband in enumerate(spectral_bands): 257 dstband = '' 258 for j, band in enumerate(bands): 259 if i + 1 == band: 260 dstband = ' dstBand="%d"' % (j + 1) 261 break 262 263 ms_relative = '0' 264 ms_name = spectral_ds[i].GetDescription() 265 if frmt.upper() == 'VRT': 266 if not os.path.isabs(ms_name): 267 ms_relative = '1' 268 ms_name = os.path.relpath(ms_name, os.path.dirname(out_name)) 269 270 vrt_xml += """ <SpectralBand%s> 271 <SourceFilename relativeToVRT="%s">%s</SourceFilename> 272 <SourceBand>%d</SourceBand> 273 </SpectralBand>\n""" % (dstband, ms_relative, ms_name, sband.GetBand()) 274 275 vrt_xml += """ </PansharpeningOptions>\n""" 276 vrt_xml += """</VRTDataset>\n""" 277 278 if frmt.upper() == 'VRT': 279 f = gdal.VSIFOpenL(out_name, 'wb') 280 if f is None: 281 print('Cannot create %s' % out_name) 282 return 1 283 gdal.VSIFWriteL(vrt_xml, 1, len(vrt_xml), f) 284 gdal.VSIFCloseL(f) 285 if verbose_vrt: 286 vrt_ds = gdal.Open(out_name, gdal.GA_Update) 287 vrt_ds.SetMetadata(vrt_ds.GetMetadata()) 288 else: 289 vrt_ds = gdal.Open(out_name) 290 if vrt_ds is None: 291 return 1 292 293 return 0 294 295 vrt_ds = gdal.Open(vrt_xml) 296 out_ds = gdal.GetDriverByName(frmt).CreateCopy(out_name, vrt_ds, 0, creation_options, callback=callback) 297 if out_ds is None: 298 return 1 299 return 0
300 301
302 -def main(argv):
303 return gdal_pansharpen(argv)
304 305 306 if __name__ == '__main__': 307 sys.exit(main(sys.argv)) 308