Package firetools :: Module GP_Util
[hide private]
[frames] | no frames]

Source Code for Module firetools.GP_Util

  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  ## Raster Conversiont Functions 
 29  ## txt(DAT) <=> Ascii <=> Raster <=> numpy 
30 -def JoinDateRasters(npylist,outfile,cleanup=False):
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
44 -def NpySeriesToRaster(study,baseseries):
45 """Convert a npy series to esri rasters""" 46 asciilist = NpyListToAscii(study.series[baseseries]) 47 rasterlist = AsciiListToRaster(asciilist,cleanup=False) 48 return rasterlist
49
50 -def RasterToNpySeries(inlist,cleanup=True):
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
66 -def NpyListToAscii(inlist):
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
78 -def DatListToAscii(inlist):
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) #(merge2?) 83 AddAsciiHeader(outfile) 84 log.info('DatListToAscii') 85 return outlist
86
87 -def AsciiListToRaster(inlist,cleanup=True,outtype='INTEGER'):
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
98 -def DatListToNpy(inlist):
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
108 -def AddAsciiHeader(infile,nrow=1024, ncol=768, xll=001, yll=001, cellsize=001, nodata=-1):
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
128 -class ProgDots:
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
147 - def step(self):
148 self.stepnow += 1 149 while self.stepnow > self.seglist[0]: 150 self.advanceBar() 151 self.calcEta() 152 self.writeBar()
153
154 - def advanceBar(self):
155 self.seglist.pop(0) 156 self.progbar.pop() 157 self.progbar.insert(0,'>')
158
159 - def calcEta(self):
160 elapsed = (time.time() - self.starttime) / 60.0 #time in min 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
167 - def writeBar(self):
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 ## shapefile functions 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
196 -def GetShpField(inshape,field):
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
205 -def PointsFromShpExtent(inshape,outpoints,cellsize):
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
216 -def GetShpVals(inshape,fields,sortexpr=""):
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 #Get Row Length 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
237 -def ExtractRasterData(points,fieldname,raster):
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
247 -def CalculateShpArea(shapefile):
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
264 -def FilterFields(shapefile,keepfields):
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
272 -def ConvertDateField(inshape,field,datetype):
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
284 -def CopyField(inshape,newfield,oldfield):
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 ## General Utilities 297 ##
298 -def MakeFolder(target):
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
313 -def ArrayToCSV(outfile,inarray,header=[]):
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
323 -def AsciiRaster2Array(inraster,ndval=-9999,nansub=True):
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#decimal date 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}#days into the year of each month 338 339 #convert input date to decimal year 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 #convert to output datetype and return 364 if outtype == "YYYY": 365 return decdate
366 367 368 # 369