On 02/03/2010 01:48 AM, Ivan wrote:
Even,
That was a little thing that I start years ago. Here is the usage:
{{{
#
=============================================================================
# Usage Message
#
=============================================================================
def Usage():
print 'Usage: mathbox.py [-formula "formula"] [outfile] [infile*]'
print
print ' -formula Formula to be applied on input pixels values'
print ' assigning the results to the output.'
print ' Use letters "a" to "j" as variables to refers'
print ' to the list of infiles in alphabetic order'
print ' outfifile Name of the output file'
print ' infile Name of the input file(s)'
print
print ' Notes: Images will be identified by the order of file
name as "a,b,c...'
print
print 'Examples:'
print
print ' mathbox.py -formula "a+b/2" result.tif dist_roads.tif
slope.tif'
print ' mathbox.py -formula "a-255" correct.tif original.tif'
print
sys.exit(1)
}}}
It is based on the "eval()" Python function, so it will actually
execute the formula as it was
entered by the user. Please let me know if that would be of any interest.
Regards,
Ivan
Some years ago I've created something similar, I still use it from time
to time. See attached python script (if that survives the list. we'll
see...). As far as I'm concerned the code is public (it's not that fancy
:)). I'm not willing to maintain this actively as part of gdal, but feel
free to adapt/adopt/use/mix'n'match it.
O by the way, it does more or less the same as Ivan's (apply a formula
or formula's on pixels on input bands/sources), but is pretty flexible
as you have all of numpy & python available in your formula. Furthermore
it evaluates per line, so no neighbourhood functions possible, but also
no memory limits :-)
Vincent.
Even Rouault wrote:
For the record, I had started some time ago to work on something
pretty close to Frank's ideas. It was based on a "generic" expression
evaluator that could accept C-like expressions ( arithmetic
operators, boolean operators, numeric constants, a few predefined
maths functions, user variables, user functions, ...). It could be
specialized for pixel operations : the specialization consists in
providing the necessary code to evaluate the variables and function
calls. I'll give more details once I find where it is burried right
now ;-).
It could potentially do more complex operations than just operations
on pixels at same location in different images.
If my memory are corrects, it could do things like :
* "pixel[ysize-1-y][x]" : to make an horizontal flip of an image
* "0.30 * source[0].pixel[y][x] + 0.59 * source[1].pixel[y][x] + 0.11
* source[2].pixel[y][x]" : to compute a grey level from RGB
* "sum(j,-1,1,sum(i,-1,1,abs(pixel[y+j][x+i] - pixel[y][x]))) / 8" :
equivalent of the TPI algorithm of gdaldem
* "abs(source[0].pixel[y][x] - source[1].pixel[y][x])" : for the
example discussed in the previous post in this thread ;-)
This worked pretty well but of course the performance of an
interpretator is way much slower than the equivalent compiled code.
Best regards,
Even
Le Tuesday 02 February 2010 17:13:54 Frank Warmerdam, vous avez écrit :
Antonio Valentino wrote:
Anyway, IMHO the point here is not to decide if we have to
*modify* the diff function so that it returns abs(b1-b2).
Greg asked to add a new one ( abs(b1-b2) ) to the *base set* of pixel
functions.
Stated that there is no technical difficulty to do this I'm wondering
what is the Frank's mind about the question.
Antonio,
I think including abs(b1-b2) just demonstrates that we haven't come
up with a convenient way of stacking more basic operations (difference
and abs).
Honestly, I'm a bit ambivalent about fleshing out the use of
pixel operations in VRT files and I've never used them myself.
I think I would prefer to wait for a more generalized expression
evaluation mechanism for VRT or possibly as a distinct driver.
However, if there are other GDAL developers keen on this functionality,
and willing to take responsibility for integration, documentation and
support then I'm willing to let them proceed.
Best regards,
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
#!/usr/bin/python
#******************************************************************************
# Copyright (c) 2004, Vincent T. Schut
#
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included
# in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
# DEALINGS IN THE SOFTWARE.
#******************************************************************************
import sys
import traceback
from osgeo import gdal, gdal_array
#from gdalnumeric import *
#from gdalconst import *
#from numpy import *
import numpy as N
#from progressbar import progressBar
def Usage():
print ' Usage: gdal_calc.py -i inputfile{,inputfile2,...] -o outputfile'
print ' -m math[;math2;...] -p outputprototype'
print ' -of outputformat -ot outputtype -ct calctype'
print
print ' For more info, type \'gdal_calc.py --help\''
def Help():
helptext = """
gdal_calc.py: applies custom math to one or more raster input files
and stores the output as a raster file. All file handling is done by
GDAL, which means that gdal_calc.py supports the same formats as your
GDAL implementation.
Explanation of options:
-i One or more input files, separated by commas.
-o Exactly one output file.
-p Prototype for the output file. Geo-information will be copied
from this file into your output file. If omitted, the first
input file will be used as prototype.
-of Output format. See the GDAL formats. The default is GeoTiff.
-ot Output type: a valid Numeric (NumPy) array typecode. Default
is 'f' (Float). After all calculations have been done, the result
will be cast to this type.
-ct Calculation type: a valid Numeric (NumPy) array typecode. All
calculations will be done using this array type.
-m One or more math expressions, separated by semicolons. Numpy is available as 'N'.
Two types of math expressions exist:
- assignments to a custom (temporary) variable. If you want to
store a calculation to use the result in a later math definition,
use this type. Syntax: var = [math], where you can choose any
name for var (as long as it is not a Python reserved word or
happens to be a variable that is used in gdal_calc...)
Example: -m "NDVI = (band2 - band1) / (band2 + band1)"
- assignments to one or more output bands. For any math line that
appears to be not an assignment, the result will be assigned to
an output band (or more if appropriate).
Example: -m "N.sin(band1)"
In the math part, there are three ways to refer to input rasters:
- band#
# is a valid input band number. Bands start to count at 1 and
counting continues through the list of input files, e.g.if you
have a 3-band file and a 2-band file, the bands of the first file
will be band1, band2, band3 and the bands of the second file will
be band4 and band5.
- file#band#
like band#, but starts counting at band1 for each input file. E.g.
band 1 of the first input file is file1band1, band 5 of the 3rd
input file is file3band5. Files start to count at 1, like bands.
- file#
maps all bands of a input file. Note that this enables you to use
NumPy functions like N.sum(file1) which will give you one output
band with all bands pixelwise summed. More tricks are easy to invent.
The math is only limited by Python and Numpy, as the (slightly
changed) math expressions are directly executed by the Python
interpreter.
For more information on Numeric typecodes and functions, see
http://www.pfdubois.com/numpy/html2/numpy.html
"""
print helptext
#band = lambda n: gdalBands[n-1].ReadAsArray(tileXStart, tileYStart, tileXSize, tileYSize)
#file = lambda n: inFile[n-1].ReadAsArray(tileXStart, tileYStart, tileXSize, tileYSize)
def doCalc():
band = lambda n: gdalBands[n-1].ReadAsArray(tileXStart, tileYStart, tileXSize, 1).astype(calcType)
file = lambda n: inFile[n-1].ReadAsArray(tileXStart, tileYStart, tileXSize, 1).astype(calcType)
outPrototype = None
inFiles = None
outFile = 'gdal_calc_out.tif'
formulae = None
outFormat = 'GTiff'
outOptions = ["TILED=YES"]
outType = 'float32'
calcType = 'float32'
args = sys.argv
argNr = 1
while argNr < len(args):
option = args[argNr]
if option == '--help':
Help()
sys.exit()
if option == '-i':
inFiles = args[argNr+1].split(',')
for inFilename in inFiles:
inFiles[inFiles.index(inFilename)] = inFilename.strip()
argNr += 2
elif option == '-o':
outFile = args[argNr+1].strip()
argNr += 2
elif option == '-p':
outPrototype = args[argNr+1].strip()
argNr += 2
elif option == '-m':
formulae = args[argNr+1].split(';')
argNr += 2
elif option == '-of':
outFormat = args[argNr+1].strip()
argNr += 2
elif option == '-ot':
outType = args[argNr+1].strip()
argNr += 2
elif option == '-ct':
calcType = args[argNr+1].strip()
argNr += 2
else:
print 'Unknow option: %s' % option
Usage()
sys.exit()
if inFiles == None:
print 'You need to specify at least one input file.'
Usage()
sys.exit()
if outFile == None:
print 'You need to specify an output file.'
Usage()
sys.exit()
if formulae == None:
print 'You need to specify at least one math directive.'
Usage()
sys.exit()
if outPrototype == None and outPrototype != 'none':
# use first input file as prototype:
outPrototype = inFiles[0]
if outFormat == None:
outFormat = 'GTiff'
# open input files.
inFile = []
xSize = None
ySize = None
numBands = 0
gdalBands = []
gdalFilebands = []
print ''
for inFilename in inFiles:
inFile.append(gdal.Open(inFilename))
if inFilename == inFiles[0]:
xSize = inFile[0].RasterXSize
ySize = inFile[0].RasterYSize
else:
if xSize != inFile[-1].RasterXSize or ySize != inFile[-1].RasterYSize:
print 'WARNING: Raster size of %s does not match other files.' % inFilename
# store band and dataset objects in lists:
for b in range(inFile[-1].RasterCount):
globalBandNr = numBands + b + 1
fileNr = len(inFile)
fileBandNr = b + 1
gdalBands.append(inFile[-1].GetRasterBand(b+1))
gdalFilebands.append(gdalBands[numBands:])
numBands += inFile[-1].RasterCount
# print mapping details:
if inFile[-1].RasterCount == 1:
print 'band%d : %s' % (numBands, inFilename)
elif inFile[-1].RasterCount == 2:
print 'band%d, band%d : %s' % (numBands-1, numBands, inFilename)
else:
print 'band%d - band%d : %s' % (numBands - inFile[-1].RasterCount + 1, numBands, inFilename)
pyFormulae = []
outBandsTotal = 0
for formula in formulae:
# rewrite band#, file# and file#band# directives in math to map to band(x) function:
newFormula = formula
globalBandNr = 0
for f in range(len(gdalFilebands))[::-1]:
for b in range(len(gdalFilebands[f]))[::-1]:
newFormula = newFormula.replace('file%dband%d' % (f+1, b+1), 'band(%d)' % globalBandNr)
globalBandNr += 1
for b in range(numBands)[::-1]:
newFormula = newFormula.replace('band%d' % (b+1), 'band(%d)' % (b+1))
for f in range(len(inFile))[::-1]:
newFormula = newFormula.replace('file%d' % (f+1), 'file(%d)' % (f+1))
# exec doesn't like nulls...
newFormula = newFormula.strip()
# prepare for testing on first tile to determine the character of the math lines:
tileXStart=0
tileYStart=0
tileXSize=xSize # currently we do processing line by line
tileYSize=1
# do the testing:
try:
# test if executing the math gives us an output band or range of output bands:
#checkResult = None
checkFormula = ('checkResult = (%s)' % (newFormula)).strip()
#checkFormula = 'checkResult = %s.astype("%s")' % (newFormula, calcType)
#exec(checkFormula) in locals()
exec(checkFormula) in globals(), locals()
#print checkResult.shape
#print checkResult.typecode()
if len(checkResult.shape) == 2 and checkResult.shape[0] == 1 and checkResult.shape[1] == tileXSize:
outBands = 1
outBandsTotal += 1
elif len(checkResult.shape) == 3 and checkResult.shape[1] == 1 and checkResult.shape[2] == tileXSize:
outBands = checkResult.shape[0]
outBandsTotal += outBands
else:
print 'Hmm... One of the maths results in and array with the following shape:', tmp.shape
print 'which does not correspond to a single or multiband array with the same'
print 'dimensions as the input files...'
print 'The offending math is:'
print formula
print 'Aborting.'
sys.exit()
#pyFormulae.append('outBands[%d:%d, tileYStart:tileYStart+tileYSize, tileXStart:tileXStart+tileXSize] = (%s).astype("%s")' % ((outBandsTotal-outBands), outBandsTotal, newFormula, calcType))
#pyFormulae.append('outBands[%d:%d, 0, tileXStart:tileXStart+tileXSize] = (%s).astype("%s")' % ((outBandsTotal-outBands), outBandsTotal, newFormula, calcType))
pyFormulae.append('outBands[%d:%d, ::] = (%s).astype("%s")' % ((outBandsTotal-outBands), outBandsTotal, newFormula, calcType))
except:
# test if formula is an assignment to a user defined var:
try:
#exec(newFormula) in locals()
exec(newFormula) in globals(), locals()
pyFormulae.append(newFormula)
except:
print 'An error occured while testing the following math line:'
print formula
print newFormula
print 'Python error report:'
traceback.print_exc()
print
print 'Aborting.'
sys.exit()
if outBandsTotal == 0:
print 'Warning: your math will not result in any output being written.'
else:
print '%d output bands will be produced.' % outBandsTotal
#outBands = zeros((outBandsTotal, ySize, xSize), calcType)
# open output file and create dataset:
driver = gdal.GetDriverByName(outFormat)
outDataset = driver.Create(outFile, xSize, ySize, outBandsTotal, gdal_array.flip_code(N.typeDict[outType]), outOptions)
# copy geo attributes from first input file
#CopyDatasetInfo(inFile[0], outDataset)
outDataset.SetGeoTransform(inFile[0].GetGeoTransform())
outDataset.SetProjection(inFile[0].GetProjection())
#outDataset.FlushCache()
# allocate output row array
outBands = N.zeros((outBandsTotal, 1, xSize), calcType)
#pbar = progressBar(0, ySize, 52)
pbarLength = 50
pbar = str(pbarLength * ".")
#oldpbar = str(pbar)
oldFloatPercentage = 0.0
print "%s %5.1f%%" % (pbar, 0.0),
for tileYStart in range(ySize):
for f in range(len(pyFormulae)):
#print pyFormulae[f]
#exec(pyFormulae[f]) in locals()
exec(pyFormulae[f]) in globals(), locals()
outBandsOutType = outBands.astype(outType)
outDataset.WriteRaster(0, tileYStart, xSize, 1, outBandsOutType.tostring());
#pbar.updateAmount(tileYStart)
floatPercentage = 100.0 * (float(tileYStart) / float(ySize-1))
#percentage = int(100.0 * (float(tileYStart) / float(ySize-1)))
pbarPosition = int(floatPercentage * (pbarLength / 100.0))
pbar = pbarPosition * "=" + (pbarLength - pbarPosition) * "-"
if floatPercentage - oldFloatPercentage >= 0.1:
print "\r%s %5.1f%%" % (pbar, floatPercentage),
sys.stdout.flush()
#oldpbar = str(pbar)
oldFloatPercentage = floatPercentage
#if outBandsTotal > 0:
# SaveArray(outBands.astype(outType), outFile, outFormat, outPrototype)
if __name__ == '__main__':
doCalc()
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
