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?
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]