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