tiling

Routines in support of tiled segmentation of very large rasters.

Main entry routine is doTiledShepherdSegmentation(). See that function for further details.

The broad idea is that the Shepherd segmentation algorithm, as implemented in the shepseg module, runs entirely in memory. For larger raster files, it is more efficient to divide the raster into tiles, segment each tile individually, and stitch the results together to create a segmentation of the whole raster.

The main caveats arise from the fact that the initial clustering is performed on a uniform subsample of the whole image, in order to give consistent segment boundaries at tile intersections. This means that for a larger raster, with a greater range of spectra, one may wish to increase the number of clusters in order to allow sufficient initial segments to characterize the variation.

Related to this, one may also consider reducing the percentile used for automatic estimation of maxSpectralDiff (see pyshepseg.shepseg.doShepherdSegmentation() and pyshepseg.shepseg.autoMaxSpectralDiff() for further details).

Because of these caveats, one should be very cautious about segmenting something like a continental-scale image. There is a lot of spectral variation across an area like a whole continent, and it may be unwise to use all the same parameters for the whole area.

Also included in this module is functionality for subsetting a large segmented image into a smaller one for checking etc. See subsetImage().

exception pyshepseg.tiling.PyShepSegSubsetError
exception pyshepseg.tiling.PyShepSegTilingError
class pyshepseg.tiling.RatPage(numIntCols, numFloatCols, startSegId, numSeg)

Hold a single page of the paged RAT

getIndexInPage(segId)

Return the index for the given segment, within the current page.

getRatVal(segId, colType, colArrayNdx)

Get the RAT entry for the given segment.

getSegmentComplete(segId)

Returns True if the segment has been flagged as complete

pageComplete()

Return True if the current page has been completed

setRatVal(segId, colType, colArrayNdx, val)

Set the RAT entry for the given segment, to be the given value.

setSegmentComplete(segId)

Flag that the given segment has had all stats calculated.

class pyshepseg.tiling.TileInfo

Class that holds the pixel coordinates of the tiles within an image.

addTile(xpos, ypos, xsize, ysize, col, row)

Add a new tile to the set

Parameters
xpos, yposint

Pixel column & row of top left pixel of tile

xsize, ysizeint

Number of pixel columns & rows in tile

col, rowint

Tile column & row

getNumTiles()

Get total number of tiles in the set

Returns
numTilesint

Total number of tiles

getTile(col, row)

Return the position and shape of the requested tile, as a single tuple of values

Parameters
col, rowint

Tile column & row

Returns
xpos, yposint

Pixel column & row of top left pixel of tile

xsize, ysizeint

Number of pixel columns & rows in tile

class pyshepseg.tiling.TiledSegmentationResult

Result of tiled segmentation

Attributes
maxSegIdshepseg.SegIdType

Largest segment ID used in final segment image

numTileRowsint

Number of rows of tiles used

numTileColsint

Number of columns of tiles used

subsamplePcntfloat

Percentage of image subsampled for clustering

maxSpectralDifffloat

The value used to limit segment merging (in all tiles)

kmeanssklearn.cluster.KMeans

The sklearn KMeans object, after fitting

hasEmptySegmentsbool

True if the segmentation contains segments with no pixels. This is an error condition, probably indicating that the merging of segments across tiles has produced inconsistent numbering. A warning message will also have been printed.

pyshepseg.tiling.calcHistogramTiled(segfile, maxSegId, writeToRat=True)

Calculate a histogram of the given segment image file.

Note that we need this function because GDAL’s GetHistogram function does not seem to work when attempting a histogram with very large numbers of entries. We want an entry for every segment, rather than an approximate count for a range of segment values, and the number of segments is very large. So, we need to write our own routine.

It works in tiles across the image, so that it can process very large images in a memory-efficient way.

For a raster which can easily fit into memory, a histogram can be calculated directly using pyshepseg.shepseg.makeSegSize().

Parameters
segfilestr or gdal.Dataset

Segmentation image file. Can be either the file name string, or an open Dataset object.

maxSegIdshepseg.SegIdType

Maximum segment ID used

writeToRatbool

If True, the completed histogram will be written to the image file’s raster attribute table. If segfile was given as a Dataset object, it would therefore need to have been opened with update access.

Returns
histint ndarray (numSegments+1, )

Histogram counts for each segment (index is segment ID number)

pyshepseg.tiling.checkForEmptySegments(outfile, maxSegId, overlapSize)

Check the final segmentation for any empty segments. These can be problematic later, and should be avoided. Prints a warning message if empty segments are found.

Parameters
outfilestr

File name of segmentation image to check

maxSegIdshepseg.SegIdType

Maximum segment ID used

overlapSizeint

Number of pixels to use in overlaps between tiles

Returns
hasEmptySegmentsbool

True if there are segment ID numbers with no pixels

pyshepseg.tiling.copyColumns(inRat, outRat)

Copy column structure from inRat to outRat. Note that this just creates the empty columns in outRat, it does not copy any data.

Parameters
inRat, outRatgdal.RasterAttributeTable

Columns found in inRat are created on outRat

Returns
numIntColsint

Number of integer columns found

numFloatColsint

Number of float columns found

pyshepseg.tiling.copySubsettedSegmentsToNew(inPage, outPagedRat, recodeDict)

Using the recodeDict, copy across the rows inPage to outPage.

inPage is processed and (taking into account of inPage.startSegId) the original input row found. This value is then looked up in recodeDict to find the row in the output RAT to copy the row from the input to.

Parameters
inPageRatPage

A page of RAT from the input file

outPagedRatnumba.typed.Dict

In-memory pages of the output RAT, as created by createPagedRat(). This is modified in-place, creating new pages as required.

recodeDictnumba.typed.Dict

Keyed by original segment ID, values are the corresponding segment IDs in the subset

pyshepseg.tiling.createPagedRat()

Create the dictionary for the paged RAT. Each element is a page of the RAT, with entries for a range of segment IDs. The key is the segment ID of the first entry in the page.

The returned dictionary is initially empty.

pyshepseg.tiling.crossesMidline(overlap, segLoc, orientation)

Check whether the given segment crosses the midline of the given overlap. If it does not, then it will lie entirely within exactly one tile, but if it does cross, then it will need to be re-coded across the midline.

Parameters
overlapshepseg.SegIdType ndarray (overlapNrows, overlapNcols)

Array of segments just for this overlap region

segLocshepseg.RowColArray

The row/col coordinates (within the overlap array) of the pixels for the segment of interest

orientation{HORIZONTAL, VERTICAL}

Indicates the orientation of the midline

Returns
crossesbool

True if the given segment crosses the midline

pyshepseg.tiling.doTiledShepherdSegmentation(infile, outfile, tileSize=4096, overlapSize=1024, minSegmentSize=50, numClusters=60, bandNumbers=None, subsamplePcnt=None, maxSpectralDiff='auto', imgNullVal=None, fixedKMeansInit=False, fourConnected=True, verbose=False, simpleTileRecode=False, outputDriver='KEA', creationOptions=[], spectDistPcntile=50, kmeansObj=None, tempfilesDriver='KEA', tempfilesExt='kea', tempfilesCreationOptions=[])

Run the Shepherd segmentation algorithm in a memory-efficient manner, suitable for large raster files. Runs the segmentation on separate (overlapping) tiles across the raster, then stitches these together into a single output segment raster.

The initial spectral clustering is performed on a sub-sample of the whole raster (using fitSpectralClustersWholeFile), to create consistent clusters. These are then used as seeds for all individual tiles. Note that subsamplePcnt is used at this stage, over the whole raster, and is not passed through to shepseg.doShepherdSegmentation() for any further sub-sampling.

Most of the arguments are passed through to shepseg.doShepherdSegmentation, and are described in the docstring for that function.

Parameters
infilestr

Filename of input raster

outfilestr

Filename of output segmentation raster

tileSizeint

Desired width & height (in pixels) of the tiles (i.e. desired tiles have shape (tileSize, tileSize). Tiles on the right and bottom edges of the input image may end up slightly larger than tileSize to ensure there are no small tiles.

overlapSizeint

Number of pixels to overlap tiles (??? should give more precise detail….)

minSegmentSizeint

Minimum number of pixels in a segment

numClustersint

Number of clusters to request in k-means clustering

bandNumberslist of int

The GDAL band numbers (i.e. start at 1) of the bands of input raster to use for segmentation

subsamplePcntfloat or None

See fitSpectralClustersWholeFile()

maxSpectralDifffloat or str

See shepseg.doShepherdSegmentation()

spectDistPcntileint

See shepseg.doShepherdSegmentation()

imgNullValfloat or None

If given, use this as the null value for the input raster. If None, use the value defined in the raster file

fixedKMeansInitbool

If True, use a fixed set of initial cluster centres for the KMeans clustering. This is good to ensure exactly reproducible results

fourConnectedbool

If True, use 4-way connectedness, otherwise use 8-way

verbosebool

If True, print informative messages during processing (to stdout)

simpleTileRecodebool

If True, use only a simple tile recoding procedure. See stitchTiles() for more detail

outputDriverstr

The short name of the GDAL format driver to use for output file

creationOptionslist of str

The GDAL output creation options to match the outputDriver

kmeansObjsklearn.cluster.KMeans

See shepseg.doShepherdSegmentation() for details

tempfilesDriverstr

Short name of GDAL driver to use for temporary raster files

tempfilesExtstr

File extension to use for temporary raster files

tempfilesCreationOptionslist of str

GDAL creation options to use for temporary raster files

Returns
tileSegResultTiledSegmentationResult
pyshepseg.tiling.doTiledShepherdSegmentation_doOne(inDs, filename, tileInfo, col, row, bandNumbers, imgNullVal, kmeansObj, minSegmentSize=50, maxSpectralDiff='auto', verbose=False, spectDistPcntile=50, fourConnected=True, tempfilesDriver='KEA', tempfilesCreationOptions=[])

Called from doTiledShepherdSegmentation(). Does a single tile, and is split out here as a separate function so it can be called from in parallel with other tiles if desired.

tileInfo is object returned from doTiledShepherdSegmentation_prepare() and col, row describe the tile that this call will process.

See doTiledShepherdSegmentation() for detailed descriptions of other parameters.

Return value is that from shepseg.doShepherdSegmentation().

pyshepseg.tiling.doTiledShepherdSegmentation_finalize(inDs, outfile, tileFilenames, tileInfo, overlapSize, tempDir, simpleTileRecode=False, outputDriver='KEA', creationOptions=[], verbose=False)

Do the stitching of tiles and check for empty segments. Call after every doTiledShepherdSegmentation_doOne() has completed for a given tiled segmentation.

Returns a tuple with (axSegId, hasEmptySegments).

pyshepseg.tiling.doTiledShepherdSegmentation_prepare(infile, tileSize=4096, overlapSize=1024, bandNumbers=None, imgNullVal=None, kmeansObj=None, verbose=False, numClusters=60, subsamplePcnt=None, fixedKMeansInit=False)

Do all the preparation for the tiled segmentation. Call this first if creating a parallel implementation, then call doTiledShepherdSegmentation_doOne() for each tile in the returned TileInfo object.

See doTiledShepherdSegmentation() for detailed parameter descriptions.

Returns a tuple with: (datasetObj, bandNumbers, kmeansObj, subsamplePcnt, imgNullVal, tileInfo)

pyshepseg.tiling.fitSpectralClustersWholeFile(inDs, bandNumbers, numClusters=60, subsamplePcnt=None, imgNullVal=None, fixedKMeansInit=False)

Given an open raster Dataset, read a selected sample of pixels and use these to fit a spectral cluster model. Uses GDAL to read the pixels, and shepseg.fitSpectralClusters() to do the fitting.

Parameters
inDsgdal.Dataset

Open GDAL Dataset object for the input raster

bandNumberslist of int (or None)

List of GDAL band numbers for the bands of interest. If None, then use all bands in the dataset. Note that GDAL band numbers start at 1.

numClustersint

Desired number of clusters

subsamplePcntfloat or None

Percentage of pixels to use in fitting. If it is None, then a suitable subsample is calculated such that around one million pixels are sampled. (Note - this would include null pixels, so if the image is dominated by nulls, this would undersample.) No further subsampling is carried out by fitSpectralClusters().

imgNullValfloat or None

Pixels with this value in the input raster are ignored. If None, the NoDataValue from the raster file is used

fixedKMeansInitbool

If True, then use a fixed estimate for the initial KMeans cluster centres. See shepseg.fitSpectralClusters() for details.

Returns
kmeansObjsklearn.cluster.KMeans

The fitted KMeans object

subsamplePcntfloat

The subsample percentage actually used

imgNullValfloat

The image null value (possibly read from the file)

pyshepseg.tiling.getImgNullValue(inDs, bandNumbers)

Return the null value for the given dataset

Parameters
inDsgdal.Dataset

Open input Dataset

bandNumberslist of int

GDAL band numbers of interest

Returns
imgNullValfloat or None

Null value from input raster, None if there is no null value

Raises
PyShepSegTilingError

If not all bands have the same null value

pyshepseg.tiling.getRatPageId(segId)

For the given segment ID, return the page ID. This is the segment ID of the first segment in the page.

pyshepseg.tiling.getTilesForFile(ds, tileSize, overlapSize)

Return a TileInfo object for a given file and input parameters.

Parameters
dsgdal.Dataset

Open GDAL Dataset object for raster to be tiles

tileSizeint

Size of tiles, in pixels. Individual tiles may end up being larger in either direction, when they meet the edge of the raster, to ensure we do not use very small tiles

overlapSizeint

Number of pixels by which tiles will overlap

Returns
tileInfoTileInfo

TileInfo object detailing the sizes and positions of all tiles across the raster

pyshepseg.tiling.overlapFilename(col, row, edge, tempDir)

Return the temporary filename used for the overlap array

Parameters
col, rowint

Tile column & row numbers

edge{right’, ‘bottom’}

Indicates from which edge of the given tile the overlap is taken

Returns
filenamestr

Temp numpy array filename for the overlap

pyshepseg.tiling.processSubsetTile(tile, recodeDict, histogramDict, maskData)

Process a tile of the subset area. Returns a new tile with the new codes. Fills in the recodeDict as it goes and also updates histogramDict.

Parameters
tileshepseg.SegIdType ndarray (tileNrows, tileNcols)

Input tile of segment IDs

recodeDictnumba.typed.Dict

Keyed by original segment ID, values are the corresponding segment IDs in the subset

histogramDictnumba.typed.Dict

Histogram counts in the subset, keyed by new segment ID

maskDataNone or int ndarray (tileNrows, tileNcols)

If not None, then is a raster mask. Only pixels which are non-zero in the mask will be included in the subset

Returns
outDatashepseg.SegIdType ndarray (tileNrows, tileNcols)

Recoded copy of the input tile.

pyshepseg.tiling.readColDataIntoPage(page, data, idx, colType, minVal)

Numba function to quickly read a column returned by rat.ReadAsArray() info a RatPage.

pyshepseg.tiling.readRATIntoPage(rat, numIntCols, numFloatCols, minVal, maxVal)

Create a new RatPage() object that represents the section of the RAT for a tile of an image. The part of the RAT between minVal and maxVal is read in and a RatPage() instance returned with the startSegId param set to minVal.

pyshepseg.tiling.readSubsampledImageBand(bandObj, subsampleProp)

Read in a sub-sampled copy of the whole of the given band.

Note that one can, in principle, do this directly using GDAL. However, if overview layers are present in the file, it will use these, and so is dependent on how these were created. Since these are often created just for display purposes, past experience has shown that they are not always to be trusted as data, so we have chosen to always go directly to the full resolution image.

Parameters
bandObjgdal.Band

An open Band object for input

subsamplePropfloat

The proportion by which to sub-sample (i.e. a value between zero and 1, applied to rows and columns separately)

Returns
imgSub<dtype> ndarray (nRowsSub, nColsSub)

A numpy array of the image subsample, equivalent to calling gdal.Band.ReadAsArray()

pyshepseg.tiling.recodeSharedSegments(tileData, overlapA, overlapB, orientation, recodeDict)

Work out a mapping which recodes segment ID numbers from the tile in tileData. Segments to be recoded are those which are in the overlap with an earlier tile, and which cross the midline of the overlap, which is where the stitchline between the tiles will fall.

Updates recodeDict, which is a dictionary keyed on the existing segment ID numbers, where the value of each entry is the segment ID number from the earlier tile, to be used to recode the segment in the current tile.

overlapA and overlapB are numpy arrays of pixels in the overlap region in question, giving the segment ID numbers in the two tiles. The values in overlapB are from the earlier tile, and those in overlapA are from the current tile.

It is critically important that the overlapping region is either at the top or the left of the current tile, as this means that the row and column numbers of pixels in the overlap arrays match the same pixels in the full tile. This cannot be used for overlaps on the right or bottom of the current tile.

Parameters
tileDatashepseg.SegIdType ndarray (tileNrows, tileNcols)

Tile subset of segment ID image

overlapA, overlapBshepseg.SegIdType ndarray (overlapNrows, overlapNcols)

Tile overlap subsets of segment ID image

orientation{HORIZONTAL, VERTICAL}

The orientation parameter defines whether we are dealing with overlap at the top (orientation == HORIZONTAL) or the left (orientation == VERTICAL).

recodeDictdict

Keys and values are both segment ID numbers. Defines the mapping which recodes segment IDs. Updated in place.

pyshepseg.tiling.recodeTile(tileData, maxSegId, tileRow, tileCol, overlapSize, tempDir, top, bottom, left, right)

Adjust the segment ID numbers in the current tile, to make them globally unique (and contiguous) across the whole mosaic.

Make use of the overlapping regions of tiles above and left, to identify shared segments, and recode those to segment IDs from the adjacent tiles (i.e. we change the current tile, not the adjacent ones). Non-shared segments are increased so they are larger than previous values.

Parameters
tileDatashepseg.SegIdType ndarray (tileNrows, tileNcols)

The array of segment IDs for a single image tile

maxSegIdshepseg.SegIdType

The current maximum segment ID for all preceding tiles.

tileRow, tileColint

The row/col numbers of this tile, within the whole-mosaic tile numbering scheme. (These are not pixel numbers, but tile grid numbers)

overlapSizeint

Number of pixels in tile overlap

tempDirstr

Name of directory for temporary files

top, bottom, left, rightint

Pixel coordinates within tile of the non-overlap region of the tile.

Returns
newTileDatashepseg.SegIdType ndarray (tileNrows, tileNcols)

A copy of tileData, with new segment ID numbers.

pyshepseg.tiling.relabelSegments(tileData, recodeDict, maxSegId, top, bottom, left, right)

Recode the segment IDs in the given tileData array.

For segment IDs which are keys in recodeDict, these are replaced with the corresponding entry. For all other segment IDs, they are replaced with sequentially increasing ID numbers, starting from one more than the previous maximum segment ID (maxSegId).

A re-coded copy of tileData is created, the original is unchanged.

Parameters
tileDatashepseg.SegIdType ndarray (tileNrows, tileNcols)

Segment IDs of tile

recodeDictdict

Keys and values are segment ID numbers. Defines mapping for segment relabelling

maxSegIdshepseg.SegIdType

Maximum segment ID number

top, bottom, left, rightint

Pixel coordinates within tile of the non-overlap region of the tile.

Returns
newTileDatashepseg.SegIdType ndarray (tileNrows, tileNcols)

Segment IDs of tile, after relabelling

newMaxSegIdshepseg.SegIdType

New maximum segment ID after relabelling

pyshepseg.tiling.setHistogramFromDictionary(dictn, histArray)

Given a dictionary of pixel counts keyed on index, write these values to the array.

pyshepseg.tiling.setSubsetRecodeFromDictionary(dictn, array)

Given the recodeDict write the original values to the array at the new indices.

pyshepseg.tiling.stitchTiles(inDs, outfile, tileFilenames, tileInfo, overlapSize, tempDir, simpleTileRecode, outputDriver, creationOptions, verbose)

Recombine individual tiles into a single segment raster output file. Segment ID values are recoded to be unique across the whole raster, and contiguous.

Parameters
inDsgdal.Dataset

Open Dataset of input raster

outfilestr

Filename of the final output raster

tileFilenamesdict

Dictionary of the individual tile filenames, keyed by a tuple of (col, row) defining which tile it is.

tileInfoTileInfo

Positions and sizes of all tiles across the raster. As returned by getTilesForFile().

overlapSizeint

The number of pixels in the overlap between tiles.

tempDirstr

Name of directory for temporary files

simpleTileRecodebool

If True, a simpler method will be used to recode segment IDs, using just a block offset to shift ID numbers. This is useful when testing, but leaves non-contiguous segment numbers. If False, then a more complicated method is used which recodes and merges segments which cross the boundary between tiles (this is the intended normal behaviour).

outputDriverstr

Short name string of the GDAL driver to use for the output file.

creationOptionslist of str

GDAL creation options for output driver

verbosebool

If True, print informative messages to stdout

Returns
maxSegIdshepseg.SegIdType

The maximum segment ID used.

pyshepseg.tiling.subsetImage(inname, outname, tlx, tly, newXsize, newYsize, outformat, creationOptions=[], origSegIdColName=None, maskImage=None)

Subset an image and “compress” the RAT so only values that are in the new image are in the RAT. Note that the image values will be recoded in the process.

gdal_translate seems to have a problem with files that have large RAT’s so while that gets fixed do the subsetting in this function.

Parameters
innamestr

Filename of input raster

outnamestr

Filename of output raster

tlx, tlyint

The x & y pixel coordinates (i.e. col & row, respectively) of the top left pixel of the image subset to extract.

newXsize, newYsizeint

The x & y size (in pixels) of the image subset to extract

outformatstr

The GDAL short name of the format driver to use for the output raster file

creationOptionslist of str

The GDAL creation options for the output file

origSegIdColNamestr or None

The name of a RAT column. If not None, this will be created in the output file, and will contain the original segment ID numbers so the new segment IDs can be linked back to the old

maskImagestr or None

If not None, then the filename of a mask raster. Only pixels which are non-zero in this mask image will be included in the subset. This image is assumed to match the shape and position of the output subset

pyshepseg.tiling.updateCounts(tileData, hist)

Fast function to increment counts for each segment ID in the given tile

pyshepseg.tiling.writeCompletedPagesForSubset(inRAT, outRAT, outPagedRat)

For the subset operation. Writes out any completed pages to outRAT using the inRAT as a template.

Parameters
inRAT, outRATgdal.RasterAttributeTable

The input and output raster attribute tables.

outPagedRatnumba.typed.Dict

The paged RAT in memory, as created by createPagedRat()