4

I have quite a big GDAL VRT with a pyramid (TIFFs in smaller zoom levels stored externally in separate folders) and the calculation of those takes quite some time.

I am also having regular updates on the VRT and I don't want to calculate the complete pyramid, because its not necessary. So my question is, do gdal_retile or gdaladdo or any other GDAL tool have the possibility to recalculate all the overviews for one image within the pyramid?

  • You could build higher res overviews on the tifs for zoomed in views and lower res overviews on the vrt for displaying large areas gis.stackexchange.com/a/148343/2856 – user2856 Jul 2 '17 at 9:42
  • What are you trying to do in the end? I had the same issue but decided to use GroServer ImageMosaic instead of VRT. – pLumo Aug 1 '17 at 14:13
  • I removed the vrt, because there was no advantage of having it. Now i am having the pyramid without a VRT file. The problem was, that as far as i know there is no script to update the overviews of this pyramid. – aleho Aug 1 '17 at 14:21
3

So, as far as I know there is not a ready-to-use solution to run incremental updates on a pyramid with gdal. So I wrote it on my own - the following script will execute gdal_retile only for the updated overviews:

    #new_file= original tif-file within the pyramid which should be updated        
    levels=[1, 2, 3, 4] #Zoom levels within the pyramid

    #1. Get the center of the tif-file as wkt-geometry        
    try:
       raster = gdal.Open(new_file)
       geoTransform = raster.GetGeoTransform()
       x_min = geoTransform[0]
       y_max = geoTransform[3]
       x_max = x_min + geoTransform[1] * raster.RasterXSize
       y_min = y_max + geoTransform[5] * raster.RasterYSize

    except RuntimeError, e:
       print 'Unable to open tif file'
       print e
       break

    x = (x_max + x_min) / 2
    y = (y_min + y_max) / 2
    wkt_point = 'POINT ({0} {1})'.format(x, y)   

    for level in levels:
       #The overviews of each level are in a separate folder (level_folder). In each folder is a shape file with the boundaries of the tiles and the path to the tif-file. 

       print 'Calculating overviews for level {0}'.format(level)

       level_folder = src_dir + '{0}/'.format(level)
       level_shape = level_folder + '{0}.shp'.format(level)

       #Load the shapefile and get the features intersecting with the center of the tif file (wkt_point) -> this will give you the names of the tifs which should be updated
       driver = ogr.GetDriverByName('ESRI Shapefile')    
       datasource = driver.Open(level_shape, 0)

       if datasource is None:
           print 'Error loading Level Shapefile'
           break

       level_layer = datasource.GetLayer()
       level_layer.SetSpatialFilter(ogr.CreateGeometryFromWkt(wkt_point))

       for feature in level_layer:

           to_be_replaced_tile = level_folder + feature.GetField('location')     

           #here i am basically writing the files which should be updated to a txt-file. The txt-file is then passed to the gdal retile command 
           return_tuple = create_txt_file(update_dir, level, level_folder, feature.GetGeometryRef())

           txt_file = return_tuple[0]
           first_file_name = return_tuple[1]

           cmd = ' '.join(['\"C:/Program Files/QGIS 2.18/bin/python.exe\" \"C:/Program Files/QGIS 2.18/bin/gdal_retile.py\"','-r', 'bilinear', '-s_srs', 'EPSG:25832', '-ps', '8000', '8000', '-pyramidOnly', '-targetDir', update_dir, '--optfile', txt_file, '-levels', '1'])

           s_info = subprocess.STARTUPINFO()
           s_info.dwFlags |= subprocess.STARTF_USESHOWWINDOW

           try:

               subprocess.check_call(cmd, startupinfo=s_info)

           except subprocess.CalledProcessError as e:
               print 'error message: {0}'.format(e.output)

           new_file = update_dir + '/1/' + first_file_name.replace('.tif', '_1_1.tif')
           new_file = new_file.replace(resource_level_folder, update_dir)
           replace_file(to_be_replaced_tile, new_file)

       print 'finished writing files for level {0}'.format(level)

    def create_txt_file(update_dir, level, resource_level_folder, geom):

       driver = ogr.GetDriverByName('ESRI Shapefile')    
       datasource = driver.Open(resource_level_folder, 0)

       if datasource is None:
          print 'Error loading Shapefile'
          return None

       resource_level_layer = datasource.GetLayer()        
       target = osr.SpatialReference()
       target.ImportFromEPSG(25832)
       bbox = geom.Buffer(-1)

       resource_level_layer.SetSpatialFilter(bbox)
       list_tiles = []
       first_file_name = ''

       for resource_feature in resource_level_layer:

          if not first_file_name:
             first_file_name = resource_feature.GetField('location')

          list_tiles.append(resource_level_folder + resource_feature.GetField('location'))

       f = open(update_dir + '/update_level{0}.txt'.format(level), 'w')
       f.write("\n".join(list_tiles))
       f.close()

       return [update_dir + '/update_level{0}.txt'.format(level), first_file_name]
|improve this answer|||||
  • could you explain in words how the script is working ? – pLumo Aug 1 '17 at 14:10
  • I added some comments and extended the code. Maybe its easier now... – aleho Aug 1 '17 at 14:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.