"""
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 . import imageio
from . import rat
[docs]class StatisticsCache(object):
"""
Allows global statistics for all the files used to be cached
Statistics are stored here associated with the filename
so if there are multiple calls to global_stats for one
file the statistics will only get calculated once.
"""
def __init__(self):
self.stats = {}
[docs] def getStats(self,fname,band,ignore):
"""
Returns tuple with (min,max,mean,stddev) if
previously cached, None otherwise.
"""
if ignore is None:
key = '%s %d' % (fname,band)
else:
key = '%s %d %f' % (fname,band,ignore)
if key in self.stats:
return self.stats[key]
else:
return None
[docs] def setStats(self,fname,band,ignore,stats):
"""
Sets tuple with (min,max,mean,stddev) in cache
"""
if ignore is None:
key = '%s %d' % (fname,band)
else:
key = '%s %d %f' % (fname,band,ignore)
self.stats[key] = stats
[docs]class AttributeTableCache(object):
"""
Allows attribute columns to be cached. This means
that an attribute column can be requested for each
block through the image, but the actual column data
will only be read the first time and cached for
subsequent calls.
"""
def __init__(self):
self.rats = {}
[docs] def getColumn(self, fname, band, colName):
"""
Returns the column if previously cached, otherwise
None.
"""
key = '%s %d %s' % (fname, band, colName)
if key in self.rats:
return self.rats[key]
else:
return None
[docs] def setColumn(self, fname, band, colName, column):
"""
Puts the column into the cache
"""
key = '%s %d %s' % (fname, band, colName)
self.rats[key] = column
[docs]class ReaderInfo(object):
"""
ReaderInfo class. Holds information about the area being
read and info on the current block
"""
def __init__(self, workingGrid, statscache, ratcache,
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))
# save the statistics cache.
self.statscache = statscache
# save the RAT cache
self.ratcache = ratcache
# 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, br) = (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 = imageio.GDALTypeToNumpyType(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
[docs] def getAttributeColumn(self, block, colName, band=1):
"""
Gets the attribute for the given block and column name
Caches columns so only first call actually extracts data
"""
(ds,fname) = self.blocklookup[id(block)]
column = self.ratcache.getColumn(fname, band, colName)
if column is None:
column = rat.readColumn(ds, colName, band)
self.ratcache.setColumn(fname, band, colName, column)
return column
[docs] def global_stats(self,block,band=1,ignore=None):
"""
Returns the (min,max,mean,stddev) for the whole band
"""
fname = self.getFilenameFor(block)
# see if we have the stats in our cache
values = self.statscache.getStats(fname,band,ignore)
if values is None:
# no, get the gdal band
bandhandle = self.getGDALBandFor(block,band)
# set ignore value if specified so that
# GDAL ignores it when calculating stats
if ignore is not None:
bandhandle.SetNoDataValue(ignore)
self.loggingstream.write("Calculating global statistics...\n")
# get GDAL to calc the stats
values = bandhandle.GetStatistics(False,True)
# set it back in our cache for next time
self.statscache.setStats(fname,band,ignore,values)
return values
[docs] def global_min(self,block,band=1,ignore=None):
"""
Returns the min for the whole band
"""
return self.global_stats(block,band,ignore)[0]
[docs] def global_max(self,block,band=1,ignore=None):
"""
Returns the max for the whole band
"""
return self.global_stats(block,band,ignore)[1]
[docs] def global_mean(self,block,band=1,ignore=None):
"""
Returns the mean for the whole band
"""
return self.global_stats(block,band,ignore)[2]
[docs] def global_stddev(self,block,band=1,ignore=None):
"""
Returns the stddev for the whole band
"""
return self.global_stats(block,band,ignore)[3]