1  """Geoprocessing Utility Functions 
  2   
  3          General Utility functions used by other modules within the tool 
  4                  Import and provide access to the geoprocessor object (gp)   
  5                  raster file format conversion functions 
  6                  generic shapefile utilities (fields, etc) 
  7                  generic utilities 
  8   
  9  """ 
 10   
 11  import os, sys, logging, time 
 12  import string, shutil, csv, pickle, copy 
 13  import numpy as num 
 14  log = logging.getLogger('firetools'.ljust(10)) 
 15   
 16  """ 
 17  try: 
 18          import arcgisscripting 
 19          gp = arcgisscripting.create(9.3) 
 20          gp.SetProduct("ArcInfo") 
 21          gp.CheckOutExtension("spatial") 
 22          gp.OverWriteOutput = 1 
 23  except: 
 24          log.info('GP import Fail') 
 25          gp = 'gp' 
 26  """ 
 27   
 28   
 29   
 31          """Load the .npy files in npylist and join them as columns in the output table""" 
 32          os.chdir(os.path.dirname(npylist[0])) 
 33          dlist = [num.load(df) for df in npylist]  
 34          dlist = [dt.flatten() for dt in dlist]   
 35          dtable = num.vstack(dlist).transpose() 
 36          num.save(outfile, dtable) 
 37           
 38          if cleanup == True:  
 39                  [os.remove(f) for f in npylist] 
 40                   
 41          return os.path.abspath(outfile+'.npy') 
  42           
 43           
 45          """Convert a npy series to esri rasters""" 
 46          asciilist = NpyListToAscii(study.series[baseseries])     
 47          rasterlist = AsciiListToRaster(asciilist,cleanup=False) 
 48          return rasterlist 
  49   
 51          outlist = [f+'.npy' for f in inlist] 
 52          asciilist = [f+'.txt' for f in inlist] 
 53          prog = ProgDots(len(inlist),title='RasterToNpy') 
 54          for (infile,asciifile,outfile) in zip(inlist,asciilist,outlist): 
 55                  gp.RasterToASCII_Conversion(infile,asciifile) 
 56                  data = num.loadtxt(asciifile,skiprows=6) 
 57                  num.save(outfile,data) 
 58                  os.remove(asciifile) 
 59                  if cleanup == True: 
 60                          gp.delete(infile) 
 61                  prog.step() 
 62           
 63           
 64          return outlist 
  65           
 67          """Covert .npy files to import ready ascii rasters""" 
 68          prog = ProgDots(len(inlist),title='NpyListToAscii') 
 69          outlist = [infile[:-4]+'.TXT' for infile in inlist] 
 70          for (infile,outfile) in zip(inlist,outlist): 
 71                  data = num.load(infile) 
 72                  num.savetxt(outfile,data,fmt='%i') 
 73                  AddAsciiHeader(outfile) 
 74                  prog.step() 
 75                   
 76          return outlist 
  77           
 79          """Convert Dat files to import ready ascii rasters""" 
 80          outlist = [infile[:-4]+'.TXT' for infile in inlist] 
 81          for (infile,outfile) in zip(inlist,outlist): 
 82                  shutil.copy(infile,outfile)  
 83                  AddAsciiHeader(outfile) 
 84          log.info('DatListToAscii') 
 85          return outlist 
  86           
 88          """Convert Ascii rasters to ESRI rasters""" 
 89          prog = ProgDots(len(inlist),title='AsciiListToRaster') 
 90          outlist = [txt[:-4] for txt in inlist] 
 91          for (infile,outfile) in zip(inlist,outlist): 
 92                  gp.ASCIIToRaster_conversion(infile, outfile, outtype ) 
 93                  if cleanup == True: 
 94                          gp.delete(infile) 
 95                  prog.step() 
 96          return outlist 
  97   
 99          """Convert DAT files to npy files""" 
100          outlist = [infile[:-4]+'.npy' for infile in inlist] 
101          prog = ProgDots(len(inlist),title='DatListToNpy') 
102          for (infile,outfile) in zip(inlist,outlist): 
103                  data = num.loadtxt(infile) 
104                  num.save(outfile,data) 
105                  prog.step() 
106          return outlist 
 107   
109          """Add the required import header for txt>raster""" 
110          header = '' 
111          header += 'NCOLS '+ str(ncol) +'\n' 
112          header += 'NROWS '+ str(nrow) +'\n' 
113          header += 'XLLCORNER '+ str(xll)+'\n' 
114          header += 'YLLCORNER '+ str(yll) +'\n' 
115          header += 'CELLSIZE '+str(cellsize)+'\n' 
116          header += 'NODATA_VALUE '+str(nodata)+'\n' 
117           
118          f = open(infile,'r') 
119          data = f.read() 
120          f.close() 
121           
122          f = open(infile,'w') 
123          f.write(header) 
124          f.write(data) 
125          f.close() 
 126   
127   
129          """Progress Bar Class with calculated averate eta 
130                  [..........] -> [>>>>......] -> [>>>>>>>>>>] 
131                  inputs: 
132                          totalsteps = number of calls to be made to step() 
133                          title = display title 
134                          segcount = number of characters in progress bar 
135                           
136                   
137          """ 
138 -        def __init__(self,totalsteps,title='',segcount=40): 
 139                  self.stepnow = 0 
140                  self.totalsteps = totalsteps 
141                  self.seglist=num.linspace(0,totalsteps,segcount+1)[1:].tolist() 
142                  self.progbar = ['.' for i in self.seglist] 
143                  self.starttime = time.time() 
144                  self.log = logging.getLogger('Progress'.ljust(10)) 
145                  self.log.info(title.ljust(50)) 
 146                   
153                   
155                  self.seglist.pop(0) 
156                  self.progbar.pop() 
157                  self.progbar.insert(0,'>') 
 158   
160                  elapsed = (time.time() - self.starttime) / 60.0  
161                  meanstep = elapsed / self.stepnow 
162                  remtime = meanstep * (self.totalsteps - self.stepnow) 
163                   
164                  self.eta = str(num.round(remtime,2)) 
165                  self.elapsed = str(num.round(elapsed,2)) 
 166   
168                  msg = string.join(self.progbar,'') 
169                  msg += 'eta/elap' 
170                  msg += '('+self.eta + '/' + self.elapsed+')' 
171                  msg += '\r' 
172                  sys.stdout.write(msg.ljust(50)) 
  173                   
174   
175   
176   
177   
178 -def ExportShapeTable(inshape="datagrid.shp", fieldlist=[], outfile='table', sortexpr="" ): 
 179          """export the table data for a shapefile 
180                  inshape = input shape file 
181                  fieldlist = list of fields to export 
182                  outfile = output table name without ".csv" 
183                  sortexpr = sorting of the table 
184                  """ 
185          log.info('Exporting table from shapefile')       
186          outfile = os.path.join(os.path.dirname(inshape), outfile+'.csv' ) 
187          data = [] 
188          rows = gp.Searchcursor(inshape,"","","",sortexpr) 
189          row = rows.Next() 
190          while row: 
191                  data.append([getattr(row,field) for field in fieldlist])  
192                  row = rows.Next() 
193                   
194          ArrayToCSV(outfile, data, fieldlist) 
 195           
197          data = [] 
198          rows = gp.Searchcursor(inshape,"","","","") 
199          row = rows.Next() 
200          while row: 
201                  data.append(getattr(row,field)) 
202                  row = rows.Next() 
203          return num.array(data) 
 204   
206          """Create a point grid based on cellsize and the extent of a shapefile"""  
207          gp.Workspace = os.path.split(outpoints)[0] 
208          gp.AddField_management(inshape, "Z", "DOUBLE") 
209          gp.PolygonToRaster_conversion(inshape,'Z','studyraster','#','#',cellsize) 
210          gp.RasterToPoint_conversion('studyraster',outpoints,'VALUE') 
211          gp.delete('studyraster') 
212          gp.DeleteField_management(inshape,"Z") 
213          FilterFields(outpoints,['POINTID']) 
214          gp.CalculateField_management(outpoints, 'POINTID', ("!FID!"), "PYTHON") 
 215   
217          """Retrieve the values from select fields in a shapefile 
218          returns dictionary with vector for each field key""" 
219          data = dict([(f,[]) for f in fields]) 
220           
221           
222          rows = gp.Searchcursor(inshape,"","","",sortexpr) 
223          row = rows.Next() 
224          while row: 
225                  for f in fields: 
226                          data[f].append(getattr(row,f)) 
227                  row = rows.Next() 
228          del row,rows 
229   
230          for f in fields: 
231                  data[f] = num.array(data[f],dtype='f') 
232   
233          return data 
 234   
235   
236           
238          """Extract data from a raster into point shapefile""" 
239          tempraster=os.path.join(os.path.dirname(raster),'extract.shp') 
240          gp.ExtractValuesToPoints_sa(points,raster,tempraster) 
241          gp.AddField_management(tempraster,fieldname,"DOUBLE") 
242          gp.CalculateField_management(tempraster,fieldname,"!RASTERVALU!","PYTHON","") 
243          gp.DeleteField(tempraster,'RASTERVALU') 
244          gp.copy(tempraster,points) 
245          gp.delete_management(tempraster) 
 246           
248          """Get total area of a shapefile""" 
249          rows = gp.Searchcursor(shapefile) 
250          row = rows.Next() 
251          totalarea = 0 
252          while row: 
253                  feat = row.shape 
254                  totalarea += float(feat.Area) 
255                  row = rows.Next() 
256          del row, rows 
257          return totalarea 
 258   
259 -def ShapeBuffer(inshape, datafield, selectvalue, buffdist, buffshape): 
 260          gp.MakeFeatureLayer(inshape,"datalyr") 
261          gp.SelectLayerByAttribute("datalyr", "NEW_SELECTION", datafield+' == '+selectvalue) 
262          gp.Buffer_Analysis('datalyr',buffshape,buffdist) 
 263             
265          """Remove fields except for FID,Shape and those specified in keepfields""" 
266          keepfields.extend(['FID','Shape']) 
267          fields = gp.ListFields(shapefile) 
268          for field in fields: 
269                  if field.Name not in keepfields: 
270                          gp.DeleteField_management(shapefile,field.Name) 
 271                           
273          """apply date conversion to all rows of a field""" 
274          rows = gp.Updatecursor(inshape) 
275          row = rows.Next() 
276          while row: 
277                  olddate = getattr(row,field) 
278                  newdate = convertdates(olddate,datetype) 
279                  setattr(row,field,newdate) 
280                  rows.UpdateRow(row) 
281                  row = rows.Next() 
282          del row,rows         
283   
285          """Copy existing field to new field""" 
286          try: 
287                  gp.AddField_management(inshape, newfield, "DOUBLE") 
288          except: 
289                  log.info('WARNING: Cant add field ' + newfield) 
290                   
291          gp.CalculateField_management(inshape, newfield, ("!" + oldfield + "!"), "PYTHON") 
 292           
293   
294   
295   
296   
297   
299          """makes or overwrites target directory and contained files""" 
300           
301          if os.path.exists(target): 
302                  for root, dirs, files in os.walk(target, topdown=False): 
303                          for name in files: 
304                                  os.remove(os.path.join(root, name)) 
305                          for name in dirs: 
306                                  os.rmdir(os.path.join(root, name)) 
307                  os.rmdir(target) 
308                   
309          os.mkdir(target) 
310          os.chdir(target) 
311          return target 
 312   
314          """Export an array""" 
315          f = open(outfile,'wb') 
316          writer = csv.writer(f) 
317          if len(header) > 0: 
318                  writer.writerow(header) 
319          for row in inarray: 
320                  writer.writerow(row) 
321          f.close() 
 322           
324          data = num.loadtxt(inraster,skiprows=6) 
325          if nansub == True: 
326                  data[data == ndval] = num.nan 
327          return data 
 328   
329 -def convertdates(datein, intype='YYYYMMDD', outtype='YYYY', defdate=0, dateprecision=3): 
 330          """convert one date format to another 
331          assign default dates as specified 
332          round to dateprecision as specified 
333          intype/outtype options: 'YYYY' 'YYYYMMDD'""" 
334   
335          decdate = 0.0 
336          monthdays = {1:31, 2:59, 3:90, 4:120, 5:151, 6:181, 7:212, 8:243, 
337          9:273, 10:304, 11:334, 12:365} 
338   
339           
340          if intype == "YYYYMMDD": 
341                  YYYYMMDD = float(datein) 
342                  year = num.floor(YYYYMMDD/10000) 
343                  MMDD = YYYYMMDD - (year*10000) 
344                  mm = num.floor(MMDD/100) 
345                  dd = MMDD - (mm*100)         
346                  if mm == 0: 
347                          decdate = year 
348                  elif mm > 12: 
349                          decdate = year     
350                  elif mm == 1: 
351                          decdate = year + (dd/365) 
352                  elif mm > 1: 
353                          m2 = mm - 1 
354                          decdate = year + (( monthdays[m2] + dd ) /365) 
355          elif intype == "YYYY": 
356                  decdate = datein 
357   
358          if num.floor(decdate) == decdate: 
359                  decdate = decdate + float(defdate) 
360   
361          decdate = num.round(decdate,dateprecision) 
362           
363           
364          if outtype == "YYYY": 
365                  return decdate 
 366   
367   
368   
369