I've been working with the gdal_retile.py script recently to retile an orthomosaic. I originally was using the solution from this post, but found that gdal_retile.py is faster with orthos > 10GB. There are a couple of features I added to that rasterio method that I want to add to the gdal_retile.py script, but I am just not sure where to add them.
Mainly, I do not want "empty" images. In my rasterio script I had if np.all(data.read(window=window) == 0: continue. These images were basically the images that were between the extent boundary and where the ortho actually contained the picture (data) - see attached image.
I've been trying to understand what the script is doing, but ultimately it's a bit beyond me. I believe I need to focus on the createTile function. I was able to find out if the data that is being written is empty, but I am unsure how to skip writing the raster when the conditions are met.
if bytearray(len(data)) == data:
print('empty')
continue
def createTile(minfo, offsetX, offsetY, width, height, tilename, OGRDS):
"""
Create tile
return name of created tile
"""
if BandType is None:
bt = minfo.band_type
else:
bt = BandType
dec = AffineTransformDecorator([minfo.ulx, minfo.scaleX, 0, minfo.uly, 0, minfo.scaleY])
s_fh = minfo.getDataSet(dec.ulx + offsetX * dec.scaleX, dec.uly + offsetY * dec.scaleY + height * dec.scaleY,
dec.ulx + offsetX * dec.scaleX + width * dec.scaleX,
dec.uly + offsetY * dec.scaleY)
if s_fh is None:
return
geotransform = [dec.ulx + offsetX * dec.scaleX, dec.scaleX, 0,
dec.uly + offsetY * dec.scaleY, 0, dec.scaleY]
if OGRDS is not None:
dec2 = AffineTransformDecorator(geotransform)
points = dec2.pointsFor(width, height)
addFeature(OGRDS, tilename, points[0], points[1])
bands = minfo.bands
if MemDriver is None:
t_fh = Driver.Create(tilename, width, height, bands, bt, CreateOptions)
else:
t_fh = MemDriver.Create(tilename, width, height, bands, bt)
if t_fh is None:
print('Creation failed, terminating gdal_tile.')
sys.exit(1)
t_fh.SetGeoTransform(geotransform)
if Source_SRS is not None:
t_fh.SetProjection(Source_SRS.ExportToWkt())
readX = min(s_fh.RasterXSize, width)
readY = min(s_fh.RasterYSize, height)
for band in range(1, bands + 1):
s_band = s_fh.GetRasterBand(band)
t_band = t_fh.GetRasterBand(band)
if minfo.ct is not None:
t_band.SetRasterColorTable(minfo.ct)
if minfo.nodata is not None:
t_band.Fill(minfo.nodata)
t_band.SetNoDataValue(minfo.nodata)
data = s_band.ReadRaster(0, 0, readX, readY, readX, readY, t_band.DataType)
if bytearray(len(data)) == data: # <---- Added these three lines
print('empty') # <-----
continue # <---------
t_band.WriteRaster(0, 0, readX, readY, data, readX, readY, t_band.DataType)
minfo.closeDataSet(s_fh)
if MemDriver is not None:
tt_fh = Driver.CreateCopy(tilename, t_fh, 0, CreateOptions)
tt_fh.FlushCache()
if Verbose:
print(tilename + " : " + str(offsetX) + "|" + str(offsetY) + "-->" + str(width) + "-" + str(height))
Blue Tiles would be excluded when saving, but the grey tile would be saved because it contains data from the ortho.
