source: grass-addons/grass7/raster/r.basin/r.basin.py

Last change on this file was 73968, checked in by neteler, 6 years ago

raster addons: Python3 compatibility changes, using 2to3 from python-tools
  • Property svn:eol-style set to native
  • Property svn:mime-type set to text/x-python
File size: 38.3 KB
Line 
1#!/usr/bin/env python
2
3############################################################################
4#
5# MODULE: r.basin
6# AUTHOR(S): Margherita Di Leo, Massimo Di Stefano
7# PURPOSE: Morphometric characterization of river basins
8# COPYRIGHT: (C) 2010-2014 by Margherita Di Leo & Massimo Di Stefano
9# dileomargherita@gmail.com
10#
11# This program is free software under the GNU General Public
12# License (>=v3.0) and comes with ABSOLUTELY NO WARRANTY.
13# See the file COPYING that comes with GRASS
14# for details.
15#
16# TODO: does r.stream.snap's snap depend on the raster resolution? hardcoded 30 below
17#
18#############################################################################
19
20#%module
21#% description: Morphometric characterization of river basins
22#% keyword: raster
23#% keyword: hydrology
24#% keyword: watershed
25#% overwrite: yes
26#%end
27
28#%option G_OPT_R_ELEV
29#% key: map
30#% description: Name of elevation raster map
31#% required: yes
32#%end
33
34#%option
35#% key: prefix
36#% type: string
37#% key_desc: prefix
38#% description: output prefix (must start with a letter)
39#% required: yes
40#%end
41
42#%option G_OPT_M_COORDS
43#% description: coordinates of the outlet (east,north)
44#% required : yes
45#%end
46
47#%option G_OPT_M_DIR
48#% key: dir
49#% description: Directory where the output will be found
50#% required : yes
51#%end
52
53#%option
54#% key: threshold
55#% type: double
56#% key_desc: threshold
57#% description: threshold
58#% required : no
59#%end
60
61#%flag
62#% key: a
63#% description: Use default threshold (1km^2)
64#%END
65
66#%flag
67#% key: c
68#% description: No maps output
69#%END
70
71import sys
72import os
73import grass.script as grass
74import math
75from numpy import zeros
76import csv
77
78# i18N
79import gettext
80gettext.install('grassmods', os.path.join(os.getenv("GISBASE"), 'locale'))
81
82# check requirements
83def check_progs():
84 found_missing = False
85 for prog in ('r.hypso', 'r.stream.basins', 'r.stream.distance', 'r.stream.extract',
86 'r.stream.order','r.stream.snap','r.stream.stats', 'r.width.funct'):
87 if not grass.find_program(prog, '--help'):
88 found_missing = True
89 grass.warning(_("'%s' required. Please install '%s' first using 'g.extension %s'") % (prog, prog, prog))
90 if found_missing:
91 grass.fatal(_("An ERROR occurred running r.basin"))
92
93def main():
94 # check dependencies
95 check_progs()
96
97 # check for unsupported locations
98 in_proj = grass.parse_command('g.proj', flags='g')
99 if in_proj['unit'].lower() == 'degree':
100 grass.fatal(_("Latitude-longitude locations are not supported"))
101 if in_proj['name'].lower() == 'xy_location_unprojected':
102 grass.fatal(_("xy-locations are not supported"))
103
104 r_elevation = options['map'].split('@')[0]
105 mapname = options['map'].replace("@"," ")
106 mapname = mapname.split()
107 mapname[0] = mapname[0].replace(".","_")
108 coordinates = options['coordinates']
109 directory = options['dir']
110 # Check if directory exists
111 if not os.path.isdir(directory):
112 os.makedirs(directory)
113 autothreshold = flags['a']
114 nomap = flags['c']
115 prefix = options['prefix']+'_'+mapname[0]
116 r_accumulation = prefix+'_accumulation'
117 r_drainage = prefix+'_drainage'
118 r_stream = prefix+'_stream'
119 r_slope = prefix+'_slope'
120 r_aspect = prefix+'_aspect'
121 r_basin = prefix+'_basin'
122 r_strahler = prefix+'_strahler'
123 r_shreve = prefix+'_shreve'
124 r_horton = prefix+'_horton'
125 r_hack = prefix+'_hack'
126 r_distance = prefix+'_dist2out'
127 r_hillslope_distance = prefix+'_hillslope_distance'
128 r_height_average = prefix+'_height_average'
129 r_aspect_mod = prefix+'_aspect_mod'
130 r_dtm_basin = prefix+'_dtm_basin'
131 r_mainchannel = prefix+'_mainchannel'
132 r_stream_e = prefix+'_stream_e'
133 r_drainage_e = prefix+'_drainage_e'
134 r_mask = prefix+'_mask'
135 r_ord_1 = prefix+'_ord_1'
136 r_average_hillslope = prefix+'_average_hillslope'
137 r_mainchannel_dim = prefix+'_mainchannel_dim'
138 r_outlet = prefix+'_r_outlet'
139 v_outlet = prefix+'_outlet'
140 v_outlet_snap = prefix+'_outlet_snap'
141 v_basin = prefix+'_basin'
142 v_mainchannel = prefix+'_mainchannel'
143 v_mainchannel_dim = prefix+'_mainchannel_dim'
144 v_network = prefix+'_network'
145 v_ord_1 = prefix+'_ord_1'
146 global tmp
147
148
149 # Save current region
150 grass.read_command('g.region', flags = 'p', save = 'original')
151
152 # Watershed SFD
153 grass.run_command('r.watershed', elevation = r_elevation,
154 accumulation = r_accumulation,
155 drainage = r_drainage,
156 convergence = 5,
157 flags = 'am')
158
159 # Managing flag
160 if autothreshold :
161 resolution = grass.region()['nsres']
162 th = 1000000 / (resolution**2)
163 grass.message( "threshold : %s" % th )
164 else :
165 th = options['threshold']
166
167 # Stream extraction
168 grass.run_command('r.stream.extract', elevation = r_elevation,
169 accumulation = r_accumulation,
170 threshold = th,
171 d8cut = 1000000000,
172 mexp = 0,
173 stream_rast = r_stream_e,
174 direction = r_drainage_e)
175
176 try:
177 # Delineation of basin
178 # Create outlet
179 grass.write_command('v.in.ascii', output = v_outlet,
180 input = "-",
181 sep = ",",
182 stdin = "%s,9999" % (coordinates))
183
184 # Snap outlet to stream network
185 # TODO: does snap depend on the raster resolution? hardcoded 30 below
186 grass.run_command('r.stream.snap', input = v_outlet,
187 output = v_outlet_snap,
188 stream_rast = r_stream_e,
189 radius = 30)
190
191 grass.run_command('v.to.rast', input = v_outlet_snap,
192 output = r_outlet,
193 use = 'cat',
194 type = 'point',
195 layer = 1,
196 value = 1)
197
198
199 grass.run_command('r.stream.basins', direction = r_drainage_e,
200 basins = r_basin,
201 points = v_outlet_snap)
202
203 grass.message( "Delineation of basin done" )
204
205 # Mask and cropping
206 elevation_name = r_elevation = r_elevation.split('@')[0]
207
208 grass.mapcalc("$r_mask = $r_basin / $r_basin",
209 r_mask = r_mask,
210 r_basin = r_basin)
211
212 grass.mapcalc("tmp = $r_accumulation / $r_mask",
213 r_accumulation = r_accumulation,
214 r_mask = r_mask)
215
216 grass.run_command('g.remove', flags='f', type='raster', name= r_accumulation, quiet = True)
217
218 grass.run_command('g.rename', raster = ('tmp',r_accumulation))
219
220 grass.mapcalc("tmp = $r_drainage / $r_mask",
221 r_drainage = r_drainage,
222 r_mask = r_mask)
223
224 grass.run_command('g.remove', flags='f', type='raster', name= r_drainage, quiet = True)
225
226 grass.run_command('g.rename', raster = ('tmp', r_drainage))
227
228 grass.mapcalc("$r_elevation_crop = $r_elevation * $r_mask",
229 r_mask = r_mask,
230 r_elevation = r_elevation,
231 r_elevation_crop = 'r_elevation_crop')
232
233 grass.mapcalc("tmp = $r_drainage_e * $r_mask",
234 r_mask = r_mask,
235 r_drainage_e = r_drainage_e)
236
237 grass.run_command('g.remove', flags='f', type='raster', name= r_drainage_e, quiet = True)
238
239 grass.run_command('g.rename', raster = ('tmp',r_drainage_e))
240
241 grass.mapcalc("tmp = $r_stream_e * $r_mask",
242 r_mask = r_mask,
243 r_stream_e = r_stream_e)
244
245 grass.run_command('g.remove', flags='f', type='raster', name= r_stream_e, quiet = True)
246 #grass.run_command('g.rename', raster = (r_stream_e,'streams'))
247
248 grass.run_command('g.rename', raster = ('tmp',r_stream_e))
249
250 grass.run_command('r.thin', input = r_stream_e,
251 output = r_stream_e+'_thin')
252
253 grass.run_command('r.to.vect', input = r_stream_e+'_thin',
254 output = v_network,
255 type = 'line')
256
257 # Creation of slope and aspect maps
258 grass.run_command('r.slope.aspect', elevation = 'r_elevation_crop',
259 slope = r_slope,
260 aspect = r_aspect)
261
262 # Basin mask (vector)
263 # Raster to vector
264 grass.run_command('r.to.vect', input = r_basin,
265 output = v_basin,
266 type = 'area',
267 flags = 'sv')
268
269 # Add two columns to the table: area and perimeter
270 grass.run_command('v.db.addcolumn', map = v_basin,
271 columns = 'area double precision')
272
273 grass.run_command('v.db.addcolumn', map = v_basin,
274 columns = 'perimeter double precision')
275
276 # Populate perimeter column
277 grass.run_command('v.to.db', map = v_basin,
278 type = 'line,boundary',
279 layer = 1,
280 qlayer = 1,
281 option = 'perimeter',
282 units = 'kilometers',
283 columns = 'perimeter')
284
285 # Read perimeter
286 tmp = grass.read_command('v.to.db', map = v_basin,
287 type = 'line,boundary',
288 layer = 1,
289 qlayer = 1,
290 option = 'perimeter',
291 units = 'kilometers',
292 qcolumn = 'perimeter',
293 flags = 'p')
294 perimeter_basin = float(tmp.split('\n')[1].split('|')[1])
295
296 # Populate area column
297 grass.run_command('v.to.db', map = v_basin,
298 type = 'line,boundary',
299 layer = 1,
300 qlayer = 1,
301 option = 'area',
302 units = 'kilometers',
303 columns = 'area')
304
305 # Read area
306 tmp = grass.read_command('v.to.db', map = v_basin,
307 type = 'line,boundary',
308 layer = 1,
309 qlayer = 1,
310 option = 'area',
311 units = 'kilometers',
312 qcolumn = 'area',
313 flags = 'p')
314 area_basin = float(tmp.split('\n')[1].split('|')[1])
315
316 # Creation of order maps: strahler, horton, hack, shreeve
317 grass.message( "Creating %s" % r_hack )
318
319 grass.run_command('r.stream.order', stream_rast = r_stream_e,
320 direction = r_drainage_e,
321 strahler = r_strahler,
322 shreve = r_shreve,
323 horton = r_horton,
324 hack = r_hack)
325
326 # Distance to outlet
327 grass.run_command('r.stream.distance', stream_rast = r_outlet,
328 direction = r_drainage_e,
329 flags = 'o',
330 distance = r_distance)
331
332
333 # hypsographic curve
334
335 grass.message( "------------------------------" )
336
337 grass.run_command('r.hypso', map = 'r_elevation_crop',
338 image = os.path.join(directory,prefix), flags = 'ab')
339
340 grass.message( "------------------------------" )
341
342 # Width Function
343
344 grass.message( "------------------------------" )
345
346 grass.run_command('r.width.funct', map = r_distance,
347 image = os.path.join(directory,prefix))
348
349 grass.message( "------------------------------" )
350
351 # Creation of map of hillslope distance to river network
352
353 grass.run_command("r.stream.distance", stream_rast = r_stream_e,
354 direction = r_drainage,
355 elevation = 'r_elevation_crop',
356 distance = r_hillslope_distance)
357
358
359 # Mean elevation
360 grass.run_command("r.stats.zonal", base = r_basin,
361 cover = "r_elevation_crop",
362 method = "average",
363 output = r_height_average)
364
365
366 grass.message("r.stats.zonal done")
367 mean_elev = float(grass.read_command('r.info', flags = 'r',
368 map = r_height_average).split('\n')[0].split('=')[1])
369 grass.message("r.info done")
370
371
372 # In Grass, aspect categories represent the number degrees of east and they increase
373 # counterclockwise: 90deg is North, 180 is West, 270 is South 360 is East.
374 # The aspect value 0 is used to indicate undefined aspect in flat areas with slope=0.
375 # We calculate the number of degree from north, increasing counterclockwise.
376 grass.mapcalc("$r_aspect_mod = if($r_aspect == 0, 0, if($r_aspect > 90, 450 - $r_aspect, 90 - $r_aspect))",
377 r_aspect = r_aspect,
378 r_aspect_mod = r_aspect_mod)
379 grass.message("r.mapcalc done")
380
381 # Centroid and mean slope
382 baricenter_slope_baricenter = grass.read_command("r.volume", input = r_slope,
383 clump = r_basin)
384
385 grass.message("r.volume done")
386
387 baricenter_slope_baricenter = baricenter_slope_baricenter.split()
388 mean_slope = baricenter_slope_baricenter[30]
389
390 # Rectangle containing basin
391 basin_east = baricenter_slope_baricenter[33]
392 basin_north = baricenter_slope_baricenter[34]
393 info_region_basin = grass.read_command("g.region",
394 vect = options['prefix']+'_'+mapname[0]+'_basin',
395 flags = 'm')
396
397 grass.message("g.region done")
398 dict_region_basin = dict(x.split('=', 1) for x in info_region_basin.split('\n') if '=' in x)
399 basin_resolution = float(dict_region_basin['nsres'])
400# x_massimo = float(dict_region_basin['n']) + (basin_resolution * 10)
401# x_minimo = float(dict_region_basin['w']) - (basin_resolution * 10)
402# y_massimo = float(dict_region_basin['e']) + (basin_resolution * 10)
403# y_minimo = float(dict_region_basin['s']) - (basin_resolution * 10)
404 nw = dict_region_basin['w'], dict_region_basin['n']
405 se = dict_region_basin['e'], dict_region_basin['s']
406 grass.message("Rectangle containing basin done")
407
408 east1,north1 = coordinates.split(',')
409 east = float(east1)
410 north = float(north1)
411
412 # Directing vector
413 delta_x = abs(float(basin_east) - east)
414 delta_y = abs(float(basin_north) - north)
415 L_orienting_vect = math.sqrt((delta_x**2)+(delta_y**2)) / 1000
416 grass.message("Directing vector done")
417
418 # Prevalent orientation
419 prevalent_orientation = math.atan(delta_y/delta_x)
420 grass.message("Prevalent orientation done")
421
422 # Compactness coefficient
423 C_comp = perimeter_basin / ( 2 * math.sqrt( area_basin / math.pi))
424 grass.message("Compactness coefficient done")
425
426 # Circularity ratio
427 R_c = ( 4 * math.pi * area_basin ) / (perimeter_basin **2)
428 grass.message("Circularity ratio done")
429
430 # Mainchannel
431 grass.mapcalc("$r_mainchannel = if($r_hack==1,1,null())",
432 r_hack = r_hack,
433 r_mainchannel = r_mainchannel)
434
435 grass.run_command("r.thin", input = r_mainchannel,
436 output = r_mainchannel+'_thin')
437 grass.run_command('r.to.vect', input = r_mainchannel+'_thin',
438 output = v_mainchannel,
439 type = 'line',
440 verbose = True)
441
442
443 # Get coordinates of the outlet (belonging to stream network)
444
445 grass.run_command('v.db.addtable', map = v_outlet_snap)
446
447 grass.run_command('v.db.addcolumn', map = v_outlet_snap,
448 columns="x double precision,y double precision" )
449
450 grass.run_command('v.to.db', map = v_outlet_snap,
451 option = "coor",
452 col = "x,y" )
453
454 namefile = os.path.join(directory, prefix + '_outlet_coors.txt')
455
456 grass.run_command('v.out.ascii', input = v_outlet_snap,
457 output = namefile,
458 cats = 1,
459 format = "point")
460
461 f = open(namefile)
462 east_o, north_o, cat = f.readline().split('|')
463
464 param_mainchannel = grass.read_command('v.what', map = v_mainchannel,
465 coordinates = '%s,%s' % (east_o,north_o),
466 distance = 5)
467 tmp = param_mainchannel.split('\n')[7]
468 mainchannel = float(tmp.split()[1]) / 1000 # km
469
470 # Topological Diameter
471 grass.mapcalc("$r_mainchannel_dim = -($r_mainchannel - $r_shreve) + 1",
472 r_mainchannel_dim = r_mainchannel_dim,
473 r_shreve = r_shreve,
474 r_mainchannel = r_mainchannel)
475 grass.run_command('r.thin', input = r_mainchannel_dim,
476 output = r_mainchannel_dim + '_thin')
477 grass.run_command('r.to.vect', input = r_mainchannel_dim + '_thin',
478 output = v_mainchannel_dim,
479 type = 'line',
480 flags = 'v',
481 verbose = True)
482 try:
483 D_topo1 = grass.read_command('v.info', map = v_mainchannel_dim,
484 layer = 1,
485 flags = 't')
486 D_topo = float(D_topo1.split('\n')[2].split('=')[1])
487 except:
488 D_topo = 1
489 grass.message( "Topological Diameter = WARNING" )
490
491 # Mean slope of mainchannel
492 grass.message("doing v.to.points")
493 grass.run_command('v.to.points',
494 input = v_mainchannel_dim,
495 output = v_mainchannel_dim+'_point',
496 type = 'line')
497 vertex = grass.read_command('v.out.ascii', verbose = True,
498 input = v_mainchannel_dim+'_point').strip().split('\n')
499 nodi = zeros((len(vertex),4),float)
500 pendenze = []
501
502 for i in range(len(vertex)):
503 x, y = float(vertex[i].split('|')[0]) , float(vertex[i].split('|')[1])
504 vertice1 = grass.read_command('r.what', verbose = True,
505 map = 'r_elevation_crop',
506 coordinates = '%s,%s' % (x,y))
507 vertice = vertice1.replace('\n','').replace('||','|').split('|')
508 nodi[i,0],nodi[i,1], nodi[i,2] = float(vertice[0]), float(vertice[1]), float(vertice[2])
509
510 for i in range(0,len(vertex)-1,2):
511 dist = math.sqrt(math.fabs((nodi[i,0] - nodi[i+1,0]))**2 + math.fabs((nodi[i,1] - nodi[i+1,1]))**2)
512 deltaz = math.fabs(nodi[i,2] - nodi[i+1,2])
513 # Control to prevent float division by zero (dist=0)
514
515 try:
516 pendenza = deltaz / dist
517 pendenze.append(pendenza)
518 mainchannel_slope = sum(pendenze) / len(pendenze) * 100
519 except :
520 pass
521
522 # Elongation Ratio
523 R_al = (2 * math.sqrt( area_basin / math.pi) ) / mainchannel
524
525 # Shape factor
526 S_f = area_basin / mainchannel
527
528 # Characteristic altitudes
529 height_basin_average = grass.read_command('r.what', map = r_height_average ,
530 cache = 500 ,
531 coordinates = '%s,%s' % (east_o , north_o ))
532 height_basin_average = height_basin_average.replace('\n','')
533 height_basin_average = float(height_basin_average.split('|')[-1])
534 minmax_height_basin = grass.read_command('r.info', flags = 'r',
535 map = 'r_elevation_crop')
536 minmax_height_basin = minmax_height_basin.strip().split('\n')
537 min_height_basin, max_height_basin = float(minmax_height_basin[0].split('=')[-1]), float(minmax_height_basin[1].split('=')[-1])
538 H1 = max_height_basin
539 H2 = min_height_basin
540 HM = H1 - H2
541
542 # Concentration time (Giandotti, 1934)
543 t_c = ((4 * math.sqrt(area_basin)) + (1.5 * mainchannel)) / (0.8 * math.sqrt(HM))
544
545 # Mean hillslope length
546 grass.run_command("r.stats.zonal", cover = r_stream_e,
547 base = r_mask,
548 method = "average",
549 output = r_average_hillslope)
550 mean_hillslope_length = float(grass.read_command('r.info', flags = 'r',
551 map = r_average_hillslope).split('\n')[0].split('=')[1])
552
553 # Magnitude
554 grass.mapcalc("$r_ord_1 = if($r_strahler==1,1,null())",
555 r_ord_1 = r_ord_1,
556 r_strahler = r_strahler)
557 grass.run_command('r.thin', input = r_ord_1,
558 output = r_ord_1+'_thin',
559 iterations = 200)
560 grass.run_command('r.to.vect', input = r_ord_1+'_thin',
561 output = v_ord_1,
562 type = 'line',
563 flags = 'v')
564 magnitudo = float(grass.read_command('v.info', map = v_ord_1,
565 layer = 1,
566 flags = 't').split('\n')[2].split('=')[1])
567
568 # First order stream frequency
569 FSF = magnitudo / area_basin
570
571 # Statistics
572
573 stream_stats = grass.read_command('r.stream.stats', stream_rast = r_strahler,
574 direction = r_drainage_e,
575 elevation = 'r_elevation_crop' )
576
577
578 print(" ------------------------------ ")
579 print("Output of r.stream.stats: ")
580 print(stream_stats)
581
582 stream_stats_summary = stream_stats.split('\n')[4].split('|')
583 stream_stats_mom = stream_stats.split('\n')[8].split('|')
584 Max_order , Num_streams , Len_streams , Stream_freq = stream_stats_summary[0] , stream_stats_summary[1] , stream_stats_summary[2] , stream_stats_summary[5]
585 Bif_ratio , Len_ratio , Area_ratio , Slope_ratio = stream_stats_mom[0] , stream_stats_mom[1] , stream_stats_mom[2] , stream_stats_mom[3]
586 drainage_density = float(Len_streams) / float(area_basin)
587
588 # Cleaning up
589 grass.run_command('g.remove', flags='f', type='raster', name= 'r_elevation_crop', quiet = True)
590 grass.run_command('g.remove', flags='f', type='raster', name= r_height_average, quiet = True)
591 grass.run_command('g.remove', flags='f', type='raster', name= r_aspect_mod, quiet = True)
592 grass.run_command('g.remove', flags='f', type='raster', name= r_mainchannel, quiet = True)
593 grass.run_command('g.remove', flags='f', type='raster', name= r_stream_e, quiet = True)
594 grass.run_command('g.remove', flags='f', type='raster', name= r_drainage_e, quiet = True)
595 grass.run_command('g.remove', flags='f', type='raster', name= r_mask, quiet = True)
596 grass.run_command('g.remove', flags='f', type='raster', name= r_ord_1, quiet = True)
597 grass.run_command('g.remove', flags='f', type='raster', name= r_average_hillslope, quiet = True)
598 grass.run_command('g.remove', flags='f', type='raster', name= r_mainchannel_dim, quiet = True)
599 grass.run_command('g.remove', flags='f', type='raster', name= r_outlet, quiet = True)
600 grass.run_command('g.remove', flags='f', type='raster', name= r_basin, quiet = True)
601 grass.run_command('g.remove', flags='f', type='raster', name= prefix+'_mainchannel_thin', quiet = True)
602 grass.run_command('g.remove', flags='f', type='raster', name= prefix+'_mainchannel_dim_thin', quiet = True)
603 grass.run_command('g.remove', flags='f', type='raster', name= prefix+'_ord_1_thin', quiet = True)
604 grass.run_command('g.remove', flags='f', type='raster', name= prefix+'_stream_e_thin', quiet = True)
605 grass.run_command('g.remove', flags='f', type='vector', name= v_mainchannel_dim+'_point', quiet = True)
606 grass.run_command('g.remove', flags='f', type='vector', name= v_mainchannel_dim, quiet = True)
607 grass.run_command('g.remove', flags='f', type='vector', name= v_ord_1, quiet = True)
608
609 if nomap :
610 grass.run_command('g.remove', flags='f', type='vector', name= v_outlet, quiet = True)
611 grass.run_command('g.remove', flags='f', type='vector', name= v_basin, quiet = True)
612 grass.run_command('g.remove', flags='f', type='vector', name= v_mainchannel, quiet = True)
613 grass.run_command('g.remove', flags='f', type='raster', name= r_accumulation, quiet = True)
614 grass.run_command('g.remove', flags='f', type='raster', name= r_drainage, quiet = True)
615 grass.run_command('g.remove', flags='f', type='raster', name= r_aspect, quiet = True)
616 grass.run_command('g.remove', flags='f', type='raster', name= r_strahler, quiet = True)
617 grass.run_command('g.remove', flags='f', type='raster', name= r_shreve, quiet = True)
618 grass.run_command('g.remove', flags='f', type='raster', name= r_horton, quiet = True)
619 grass.run_command('g.remove', flags='f', type='raster', name= r_hack, quiet = True)
620 grass.run_command('g.remove', flags='f', type='raster', name= r_distance, quiet = True)
621 grass.run_command('g.remove', flags='f', type='raster', name= r_hillslope_distance, quiet = True)
622 grass.run_command('g.remove', flags='f', type='raster', name= r_slope, quiet = True)
623
624 ####################################################
625
626 parametri_bacino = {}
627 parametri_bacino["mean_slope"] = float(mean_slope)
628 parametri_bacino["mean_elev"] = float(mean_elev)
629 parametri_bacino["basin_east"] = float(basin_east)
630 parametri_bacino["basin_north"] = float(basin_north)
631 parametri_bacino["basin_resolution"] = float(basin_resolution)
632 parametri_bacino["nw"] = nw
633 parametri_bacino["se"] = se
634 parametri_bacino["area_basin"] = float(area_basin)
635 parametri_bacino["perimeter_basin"] = float(perimeter_basin)
636 parametri_bacino["L_orienting_vect"] = float(L_orienting_vect)
637 parametri_bacino["prevalent_orientation"] = float(prevalent_orientation)
638 parametri_bacino["C_comp"] = float(C_comp)
639 parametri_bacino["R_c"] = float(R_c)
640 parametri_bacino["mainchannel"] = float(mainchannel)
641 parametri_bacino["D_topo"] = float(D_topo)
642 parametri_bacino["mainchannel_slope" ] = float(mainchannel_slope)
643 parametri_bacino["R_al"] = float(R_al)
644 parametri_bacino["S_f"] = float(S_f)
645 parametri_bacino["H1"] = float(H1)
646 parametri_bacino["H2"] = float(H2)
647 parametri_bacino["HM"] = float(HM)
648 parametri_bacino["t_c"] = float(t_c)
649 parametri_bacino["mean_hillslope_length"] = float(mean_hillslope_length)
650 parametri_bacino["magnitudo"] = float(magnitudo)
651 parametri_bacino["Max_order"] = float(Max_order)
652 parametri_bacino["Num_streams"] = float(Num_streams)
653 parametri_bacino["Len_streams"] = float(Len_streams)
654 parametri_bacino["Stream_freq"] = float(Stream_freq)
655 parametri_bacino["Bif_ratio"] = float(Bif_ratio)
656 parametri_bacino["Len_ratio"] = float(Len_ratio)
657 parametri_bacino["Area_ratio"] = float(Area_ratio)
658 parametri_bacino["Slope_ratio"] = float(Slope_ratio)
659 parametri_bacino["drainage_density"] = float(drainage_density)
660 parametri_bacino["FSF"] = float(FSF)
661
662
663 # create .csv file
664 csvfile = os.path.join( directory, prefix + '_parameters.csv' )
665 with open(csvfile, 'w') as f:
666 writer = csv.writer(f)
667 writer.writerow(['Morphometric parameters of basin:'])
668 writer.writerow([' '])
669 writer.writerow(['Easting Centroid of basin'] + [basin_east])
670 writer.writerow(['Northing Centroid of basin'] + [basin_north])
671 writer.writerow(['Rectangle containing basin N-W'] + [nw])
672 writer.writerow(['Rectangle containing basin S-E'] + [se])
673 writer.writerow(['Area of basin [km^2]'] + [area_basin])
674 writer.writerow(['Perimeter of basin [km]'] + [perimeter_basin])
675 writer.writerow(['Max Elevation [m s.l.m.]'] + [H1])
676 writer.writerow(['Min Elevation [m s.l.m.]'] + [H2])
677 writer.writerow(['Elevation Difference [m]'] + [HM])
678 writer.writerow(['Mean Elevation'] + [mean_elev])
679 writer.writerow(['Mean Slope'] + [mean_slope])
680 writer.writerow(['Length of Directing Vector [km]'] + [L_orienting_vect])
681 writer.writerow(['Prevalent Orientation [degree from north, counterclockwise]'] + [prevalent_orientation])
682 writer.writerow(['Compactness Coefficient'] + [C_comp])
683 writer.writerow(['Circularity Ratio'] + [R_c])
684 writer.writerow(['Topological Diameter'] + [D_topo])
685 writer.writerow(['Elongation Ratio'] + [R_al])
686 writer.writerow(['Shape Factor'] + [S_f])
687 writer.writerow(['Concentration Time (Giandotti, 1934) [hr]'] + [t_c])
688 writer.writerow(['Length of Mainchannel [km]'] + [mainchannel])
689 writer.writerow(['Mean slope of mainchannel [percent]'] + [mainchannel_slope])
690 writer.writerow(['Mean hillslope length [m]'] + [mean_hillslope_length])
691 writer.writerow(['Magnitudo'] + [magnitudo])
692 writer.writerow(['Max order (Strahler)'] + [Max_order])
693 writer.writerow(['Number of streams'] + [Num_streams])
694 writer.writerow(['Total Stream Length [km]'] + [Len_streams])
695 writer.writerow(['First order stream frequency'] + [FSF])
696 writer.writerow(['Drainage Density [km/km^2]'] + [drainage_density])
697 writer.writerow(['Bifurcation Ratio (Horton)'] + [Bif_ratio])
698 writer.writerow(['Length Ratio (Horton)'] + [Len_ratio])
699 writer.writerow(['Area ratio (Horton)'] + [Area_ratio])
700 writer.writerow(['Slope ratio (Horton)'] + [Slope_ratio])
701
702 # Create summary (transposed)
703 csvfileT = os.path.join( directory, prefix + '_parametersT.csv' ) # transposed
704 with open(csvfileT, 'w') as f:
705 writer = csv.writer(f)
706 writer.writerow(['x'] +
707 ['y'] +
708 ['Easting_Centroid_basin'] +
709 ['Northing_Centroid_basin'] +
710 ['Rectangle_containing_basin_N_W'] +
711 ['Rectangle_containing_basin_S_E'] +
712 ['Area_of_basin_km2'] +
713 ['Perimeter_of_basin_km'] +
714 ['Max_Elevation'] +
715 ['Min_Elevation'] +
716 ['Elevation_Difference'] +
717 ['Mean_Elevation'] +
718 ['Mean_Slope'] +
719 ['Length_of_Directing_Vector_km'] +
720 ['Prevalent_Orientation_deg_from_north_ccw'] +
721 ['Compactness_Coefficient'] +
722 ['Circularity_Ratio'] +
723 ['Topological_Diameter'] +
724 ['Elongation_Ratio'] +
725 ['Shape_Factor'] +
726 ['Concentration_Time_hr'] +
727 ['Length_of_Mainchannel_km'] +
728 ['Mean_slope_of_mainchannel_percent'] +
729 ['Mean_hillslope_length_m'] +
730 ['Magnitudo'] +
731 ['Max_order_Strahler'] +
732 ['Number_of_streams'] +
733 ['Total_Stream_Length_km'] +
734 ['First_order_stream_frequency'] +
735 ['Drainage_Density_km_over_km2'] +
736 ['Bifurcation_Ratio_Horton'] +
737 ['Length_Ratio_Horton'] +
738 ['Area_ratio_Horton'] +
739 ['Slope_ratio_Horton'] )
740 writer.writerow([east_o]
741 + [north_o]
742 + [basin_east]
743 + [basin_north]
744 + [nw]
745 + [se]
746 + [area_basin]
747 + [perimeter_basin]
748 + [H1]
749 + [H2]
750 + [HM]
751 + [mean_elev]
752 + [mean_slope]
753 + [L_orienting_vect]
754 + [prevalent_orientation]
755 + [C_comp]
756 + [R_c]
757 + [D_topo]
758 + [R_al]
759 + [S_f]
760 + [t_c]
761 + [mainchannel]
762 + [mainchannel_slope]
763 + [mean_hillslope_length]
764 + [magnitudo]
765 + [Max_order]
766 + [Num_streams]
767 + [Len_streams]
768 + [FSF]
769 + [drainage_density]
770 + [Bif_ratio]
771 + [Len_ratio]
772 + [Area_ratio]
773 + [Slope_ratio])
774
775
776 # Import table "rbasin_summary", joins it to "outlet_snap", then drops it
777 grass.message("db.in.ogr: importing CSV table <%s>..." % csvfileT)
778 grass.run_command("db.in.ogr", input = csvfileT,
779 output = "rbasin_summary")
780
781 grass.run_command("v.db.join", map = v_outlet_snap,
782 otable = "rbasin_summary",
783 column = "y",
784 ocolumn = "y")
785 grass.run_command("db.droptable", table = "rbasin_summary", flags = 'f')
786
787 grass.message( "\n" )
788 grass.message( "----------------------------------" )
789 grass.message( "Morphometric parameters of basin :" )
790 grass.message( "----------------------------------\n" )
791 grass.message( "Easting Centroid of basin : %s " % basin_east )
792 grass.message( "Northing Centroid of Basin : %s " % basin_north )
793 grass.message( "Rectangle containing basin N-W : %s , %s " % nw )
794 grass.message( "Rectangle containing basin S-E : %s , %s " % se )
795 grass.message( "Area of basin [km^2] : %s " % area_basin )
796 grass.message( "Perimeter of basin [km] : %s " % perimeter_basin )
797 grass.message( "Max Elevation [m s.l.m.] : %s " % H1 )
798 grass.message( "Min Elevation [m s.l.m.]: %s " % H2 )
799 grass.message( "Elevation Difference [m]: %s " % HM )
800 grass.message( "Mean Elevation [m s.l.m.]: %s " % mean_elev )
801 grass.message( "Mean Slope : %s " % mean_slope )
802 grass.message( "Length of Directing Vector [km] : %s " % L_orienting_vect )
803 grass.message( "Prevalent Orientation [degree from north, counterclockwise] : %s " % prevalent_orientation )
804 grass.message( "Compactness Coefficient : %s " % C_comp )
805 grass.message( "Circularity Ratio : %s " % R_c )
806 grass.message( "Topological Diameter : %s " % D_topo )
807 grass.message( "Elongation Ratio : %s " % R_al )
808 grass.message( "Shape Factor : %s " % S_f )
809 grass.message( "Concentration Time (Giandotti, 1934) [hr] : %s " % t_c )
810 grass.message( "Length of Mainchannel [km] : %s " % mainchannel )
811 grass.message( "Mean slope of mainchannel [percent] : %f " % mainchannel_slope )
812 grass.message( "Mean hillslope length [m] : %s " % mean_hillslope_length )
813 grass.message( "Magnitudo : %s " % magnitudo )
814 grass.message( "Max order (Strahler) : %s " % Max_order )
815 grass.message( "Number of streams : %s " % Num_streams )
816 grass.message( "Total Stream Length [km] : %s " % Len_streams )
817 grass.message( "First order stream frequency : %s " % FSF )
818 grass.message( "Drainage Density [km/km^2] : %s " % drainage_density )
819 grass.message( "Bifurcation Ratio (Horton) : %s " % Bif_ratio )
820 grass.message( "Length Ratio (Horton) : %s " % Len_ratio )
821 grass.message( "Area ratio (Horton) : %s " % Area_ratio )
822 grass.message( "Slope ratio (Horton): %s " % Slope_ratio )
823 grass.message( "------------------------------" )
824 grass.message( "\n" )
825 grass.message( "Done!")
826
827 except:
828 grass.message( "\n" )
829 grass.message( "------------------------------" )
830 grass.message( "\n" )
831 grass.message( "An ERROR occurred running r.basin" )
832 grass.message( "Please check for error messages above or try with another pairs of outlet coordinates" )
833
834
835 # Set region to original
836 grass.read_command('g.region', flags = 'p', region = 'original')
837 grass.run_command('g.remove', flags = 'f', type = 'region', name = 'original')
838
839if __name__ == "__main__":
840 options, flags = grass.parser()
841 sys.exit(main())
Note: See TracBrowser for help on using the repository browser.