1 import re
2 import sys
3 import os
4 import time
5 import shutil
6 import subprocess
7 from madrona.common.utils import get_logger
8 log = get_logger()
9
10 DJANGO = True
11 try:
12 from django.conf import settings
13 except ImportError:
14
15 DJANGO = False
16
17
18 '''
19 Sets up a grass environment in the current running process. Allows grass commands
20 to be run in a shell environment. Optimally grass python bindings would be used in the
21 future once they become more stable and full-featured
22
23 Grass is not thread-safe, see http://grass.osgeo.org/wiki/Parallel_GRASS_jobs
24 Because of this there can be issues with race conditions while working within
25 the same mapset with more than one process. The solution is to create a new
26 mapset for each analysis run, copy the necessary maps into the new mapset,
27 run the analysis, and then remove the mapset when you're done.
28
29 @author: Tim Welch, Ecotrust 2009
30 @author: Matthew Perry, Ecotrust 2011
31 Concepts taken from the PyWPS project.
32
33 Example usage (django):
34 g = Grass('mlpa_sc')
35 result = g.run('r.info soilmoisture')
36
37 # Django Settings:
38 # GRASS_GISBASE = '/usr/local/grass-6.4.1RC2'
39 # GRASS_GISDBASE = '/home/grass'
40
41 Example usage (outside django):
42 g = Grass('mlpa_sc', gisdbase="/home/grass", gisbase="/usr/local/grass-6.4.1RC2")
43 result = g.run('r.info soilmoisture')
44
45 '''
49 '''
50 If being used to load the source mapset then tmpMapset should be None or False.
51 If you're running any commands with outputs, use tmpMapset=True (autogenerates
52 tmpMapset name) or specify a unique name.
53 autoclean=True/False determines if tmpMapset gets deleted at Grass.__del__()
54 gisdbase and gisbase should be set in django settings or overriden with kwargs
55 '''
56 - def __init__(self, location, srcMapset='PERMANENT', tmpMapset=True,
57 gisdbase=None, gisbase=None, autoclean=True):
58 self.autoclean = autoclean
59 self.verbose = False
60
61 self.location = location
62 if tmpMapset == True:
63
64 self.tmpMapset = 'tmpmapset_%s_%s' % (srcMapset, str(time.time()))
65 elif tmpMapset:
66 self.tmpMapset = tmpMapset
67 self.srcMapset = srcMapset
68 self.templateMapset = 'mm_template_mapset'
69
70 if DJANGO and not gisdbase:
71 self.gisdbase = settings.GRASS_GISDBASE
72 else:
73 self.gisdbase = gisdbase
74 if DJANGO and not gisbase:
75 self.gisbase = settings.GRASS_GISBASE
76 else:
77 self.gisbase = gisbase
78
79 self.locationPath = os.path.join(self.gisdbase, self.location)
80 self.tmpMapsetPath = os.path.join(self.locationPath, self.tmpMapset)
81 self.srcMapsetPath = os.path.join(self.locationPath, self.srcMapset)
82 self.templatePath = os.path.join(self.locationPath, self.templateMapset)
83
84 if not tmpMapset:
85 self.curMapset = srcMapset
86 self.grassRcFile = os.path.join(gisdbase, "mm_grass6rc")
87 self.setupTmpGrassEnv(mktmp=False)
88 else:
89 self.curMapset = self.tmpMapset
90 self.grassRcFile = os.path.join(self.tmpMapsetPath, "mm_grass6rc")
91 self.setupTmpGrassEnv()
92
93
94
95
97 if self.autoclean:
98 self.cleanup()
99
101 if self.tmpMapset:
102 return u'Grass instance at temporary mapset %s (%s)' % (self.tmpMapset, self.gisdbase)
103 else:
104 return u'Grass instance at permanent mapset %s' % self.srcMapsetPath
105
106 @property
108 rev = self.run('g.version -r')
109 for line in rev.split("\n"):
110 parts = line.split(":")
111 parts = [p.strip() for p in parts]
112 if parts[0].lower() == "revision":
113 return int(parts[1])
114 return None
115
116 @property
118 try:
119 return settings.GRASS_TMP
120 except:
121 return '/tmp'
122
123 @property
125 try:
126 return settings.GRASS_PATH
127 except:
128 return '%(GRASS)s/bin:%(GRASS)s/scripts:/usr/local/sbin:/usr/local/bin:' \
129 '/usr/sbin:/usr/bin:/sbin:/bin' % {'GRASS': self.gisbase}
130
131 @property
133 try:
134 return settings.GRASS_VERSION
135 except:
136 try:
137 return self.gisbase.split('-')[1].replace('/','')
138 except:
139 return "6.4.1"
140
141 @property
143 try:
144 return settings.GRASS_LIB_PATH
145 except:
146 return '%s/lib' % self.gisbase
147
148 '''
149 Setup grass environment to use a temporary mapset as the current mapset.
150 Base maps can be used from the src mapset as it will be in your
151 mapset path, but results will go to the temporary mapset.
152 grassenv settings are system specific. To get these values start a
153 grass shell and echo the environment variables to get the necessary
154 values. eg. >grass64 ; >echo $GISBASE
155 '''
157 grassenv = {
158 'HOME': self.grass_tmp,
159 'PATH': self.grass_path,
160 'GRASS_VERSION': self.grass_version,
161 'GISBASE': self.gisbase,
162 'LD_LIBRARY_PATH': self.grass_lib_path,
163 'LOCATION_NAME': self.location,
164 'MAPSET': self.srcMapset,
165 'GISDBASE': self.gisdbase,
166 'GRASS_GUI': 'text',
167 'GIS_LOCK': str(os.getpid()),
168 'GISRC': self.grassRcFile
169 }
170
171 for key in grassenv.keys():
172 self.setEnv(key, grassenv[key])
173
174
175 if mktmp:
176 self.createTempMapset()
177
178 os.chdir(self.locationPath)
179
180 self.writeGrassRc()
181 return grassenv
182
183 '''
184 Create a temporary mapset
185 '''
187 if not os.path.exists(self.tmpMapsetPath):
188 if os.path.exists(self.templatePath):
189
190 shutil.copytree(self.templatePath, self.tmpMapsetPath)
191 else:
192
193 os.mkdir(self.tmpMapsetPath)
194 shutil.copyfile(
195 os.path.join(self.srcMapsetPath, 'DEFAULT_WIND'),
196 os.path.join(self.tmpMapsetPath, 'DEFAULT_WIND'))
197 shutil.copyfile(
198 os.path.join(self.srcMapsetPath, 'WIND'),
199 os.path.join(self.tmpMapsetPath, 'WIND'))
200 myname = open(os.path.join(self.tmpMapsetPath, 'MYNAME'),'w')
201 myname.write(self.tmpMapset)
202 myname.close()
203
204 '''
205 Cleanup disk such as temporary mapset
206 '''
208 if self.tmpMapset:
209 shutil.rmtree(self.tmpMapsetPath)
210
211 '''
212 Output grassrc file used by Grass. Set to use the current mapset
213 '''
215 gisrc = open(self.grassRcFile,"w")
216 gisrc.write("LOCATION_NAME: %s\n" % (self.location))
217 gisrc.write("MAPSET: %s\n" % (self.curMapset))
218 gisrc.write("DIGITIZER: none\n")
219 gisrc.write("GISDBASE: %s\n" % (self.gisdbase))
220 gisrc.write("OVERWRITE: 1\n")
221 gisrc.write("GRASS_GUI: text\n")
222 gisrc.close()
223 os.environ['GISRC'] = self.grassRcFile
224
225 '''
226 Set a single environment variable to the given value
227 '''
229 origValue = os.getenv(key)
230
231
232 os.putenv(key,value)
233 os.environ[key] = value
234 return
235
236 '''
237 Returns the raster path for the current mapset
238 '''
240 return self.gisdbase + self.location + '/' + self.curMapset + '/cell/'
241
242 - def run(self, cmd, nice=None):
243 if nice:
244 cmd = 'nice -n %d %s' % (nice,cmd)
245 if self.verbose:
246 log.debug(cmd)
247 proc = subprocess.Popen(cmd, shell=True,
248 stdout=subprocess.PIPE,
249 stderr=subprocess.PIPE,
250 env=self.setupTmpGrassEnv(),)
251 out, err = proc.communicate()
252 returncode = proc.returncode
253
254 if returncode != 0:
255 raise Exception("\nCommand failed with return code %s: \n %s \n %s" % (returncode, cmd, err))
256 elif err and self.verbose:
257 log.debug(err)
258 return out
259
260
261
262 '''
263 Copy a vector map from the base mapset to the temporary mapset
264 '''
266 command = "g.copy %s=%s@%s,%s" % (type, mapName, self.srcMapset, mapName)
267 self.run(command)
268
269 '''
270 Input a vector ogr datasource (input) and save as a grass vector map (output)
271 '''
273 command = 'v.in.ogr -o dsn=%s output=%s' % (input, output)
274 self.run(command)
275
276 '''
277 Input a raster ogr datasource (input) and save as a grass raster map (output)
278 '''
280 command = 'r.in.gdal -o input=%s output=%s' % (input, output)
281 self.run(command)
282
283 '''
284 Calculates the sum of all cell values where input is the map name
285 '''
287 value = 0.0
288 valueResult = self.run('r.sum ' + input)
289 valueList = re.split('SUM = ', valueResult)
290 if len(valueList) > 0:
291 value = float(valueList[1])
292 return value
293 else:
294 return None
295
296 '''
297 Calculate the area of a given raster (positive value cells) where input is the name
298 of the map layer and cell_size is the width of each square cell in the raster
299 '''
300 - def r_area(self, input, cell_size):
301
302 maskCmd = 'r.mapcalc "maskMap=if(%s, 1)"' % (input)
303 self.run(maskCmd, nice=1)
304
305
306 area = 0.0
307 cells = 0
308 areaResult = self.run('r.sum maskMap')
309 areaList = re.split('SUM = ', areaResult)
310 if len(areaList) > 0:
311 cells = float(areaList[1])
312 area = self.__areaOfCells(cells, cell_size)
313 return (cells, area)
314 else:
315 return None
316
317 '''
318 Intersect raster map m1 with m2, storing the cell value of m2 in the resulting raster
319 '''
321 intersectCmd = 'r.mapcalc "%s=if(%s,%s)"' % (result, m1, m2)
322 self.run(intersectCmd, nice=1)
323
324 '''
325 Converts vMap a vector map, to rMap a raster map. Where vectors overlap, raster cells
326 are given the value of val
327 '''
328 - def v_to_r(self, vMap, rMap, val):
329 command = 'v.to.rast use=val value=%s in=%s out=%s' % (val, vMap, rMap)
330 self.run(command, nice=1)
331
332 '''
333 Calculates area where numCells is the number of square cells composing the area
334 and cell_size is the width of the square cell
335 '''
337 return numCells * (cell_size * cell_size)
338
339 '''
340 Removes a grass map layer where type is the layer type 'vect' or 'rast' and layers
341 is the name of the map layer
342 '''
344 command = 'g.remove %s=%s' % (type,layer)
345 self.run(command)
346
347 '''
348 returns a dict of available raster and vector datasets as a python list
349 '''
351 cmd = 'g.mlist type=rast'
352 rasts = self.run(cmd)
353 cmd = 'g.mlist type=vect'
354 vects = self.run(cmd)
355 return {'rast': rasts.split(), 'vect': vects.split()}
356