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