kealib.build_neighbours

  1"""
  2 *  build_neighbours.py
  3 *  LibKEA
  4 *
  5 *  Created by Sam Gillingham on 09/02/2023.
  6 *  Copyright 2012 LibKEA. All rights reserved.
  7 *
  8 *  This file is part of LibKEA.
  9 *
 10 *  Permission is hereby granted, free of charge, to any person 
 11 *  obtaining a copy of this software and associated documentation 
 12 *  files (the "Software"), to deal in the Software without restriction, 
 13 *  including without limitation the rights to use, copy, modify, 
 14 *  merge, publish, distribute, sublicense, and/or sell copies of the 
 15 *  Software, and to permit persons to whom the Software is furnished 
 16 *  to do so, subject to the following conditions:
 17 *
 18 *  The above copyright notice and this permission notice shall be 
 19 *  included in all copies or substantial portions of the Software.
 20 *
 21 *  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 
 22 *  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES 
 23 *  OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. 
 24 *  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR 
 25 *  ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF 
 26 *  CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION 
 27 *  WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 28 *
 29"""
 30
 31import os
 32import argparse
 33import numpy
 34from kealib import extrat
 35
 36DFLT_TILESIZE = 2048
 37
 38KEA_TO_NUMPY = {extrat.KEADataType.t8int: numpy.int8,
 39    extrat.KEADataType.t16int: numpy.int16,
 40    extrat.KEADataType.t32int: numpy.int32,
 41    extrat.KEADataType.t64int: numpy.int64,
 42    extrat.KEADataType.t8uint: numpy.uint8,
 43    extrat.KEADataType.t16uint: numpy.uint16,
 44    extrat.KEADataType.t32uint: numpy.uint32,
 45    extrat.KEADataType.t64uint: numpy.uint64,
 46    extrat.KEADataType.t32float: numpy.float32,
 47    extrat.KEADataType.t64float: float}
 48
 49def getCmdargs():
 50    """     
 51    Get the command line arguments.
 52    """
 53    p = argparse.ArgumentParser()
 54    p.add_argument("-i", "--infile", help=("Input Raster file. Neighbours will "
 55                + "be written back to this file. This file should have a Raster"
 56                + "Attribute Table"), required=True)
 57    p.add_argument("-t", "--tilesize", default=DFLT_TILESIZE,
 58        help="Size (in pixels) of tiles to chop input image into for processing."+
 59                " (default=%(default)s)", type=int)
 60    p.add_argument("--eightway", default=False, action="store_true",
 61        help="Use 8-way instead of 4-way")
 62    p.add_argument("--band", "-b", default=1, type=int, help="Band to use in "+
 63                "input dataset (1-based) (default=%(default)s)")
 64    cmdargs = p.parse_args()
 65    
 66    return cmdargs
 67    
 68def riosAccumulate(info, inputs, outputs, otherargs):
 69    """
 70    Calls NeighbourAccumulator.addArray() with each tile
 71    """
 72    # convert to 0 baed for numpy
 73    data = inputs.input[otherargs.band - 1]
 74    otherargs.accumulator.addArray(data)
 75
 76def main():
 77    """
 78    Main routine. For calling from the command line 
 79    """
 80    cmdargs = getCmdargs()
 81
 82    # we need the dataset for NeighbourAccumulator
 83    # Open using GDAL here, test routines can supply a fake
 84    # GDAL dataset.
 85    # only import here so we don't need GDAL for testing
 86    from osgeo import gdal
 87    gdal.UseExceptions()
 88    ds = gdal.Open(cmdargs.infile, gdal.GA_Update)
 89    
 90    buildNeighbours(ds, cmdargs.band, cmdargs.tilesize, cmdargs.eightway)
 91    
 92def buildNeighbours(ds, band, tilesize=DFLT_TILESIZE, eightway=False):
 93    """
 94    Does the actual work or building the neighbours. 
 95    
 96    ds should be a GDAL dataset object
 97    """
 98    fieldInfo = extrat.getFieldByName(ds, band, "Histogram")
 99    size = extrat.getSize(ds, band)
100    histogram = extrat.getField(ds, band, fieldInfo, 0, size)
101    
102    # will through exception if not set, but we need a nodata otherwise
103    # we can't do a read with margin
104    nodata = extrat.getNoDataValue(ds, band)
105    
106    keaDataType = extrat.getImageBandDataType(ds, band)
107    dataType = KEA_TO_NUMPY[keaDataType]
108
109    imageInfo = extrat.getSpatialInfo(ds)
110
111    accumulator = extrat.NeighbourAccumulator(histogram, ds, 
112                        band, not eightway)    
113    
114    numXtiles = int(numpy.ceil(imageInfo.xSize / tilesize))
115    numYtiles = int(numpy.ceil(imageInfo.ySize / tilesize))
116
117    for tileRow in range(numYtiles):
118        ypos = tileRow * tilesize
119        ysize = min(tilesize, (imageInfo.ySize - ypos))
120        for tileCol in range(numXtiles):
121            xpos = tileCol * tilesize
122            xsize = min(tilesize, (imageInfo.xSize - xpos))
123
124            data = readBlockWithMargin(ds, band, xpos, ypos, 
125                        xsize, ysize, 1, dataType, nodata, 
126                        imageInfo.xSize, imageInfo.ySize)
127            accumulator.addArray(data)
128                        
129    
130def readBlockWithMargin(ds, band, xoff, yoff, xsize, ysize, margin, datatype, 
131        nodata, RasterXSize, RasterYSize):
132        """
133        Stolen from RIOS
134        
135        A 'drop-in' look-alike for the ReadAsArray function in GDAL,
136        but with the option of specifying a margin width, such that
137        the block actually read and returned will be larger by that many pixels. 
138        The returned array will ALWAYS contain these extra rows/cols, and 
139        if they do not exist in the file (e.g. because the margin would push off 
140        the edge of the file) then they will be filled with the given nullVal. 
141        Otherwise they will be read from the file along with the rest of the block. 
142        
143        Variables within this function which have _margin as suffix are intended to 
144        designate variables which include the margin, as opposed to those without. 
145        
146        This routine will cope with any specified region, even if it is entirely outside
147        the given raster. The returned block would, in that case, be filled
148        entirely with the null value. 
149        
150        """
151        
152        # Create the final array, with margin, but filled with the null value. 
153        xSize_margin = xsize + 2 * margin
154        ySize_margin = ysize + 2 * margin
155        outBlockShape = (ySize_margin, xSize_margin)
156        
157        # Create the empty output array, filled with the appropriate null value. 
158        block_margin = numpy.full(outBlockShape, nodata, dtype=datatype)
159
160        # Calculate the bounds of the block which we will actually read from the file,
161        # based on what we have been asked for, what margin size, and how close we
162        # are to the edge of the file. 
163        
164        # The bounds of the whole image in the file
165        imgLeftBound = 0
166        imgTopBound = 0
167        imgRightBound = RasterXSize
168        imgBottomBound = RasterYSize
169        
170        # The region we will, in principle, read from the file. Note that xSize_margin 
171        # and ySize_margin are already calculated above
172        xoff_margin = xoff - margin
173        yoff_margin = yoff - margin
174        
175        # Restrict this to what is available in the file
176        xoff_margin_file = max(xoff_margin, imgLeftBound)
177        xoff_margin_file = min(xoff_margin_file, imgRightBound)
178        xright_margin_file = xoff_margin + xSize_margin
179        xright_margin_file = min(xright_margin_file, imgRightBound)
180        xSize_margin_file = xright_margin_file - xoff_margin_file
181
182        yoff_margin_file = max(yoff_margin, imgTopBound)
183        yoff_margin_file = min(yoff_margin_file, imgBottomBound)
184        ySize_margin_file = min(ySize_margin, imgBottomBound - yoff_margin_file)
185        ybottom_margin_file = yoff_margin + ySize_margin
186        ybottom_margin_file = min(ybottom_margin_file, imgBottomBound)
187        ySize_margin_file = ybottom_margin_file - yoff_margin_file
188        
189        # How many pixels on each edge of the block we end up NOT reading from 
190        # the file, and thus have to leave as null in the array
191        notRead_left = xoff_margin_file - xoff_margin
192        notRead_right = xSize_margin - (notRead_left + xSize_margin_file)
193        notRead_top = yoff_margin_file - yoff_margin
194        notRead_bottom = ySize_margin - (notRead_top + ySize_margin_file)
195        
196        # The upper bounds on the slices specified to receive the data
197        slice_right = xSize_margin - notRead_right
198        slice_bottom = ySize_margin - notRead_bottom
199        
200        if xSize_margin_file > 0 and ySize_margin_file > 0:
201            # Now read in the part of the array which we can actually read from the file.
202            # Read each layer separately, to honour the layerselection
203            
204            # The part of the final array we are filling
205            imageSlice = (slice(notRead_top, slice_bottom), slice(notRead_left, slice_right))
206            
207            block_margin[imageSlice] = extrat.getImageBlock(ds, band, xoff_margin_file, 
208                    yoff_margin_file, xSize_margin_file, ySize_margin_file) 
209            
210        return block_margin
211
212