Just fyi: I've been using my own brewn version of a band-math like
command line gdal python script for some years now (have sent it to the
gdal list once, like a year ago or something), you might be interested
to take a look at the source and compare approaches. I've attached the
script with this email. Afaic it's public domain, do with it what you
want :-)
Cheers,
Vincent.
On 09/17/2010 07:38 PM, Ivan Lucena wrote:
Hi Chris,
That sounds good.
Is that by any chance based on this idea here:
[http://trac.osgeo.org/gdal/wiki/SummerOfCode]?
Regards,
Ivan
-------Original Message-------
From: Chris Yesson<[email protected]>
To: [email protected]
Subject: [gdal-dev] gdal_calc.py
Sent: Sep 17 '10 09:29
Hi all,
I have been playing around with a python/gdalnumeric-based raster
calculator for a while now (http://static.zsl.org/files/gdal-calc-1218.py).
I thought it was time to try to share this properly via GDAL.
Could this go on the sample python scripts directory?
Could I get advice on what I would need to do to get it into the main
scripts?
Many thanks
- Chris
_Dr Chris Yesson_
_Postdoctoral Researcher_
_Institute__ of Zoology_
_Zoological Society of London_
_Regent's Park_
_London__ NW1 4RY_
_020 7449 6267_
[LINK: http://www.zsl.org/chrisyesson/] _http://www.zsl.org/chrisyesson/_
[LINK: http://eu-fp7-coralfish.net/] _http://eu-fp7-coralfish.net/_
The Zoological Society of London is incorporated by Royal Charter
Principal Office England. Company Number RC000749
Registered address:
Regent's Park, London, England NW1 4RY
Registered Charity in England and Wales no. 208728
_________________________________________________________________________
This e-mail has been sent in confidence to the named addressee(s).
If you are not the intended recipient, you must not disclose or distribute
it in any form, and you are asked to contact the sender immediately.
Views or opinions expressed in this communication may not be those
of The Zoological Society of London and, therefore, The Zoological
Society of London does not accept legal responsibility for the contents
of this message. The recipient(s) must be aware that e-mail is not a
secure communication medium and that the contents of this mail may
have been altered by a third party in transit.
If you have any issues regarding this mail please contact:
[LINK: mailto:[email protected]] [email protected].
___________________________________________________________________________
This message has been scanned for viruses by [LINK:
http://www.mailcontrol.com/] MailControl, a service from [LINK:
http://www.blackspider.com/] BlackSpider Technologies.
--------------------
_______________________________________________
gdal-dev mailing list
[LINK: [email protected]] [email protected]
[LINK: http://lists.osgeo.org/mailman/listinfo/gdal-dev]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
_______________________________________________
gdal-dev mailing list
[email protected]
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
import scipy
import scipy.stats
#from progressbar import progressBar
N.seterr(all='ignore')
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.
-y Number of lines to process per iteration (defaults to input file blocksize)
-nd NoData value to set in metadata of output file
-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, tileHeight).astype(calcType)
file = lambda n: inFile[n-1].ReadAsArray(tileXStart, tileYStart, tileXSize, tileHeight).astype(calcType)
outPrototype = None
inFiles = None
outFile = 'gdal_calc_out.tif'
formulae = None
outFormat = 'GTiff'
outOptions = ['INTERLEAVE=BAND']
outType = 'float32'
calcType = 'float32'
tileYSize = 0
ndv = None
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
elif option == '-y':
tileYSize = int(args[argNr+1].strip())
argNr += 2
elif option == '-nd':
ndv = float(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
blockSize = inFile[0].GetRasterBand(1).GetBlockSize()
print blockSize
if blockSize[0] != xSize:
outOptions.append('TILED=YES')
outOptions.append('BLOCKXSIZE=%i' % blockSize[0])
outOptions.append('BLOCKYSIZE=%i' % blockSize[1])
if tileYSize == 0:
tileYSize = blockSize[1]
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 = []
pyFormulaeBands = []
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()
tileHeight = 1
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('result[...] = (%s).astype("%s")' % (newFormula, calcType))
pyFormulaeBands.append(outBands)
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)
outOptions.append('BLOCKXSIZE=%i' % blockSize[0])
outOptions.append('BLOCKYSIZE=%i' % blockSize[1])
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())
if ndv is not None:
for b in range(outBandsTotal):
outDataset.GetRasterBand(b+1).SetNoDataValue(ndv)
#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(0, ySize, tileYSize):
tileHeight = min(tileYSize, ySize-tileYStart)
outBand = 1
for f in range(len(pyFormulae)):
if pyFormulae[f].startswith('result'):
if pyFormulaeBands[f]==1:
result = N.zeros((tileHeight, xSize), dtype=outType)
else:
result = N.zeros((pyFormulaeBands[f], tileHeight, xSize), dtype=outType)
else:
result = None
#print pyFormulae[f]
#exec(pyFormulae[f]) in locals()
exec(pyFormulae[f]) in globals(), locals()
#print f, pyFormulae[f]
#print result.min(), result.max()
if result is not None:
if pyFormulaeBands[f]==1:
gdal_array.BandWriteArray(outDataset.GetRasterBand(outBand), result, 0, tileYStart)
else:
for bi in range(pyFormulaeBands[f]):
gdal_array.BandWriteArray(outDataset.GetRasterBand(outBand + bi), result[bi], 0, tileYStart)
outBand += pyFormulaeBands[f]
floatPercentage = 100.0 * (float(tileYStart + (tileHeight / len(pyFormulae) * (f+1))) / 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%% || %f %f %f" % (pbar, floatPercentage, N.min(result), N.mean(result), N.max(result)),
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
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
