Source Code for Module pgeo.gis.scripts.gdal_calculations

  1  #! /usr/bin/python 
  2  # -*- coding: utf-8 -*- 
  3  #****************************************************************************** 
  4  # 
  5  #  Project:  GDAL 
  6  #  Purpose:  Command line raster calculator for gdal supported files 
  7  #  Author:   Chris Yesson, chris.yesson@ioz.ac.uk 
  8  # 
  9  #****************************************************************************** 
 10  #  Copyright (c) 2010, Chris Yesson <chris.yesson@ioz.ac.uk> 
 11  # 
 12  #  Permission is hereby granted, free of charge, to any person obtaining a 
 13  #  copy of this software and associated documentation files (the "Software"), 
 14  #  to deal in the Software without restriction, including without limitation 
 15  #  the rights to use, copy, modify, merge, publish, distribute, sublicense, 
 16  #  and/or sell copies of the Software, and to permit persons to whom the 
 17  #  Software is furnished to do so, subject to the following conditions: 
 18  # 
 19  #  The above copyright notice and this permission notice shall be included 
 20  #  in all copies or substantial portions of the Software. 
 21  # 
 22  #  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 
 23  #  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 
 24  #  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 
 25  #  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 
 26  #  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 
 27  #  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 
 28  #  DEALINGS IN THE SOFTWARE. 
 29  #****************************************************************************** 
 30   
 31  ################################################################ 
 32  # Command line raster calculator for gdal supported files. Use any basic arithmetic supported by numpy arrays such as +-*\ along with logical operators such as >.  Note that all files must be the same dimensions, but no projection checking is performed.  Use gdal_calc.py --help for list of options. 
 33   
 34  # example 1 - add two files together 
 35  # gdal_calc.py -A input1.tif -B input2.tif --outfile=result.tif --calc="A+B" 
 36   
 37  # example 2 - average of two layers 
 38  # gdal_calc.py -A input.tif -B input2.tif --outfile=result.tif --calc="(A+B)/2" 
 39   
 40  # example 3 - set values of zero and below to null 
 41  # gdal_calc.py -A input.tif --outfile=result.tif --calc="A*(A>0)" --NoDataValue=0 
 42  ################################################################ 
 43   
 44  try: 
 45      from osgeo import gdal 
 46      from osgeo.gdalnumeric import * 
 47  except ImportError: 
 48      import gdal 
 49      from gdalnumeric import * 
 50   
 51  from optparse import OptionParser 
 52  import os 
 53  import sys 
 54   
 55  # set up some default nodatavalues for each datatype 
 56  DefaultNDVLookup={'Byte':255, 'UInt16':65535, 'Int16':-32767, 'UInt32':4294967293, 'Int32':-2147483647, 'Float32':1.175494351E-38, 'Float64':1.7976931348623158E+308} 
 57   
 58  #3.402823466E+38 
 59   
 60  ################################################################ 
 61 -def doit(opts, AlphaList, len_argv): 
 62   
 63      if opts.debug: 
 64          print("gdal_calc.py starting calculation %s" %(opts.calc)) 
 65   
 66      ################################################################ 
 67      # fetch details of input layers 
 68      ################################################################ 
 69   
 70      # set up some lists to store data for each band 
 71      myFiles=[] 
 72      myBands=[] 
 73      myAlphaList=[] 
 74      myDataType=[] 
 75      myDataTypeNum=[] 
 76      myNDV=[] 
 77      DimensionsCheck=None 
 78   
 79      # loop through input files - checking dimensions 
 80      for i,myI in enumerate(AlphaList[0:len_argv-1]): 
 81          myF = eval("opts.%s" %(myI)) 
 82          myBand = eval("opts.%s_band" %(myI)) 
 83   
 84          if myF: 
 85              myFiles.append(gdal.Open(myF, gdal.GA_ReadOnly)) 
 86              # check if we have asked for a specific band... 
 87              if myBand: 
 88                  myBands.append(myBand) 
 89              else: 
 90                  myBands.append(1) 
 91              myAlphaList.append(myI) 
 92              myDataType.append(gdal.GetDataTypeName(myFiles[i].GetRasterBand(myBands[i]).DataType)) 
 93              myDataTypeNum.append(myFiles[i].GetRasterBand(myBands[i]).DataType) 
 94              myNDV.append(myFiles[i].GetRasterBand(myBands[i]).GetNoDataValue()) 
 95              # check that the dimensions of each layer are the same 
 96              if DimensionsCheck: 
 97                  if DimensionsCheck!=[myFiles[i].RasterXSize, myFiles[i].RasterYSize]: 
 98                      print("Error! Dimensions of file %s (%i, %i) are different from other files (%i, %i).  Cannot proceed" % \ 
 99                            (myF,myFiles[i].RasterXSize, myFiles[i].RasterYSize,DimensionsCheck[0],DimensionsCheck[1])) 
100                      return 
101              else: 
102                  DimensionsCheck=[myFiles[i].RasterXSize, myFiles[i].RasterYSize] 
103   
104              if opts.debug: 
105                  print("file %s: %s, dimensions: %s, %s, type: %s" %(myI,myF,DimensionsCheck[0],DimensionsCheck[1],myDataType[i])) 
106   
107      ################################################################ 
108      # set up output file 
109      ################################################################ 
110   
111      # open output file exists 
112      if os.path.isfile(opts.outF) and not opts.overwrite: 
113          if opts.debug: 
114              print("Output file %s exists - filling in results into file" %(opts.outF)) 
115          myOut=gdal.Open(opts.outF, gdal.GA_Update) 
116          if [myOut.RasterXSize,myOut.RasterYSize] != DimensionsCheck: 
117              print("Error! Output exists, but is the wrong size.  Use the --overwrite option to automatically overwrite the existing file") 
118              return 
119          myOutB=myOut.GetRasterBand(1) 
120          myOutNDV=myOutB.GetNoDataValue() 
121          myOutType=gdal.GetDataTypeName(myOutB.DataType) 
122   
123      else: 
124          # remove existing file and regenerate 
125          if os.path.isfile(opts.outF): 
126              os.remove(opts.outF) 
127          # create a new file 
128          if opts.debug: 
129              print("Generating output file %s" %(opts.outF)) 
130   
131          # find data type to use 
132          if not opts.type: 
133              # use the largest type of the input files 
134              myOutType=gdal.GetDataTypeName(max(myDataTypeNum)) 
135          else: 
136              myOutType=opts.type 
137   
138          # create file 
139          myOutDrv = gdal.GetDriverByName(opts.format) 
140          myOut = myOutDrv.Create( 
141              opts.outF, DimensionsCheck[0], DimensionsCheck[1], 1, 
142              gdal.GetDataTypeByName(myOutType), opts.creation_options) 
143   
144          # set output geo info based on first input layer 
145          myOut.SetGeoTransform(myFiles[0].GetGeoTransform()) 
146          myOut.SetProjection(myFiles[0].GetProjection()) 
147   
148          myOutB=myOut.GetRasterBand(1) 
149          if opts.NoDataValue!=None: 
150              myOutNDV=opts.NoDataValue 
151          else: 
152              myOutNDV=DefaultNDVLookup[myOutType] 
153          myOutB.SetNoDataValue(myOutNDV) 
154          # write to band 
155          myOutB=None 
156          # refetch band 
157          myOutB=myOut.GetRasterBand(1) 
158   
159      if opts.debug: 
160          print("output file: %s, dimensions: %s, %s, type: %s" %(opts.outF,myOut.RasterXSize,myOut.RasterYSize,gdal.GetDataTypeName(myOutB.DataType))) 
161   
162      ################################################################ 
163      # find block size to chop grids into bite-sized chunks  
164      ################################################################ 
165   
166      # use the block size of the first layer to read efficiently 
167      myBlockSize=myFiles[0].GetRasterBand(myBands[0]).GetBlockSize(); 
168      # store these numbers in variables that may change later 
169      nXValid = myBlockSize[0] 
170      nYValid = myBlockSize[1] 
171      # find total x and y blocks to be read 
172      nXBlocks = (int)((DimensionsCheck[0] + myBlockSize[0] - 1) / myBlockSize[0]); 
173      nYBlocks = (int)((DimensionsCheck[1] + myBlockSize[1] - 1) / myBlockSize[1]); 
174      myBufSize = myBlockSize[0]*myBlockSize[1] 
175   
176      if opts.debug: 
177          print("using blocksize %s x %s" %(myBlockSize[0], myBlockSize[1])) 
178   
179      # variables for displaying progress 
180      ProgressCt=-1 
181      ProgressMk=-1 
182      ProgressEnd=nXBlocks*nYBlocks 
183   
184      ################################################################ 
185      # start looping through blocks of data 
186      ################################################################ 
187   
188      # loop through X-lines 
189      for X in range(0,nXBlocks): 
190   
191          # in the rare (impossible?) case that the blocks don't fit perfectly 
192          # change the block size of the final piece 
193          if X==nXBlocks-1: 
194              nXValid = DimensionsCheck[0] - X * myBlockSize[0] 
195              myBufSize = nXValid*nYValid 
196   
197          # find X offset 
198          myX=X*myBlockSize[0] 
199   
200          # reset buffer size for start of Y loop 
201          nYValid = myBlockSize[1] 
202          myBufSize = nXValid*nYValid 
203   
204          # loop through Y lines 
205          for Y in range(0,nYBlocks): 
206              ProgressCt+=1 
207              if 10*ProgressCt/ProgressEnd%10!=ProgressMk: 
208                  ProgressMk=10*ProgressCt/ProgressEnd%10 
209                  from sys import version_info 
210                  if version_info >= (3,0,0): 
211                      exec('print("%d.." % (10*ProgressMk), end=" ")') 
212                  else: 
213                      exec('print 10*ProgressMk, "..",') 
214   
215              # change the block size of the final piece 
216              if Y==nYBlocks-1: 
217                  nYValid = DimensionsCheck[1] - Y * myBlockSize[1] 
218                  myBufSize = nXValid*nYValid 
219   
220              # find Y offset 
221              myY=Y*myBlockSize[1] 
222   
223              # create empty buffer to mark where nodata occurs 
224              myNDVs=numpy.zeros(myBufSize) 
225              myNDVs.shape=(nYValid,nXValid) 
226   
227              # fetch data for each input layer 
228              for i,Alpha in enumerate(myAlphaList): 
229   
230                  # populate lettered arrays with values 
231                  myval=BandReadAsArray(myFiles[i].GetRasterBand(myBands[i]), 
232                                        xoff=myX, yoff=myY, 
233                                        win_xsize=nXValid, win_ysize=nYValid) 
234   
235                  # fill in nodata values 
236                  myNDVs=1*numpy.logical_or(myNDVs==1, myval==myNDV[i]) 
237   
238                  # create an array of values for this block 
239                  exec("%s=myval" %Alpha) 
240                  myval=None 
241   
242   
243              # try the calculation on the array blocks 
244              try: 
245                  myResult = eval(opts.calc) 
246              except: 
247                  print("evaluation of calculation %s failed" %(opts.calc)) 
248                  raise 
249   
250              # propogate nodata values  
251              # (set nodata cells to zero then add nodata value to these cells) 
252              myResult = ((1*(myNDVs==0))*myResult) + (myOutNDV*myNDVs) 
253   
254              # write data block to the output file 
255              BandWriteArray(myOutB, myResult, xoff=myX, yoff=myY) 
256   
257      #print("Finished - Results written to %s" %opts.outF) 
258   
259      return 
260   
261  ################################################################ 
262 -def main(): 
263      print "gdal_calculations " + str(sys.argv) 
264      index = 0 
265      for args in sys.argv: 
266          if "alphalist" in args: 
267              AlphaList = sys.argv[index+1] 
268          index +=1 
269   
270      # getting the Alphalist (list of the layer's alias) 
271      AlphaList = eval(AlphaList) 
272   
273      parser = OptionParser() 
274      # define options 
275      parser.add_option("--calc", dest="calc", help="calculation in gdalnumeric syntax using +-/* or any numpy array functions (i.e. logical_and())") 
276      # hack to limit the number of input file options close to required number 
277   
278      for myAlpha in AlphaList[0:len(sys.argv)-1]: 
279          eval('parser.add_option("--%s", dest="%s", help="input gdal raster file, note you can use any letter A-Z")' %(myAlpha, myAlpha)) 
280          eval('parser.add_option("--%s_band", dest="%s_band", default=0, type=int, help="number of raster band for file %s")' %(myAlpha, myAlpha, myAlpha)) 
281   
282      parser.add_option("--alphalist", dest="",  default='', help="") 
283      parser.add_option("--outfile", dest="outF", default='gdal_calc.tif', help="output file to generate or fill.") 
284      parser.add_option("--NoDataValue", dest="NoDataValue", type=float, help="set output nodatavalue (Defaults to datatype specific values)") 
285      parser.add_option("--type", dest="type", help="set datatype must be one of %s" % list(DefaultNDVLookup.keys())) 
286      parser.add_option("--format", dest="format", default="GTiff", help="GDAL format for output file (default 'GTiff')") 
287      parser.add_option("--creation-option", "--co", dest="creation_options", default=[], action="append", 
288                        help="Passes a creation option to the output format driver. Multiple" 
289                             "options may be listed. See format specific documentation for legal" 
290                             "creation options for each format.") 
291      parser.add_option("--overwrite", dest="overwrite", action="store_true", help="overwrite output file if it already exists") 
292      parser.add_option("--debug", dest="debug", action="store_true", help="print debugging information") 
293   
294      (opts, args) = parser.parse_args() 
295      if len(sys.argv) == 1: 
296          parser.print_help() 
297      elif not opts.calc: 
298          print("No calculation provided.  Nothing to do!") 
299          parser.print_help() 
300      else: 
301          doit(opts, AlphaList, len(sys.argv)) 
302   
303  if __name__ == "__main__": 
304      main() 
305