osgeo.utils.gdal

1 #!/usr/bin/env python3 2 ############################################################################### 3 # $Id$ 4 # 5 # Project: InSAR Peppers 6 # Purpose: Module to extract data from many rasters into one output. 7 # Author: Frank Warmerdam, warmerdam@pobox.com 8 # 9 ############################################################################### 10 # Copyright (c) 2000, Atlantis Scientific Inc. (www.atlsci.com) 11 # Copyright (c) 2009-2011, Even Rouault <even dot rouault at spatialys.com> 12 # 13 # This library is free software; you can redistribute it and/or 14 # modify it under the terms of the GNU Library General Public 15 # License as published by the Free Software Foundation; either 16 # version 2 of the License, or (at your option) any later version. 17 # 18 # This library is distributed in the hope that it will be useful, 19 # but WITHOUT ANY WARRANTY; without even the implied warranty of 20 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 21 # Library General Public License for more details. 22 # 23 # You should have received a copy of the GNU Library General Public 24 # License along with this library; if not, write to the 25 # Free Software Foundation, Inc., 59 Temple Place - Suite 330, 26 # Boston, MA 02111-1307, USA. 27 ############################################################################### 28 # changes 29Apr2011 29 # If the input image is a multi-band one, use all the channels in 30 # building the stack. 31 # anssi.pekkarinen@fao.org 32 33 import math 34 import os.path 35 import sys 36 import time 37 38 from osgeo import gdal 39 40 progress = gdal.TermProgress_nocb 41 42 __version__ = '$id$'[5:-1] 43 44

45 -def DoesDriverHandleExtension(drv, ext):

46 exts = drv.GetMetadataItem(gdal.DMD_EXTENSIONS) 47 return exts is not None and exts.lower().find(ext.lower()) >= 0

48 49

50 -def GetExtension(filename):

51 ext = os.path.splitext(filename)[1] 52 if ext.startswith('.'): 53 ext = ext[1:] 54 return ext

55 56

57 -def GetOutputDriversFor(filename):

58 drv_list = [] 59 ext = GetExtension(filename) 60 for i in range(gdal.GetDriverCount()): 61 drv = gdal.GetDriver(i) 62 if (drv.GetMetadataItem(gdal.DCAP_CREATE) is not None or 63 drv.GetMetadataItem(gdal.DCAP_CREATECOPY) is not None) and \ 64 drv.GetMetadataItem(gdal.DCAP_RASTER) is not None: 65 if ext and DoesDriverHandleExtension(drv, ext): 66 drv_list.append(drv.ShortName) 67 else: 68 prefix = drv.GetMetadataItem(gdal.DMD_CONNECTION_PREFIX) 69 if prefix is not None and filename.lower().startswith(prefix.lower()): 70 drv_list.append(drv.ShortName) 71 72 # GMT is registered before netCDF for opening reasons, but we want 73 # netCDF to be used by default for output. 74 if ext.lower() == 'nc' and not drv_list and \ 75 drv_list[0].upper() == 'GMT' and drv_list[1].upper() == 'NETCDF': 76 drv_list = ['NETCDF', 'GMT'] 77 78 return drv_list

79 80

81 -def GetOutputDriverFor(filename):

82 drv_list = GetOutputDriversFor(filename) 83 ext = GetExtension(filename) 84 if not drv_list: 85 if not ext: 86 return 'GTiff' 87 else: 88 raise Exception("Cannot guess driver for %s" % filename) 89 elif len(drv_list) > 1: 90 print("Several drivers matching %s extension. Using %s" % (ext if ext else '', drv_list[0])) 91 return drv_list[0]

92 93 94 # =============================================================================

95 -def raster_copy(s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n, 96 t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n, 97 nodata=None, verbose=0):

98 99 if verbose != 0: 100 print('Copy %d,%d,%d,%d to %d,%d,%d,%d.' 101 % (s_xoff, s_yoff, s_xsize, s_ysize, 102 t_xoff, t_yoff, t_xsize, t_ysize)) 103 104 if nodata is not None: 105 return raster_copy_with_nodata( 106 s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n, 107 t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n, 108 nodata) 109 110 s_band = s_fh.GetRasterBand(s_band_n) 111 m_band = None 112 # Works only in binary mode and doesn't take into account 113 # intermediate transparency values for compositing. 114 if s_band.GetMaskFlags() != gdal.GMF_ALL_VALID: 115 m_band = s_band.GetMaskBand() 116 elif s_band.GetColorInterpretation() == gdal.GCI_AlphaBand: 117 m_band = s_band 118 if m_band is not None: 119 return raster_copy_with_mask( 120 s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n, 121 t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n, 122 m_band) 123 124 s_band = s_fh.GetRasterBand(s_band_n) 125 t_band = t_fh.GetRasterBand(t_band_n) 126 127 data = s_band.ReadRaster(s_xoff, s_yoff, s_xsize, s_ysize, 128 t_xsize, t_ysize, t_band.DataType) 129 t_band.WriteRaster(t_xoff, t_yoff, t_xsize, t_ysize, 130 data, t_xsize, t_ysize, t_band.DataType) 131 132 return 0

133 134 # ============================================================================= 135 136

137 -def raster_copy_with_nodata(s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n, 138 t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n, 139 nodata):

140 try: 141 import numpy as Numeric 142 except ImportError: 143 import Numeric 144 145 s_band = s_fh.GetRasterBand(s_band_n) 146 t_band = t_fh.GetRasterBand(t_band_n) 147 148 data_src = s_band.ReadAsArray(s_xoff, s_yoff, s_xsize, s_ysize, 149 t_xsize, t_ysize) 150 data_dst = t_band.ReadAsArray(t_xoff, t_yoff, t_xsize, t_ysize) 151 152 if not Numeric.isnan(nodata): 153 nodata_test = Numeric.equal(data_src, nodata) 154 else: 155 nodata_test = Numeric.isnan(data_src) 156 157 to_write = Numeric.choose(nodata_test, (data_src, data_dst)) 158 159 t_band.WriteArray(to_write, t_xoff, t_yoff) 160 161 return 0

162 163 # ============================================================================= 164 165

166 -def raster_copy_with_mask(s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n, 167 t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n, 168 m_band):

169 try: 170 import numpy as Numeric 171 except ImportError: 172 import Numeric 173 174 s_band = s_fh.GetRasterBand(s_band_n) 175 t_band = t_fh.GetRasterBand(t_band_n) 176 177 data_src = s_band.ReadAsArray(s_xoff, s_yoff, s_xsize, s_ysize, 178 t_xsize, t_ysize) 179 data_mask = m_band.ReadAsArray(s_xoff, s_yoff, s_xsize, s_ysize, 180 t_xsize, t_ysize) 181 data_dst = t_band.ReadAsArray(t_xoff, t_yoff, t_xsize, t_ysize) 182 183 mask_test = Numeric.equal(data_mask, 0) 184 to_write = Numeric.choose(mask_test, (data_src, data_dst)) 185 186 t_band.WriteArray(to_write, t_xoff, t_yoff) 187 188 return 0

189 190 # ============================================================================= 191 192

193 -def names_to_fileinfos(names):

194 """ 195 Translate a list of GDAL filenames, into file_info objects. 196 197 names -- list of valid GDAL dataset names. 198 199 Returns a list of file_info objects. There may be less file_info objects 200 than names if some of the names could not be opened as GDAL files. 201 """ 202 203 file_infos = [] 204 for name in names: 205 fi = file_info() 206 if fi.init_from_name(name) == 1: 207 file_infos.append(fi) 208 209 return file_infos

210 211 # ***************************************************************************** 212 213

214 -class file_info(object):

215 """A class holding information about a GDAL file.""" 216

217 - def __init__(self):

218 self.band_type = None 219 self.bands = None 220 self.ct = None 221 self.filename = None 222 self.geotransform = None 223 self.lrx = None 224 self.lry = None 225 self.projection = None 226 self.ulx = None 227 self.uly = None 228 self.xsize = None 229 self.ysize = None

230

231 - def init_from_name(self, filename):

232 """ 233 Initialize file_info from filename 234 235 filename -- Name of file to read. 236 237 Returns 1 on success or 0 if the file can't be opened. 238 """ 239 fh = gdal.Open(filename) 240 if fh is None: 241 return 0 242 243 self.filename = filename 244 self.bands = fh.RasterCount 245 self.xsize = fh.RasterXSize 246 self.ysize = fh.RasterYSize 247 self.band_type = fh.GetRasterBand(1).DataType 248 self.projection = fh.GetProjection() 249 self.geotransform = fh.GetGeoTransform() 250 self.ulx = self.geotransform[0] 251 self.uly = self.geotransform[3] 252 self.lrx = self.ulx + self.geotransform[1] * self.xsize 253 self.lry = self.uly + self.geotransform[5] * self.ysize 254 255 ct = fh.GetRasterBand(1).GetRasterColorTable() 256 if ct is not None: 257 self.ct = ct.Clone() 258 else: 259 self.ct = None 260 261 return 1

262

263 - def report(self):

264 print('Filename: ' + self.filename) 265 print('File Size: %dx%dx%d' 266 % (self.xsize, self.ysize, self.bands)) 267 print('Pixel Size: %f x %f' 268 % (self.geotransform[1], self.geotransform[5])) 269 print('UL:(%f,%f) LR:(%f,%f)' 270 % (self.ulx, self.uly, self.lrx, self.lry))

271

272 - def copy_into(self, t_fh, s_band=1, t_band=1, nodata_arg=None, verbose=0):

273 """ 274 Copy this files image into target file. 275 276 This method will compute the overlap area of the file_info objects 277 file, and the target gdal.Dataset object, and copy the image data 278 for the common window area. It is assumed that the files are in 279 a compatible projection ... no checking or warping is done. However, 280 if the destination file is a different resolution, or different 281 image pixel type, the appropriate resampling and conversions will 282 be done (using normal GDAL promotion/demotion rules). 283 284 t_fh -- gdal.Dataset object for the file into which some or all 285 of this file may be copied. 286 287 Returns 1 on success (or if nothing needs to be copied), and zero one 288 failure. 289 """ 290 t_geotransform = t_fh.GetGeoTransform() 291 t_ulx = t_geotransform[0] 292 t_uly = t_geotransform[3] 293 t_lrx = t_geotransform[0] + t_fh.RasterXSize * t_geotransform[1] 294 t_lry = t_geotransform[3] + t_fh.RasterYSize * t_geotransform[5] 295 296 # figure out intersection region 297 tgw_ulx = max(t_ulx, self.ulx) 298 tgw_lrx = min(t_lrx, self.lrx) 299 if t_geotransform[5] < 0: 300 tgw_uly = min(t_uly, self.uly) 301 tgw_lry = max(t_lry, self.lry) 302 else: 303 tgw_uly = max(t_uly, self.uly) 304 tgw_lry = min(t_lry, self.lry) 305 306 # do they even intersect? 307 if tgw_ulx >= tgw_lrx: 308 return 1 309 if t_geotransform[5] < 0 and tgw_uly <= tgw_lry: 310 return 1 311 if t_geotransform[5] > 0 and tgw_uly >= tgw_lry: 312 return 1 313 314 # compute target window in pixel coordinates. 315 tw_xoff = int((tgw_ulx - t_geotransform[0]) / t_geotransform[1] + 0.1) 316 tw_yoff = int((tgw_uly - t_geotransform[3]) / t_geotransform[5] + 0.1) 317 tw_xsize = int((tgw_lrx - t_geotransform[0]) / t_geotransform[1] + 0.5) \ 318 - tw_xoff 319 tw_ysize = int((tgw_lry - t_geotransform[3]) / t_geotransform[5] + 0.5) \ 320 - tw_yoff 321 322 if tw_xsize < 1 or tw_ysize < 1: 323 return 1 324 325 # Compute source window in pixel coordinates. 326 sw_xoff = int((tgw_ulx - self.geotransform[0]) / self.geotransform[1] + 0.1) 327 sw_yoff = int((tgw_uly - self.geotransform[3]) / self.geotransform[5] + 0.1) 328 sw_xsize = int((tgw_lrx - self.geotransform[0]) / 329 self.geotransform[1] + 0.5) - sw_xoff 330 sw_ysize = int((tgw_lry - self.geotransform[3]) / 331 self.geotransform[5] + 0.5) - sw_yoff 332 333 if sw_xsize < 1 or sw_ysize < 1: 334 return 1 335 336 # Open the source file, and copy the selected region. 337 s_fh = gdal.Open(self.filename) 338 339 return raster_copy(s_fh, sw_xoff, sw_yoff, sw_xsize, sw_ysize, s_band, 340 t_fh, tw_xoff, tw_yoff, tw_xsize, tw_ysize, t_band, 341 nodata_arg, verbose)

342 343 344 # =============================================================================

345 -def Usage():

346 print('Usage: gdal_merge.py [-o out_filename] [-of out_format] [-co NAME=VALUE]*') 347 print(' [-ps pixelsize_x pixelsize_y] [-tap] [-separate] [-q] [-v] [-pct]') 348 print(' [-ul_lr ulx uly lrx lry] [-init "value [value...]"]') 349 print(' [-n nodata_value] [-a_nodata output_nodata_value]') 350 print(' [-ot datatype] [-createonly] input_files') 351 print(' [--help-general]') 352 print('')

353 354

355 -def main(argv=None):

356 verbose = 0 357 quiet = 0 358 names = [] 359 frmt = None 360 out_file = 'out.tif' 361 362 ulx = None 363 psize_x = None 364 separate = 0 365 copy_pct = 0 366 nodata = None 367 a_nodata = None 368 create_options = [] 369 pre_init = [] 370 band_type = None 371 createonly = 0 372 bTargetAlignedPixels = False 373 start_time = time.time() 374 375 gdal.AllRegister() 376 if argv is None: 377 argv = argv 378 argv = gdal.GeneralCmdLineProcessor(argv) 379 if argv is None: 380 sys.exit(0) 381 382 # Parse command line arguments. 383 i = 1 384 while i < len(argv): 385 arg = argv[i] 386 387 if arg == '-o': 388 i = i + 1 389 out_file = argv[i] 390 391 elif arg == '-v': 392 verbose = 1 393 394 elif arg == '-q' or arg == '-quiet': 395 quiet = 1 396 397 elif arg == '-createonly': 398 createonly = 1 399 400 elif arg == '-separate': 401 separate = 1 402 403 elif arg == '-seperate': 404 separate = 1 405 406 elif arg == '-pct': 407 copy_pct = 1 408 409 elif arg == '-ot': 410 i = i + 1 411 band_type = gdal.GetDataTypeByName(argv[i]) 412 if band_type == gdal.GDT_Unknown: 413 print('Unknown GDAL data type: %s' % argv[i]) 414 sys.exit(1) 415 416 elif arg == '-init': 417 i = i + 1 418 str_pre_init = argv[i].split() 419 for x in str_pre_init: 420 pre_init.append(float(x)) 421 422 elif arg == '-n': 423 i = i + 1 424 nodata = float(argv[i]) 425 426 elif arg == '-a_nodata': 427 i = i + 1 428 a_nodata = float(argv[i]) 429 430 elif arg == '-f' or arg == '-of': 431 i = i + 1 432 frmt = argv[i] 433 434 elif arg == '-co': 435 i = i + 1 436 create_options.append(argv[i]) 437 438 elif arg == '-ps': 439 psize_x = float(argv[i + 1]) 440 psize_y = -1 * abs(float(argv[i + 2])) 441 i = i + 2 442 443 elif arg == '-tap': 444 bTargetAlignedPixels = True 445 446 elif arg == '-ul_lr': 447 ulx = float(argv[i + 1]) 448 uly = float(argv[i + 2]) 449 lrx = float(argv[i + 3]) 450 lry = float(argv[i + 4]) 451 i = i + 4 452 453 elif arg[:1] == '-': 454 print('Unrecognized command option: %s' % arg) 455 Usage() 456 sys.exit(1) 457 458 else: 459 names.append(arg) 460 461 i = i + 1 462 463 if not names: 464 print('No input files selected.') 465 Usage() 466 sys.exit(1) 467 468 if frmt is None: 469 frmt = GetOutputDriverFor(out_file) 470 471 Driver = gdal.GetDriverByName(frmt) 472 if Driver is None: 473 print('Format driver %s not found, pick a supported driver.' % frmt) 474 sys.exit(1) 475 476 DriverMD = Driver.GetMetadata() 477 if 'DCAP_CREATE' not in DriverMD: 478 print('Format driver %s does not support creation and piecewise writing.\nPlease select a format that does, such as GTiff (the default) or HFA (Erdas Imagine).' % frmt) 479 sys.exit(1) 480 481 # Collect information on all the source files. 482 file_infos = names_to_fileinfos(names) 483 484 if ulx is None: 485 ulx = file_infos[0].ulx 486 uly = file_infos[0].uly 487 lrx = file_infos[0].lrx 488 lry = file_infos[0].lry 489 490 for fi in file_infos: 491 ulx = min(ulx, fi.ulx) 492 uly = max(uly, fi.uly) 493 lrx = max(lrx, fi.lrx) 494 lry = min(lry, fi.lry) 495 496 if psize_x is None: 497 psize_x = file_infos[0].geotransform[1] 498 psize_y = file_infos[0].geotransform[5] 499 500 if band_type is None: 501 band_type = file_infos[0].band_type 502 503 # Try opening as an existing file. 504 gdal.PushErrorHandler('CPLQuietErrorHandler') 505 t_fh = gdal.Open(out_file, gdal.GA_Update) 506 gdal.PopErrorHandler() 507 508 # Create output file if it does not already exist. 509 if t_fh is None: 510 511 if bTargetAlignedPixels: 512 ulx = math.floor(ulx / psize_x) * psize_x 513 lrx = math.ceil(lrx / psize_x) * psize_x 514 lry = math.floor(lry / -psize_y) * -psize_y 515 uly = math.ceil(uly / -psize_y) * -psize_y 516 517 geotransform = [ulx, psize_x, 0, uly, 0, psize_y] 518 519 xsize = int((lrx - ulx) / geotransform[1] + 0.5) 520 ysize = int((lry - uly) / geotransform[5] + 0.5) 521 522 if separate != 0: 523 bands = 0 524 525 for fi in file_infos: 526 bands = bands + fi.bands 527 else: 528 bands = file_infos[0].bands 529 530 t_fh = Driver.Create(out_file, xsize, ysize, bands, 531 band_type, create_options) 532 if t_fh is None: 533 print('Creation failed, terminating gdal_merge.') 534 sys.exit(1) 535 536 t_fh.SetGeoTransform(geotransform) 537 t_fh.SetProjection(file_infos[0].projection) 538 539 if copy_pct: 540 t_fh.GetRasterBand(1).SetRasterColorTable(file_infos[0].ct) 541 else: 542 if separate != 0: 543 bands = 0 544 for fi in file_infos: 545 bands = bands + fi.bands 546 if t_fh.RasterCount < bands: 547 print('Existing output file has less bands than the input files. You should delete it before. Terminating gdal_merge.') 548 sys.exit(1) 549 else: 550 bands = min(file_infos[0].bands, t_fh.RasterCount) 551 552 # Do we need to set nodata value ? 553 if a_nodata is not None: 554 for i in range(t_fh.RasterCount): 555 t_fh.GetRasterBand(i + 1).SetNoDataValue(a_nodata) 556 557 # Do we need to pre-initialize the whole mosaic file to some value? 558 if pre_init is not None: 559 if t_fh.RasterCount <= len(pre_init): 560 for i in range(t_fh.RasterCount): 561 t_fh.GetRasterBand(i + 1).Fill(pre_init[i]) 562 elif len(pre_init) == 1: 563 for i in range(t_fh.RasterCount): 564 t_fh.GetRasterBand(i + 1).Fill(pre_init[0]) 565 566 # Copy data from source files into output file. 567 t_band = 1 568 569 if quiet == 0 and verbose == 0: 570 progress(0.0) 571 fi_processed = 0 572 573 for fi in file_infos: 574 if createonly != 0: 575 continue 576 577 if verbose != 0: 578 print("") 579 print("Processing file %5d of %5d, %6.3f%% completed in %d minutes." 580 % (fi_processed + 1, len(file_infos), 581 fi_processed * 100.0 / len(file_infos), 582 int(round((time.time() - start_time) / 60.0)))) 583 fi.report() 584 585 if separate == 0: 586 for band in range(1, bands + 1): 587 fi.copy_into(t_fh, band, band, nodata, verbose) 588 else: 589 for band in range(1, fi.bands + 1): 590 fi.copy_into(t_fh, band, t_band, nodata, verbose) 591 t_band = t_band + 1 592 593 fi_processed = fi_processed + 1 594 if quiet == 0 and verbose == 0: 595 progress(fi_processed / float(len(file_infos))) 596 597 # Force file to be closed. 598 t_fh = None

599 600 601 if __name__ == '__main__': 602 sys.exit(main(sys.argv)) 603

Source Code for Module osgeo.utils.gdal_merge