"""
This module contains the ReaderInfo class
which holds information about the area being
read and info on the current block
"""
# This file is part of RIOS - Raster I/O Simplification
# Copyright (C) 2012 Sam Gillingham, Neil Flood
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
import math
import numpy
from osgeo import gdal_array
from . import imageio
[docs]class ReaderInfo(object):
"""
ReaderInfo class. Holds information about the area being
read and info on the current block
"""
def __init__(self, workingGrid,
windowxsize, windowysize, overlap, loggingstream):
self.loggingstream = loggingstream
# grab the working grid
self.workingGrid = workingGrid
# save the window size and overlap
self.windowxsize = windowxsize
self.windowysize = windowysize
self.overlap = overlap
# work out the area being read
self.xsize = int(round((self.workingGrid.xMax -
self.workingGrid.xMin) / self.workingGrid.xRes))
self.ysize = int(round((self.workingGrid.yMax -
self.workingGrid.yMin) / self.workingGrid.yRes))
# total number of blocks
self.xtotalblocks = int(math.ceil(float(self.xsize) / self.windowxsize))
self.ytotalblocks = int(math.ceil(float(self.ysize) / self.windowysize))
# The feilds below apply to a particular block
# and are filled in after this object is copied
# to make it specific fir each block
self.blockwidth = None
self.blockheight = None
self.blocktl = None
self.blockbr = None
self.xblock = None
self.yblock = None
# dictionary keyed by id() of the number array
# value is a tuple with the GDAL dataset object
# that corresponds to it, and the original filename
self.blocklookup = {}
[docs] def setBlockDataset(self, block, dataset, filename):
"""
Saves a match between the numpy block read
and it's GDAL dataset. So we can look up the
dataset later given a block.
This routine is for internal use by RIOS. Its use in any other
context is not sensible.
"""
self.blocklookup[id(block)] = (dataset, filename)
[docs] def getWindowSize(self):
"""
Returns the size of the current window. Returns a
tuple (numCols, numRows)
"""
return (self.windowxsize, self.windowysize)
[docs] def getOverlapSize(self):
"""
Returns the size of the pixel overlap between
each window. This is the number of pixels added as
margin around each block
"""
return self.overlap
[docs] def getTotalSize(self):
"""
Returns the total size (in pixels) of the dataset
being processed
"""
return (self.xsize, self.ysize)
[docs] def getProjection(self):
"""
Return the WKT describing the current
projection system
"""
return self.workingGrid.projection
[docs] def getTotalBlocks(self):
"""
Returns the total number of blocks the dataset
has been split up into for processing
"""
return (self.xtotalblocks, self.ytotalblocks)
[docs] def setBlockSize(self, blockwidth, blockheight):
"""
Sets the size of the current block
This routine is for internal use by RIOS. Its use in any other
context is not sensible.
"""
self.blockwidth = blockwidth
self.blockheight = blockheight
[docs] def getBlockSize(self):
"""
Get the size of the current block. Returns a tuple::
(numCols, numRows)
for the current block. Mostly the same as the window size,
except on the edge of the raster.
"""
return (self.blockwidth, self.blockheight)
[docs] def setBlockBounds(self, blocktl, blockbr):
"""
Sets the coordinate bounds of the current block
This routine is for internal use by RIOS. Its use in any other
context is not sensible.
"""
self.blocktl = blocktl
self.blockbr = blockbr
[docs] def getBlockCoordArrays(self):
"""
Return a tuple of the world coordinates for every pixel
in the current block. Each array has the same shape as the
current block. Return value is a tuple::
(xBlock, yBlock)
where the values in xBlock are the X coordinates of the centre
of each pixel, and similarly for yBlock.
The coordinates returned are for the pixel centres. This is
slightly inconsistent with usual GDAL usage, but more likely to
be what one wants.
"""
(tl, _) = (self.blocktl, self.blockbr)
(nCols, nRows) = self.getBlockSize()
nCols += 2 * self.overlap
nRows += 2 * self.overlap
(xRes, yRes) = self.getPixelSize()
(rowNdx, colNdx) = numpy.mgrid[0:nRows, 0:nCols]
xBlock = tl.x - self.overlap * xRes + xRes / 2.0 + colNdx * xRes
yBlock = tl.y + self.overlap * yRes - yRes / 2.0 - rowNdx * yRes
return (xBlock, yBlock)
[docs] def setBlockCount(self, xblock, yblock):
"""
Sets the count of the current block
This routine is for internal use by RIOS. Its use in any other
context is not sensible.
"""
self.xblock = xblock
self.yblock = yblock
[docs] def getBlockCount(self):
"""
Gets the count of the current block
"""
return (self.xblock, self.yblock)
[docs] def getPixelSize(self):
"""
Gets the current pixel size and returns it as a tuple (x and y)
"""
return (self.workingGrid.xRes, self.workingGrid.yRes)
[docs] def getPixRowColBlock(self, x, y):
"""
Return the row/column numbers, within the current block,
for the pixel which contains the given (x, y) coordinate.
The coordinates of (x, y) are in the world coordinate
system of the reference grid. The row/col numbers are
suitable to use as array indices in the array(s) for the
current block. If the nominated pixel is not contained
within the current block, the row and column numbers are
both None (hence this should be checked).
Return value is a tuple of 2 int values
(row, col)
"""
transform = self.workingGrid.makeGeoTransform()
imgRowCol = imageio.wld2pix(transform, x, y)
imgRow = imgRowCol.y
imgCol = imgRowCol.x
blockStartRow = self.yblock * self.windowysize - self.overlap
blockStartCol = self.xblock * self.windowxsize - self.overlap
blockRow = int(imgRow - blockStartRow)
blockCol = int(imgCol - blockStartCol)
if ((blockRow < 0 or blockRow > (self.windowysize + 2 * self.overlap)) or
(blockCol < 0 or blockCol > (self.windowxsize + 2 * self.overlap))):
blockRow = None
blockCol = None
return (blockRow, blockCol)
[docs] def getPixColRow(self, x, y):
"""
This function is for internal use only. The user should
be looking at getBlockCoordArrays() or getPixRowColBlock()
for dealing with blocks and coordinates.
Get the (col, row) relative to the current image grid,
for the nominated pixel within the current block. The
given (x, y) are column/row numbers (starting at zero),
and the return is a tuple::
(column, row)
where these are relative to the whole of the current
working grid. If working with a single raster, this is the same
as for that raster, but if working with multiple rasters,
the working grid is the intersection or union of them.
Note that this function will give incorrect/misleading results
if used in conjunction with a block overlap.
"""
col = self.xblock * self.windowxsize + x
row = self.yblock * self.windowysize + y
return (col, row)
[docs] def isFirstBlock(self):
"""
Returns True if this is the first block to be processed
"""
return self.xblock == 0 and self.yblock == 0
[docs] def isLastBlock(self):
"""
Returns True if this is the last block to be processed
"""
xtotalblocksminus1 = self.xtotalblocks - 1
ytotalblocksminus1 = self.ytotalblocks - 1
return self.xblock == xtotalblocksminus1 and self.yblock == ytotalblocksminus1
[docs] def getFilenameFor(self, block):
"""
Get the input filename of a dataset
"""
# can't use ds.GetDescription() as may have been resampled
(ds, fname) = self.blocklookup[id(block)]
return fname
[docs] def getGDALDatasetFor(self, block):
"""
Get the underlying GDAL handle of a dataset
"""
(ds, fname) = self.blocklookup[id(block)]
return ds
[docs] def getGDALBandFor(self, block, band):
"""
Get the underlying GDAL handle for a band of a dataset
"""
ds = self.getGDALDatasetFor(block)
return ds.GetRasterBand(band)
[docs] def getNoDataValueFor(self, block, band=1):
"""
Returns the 'no data' value for the dataset
underlying the block. This should be the
same as what was set for the stats ignore value
when that dataset was created.
The value is cast to the same data type as the
dataset.
"""
ds = self.getGDALDatasetFor(block)
band = ds.GetRasterBand(band)
novalue = band.GetNoDataValue()
# if there is a valid novalue, cast it to the type
# of the dataset. Note this creates a numpy 0-d array
if novalue is not None:
numpytype = gdal_array.GDALTypeCodeToNumericTypeCode(band.DataType)
novalue = numpy.cast[numpytype](novalue)
return novalue
[docs] def getPercent(self):
"""
Returns the percent complete.
"""
percent = int(float(self.yblock * self.xtotalblocks + self.xblock) /
float(self.xtotalblocks * self.ytotalblocks) * 100)
return percent