-
Notifications
You must be signed in to change notification settings - Fork 10
/
Copy pathpygeonet_grass_py2.py
152 lines (148 loc) · 6.72 KB
/
pygeonet_grass_py2.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
import os
import sys
import shutil
import subprocess
from time import clock
import ConfigParser
#from pygeonet_rasterio import *
import grass.script as g
import grass.script.setup as gsetup
def grass(demFileName, geonetResultsDir, pmGrassGISfileName):
grass7bin = 'grass76'
if sys.platform.startswith('win'):
# MS Windows
grass7bin = r'C:\Program Files\GRASS GIS 7.6\grass76.bat'
# uncomment when using standalone WinGRASS installer
# grass7bin = r'C:\software\GRASS GIS 7.6\grass76.bat'
# this can be avoided if GRASS executable is added to PATH
elif sys.platform == 'darwin':
# Mac OS X
# TODO: this have to be checked, maybe unix way is good enough
grass7bin = r'C:\software\GRASS_GIS_7.6\grass76.bat'
startcmd = [grass7bin, '--config', 'path']
p = subprocess.Popen(startcmd, shell=False,
stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out, err = p.communicate()
if p.returncode != 0:
print >>sys.stderr, "ERROR: Cannot find GRASS GIS 7 " \
"start script (%s)" % startcmd
sys.exit(-1)
gisbase = out.strip('\n\r')
gisdb = os.path.join(os.path.expanduser("~"), "grassdata")
mswin = sys.platform.startswith('win')
if mswin:
gisdbdir = os.path.join(os.path.expanduser("~"), "Personal\grassdata")
else:
gisdbdir = os.path.join(os.path.expanduser("~"), "grassdata")
gisdbdir = r"C:\Users\\alecc\Documents\grassdata1"
geotiff = pmGrassGISfileName
locationGeonet = 'geonet'
grassGISlocation = os.path.join(gisdbdir, locationGeonet)
if os.path.exists(grassGISlocation):
print "Cleaning existing Grass location"
shutil.rmtree(grassGISlocation)
gsetup.init(gisbase, gisdbdir, locationGeonet, 'PERMANENT')
mapsetGeonet = 'geonetuser'
print 'Making the geonet location'
g.run_command('g.proj', georef=geotiff, location = locationGeonet)
location = locationGeonet
mapset = mapsetGeonet
print 'Existing Mapsets after making locations:'
g.read_command('g.mapsets', flags = 'l')
print 'Setting GRASSGIS environ'
gsetup.init(gisbase, gisdbdir, locationGeonet, 'PERMANENT')
## g.gisenv()
print 'Making mapset now'
g.run_command('g.mapset', flags = 'c', mapset = mapsetGeonet,\
location = locationGeonet, dbase = gisdbdir)
# gsetup initialization
gsetup.init(gisbase, gisdbdir, locationGeonet, mapsetGeonet)
# Read the filtered DEM
print 'Import filtered DEM into GRASSGIS and '\
'name the new layer with the DEM name'
tmpfile = demFileName
geotiffmapraster = tmpfile.split('.')[0]
print 'GRASSGIS layer name: ',geotiffmapraster
g.run_command('r.in.gdal', input=geotiff, \
output=geotiffmapraster,overwrite=True)
#Flow computation for massive grids (float version)
print "Calling the r.watershed command from GRASS GIS"
subbasinThreshold = 1500
print ('using swap memory option for large size DEM')
g.run_command('r.watershed',flags ='am',overwrite=True,\
elevation=geotiffmapraster, \
threshold=subbasinThreshold, \
drainage = 'dra1v23')
g.run_command('r.watershed',flags ='am',overwrite=True,\
elevation=geotiffmapraster, \
threshold=subbasinThreshold, \
accumulation='acc1v23')
print 'Identify outlets by negative flow direction'
g.run_command('r.mapcalc',overwrite=True,\
expression='outletmap = if(dra1v23 >= 0,null(),1)')
print 'Convert outlet raster to vector'
g.run_command('r.to.vect',overwrite=True,\
input = 'outletmap', output = 'outletsmapvec',\
type='point')
print 'Delineate basins according to outlets'
g.run_command('r.stream.basins',overwrite=True,\
direction='dra1v23',points='outletsmapvec',\
basins = 'outletbasins')
# Save the outputs as TIFs
outlet_filename = geotiffmapraster + '_outlets.tif'
g.run_command('r.out.gdal',overwrite=True,\
input='outletmap', type='Float32',\
output=os.path.join(geonetResultsDir,
outlet_filename),\
format='GTiff')
outputFAC_filename = geotiffmapraster + '_fac.tif'
g.run_command('r.out.gdal',overwrite=True,\
input='acc1v23', type='Float64',\
output=os.path.join(geonetResultsDir,
outputFAC_filename),\
format='GTiff')
outputFDR_filename = geotiffmapraster + '_fdr.tif'
g.run_command('r.out.gdal',overwrite=True,\
input = "dra1v23", type='Int32',\
output=os.path.join(geonetResultsDir,
outputFDR_filename),\
format='GTiff')
outputBAS_filename = geotiffmapraster + '_basins.tif'
g.run_command('r.out.gdal',overwrite=True,\
input = "outletbasins", type='Int16',\
output=os.path.join(geonetResultsDir,
outputBAS_filename),\
format='GTiff')
#def path_finder():
# config = ConfigParser.RawConfigParser()
# cfg =[]
# for dirpath, subdirs, files in os.walk(os.getcwd()):
# for y in subdirs:
# if y.endswith
# for x in files:
# if x.endswith("_blah.cfg"):
# cfg.append(os.path.join(dirpath, x))
# print(cfg)
# config.read('GeoFlood.cfg')
# demFileName=config.get('Section','dem_name')
# geoNetHomeDir = config.get('Section', 'geofloodhomedir')
# geonetResultsDir = os.path.join(geoNetHomeDir,"Outputs","GIS")
# pmGrassGISfileName = os.path.join(geonetResultsDir, "PM_filtered_grassgis.tif")
# return demFileName, geonetResultsDir, pmGrassGISfileName
def main():
abspath=os.path.abspath(__file__)
geonet_dir = os.path.dirname(os.path.dirname(abspath))
config = ConfigParser.RawConfigParser()
config.read(os.path.join(geonet_dir,"GeoFlood.cfg"))
projectname=config.get('Section','projectname')
dem_name=config.get('Section','dem_name')
demFileName=dem_name+'.tif'
geonetResultsDir=os.path.join(geonet_dir,"Outputs","GIS",projectname)
pmGrassGISfileName = os.path.join(geonetResultsDir, "PM_filtered_grassgis.tif")
grass(demFileName, geonetResultsDir, pmGrassGISfileName)
if __name__ == '__main__':
#demFileName, geonetResultsDir, pmGrassGISfileName=path_finder()
t0 = clock()
main()
t1 = clock()
print "time taken to complete flow accumulation:", t1-t0, " seconds"