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

Source Code for Module osgeo.utils.gdal_retile

  1  #!/usr/bin/env python3 
  2  ############################################################################### 
  3  #  $Id$ 
  4  # 
  5  # Purpose:  Module for retiling (merging) tiles and building tiled pyramids 
  6  # Author:   Christian Meuller, christian.mueller@nvoe.at 
  7  # UseDirForEachRow support by Chris Giesey & Elijah Robison 
  8  # 
  9  ############################################################################### 
 10  # Copyright (c) 2007, Christian Mueller 
 11  # Copyright (c) 2009-2012, 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  from __future__ import print_function 
 32  import os 
 33  import sys 
 34  import shutil 
 35   
 36  from osgeo import gdal 
 37  from osgeo import ogr 
 38  from osgeo import osr 
 39   
 40  progress = gdal.TermProgress_nocb 
 41   
 42   
43 -class AffineTransformDecorator(object):
44 """ A class providing some useful methods for affine Transformations """ 45
46 - def __init__(self, transform):
47 self.lrx = None 48 self.lry = None 49 self.geotransform = transform 50 self.scaleX = self.geotransform[1] 51 self.scaleY = self.geotransform[5] 52 if self.scaleY > 0: 53 self.scaleY *= -1 54 self.ulx = self.geotransform[0] 55 self.uly = self.geotransform[3]
56
57 - def pointsFor(self, width, height):
58 xlist = [] 59 ylist = [] 60 w = self.scaleX * width 61 h = self.scaleY * height 62 63 xlist.append(self.ulx) 64 ylist.append(self.uly) 65 xlist.append(self.ulx + w) 66 ylist.append(self.uly) 67 xlist.append(self.ulx + w) 68 ylist.append(self.uly + h) 69 xlist.append(self.ulx) 70 ylist.append(self.uly + h) 71 return [xlist, ylist]
72 73
74 -class DataSetCache(object):
75 """ A class for caching source tiles """ 76
77 - def __init__(self):
78 self.cacheSize = 8 79 self.queue = [] 80 self.dict = {}
81
82 - def get(self, name):
83 84 if name in self.dict: 85 return self.dict[name] 86 result = gdal.Open(name) 87 if result is None: 88 print("Error opening: %s" % NameError) 89 sys.exit(1) 90 if len(self.queue) == self.cacheSize: 91 toRemove = self.queue.pop(0) 92 del self.dict[toRemove] 93 self.queue.append(name) 94 self.dict[name] = result 95 return result
96
97 - def __del__(self):
98 for dataset in self.dict.values(): 99 del dataset 100 del self.queue 101 del self.dict
102 103
104 -class tile_info(object):
105 """ A class holding info how to tile """ 106
107 - def __init__(self, xsize, ysize, tileWidth, tileHeight, overlap):
108 self.width = xsize 109 self.height = ysize 110 self.tileWidth = tileWidth 111 self.tileHeight = tileHeight 112 self.countTilesX = 1 113 if xsize > tileWidth: 114 self.countTilesX += int((xsize - tileWidth + (tileWidth - overlap) - 1) / (tileWidth - overlap)) 115 self.countTilesY = 1 116 if ysize > tileHeight: 117 self.countTilesY += int((ysize - tileHeight + (tileHeight - overlap) - 1) / (tileHeight - overlap)) 118 self.overlap = overlap
119
120 - def report(self):
121 print('tileWidth: %d' % self.tileWidth) 122 print('tileHeight: %d' % self.tileHeight) 123 print('countTilesX: %d' % self.countTilesX) 124 print('countTilesY: %d' % self.countTilesY) 125 print('overlap: %d' % self.overlap)
126 127
128 -class mosaic_info(object):
129 """A class holding information about a GDAL file or a GDAL fileset""" 130
131 - def __init__(self, filename, inputDS):
132 """ 133 Initialize mosaic_info from filename 134 135 filename -- Name of file to read. 136 137 """ 138 self.TempDriver = gdal.GetDriverByName("MEM") 139 self.filename = filename 140 self.cache = DataSetCache() 141 self.ogrTileIndexDS = inputDS 142 143 self.ogrTileIndexDS.GetLayer().ResetReading() 144 feature = self.ogrTileIndexDS.GetLayer().GetNextFeature() 145 imgLocation = feature.GetField(0) 146 147 fhInputTile = self.cache.get(imgLocation) 148 149 self.bands = fhInputTile.RasterCount 150 self.band_type = fhInputTile.GetRasterBand(1).DataType 151 self.projection = fhInputTile.GetProjection() 152 self.nodata = fhInputTile.GetRasterBand(1).GetNoDataValue() 153 154 dec = AffineTransformDecorator(fhInputTile.GetGeoTransform()) 155 self.scaleX = dec.scaleX 156 self.scaleY = dec.scaleY 157 ct = fhInputTile.GetRasterBand(1).GetRasterColorTable() 158 if ct is not None: 159 self.ct = ct.Clone() 160 else: 161 self.ct = None 162 self.ci = [0] * self.bands 163 for iband in range(self.bands): 164 self.ci[iband] = fhInputTile.GetRasterBand(iband + 1).GetRasterColorInterpretation() 165 166 extent = self.ogrTileIndexDS.GetLayer().GetExtent() 167 self.ulx = extent[0] 168 self.uly = extent[3] 169 self.lrx = extent[1] 170 self.lry = extent[2] 171 172 self.xsize = int(round((self.lrx - self.ulx) / self.scaleX)) 173 self.ysize = abs(int(round((self.uly - self.lry) / self.scaleY)))
174
175 - def __del__(self):
176 del self.cache 177 del self.ogrTileIndexDS
178
179 - def getDataSet(self, minx, miny, maxx, maxy):
180 181 self.ogrTileIndexDS.GetLayer().ResetReading() 182 self.ogrTileIndexDS.GetLayer().SetSpatialFilterRect(minx, miny, maxx, maxy) 183 features = [] 184 envelope = None 185 while True: 186 feature = self.ogrTileIndexDS.GetLayer().GetNextFeature() 187 if feature is None: 188 break 189 features.append(feature) 190 if envelope is None: 191 envelope = feature.GetGeometryRef().GetEnvelope() 192 else: 193 featureEnv = feature.GetGeometryRef().GetEnvelope() 194 envelope = (min(featureEnv[0], envelope[0]), max(featureEnv[1], envelope[1]), 195 min(featureEnv[2], envelope[2]), max(featureEnv[3], envelope[3])) 196 197 if envelope is None: 198 return None 199 200 # enlarge to query rect if necessary 201 envelope = (min(minx, envelope[0]), max(maxx, envelope[1]), 202 min(miny, envelope[2]), max(maxy, envelope[3])) 203 204 self.ogrTileIndexDS.GetLayer().SetSpatialFilter(None) 205 206 # merge tiles 207 208 resultSizeX = int((maxx - minx) / self.scaleX + 0.5) 209 resultSizeY = int((miny - maxy) / self.scaleY + 0.5) 210 211 resultDS = self.TempDriver.Create("TEMP", resultSizeX, resultSizeY, self.bands, self.band_type, []) 212 resultDS.SetGeoTransform([minx, self.scaleX, 0, maxy, 0, self.scaleY]) 213 214 for bandNr in range(1, self.bands + 1): 215 t_band = resultDS.GetRasterBand(bandNr) 216 if self.nodata is not None: 217 t_band.Fill(self.nodata) 218 t_band.SetNoDataValue(self.nodata) 219 220 for feature in features: 221 featureName = feature.GetField(0) 222 sourceDS = self.cache.get(featureName) 223 dec = AffineTransformDecorator(sourceDS.GetGeoTransform()) 224 225 dec.lrx = dec.ulx + sourceDS.RasterXSize * dec.scaleX 226 dec.lry = dec.uly + sourceDS.RasterYSize * dec.scaleY 227 228 # Find the intersection region 229 tgw_ulx = max(dec.ulx, minx) 230 tgw_lrx = min(dec.lrx, maxx) 231 if self.scaleY < 0: 232 tgw_uly = min(dec.uly, maxy) 233 tgw_lry = max(dec.lry, miny) 234 else: 235 tgw_uly = max(dec.uly, maxy) 236 tgw_lry = min(dec.lry, miny) 237 238 # Compute source window in pixel coordinates. 239 sw_xoff = int((tgw_ulx - dec.ulx) / dec.scaleX + 0.5) 240 sw_yoff = int((tgw_uly - dec.uly) / dec.scaleY + 0.5) 241 sw_xsize = min(sourceDS.RasterXSize, int((tgw_lrx - dec.ulx) / dec.scaleX + 0.5)) - sw_xoff 242 sw_ysize = min(sourceDS.RasterYSize, int((tgw_lry - dec.uly) / dec.scaleY + 0.5)) - sw_yoff 243 if sw_xsize <= 0 or sw_ysize <= 0: 244 continue 245 246 # Compute target window in pixel coordinates 247 tw_xoff = int((tgw_ulx - minx) / self.scaleX + 0.5) 248 tw_yoff = int((tgw_uly - maxy) / self.scaleY + 0.5) 249 tw_xsize = min(resultDS.RasterXSize, int((tgw_lrx - minx) / self.scaleX + 0.5)) - tw_xoff 250 tw_ysize = min(resultDS.RasterYSize, int((tgw_lry - maxy) / self.scaleY + 0.5)) - tw_yoff 251 if tw_xsize <= 0 or tw_ysize <= 0: 252 continue 253 254 assert tw_xoff >= 0 255 assert tw_yoff >= 0 256 assert sw_xoff >= 0 257 assert sw_yoff >= 0 258 259 for bandNr in range(1, self.bands + 1): 260 s_band = sourceDS.GetRasterBand(bandNr) 261 t_band = resultDS.GetRasterBand(bandNr) 262 if self.ct is not None: 263 t_band.SetRasterColorTable(self.ct) 264 t_band.SetRasterColorInterpretation(self.ci[bandNr - 1]) 265 266 data = s_band.ReadRaster(sw_xoff, sw_yoff, sw_xsize, sw_ysize, tw_xsize, tw_ysize, self.band_type) 267 if data is None: 268 print(gdal.GetLastErrorMsg()) 269 270 t_band.WriteRaster(tw_xoff, tw_yoff, tw_xsize, tw_ysize, data) 271 272 return resultDS
273
274 - def closeDataSet(self, memDS):
275 del memDS
276 # self.TempDriver.Delete("TEMP") 277
278 - def report(self):
279 print('Filename: ' + self.filename) 280 print('File Size: %dx%dx%d' 281 % (self.xsize, self.ysize, self.bands)) 282 print('Pixel Size: %f x %f' 283 % (self.scaleX, self.scaleY)) 284 print('UL:(%f,%f) LR:(%f,%f)' 285 % (self.ulx, self.uly, self.lrx, self.lry))
286 287
288 -def getTileIndexFromFiles(g):
289 if g.Verbose: 290 print("Building internal Index for %d tile(s) ..." % len(g.Names), end=" ") 291 292 ogrTileIndexDS = createTileIndex(g.Verbose, "TileIndex", g.TileIndexFieldName, None, g.TileIndexDriverTyp) 293 for inputTile in g.Names: 294 295 fhInputTile = gdal.Open(inputTile) 296 if fhInputTile is None: 297 return None 298 299 dec = AffineTransformDecorator(fhInputTile.GetGeoTransform()) 300 points = dec.pointsFor(fhInputTile.RasterXSize, fhInputTile.RasterYSize) 301 302 addFeature(g.TileIndexFieldName, ogrTileIndexDS, inputTile, points[0], points[1]) 303 del fhInputTile 304 305 if g.Verbose: 306 print("finished") 307 # ogrTileIndexDS.GetLayer().SyncToDisk() 308 return ogrTileIndexDS
309 310
311 -def getTargetDir(g, level=-1):
312 if level == -1: 313 return g.TargetDir 314 return g.TargetDir + str(level) + os.sep
315 316
317 -def tileImage(g, minfo, ti):
318 """ 319 320 Tile image in mosaicinfo minfo based on tileinfo ti 321 322 returns list of created tiles 323 324 """ 325 326 g.LastRowIndx = -1 327 OGRDS = createTileIndex(g.Verbose, "TileResult_0", g.TileIndexFieldName, g.Source_SRS, g.TileIndexDriverTyp) 328 329 yRange = list(range(1, ti.countTilesY + 1)) 330 xRange = list(range(1, ti.countTilesX + 1)) 331 332 if not g.Quiet and not g.Verbose: 333 progress(0.0) 334 processed = 0 335 total = len(xRange) * len(yRange) 336 337 for yIndex in yRange: 338 for xIndex in xRange: 339 offsetY = (yIndex - 1) * (ti.tileHeight - ti.overlap) 340 offsetX = (xIndex - 1) * (ti.tileWidth - ti.overlap) 341 height = ti.tileHeight 342 width = ti.tileWidth 343 if g.UseDirForEachRow: 344 tilename = getTileName(g, minfo, ti, xIndex, yIndex, 0) 345 else: 346 tilename = getTileName(g, minfo, ti, xIndex, yIndex) 347 348 if offsetX + width > ti.width: 349 width = ti.width - offsetX 350 if offsetY + height > ti.height: 351 height = ti.height - offsetY 352 353 feature_only = g.Resume and os.path.exists(tilename) 354 createTile(g, minfo, offsetX, offsetY, width, height, tilename, OGRDS, feature_only) 355 356 if not g.Quiet and not g.Verbose: 357 processed += 1 358 progress(processed / float(total)) 359 360 if g.TileIndexName is not None: 361 if g.UseDirForEachRow and not g.PyramidOnly: 362 shapeName = getTargetDir(g, 0) + g.TileIndexName 363 else: 364 shapeName = getTargetDir(g) + g.TileIndexName 365 copyTileIndexToDisk(g, OGRDS, shapeName) 366 367 if g.CsvFileName is not None: 368 if g.UseDirForEachRow and not g.PyramidOnly: 369 csvName = getTargetDir(g, 0) + g.CsvFileName 370 else: 371 csvName = getTargetDir(g) + g.CsvFileName 372 copyTileIndexToCSV(g, OGRDS, csvName) 373 374 return OGRDS
375 376
377 -def copyTileIndexToDisk(g, OGRDS, fileName):
378 SHAPEDS = createTileIndex(g.Verbose, fileName, g.TileIndexFieldName, OGRDS.GetLayer().GetSpatialRef(), "ESRI Shapefile") 379 OGRDS.GetLayer().ResetReading() 380 while True: 381 feature = OGRDS.GetLayer().GetNextFeature() 382 if feature is None: 383 break 384 newFeature = feature.Clone() 385 basename = os.path.basename(feature.GetField(0)) 386 if g.UseDirForEachRow: 387 t = os.path.split(os.path.dirname(feature.GetField(0))) 388 basename = t[1] + "/" + basename 389 newFeature.SetField(0, basename) 390 SHAPEDS.GetLayer().CreateFeature(newFeature) 391 closeTileIndex(SHAPEDS)
392 393
394 -def copyTileIndexToCSV(g, OGRDS, fileName):
395 csvfile = open(fileName, 'w') 396 OGRDS.GetLayer().ResetReading() 397 while True: 398 feature = OGRDS.GetLayer().GetNextFeature() 399 if feature is None: 400 break 401 basename = os.path.basename(feature.GetField(0)) 402 if g.UseDirForEachRow: 403 t = os.path.split(os.path.dirname(feature.GetField(0))) 404 basename = t[1] + "/" + basename 405 csvfile.write(basename) 406 geom = feature.GetGeometryRef() 407 coords = geom.GetEnvelope() 408 409 for coord in coords: 410 csvfile.write(g.CsvDelimiter) 411 csvfile.write("%f" % coord) 412 csvfile.write("\n") 413 414 csvfile.close()
415 416
417 -def createPyramidTile(g, levelMosaicInfo, offsetX, offsetY, width, height, tileName, OGRDS, feature_only):
418 419 temp_tilename = tileName + '.tmp' 420 421 sx = levelMosaicInfo.scaleX * 2 422 sy = levelMosaicInfo.scaleY * 2 423 424 dec = AffineTransformDecorator([levelMosaicInfo.ulx + offsetX * sx, sx, 0, 425 levelMosaicInfo.uly + offsetY * sy, 0, sy]) 426 427 if OGRDS is not None: 428 points = dec.pointsFor(width, height) 429 addFeature(g.TileIndexFieldName, OGRDS, tileName, points[0], points[1]) 430 431 if feature_only: 432 return 433 434 s_fh = levelMosaicInfo.getDataSet(dec.ulx, dec.uly + height * dec.scaleY, 435 dec.ulx + width * dec.scaleX, dec.uly) 436 if s_fh is None: 437 return 438 439 if g.BandType is None: 440 bt = levelMosaicInfo.band_type 441 else: 442 bt = g.BandType 443 444 geotransform = [dec.ulx, dec.scaleX, 0, dec.uly, 0, dec.scaleY] 445 446 bands = levelMosaicInfo.bands 447 448 if g.MemDriver is None: 449 t_fh = g.Driver.Create(temp_tilename, width, height, bands, bt, g.CreateOptions) 450 else: 451 t_fh = g.MemDriver.Create('', width, height, bands, bt) 452 453 if t_fh is None: 454 print('Creation failed, terminating gdal_tile.') 455 sys.exit(1) 456 457 t_fh.SetGeoTransform(geotransform) 458 t_fh.SetProjection(levelMosaicInfo.projection) 459 for band in range(1, bands + 1): 460 t_band = t_fh.GetRasterBand(band) 461 if levelMosaicInfo.ct is not None: 462 t_band.SetRasterColorTable(levelMosaicInfo.ct) 463 t_band.SetRasterColorInterpretation(levelMosaicInfo.ci[band - 1]) 464 465 if levelMosaicInfo.nodata is not None: 466 for band in range(1, bands + 1): 467 t_band = t_fh.GetRasterBand(band) 468 t_band.Fill(levelMosaicInfo.nodata) 469 t_band.SetNoDataValue(levelMosaicInfo.nodata) 470 471 res = gdal.ReprojectImage(s_fh, t_fh, None, None, g.ResamplingMethod) 472 if res != 0: 473 print("Reprojection failed for %s, error %d" % (temp_tilename, res)) 474 sys.exit(1) 475 476 levelMosaicInfo.closeDataSet(s_fh) 477 478 if g.MemDriver is None: 479 t_fh.FlushCache() 480 else: 481 tt_fh = g.Driver.CreateCopy(temp_tilename, t_fh, 0, g.CreateOptions) 482 tt_fh.FlushCache() 483 tt_fh = None 484 485 t_fh = None 486 487 if os.path.exists(tileName): 488 os.remove(tileName) 489 shutil.move(temp_tilename, tileName) 490 491 if g.Verbose: 492 print(tileName + " : " + str(offsetX) + "|" + str(offsetY) + "-->" + str(width) + "-" + str(height))
493 494
495 -def createTile(g, minfo, offsetX, offsetY, width, height, tilename, OGRDS, feature_only):
496 """ 497 498 Create tile 499 500 """ 501 temp_tilename = tilename + '.tmp' 502 503 if g.BandType is None: 504 bt = minfo.band_type 505 else: 506 bt = g.BandType 507 508 dec = AffineTransformDecorator([minfo.ulx, minfo.scaleX, 0, minfo.uly, 0, minfo.scaleY]) 509 510 geotransform = [dec.ulx + offsetX * dec.scaleX, dec.scaleX, 0, 511 dec.uly + offsetY * dec.scaleY, 0, dec.scaleY] 512 513 if OGRDS is not None: 514 dec2 = AffineTransformDecorator(geotransform) 515 points = dec2.pointsFor(width, height) 516 addFeature(g.TileIndexFieldName, OGRDS, tilename, points[0], points[1]) 517 518 if feature_only: 519 return 520 521 s_fh = minfo.getDataSet(dec.ulx + offsetX * dec.scaleX, dec.uly + offsetY * dec.scaleY + height * dec.scaleY, 522 dec.ulx + offsetX * dec.scaleX + width * dec.scaleX, 523 dec.uly + offsetY * dec.scaleY) 524 if s_fh is None: 525 return 526 527 bands = minfo.bands 528 529 if g.MemDriver is None: 530 t_fh = g.Driver.Create(temp_tilename, width, height, bands, bt, g.CreateOptions) 531 else: 532 t_fh = g.MemDriver.Create('', width, height, bands, bt) 533 534 if t_fh is None: 535 print('Creation failed, terminating gdal_tile.') 536 sys.exit(1) 537 538 t_fh.SetGeoTransform(geotransform) 539 if g.Source_SRS is not None: 540 t_fh.SetProjection(g.Source_SRS.ExportToWkt()) 541 542 readX = min(s_fh.RasterXSize, width) 543 readY = min(s_fh.RasterYSize, height) 544 for band in range(1, bands + 1): 545 s_band = s_fh.GetRasterBand(band) 546 t_band = t_fh.GetRasterBand(band) 547 if minfo.ct is not None: 548 t_band.SetRasterColorTable(minfo.ct) 549 if minfo.nodata is not None: 550 t_band.Fill(minfo.nodata) 551 t_band.SetNoDataValue(minfo.nodata) 552 553 data = s_band.ReadRaster(0, 0, readX, readY, readX, readY, t_band.DataType) 554 t_band.WriteRaster(0, 0, readX, readY, data, readX, readY, t_band.DataType) 555 556 minfo.closeDataSet(s_fh) 557 558 if g.MemDriver is None: 559 t_fh.FlushCache() 560 else: 561 tt_fh = g.Driver.CreateCopy(temp_tilename, t_fh, 0, g.CreateOptions) 562 tt_fh.FlushCache() 563 tt_fh = None 564 565 t_fh = None 566 567 if os.path.exists(tilename): 568 os.remove(tilename) 569 shutil.move(temp_tilename, tilename) 570 571 if g.Verbose: 572 print(tilename + " : " + str(offsetX) + "|" + str(offsetY) + "-->" + str(width) + "-" + str(height))
573 574
575 -def createTileIndex(Verbose, dsName, fieldName, srs, driverName):
576 OGRDriver = ogr.GetDriverByName(driverName) 577 if OGRDriver is None: 578 print('ESRI Shapefile driver not found') 579 sys.exit(1) 580 581 OGRDataSource = OGRDriver.Open(dsName) 582 if OGRDataSource is not None: 583 OGRDataSource.Destroy() 584 OGRDriver.DeleteDataSource(dsName) 585 if Verbose: 586 print('truncating index ' + dsName) 587 588 OGRDataSource = OGRDriver.CreateDataSource(dsName) 589 if OGRDataSource is None: 590 print('Could not open datasource ' + dsName) 591 sys.exit(1) 592 593 OGRLayer = OGRDataSource.CreateLayer("index", srs, ogr.wkbPolygon) 594 if OGRLayer is None: 595 print('Could not create Layer') 596 sys.exit(1) 597 598 OGRFieldDefn = ogr.FieldDefn(fieldName, ogr.OFTString) 599 if OGRFieldDefn is None: 600 print('Could not create FieldDefn for ' + fieldName) 601 sys.exit(1) 602 603 OGRFieldDefn.SetWidth(256) 604 if OGRLayer.CreateField(OGRFieldDefn) != 0: 605 print('Could not create Field for ' + fieldName) 606 sys.exit(1) 607 608 return OGRDataSource
609 610
611 -def addFeature(TileIndexFieldName, OGRDataSource, location, xlist, ylist):
612 OGRLayer = OGRDataSource.GetLayer() 613 OGRFeature = ogr.Feature(OGRLayer.GetLayerDefn()) 614 if OGRFeature is None: 615 print('Could not create Feature') 616 sys.exit(1) 617 618 OGRFeature.SetField(TileIndexFieldName, location) 619 wkt = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f ))' % (xlist[0], ylist[0], 620 xlist[1], ylist[1], xlist[2], ylist[2], xlist[3], ylist[3], xlist[0], ylist[0]) 621 OGRGeometry = ogr.CreateGeometryFromWkt(wkt, OGRLayer.GetSpatialRef()) 622 if OGRGeometry is None: 623 print('Could not create Geometry') 624 sys.exit(1) 625 626 OGRFeature.SetGeometryDirectly(OGRGeometry) 627 628 OGRLayer.CreateFeature(OGRFeature) 629 OGRFeature.Destroy()
630 631
632 -def closeTileIndex(OGRDataSource):
633 OGRDataSource.Destroy()
634 635
636 -def buildPyramid(g, minfo, createdTileIndexDS, tileWidth, tileHeight, overlap):
637 inputDS = createdTileIndexDS 638 for level in range(1, g.Levels + 1): 639 g.LastRowIndx = -1 640 levelMosaicInfo = mosaic_info(minfo.filename, inputDS) 641 levelOutputTileInfo = tile_info(int(levelMosaicInfo.xsize / 2), int(levelMosaicInfo.ysize / 2), tileWidth, tileHeight, overlap) 642 inputDS = buildPyramidLevel(g, levelMosaicInfo, levelOutputTileInfo, level)
643 644
645 -def buildPyramidLevel(g, levelMosaicInfo, levelOutputTileInfo, level):
646 yRange = list(range(1, levelOutputTileInfo.countTilesY + 1)) 647 xRange = list(range(1, levelOutputTileInfo.countTilesX + 1)) 648 649 OGRDS = createTileIndex(g.Verbose, "TileResult_" + str(level), g.TileIndexFieldName, g.Source_SRS, g.TileIndexDriverTyp) 650 651 for yIndex in yRange: 652 for xIndex in xRange: 653 offsetY = (yIndex - 1) * (levelOutputTileInfo.tileHeight - levelOutputTileInfo.overlap) 654 offsetX = (xIndex - 1) * (levelOutputTileInfo.tileWidth - levelOutputTileInfo.overlap) 655 height = levelOutputTileInfo.tileHeight 656 width = levelOutputTileInfo.tileWidth 657 658 if offsetX + width > levelOutputTileInfo.width: 659 width = levelOutputTileInfo.width - offsetX 660 if offsetY + height > levelOutputTileInfo.height: 661 height = levelOutputTileInfo.height - offsetY 662 663 tilename = getTileName(g, levelMosaicInfo, levelOutputTileInfo, xIndex, yIndex, level) 664 665 feature_only = g.Resume and os.path.exists(tilename) 666 createPyramidTile(g, levelMosaicInfo, offsetX, offsetY, width, height, tilename, OGRDS, feature_only) 667 668 if g.TileIndexName is not None: 669 shapeName = getTargetDir(g, level) + g.TileIndexName 670 copyTileIndexToDisk(g, OGRDS, shapeName) 671 672 if g.CsvFileName is not None: 673 csvName = getTargetDir(g, level) + g.CsvFileName 674 copyTileIndexToCSV(g, OGRDS, csvName) 675 676 return OGRDS
677 678
679 -def getTileName(g, minfo, ti, xIndex, yIndex, level=-1):
680 """ 681 creates the tile file name 682 """ 683 maxim = ti.countTilesX 684 if ti.countTilesY > maxim: 685 maxim = ti.countTilesY 686 countDigits = len(str(maxim)) 687 parts = os.path.splitext(os.path.basename(minfo.filename)) 688 if parts[0][0] == "@": # remove possible leading "@" 689 parts = (parts[0][1:len(parts[0])], parts[1]) 690 691 yIndex_str = ("%0" + str(countDigits) + "i") % (yIndex,) 692 xIndex_str = ("%0" + str(countDigits) + "i") % (xIndex,) 693 694 if g.UseDirForEachRow: 695 frmt = getTargetDir(g, level) + str(yIndex) + os.sep + parts[0] + "_" + yIndex_str + "_" + xIndex_str 696 # See if there was a switch in the row, if so then create new dir for row. 697 if g.LastRowIndx < yIndex: 698 g.LastRowIndx = yIndex 699 if not os.path.exists(getTargetDir(g, level) + str(yIndex)): 700 os.mkdir(getTargetDir(g, level) + str(yIndex)) 701 else: 702 frmt = getTargetDir(g, level) + parts[0] + "_" + yIndex_str + "_" + xIndex_str 703 # Check for the extension that should be used. 704 if g.Extension is None: 705 frmt += parts[1] 706 else: 707 frmt += "." + g.Extension 708 return frmt
709 710
711 -def UsageFormat():
712 print('Valid formats:') 713 count = gdal.GetDriverCount() 714 for index in range(count): 715 driver = gdal.GetDriver(index) 716 print(driver.ShortName)
717 718 # ============================================================================= 719 720
721 -def Usage():
722 print('Usage: gdal_retile.py ') 723 print(' [-v] [-q] [-co NAME=VALUE]* [-of out_format]') 724 print(' [-ps pixelWidth pixelHeight]') 725 print(' [-overlap val_in_pixel]') 726 print(' [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/') 727 print(' CInt16/CInt32/CFloat32/CFloat64}]') 728 print(' [ -tileIndex tileIndexName [-tileIndexField fieldName]]') 729 print(' [ -csv fileName [-csvDelim delimiter]]') 730 print(' [-s_srs srs_def] [-pyramidOnly] -levels numberoflevels') 731 print(' [-r {near/bilinear/cubic/cubicspline/lanczos}]') 732 print(' [-useDirForEachRow] [-resume]') 733 print(' -targetDir TileDirectory input_files')
734 735 # ============================================================================= 736
737 -def main(args=None, g=None):
738 if g is None: 739 g = RetileGlobals() 740 741 gdal.AllRegister() 742 743 if args is None: 744 args = sys.argv 745 argv = gdal.GeneralCmdLineProcessor(args) 746 if argv is None: 747 return 1 748 749 # Parse command line arguments. 750 i = 1 751 while i < len(argv): 752 arg = argv[i] 753 754 if arg == '-of' or arg == '-f': 755 i += 1 756 g.Format = argv[i] 757 elif arg == '-ot': 758 i += 1 759 g.BandType = gdal.GetDataTypeByName(argv[i]) 760 if g.BandType == gdal.GDT_Unknown: 761 print('Unknown GDAL data type: %s' % argv[i]) 762 return 1 763 elif arg == '-co': 764 i += 1 765 g.CreateOptions.append(argv[i]) 766 767 elif arg == '-v': 768 g.Verbose = True 769 elif arg == '-q': 770 g.Quiet = True 771 772 elif arg == '-targetDir': 773 i += 1 774 g.TargetDir = argv[i] 775 776 if not os.path.exists(g.TargetDir): 777 print("TargetDir " + g.TargetDir + " does not exist") 778 return 1 779 if g.TargetDir[len(g.TargetDir) - 1:] != os.sep: 780 g.TargetDir = g.TargetDir + os.sep 781 782 elif arg == '-ps': 783 i += 1 784 g.TileWidth = int(argv[i]) 785 i += 1 786 g.TileHeight = int(argv[i]) 787 788 elif arg == '-overlap': 789 i += 1 790 g.Overlap = int(argv[i]) 791 792 elif arg == '-r': 793 i += 1 794 ResamplingMethodString = argv[i] 795 if ResamplingMethodString == "near": 796 g.ResamplingMethod = gdal.GRA_NearestNeighbour 797 elif ResamplingMethodString == "bilinear": 798 g.ResamplingMethod = gdal.GRA_Bilinear 799 elif ResamplingMethodString == "cubic": 800 g.ResamplingMethod = gdal.GRA_Cubic 801 elif ResamplingMethodString == "cubicspline": 802 g.ResamplingMethod = gdal.GRA_CubicSpline 803 elif ResamplingMethodString == "lanczos": 804 g.ResamplingMethod = gdal.GRA_Lanczos 805 else: 806 print("Unknown resampling method: %s" % ResamplingMethodString) 807 return 1 808 elif arg == '-levels': 809 i += 1 810 g.Levels = int(argv[i]) 811 if g.Levels < 1: 812 print("Invalid number of levels : %d" % g.Levels) 813 return 1 814 elif arg == '-s_srs': 815 i += 1 816 g.Source_SRS = osr.SpatialReference() 817 if g.Source_SRS.SetFromUserInput(argv[i]) != 0: 818 print('invalid -s_srs: ' + argv[i]) 819 return 1 820 821 elif arg == "-pyramidOnly": 822 g.PyramidOnly = True 823 elif arg == '-tileIndex': 824 i += 1 825 g.TileIndexName = argv[i] 826 parts = os.path.splitext(g.TileIndexName) 827 if not parts[1]: 828 g.TileIndexName += ".shp" 829 830 elif arg == '-tileIndexField': 831 i += 1 832 g.TileIndexFieldName = argv[i] 833 elif arg == '-csv': 834 i += 1 835 g.CsvFileName = argv[i] 836 parts = os.path.splitext(g.CsvFileName) 837 if not parts[1]: 838 g.CsvFileName += ".csv" 839 elif arg == '-csvDelim': 840 i += 1 841 g.CsvDelimiter = argv[i] 842 elif arg == '-useDirForEachRow': 843 g.UseDirForEachRow = True 844 elif arg == "-resume": 845 g.Resume = True 846 elif arg[:1] == '-': 847 print('Unrecognized command option: %s' % arg) 848 Usage() 849 return 1 850 851 else: 852 g.Names.append(arg) 853 i += 1 854 855 if not g.Names: 856 print('No input files selected.') 857 Usage() 858 return 1 859 860 if (g.TileWidth == 0 or g.TileHeight == 0): 861 print("Invalid tile dimension %d,%d" % (g.TileWidth, g.TileHeight)) 862 return 1 863 if (g.TileWidth - g.Overlap <= 0 or g.TileHeight - g.Overlap <= 0): 864 print("Overlap too big w.r.t tile height/width") 865 return 1 866 867 if g.TargetDir is None: 868 print("Missing Directory for Tiles -targetDir") 869 Usage() 870 return 1 871 872 # create level 0 directory if needed 873 if g.UseDirForEachRow and not g.PyramidOnly: 874 leveldir = g.TargetDir + str(0) + os.sep 875 if not os.path.exists(leveldir): 876 os.mkdir(leveldir) 877 878 if g.Levels > 0: # prepare Dirs for pyramid 879 startIndx = 1 880 for levelIndx in range(startIndx, g.Levels + 1): 881 leveldir = g.TargetDir + str(levelIndx) + os.sep 882 if os.path.exists(leveldir): 883 continue 884 os.mkdir(leveldir) 885 if not os.path.exists(leveldir): 886 print("Cannot create level dir: %s" % leveldir) 887 return 1 888 if g.Verbose: 889 print("Created level dir: %s" % leveldir) 890 891 g.Driver = gdal.GetDriverByName(g.Format) 892 if g.Driver is None: 893 print('Format driver %s not found, pick a supported driver.' % g.Format) 894 UsageFormat() 895 return 1 896 897 DriverMD = g.Driver.GetMetadata() 898 g.Extension = DriverMD.get(gdal.DMD_EXTENSION) 899 if 'DCAP_CREATE' not in DriverMD: 900 g.MemDriver = gdal.GetDriverByName("MEM") 901 902 tileIndexDS = getTileIndexFromFiles(g) 903 if tileIndexDS is None: 904 print("Error building tile index") 905 return 1 906 minfo = mosaic_info(g.Names[0], tileIndexDS) 907 ti = tile_info(minfo.xsize, minfo.ysize, g.TileWidth, g.TileHeight, g.Overlap) 908 909 if g.Source_SRS is None and minfo.projection: 910 g.Source_SRS = osr.SpatialReference() 911 if g.Source_SRS.SetFromUserInput(minfo.projection) != 0: 912 print('invalid projection ' + minfo.projection) 913 return 1 914 915 if g.Verbose: 916 minfo.report() 917 ti.report() 918 919 if not g.PyramidOnly: 920 dsCreatedTileIndex = tileImage(g, minfo, ti) 921 tileIndexDS.Destroy() 922 else: 923 dsCreatedTileIndex = tileIndexDS 924 925 if g.Levels > 0: 926 buildPyramid(g, minfo, dsCreatedTileIndex, g.TileWidth, g.TileHeight, g.Overlap) 927 928 if g.Verbose: 929 print("FINISHED") 930 return 0
931 932
933 -class RetileGlobals():
934 __slots__ = ['Verbose', 935 'Quiet', 936 'CreateOptions', 937 'Names', 938 'TileWidth', 939 'TileHeight', 940 'Overlap', 941 'Format', 942 'BandType', 943 'Driver', 944 'Extension', 945 'MemDriver', 946 'TileIndexFieldName', 947 'TileIndexName', 948 'TileIndexDriverTyp', 949 'CsvDelimiter', 950 'CsvFileName', 951 'Source_SRS', 952 'TargetDir', 953 'ResamplingMethod', 954 'Levels', 955 'PyramidOnly', 956 'LastRowIndx', 957 'UseDirForEachRow', 958 'Resume'] 959
960 - def __init__(self):
961 """ Only used for unit tests """ 962 self.Verbose = False 963 self.Quiet = False 964 self.CreateOptions = [] 965 self.Names = [] 966 self.TileWidth = 256 967 self.TileHeight = 256 968 self.Overlap = 0 969 self.Format = 'GTiff' 970 self.BandType = None 971 self.Driver = None 972 self.Extension = None 973 self.MemDriver = None 974 self.TileIndexFieldName = 'location' 975 self.TileIndexName = None 976 self.TileIndexDriverTyp = "Memory" 977 self.CsvDelimiter = ";" 978 self.CsvFileName = None 979 980 self.Source_SRS = None 981 self.TargetDir = None 982 self.ResamplingMethod = gdal.GRA_NearestNeighbour 983 self.Levels = 0 984 self.PyramidOnly = False 985 self.LastRowIndx = -1 986 self.UseDirForEachRow = False 987 self.Resume = False
988 989 990 if __name__ == '__main__': 991 sys.exit(main(sys.argv)) 992