Index: swig/python/scripts/gdal_calc.py =================================================================== --- swig/python/scripts/gdal_calc.py (revision 26841) +++ swig/python/scripts/gdal_calc.py (working copy) @@ -1,4 +1,4 @@ -#!/usr/bin/env python +#! /usr/bin/python # -*- coding: utf-8 -*- #****************************************************************************** # @@ -107,12 +107,28 @@ if opts.debug: print("file %s: %s, dimensions: %s, %s, type: %s" %(myI,myF,DimensionsCheck[0],DimensionsCheck[1],myDataType[i])) + # process allBands option + allBandsIndex=None + allBandsCount=1 + if opts.allBands: + try: + allBandsIndex=myAlphaList.index(opts.allBands) + except ValueError: + print("Error! allBands option was given but Band %s not found. Cannot proceed" % (opts.allBands)) + return + allBandsCount=myFiles[allBandsIndex].RasterCount + if allBandsCount <= 1: + allBandsIndex=None + ################################################################ # set up output file ################################################################ # open output file exists if os.path.isfile(opts.outF) and not opts.overwrite: + if allBandsIndex is not None: + print("Error! allBands option was given but Output file exists, must use --overwrite option!") + return if opts.debug: print("Output file %s exists - filling in results into file" %(opts.outF)) myOut=gdal.Open(opts.outF, gdal.GA_Update) @@ -141,24 +157,24 @@ # create file myOutDrv = gdal.GetDriverByName(opts.format) myOut = myOutDrv.Create( - opts.outF, DimensionsCheck[0], DimensionsCheck[1], 1, + opts.outF, DimensionsCheck[0], DimensionsCheck[1], allBandsCount, gdal.GetDataTypeByName(myOutType), opts.creation_options) # set output geo info based on first input layer myOut.SetGeoTransform(myFiles[0].GetGeoTransform()) myOut.SetProjection(myFiles[0].GetProjection()) - myOutB=myOut.GetRasterBand(1) if opts.NoDataValue!=None: myOutNDV=opts.NoDataValue else: myOutNDV=DefaultNDVLookup[myOutType] - myOutB.SetNoDataValue(myOutNDV) - # write to band - myOutB=None - # refetch band - myOutB=myOut.GetRasterBand(1) + for i in range(1,allBandsCount+1): + myOutB=myOut.GetRasterBand(i) + myOutB.SetNoDataValue(myOutNDV) + # write to band + myOutB=None + if opts.debug: print("output file: %s, dimensions: %s, %s, type: %s" %(opts.outF,myOut.RasterXSize,myOut.RasterYSize,gdal.GetDataTypeName(myOutB.DataType))) @@ -182,81 +198,92 @@ # variables for displaying progress ProgressCt=-1 ProgressMk=-1 - ProgressEnd=nXBlocks*nYBlocks + ProgressEnd=nXBlocks*nYBlocks*allBandsCount ################################################################ - # start looping through blocks of data + # start looping through each band in allBandsCount ################################################################ - # loop through X-lines - for X in range(0,nXBlocks): + for bandNo in range(1,allBandsCount+1): - # in the rare (impossible?) case that the blocks don't fit perfectly - # change the block size of the final piece - if X==nXBlocks-1: - nXValid = DimensionsCheck[0] - X * myBlockSize[0] - myBufSize = nXValid*nYValid + ################################################################ + # start looping through blocks of data + ################################################################ - # find X offset - myX=X*myBlockSize[0] + # loop through X-lines + for X in range(0,nXBlocks): - # reset buffer size for start of Y loop - nYValid = myBlockSize[1] - myBufSize = nXValid*nYValid - - # loop through Y lines - for Y in range(0,nYBlocks): - ProgressCt+=1 - if 10*ProgressCt/ProgressEnd%10!=ProgressMk: - ProgressMk=10*ProgressCt/ProgressEnd%10 - from sys import version_info - if version_info >= (3,0,0): - exec('print("%d.." % (10*ProgressMk), end=" ")') - else: - exec('print 10*ProgressMk, "..",') - + # in the rare (impossible?) case that the blocks don't fit perfectly # change the block size of the final piece - if Y==nYBlocks-1: - nYValid = DimensionsCheck[1] - Y * myBlockSize[1] + if X==nXBlocks-1: + nXValid = DimensionsCheck[0] - X * myBlockSize[0] myBufSize = nXValid*nYValid - # find Y offset - myY=Y*myBlockSize[1] + # find X offset + myX=X*myBlockSize[0] - # create empty buffer to mark where nodata occurs - myNDVs=numpy.zeros(myBufSize) - myNDVs.shape=(nYValid,nXValid) + # reset buffer size for start of Y loop + nYValid = myBlockSize[1] + myBufSize = nXValid*nYValid - # fetch data for each input layer - for i,Alpha in enumerate(myAlphaList): + # loop through Y lines + for Y in range(0,nYBlocks): + ProgressCt+=1 + if 10*ProgressCt/ProgressEnd%10!=ProgressMk: + ProgressMk=10*ProgressCt/ProgressEnd%10 + from sys import version_info + if version_info >= (3,0,0): + exec('print("%d.." % (10*ProgressMk), end=" ")') + else: + exec('print 10*ProgressMk, "..",') - # populate lettered arrays with values - myval=BandReadAsArray(myFiles[i].GetRasterBand(myBands[i]), - xoff=myX, yoff=myY, - win_xsize=nXValid, win_ysize=nYValid) + # change the block size of the final piece + if Y==nYBlocks-1: + nYValid = DimensionsCheck[1] - Y * myBlockSize[1] + myBufSize = nXValid*nYValid - # fill in nodata values - myNDVs=1*numpy.logical_or(myNDVs==1, myval==myNDV[i]) + # find Y offset + myY=Y*myBlockSize[1] - # create an array of values for this block - exec("%s=myval" %Alpha) - myval=None + # create empty buffer to mark where nodata occurs + myNDVs=numpy.zeros(myBufSize) + myNDVs.shape=(nYValid,nXValid) + # fetch data for each input layer + for i,Alpha in enumerate(myAlphaList): - # try the calculation on the array blocks - try: - myResult = eval(opts.calc) - except: - print("evaluation of calculation %s failed" %(opts.calc)) - raise + # populate lettered arrays with values + if allBandsIndex is not None and allBandsIndex==i: + myBandNo=bandNo + else: + myBandNo=myBands[i] + myval=BandReadAsArray(myFiles[i].GetRasterBand(myBandNo), + xoff=myX, yoff=myY, + win_xsize=nXValid, win_ysize=nYValid) - # propogate nodata values - # (set nodata cells to zero then add nodata value to these cells) - myResult = ((1*(myNDVs==0))*myResult) + (myOutNDV*myNDVs) + # fill in nodata values + myNDVs=1*numpy.logical_or(myNDVs==1, myval==myNDV[i]) - # write data block to the output file - BandWriteArray(myOutB, myResult, xoff=myX, yoff=myY) + # create an array of values for this block + exec("%s=myval" %Alpha) + myval=None + + # try the calculation on the array blocks + try: + myResult = eval(opts.calc) + except: + print("evaluation of calculation %s failed" %(opts.calc)) + raise + + # propogate nodata values + # (set nodata cells to zero then add nodata value to these cells) + myResult = ((1*(myNDVs==0))*myResult) + (myOutNDV*myNDVs) + + # write data block to the output file + myOutB=myOut.GetRasterBand(bandNo) + BandWriteArray(myOutB, myResult, xoff=myX, yoff=myY) + print("100 - Done") #print("Finished - Results written to %s" %opts.outF) @@ -283,6 +310,7 @@ help="Passes a creation option to the output format driver. Multiple" "options may be listed. See format specific documentation for legal" "creation options for each format.") + parser.add_option("--allBands", dest="allBands", default="", help="process all Bands for given raster (A-Z)") parser.add_option("--overwrite", dest="overwrite", action="store_true", help="overwrite output file if it already exists") parser.add_option("--debug", dest="debug", action="store_true", help="print debugging information")