1

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.

enter image description here

1 Answer 1

1

I found that the .tif was actually being created on the lines if MemDriver is None: and not where t_band.WriteRaster() is located. So checking for an empty .tif should go right above that. See the new code below:

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


    #### ADDED HERE ####
    data = s_fh.GetRasterBand(1).ReadRaster()
    if bytearray(len(data)) == data:
        return
    #### ENDED HERE ####


    if MemDriver is None:
        t_fh = Driver.Create(tilename, width, height, bands, bt, CreateOptions)
    else:
        t_fh = MemDriver.Create(tilename, width, height, bands, bt)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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