hrudelin_1_init.py 16.8 KB
Newer Older
1
# coding: utf-8
christine.barachet's avatar
   
christine.barachet committed
2
3
4
5
6
'''

 MODULE:       hru-delin_init

 AUTHOR(S):    adapted from GRASS-HRU (ILMS) - JENA University
7
               by IRSTEA - Christine Barachet
christine.barachet's avatar
   
christine.barachet committed
8
9
10
11
12
13
14
15

 PURPOSE:      Prepare files for the next steps of HRU Delineation
                   shapefile : selected gauges
                   rasters :   bounded DEM and others (geology, soils, landuse)
                               DEM, slope and aspect with values of reclassification
                               calculates oriented flows and subbasins

'''
christine.barachet's avatar
   
christine.barachet committed
16

Julien Veyssier's avatar
Julien Veyssier committed
17
18
# to keep python2 compatibility
from __future__ import print_function
19
import string, os, shutil, sys, glob, types
christine.barachet's avatar
   
christine.barachet committed
20
import math
Julien Veyssier's avatar
Julien Veyssier committed
21
22
try:
    import ConfigParser
23
except Exception as e:
Julien Veyssier's avatar
Julien Veyssier committed
24
    import configparser as ConfigParser
christine.barachet's avatar
   
christine.barachet committed
25
import grass.script as grass
26
from grass.script.utils import decode, encode
christine.barachet's avatar
   
christine.barachet committed
27

28
29
30
try:
    from osgeo import gdal
    from osgeo import ogr
31
except Exception as e:
32
33
34
35
36
37
38
39
40
    print('!!! %spython3-gdal was not found on your system, did you install it?%s' % (COLOR_RED, COLOR_RESET))
    print('On Debian/Ubuntu/Linux Mint you can install it with:\n')
    print('%ssudo apt install python3-gdal%s\n' % (COLOR_GREEN, COLOR_RESET))
    print('or if you don\'t have superuser access:\n')
    print('%spip3 install GDAL%s\n' % (COLOR_GREEN, COLOR_RESET))
    sys.exit(1)

try:
    import dbf
41
except Exception as e:
42
43
44
45
46
47
    print('!!! %sdbf python3 module was not found on your system, no .dbf file will be generated.%s' % (COLOR_YELLOW, COLOR_RESET))
    print('On Debian/Ubuntu/Linux Mint you can install it with:\n')
    print('%ssudo apt install python3-dbf%s\n' % (COLOR_GREEN, COLOR_RESET))
    print('or if you don\'t have superuser access:\n')
    print('%spip3 install dbf%s\n' % (COLOR_GREEN, COLOR_RESET))

48
MY_ABS_PATH=os.path.abspath(__file__)
49
50
MY_DIR=os.path.dirname(MY_ABS_PATH)
DATA_PATH=os.path.join(MY_DIR, '..', 'data')
christine.barachet's avatar
   
christine.barachet committed
51

52
# cutting the rasters according to user selection
53
def cutting_raster(raster,raster_out, xmin_slect, ymax_slect, xmax_slect, ymin_slect):
Julien Veyssier's avatar
Julien Veyssier committed
54
55
56
    ds = gdal.Open(raster)
    ds = gdal.Translate(raster_out, ds, projWin=[xmin_slect, ymax_slect, xmax_slect, ymin_slect])
    ds = None
christine.barachet's avatar
   
christine.barachet committed
57

58
59
60
61
62
63
64
65
66
67
68
69
def get_raster_bounds(layer, xmin_slect, ymax_slect, xmax_slect, ymin_slect):
    # this is kind of sketchy and should be done right
    ds = gdal.Open(layer)
    infoLines = gdal.Info(ds).split('\n')
    for l in infoLines:
        if 'Upper Left' in l:
            xmin_slect = max(xmin_slect, float(l.split(')')[0].split('(')[1].split(',')[0]))
            ymax_slect = min(ymax_slect, float(l.split(')')[0].split('(')[1].split(',')[1]))
        elif 'Lower Right' in l:
            xmax_slect = min(xmax_slect, float(l.split(')')[0].split('(')[1].split(',')[0]))
            ymin_slect = max(ymin_slect, float(l.split(')')[0].split('(')[1].split(',')[1]))
    ds = None
70

71
    return xmin_slect, ymax_slect, xmax_slect, ymin_slect
72

73
def get_polygon_bounds(shape, xmin_slect, ymax_slect, xmax_slect, ymin_slect):
74
75
76
77
78
79
80
    ds = ogr.Open(shape)
    layer = ds.GetLayer()
    xmin, xmax, ymin, ymax = layer.GetExtent()
    xmin_slect = max(xmin_slect, xmin)
    xmax_slect = min(xmax_slect, xmax)
    ymin_slect = max(ymin_slect, ymin)
    ymax_slect = min(ymax_slect, ymax)
81

82
83
84
    return xmin_slect, ymax_slect, xmax_slect, ymin_slect

def get_coords(parms, xmin_slect, ymax_slect, xmax_slect, ymin_slect):
Julien Veyssier's avatar
Julien Veyssier committed
85
86
87
88
    xmin_slect = max(xmin_slect, float(parms.get('surface', 'west')))
    ymax_slect = min(ymax_slect, float(parms.get('surface', 'north')))
    xmax_slect = min(xmax_slect, float(parms.get('surface', 'east')))
    ymin_slect = max(ymin_slect, float(parms.get('surface', 'south')))
89
    return xmin_slect, ymax_slect, xmax_slect, ymin_slect
christine.barachet's avatar
   
christine.barachet committed
90

christine.barachet's avatar
   
christine.barachet committed
91
try:
92
93
    # Python 3
    from subprocess import DEVNULL
christine.barachet's avatar
   
christine.barachet committed
94
except ImportError:
christine.barachet's avatar
   
christine.barachet committed
95
    #DEVNULL = open('cbtrace_step1', 'wb')
96
    DEVNULL = open(os.devnull, 'wb')
christine.barachet's avatar
   
christine.barachet committed
97

98
99
100
101
102
103
104
105
'''
 MAIN
'''
def main(parms_file):
    print('---------- HRU-delin Step 1 started ---------------------------------------------')

    configFileDir = os.path.dirname(parms_file)
    # create main env
Julien Veyssier's avatar
Julien Veyssier committed
106
    buildGrassEnv(os.path.join(configFileDir, 'grass_db'), 'hru-delin')
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
    os.environ['GISRC'] = os.path.join(configFileDir, 'grass_db', 'grassdata', 'hru-delin', '.grassrc')
    # Get parameters from configuration file
    parms = ConfigParser.ConfigParser(allow_no_value=True)
    parms.read(parms_file)
    directory = parms.get('dir_in', 'dir')
    # manage absolute and relative paths
    if not os.path.isabs(directory):
        directory = os.path.join(configFileDir, directory)
    directory_out = parms.get('dir_out', 'files')
    if not os.path.isabs(directory_out):
        directory_out = os.path.join(configFileDir, directory_out)

    # Initializing the variables for cutting up of layers
    xmin_slect = ymin_slect = 0.0
    ymax_slect = xmax_slect = 9999999999.9

    dem = os.path.join(directory, str(parms.get('files_in', 'dem')))
    dem_name = os.path.basename(parms.get('files_in', 'dem')).split('.')[0]

    selection = (parms.get('surface', 'selection'))
127
    xmin_slect, ymax_slect, xmax_slect, ymin_slect = get_raster_bounds(dem, xmin_slect, ymax_slect, xmax_slect, ymin_slect)
128
129
130
131

    if selection == 'total':
        pass
    elif selection == 'polygon':
132
        xmin_slect, ymax_slect, xmax_slect, ymin_slect = get_polygon_bounds(os.path.join(directory, parms.get('surface', 'polygon')), xmin_slect, ymax_slect, xmax_slect, ymin_slect)
133
    elif selection == 'coords':
134
        xmin_slect, ymax_slect, xmax_slect, ymin_slect = get_coords(parms, xmin_slect, ymax_slect, xmax_slect, ymin_slect)
135
136
137
138
139
140
141
142
    else:
        print('==============================================================================')
        print('No surface specified in hrudelin_config or invalid parameter, full DEM applied')
        print('==============================================================================')

    print('---------- HRU-delin Step 1 : Preparing the layers')
    # Cutting up of the DEM
    dem_cut = os.path.join(directory_out, 'step1_dem_cut.tif')
143
    cutting_raster(dem, dem_cut, xmin_slect, ymax_slect, xmax_slect, ymin_slect)
144
145
146
147
148

    irr_rast = str(parms.get('irrigation', 'irrig_rast'))
    if irr_rast != '':
        irr_rast = os.path.join(directory, irr_rast)
        layer_cut = os.path.join(directory_out, 'step1_irrigation.tif')
149
        cutting_raster(irr_rast, layer_cut, xmin_slect, ymax_slect, xmax_slect, ymin_slect)
150
151
152
153

    # Cutting up of the gauges shapefile
    gauges_shp = os.path.join(directory, parms.get('files_in', 'gauges'))
    gauges_selected = os.path.join(directory_out, 'gauges_selected.shp')
154
155
    ## this ensures intersection with DEM
    ogr2ogrMain(['', '-spat', str(xmin_slect), str(ymin_slect), str(xmax_slect), str(ymax_slect), gauges_selected, gauges_shp])
156
157
158

    # Cutting up and save other rasters (landuse, soils, geology, ...)
    for data in parms.items('data'):
159
        xmin_slect, ymax_slect, xmax_slect, ymin_slect = get_raster_bounds(os.path.join(directory, data[1]), xmin_slect, ymax_slect, xmax_slect, ymin_slect)
160
161
162
    for data in parms.items('data'):
        layer = os.path.join(directory, data[1])
        layer_cut = os.path.join(directory_out, 'step1_' + data[1])
163
        cutting_raster(layer, layer_cut, xmin_slect, ymax_slect, xmax_slect, ymin_slect)
164
165

    dem_wk = 'dem_wk'
166
167
168
    grass_run_command('g.proj', flags='c', georef=dem_cut, stdout=DEVNULL, stderr=DEVNULL)
    grass_run_command('r.in.gdal', flags='o', input=dem_cut, output=dem_wk, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
    grass_run_command('g.region', flags='sp', raster=dem_wk, stdout=DEVNULL, stderr=DEVNULL)
169
170
171
172
173
174
175
176
177
178
179
180

    # Fill the DEM (depressionless)
    dem_filled = elev = dem_wk
    if (parms.get('demfill', 'demfill') == 'yes'):
        print('---------- HRU-delin Step 1 : Filling the DEM')
        dem_filled = 'dem_filled'
        dem_temp = 'dem_temp'
        dir_temp = 'dir_temp'
        unfilled_areas = 'unfilled_areas'

        areas = 1
        while (areas > 0):
181
182
183
            grass_run_command('r.fill.dir', input=elev, output=dem_filled, direction=dir_temp, areas=unfilled_areas,
                overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
            areas = int(grass_parse_command('r.info', flags='r', map=unfilled_areas, stdout=DEVNULL, stderr=DEVNULL)['max'])
184
185
186
187
188
189
            elev = dem_filled

    # Slope and aspect derivation
    print('---------- HRU-delin Step 1 : Derivating layers from MNT')
    dem_slope = 'dem_slope'
    dem_aspect = 'dem_aspect'
190
191
    grass_run_command('r.slope.aspect', elevation=dem_filled, slope=dem_slope, aspect=dem_aspect,
        overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
192
193
194

    # Reclassifying rasters
    # DEM
195
    dem_rcl = 'dem_rcl'
196
197

    if (parms.get('reclass_dem', 'rules_auto_dem')) == 'yes':
198
        pRecode = grass_feed_command('r.recode', input=dem_filled, output=dem_rcl, rules='-', quiet=True)
199

200
        step_reclass = int(parms.get('reclass_dem', 'step_dem'))
201
202
        alti_min = float(grass_parse_command('r.info', flags='r', map=dem_filled, stdout=DEVNULL, stderr=DEVNULL)['min'])
        alti_max = float(grass_parse_command('r.info', flags='r', map=dem_filled, stdout=DEVNULL, stderr=DEVNULL)['max'])
203

204
205
        # just a copy, still do the pipe
        recl_dem_file = os.path.join(directory_out, 'reclass_rules_dem')
206
207
208
209
        with open(recl_dem_file, 'w') as dem_rules:
            new_max = math.floor(alti_max / step_reclass) * step_reclass + step_reclass
            if alti_min < 0:
                x = 1
210
211
212
                rule = '%.1f:0:1:%d\n' % (alti_min, x)
                pRecode.stdin.write(encode(rule))
                dem_rules.write(rule)
213
214
215
216
217
218
219
                new_min = 0
            else:
                x = 0
                new_min = math.floor(alti_min / step_reclass) * step_reclass
            while new_min < new_max:
                x += 1
                # keep the decimal to make sure r.recode does not round things
220
221
222
                rule = '%.1f:%.1f:%d:%d\n' % (new_min, new_min+step_reclass, x, x)
                pRecode.stdin.write(encode(rule))
                dem_rules.write(rule)
223
                new_min += step_reclass
224
225
        pRecode.stdin.close()
        pRecode.wait()
226
    else:
227
        recl_dem_file = os.path.join(configFileDir, 'reclass_rules_dem')
228
229
        if not os.path.isfile(recl_dem_file):
            sys.exit('------------> File reclass_rules_dem is missing')
230
        # r.recode does a better job with grass7, values close to interval limits are now correctly reclassified
231
        grass_run_command('r.recode', input=dem_filled, output=dem_rcl, rules=recl_dem_file, stdout=DEVNULL, stderr=DEVNULL)
232
233
234
        # Save the rules to directory_out
        recl_dem_name = os.path.basename(recl_dem_file)
        shutil.copyfile(recl_dem_file, os.path.join(directory_out, recl_dem_name))
235

236
    # Slope
237
238
    slp_rcl = 'slp_rcl'

239
    if (parms.get('reclass_slope', 'rules_auto_slope')) == 'yes':
240
        pRecode = grass_feed_command('r.recode', input=dem_slope, output=slp_rcl, rules='-', quiet=True)
241
242
243
        recl_slope_file = os.path.join(directory_out, 'reclass_rules_slope')
        default_recl_slope_file = os.path.join(DATA_PATH, 'reclass_default_rules_slope')

244
245
        slp_min = float(grass_parse_command('r.info', flags='r', map=dem_slope, stdout=DEVNULL, stderr=DEVNULL)['min'])
        slp_max = float(grass_parse_command('r.info', flags='r', map=dem_slope, stdout=DEVNULL, stderr=DEVNULL)['max'])
246
247
        with open(default_recl_slope_file, 'r') as default_slp_rules:
            slp_class = default_slp_rules.read()
248
249
250
251
        max_slp = str(slp_max)
        slp_class_new = slp_class.replace('max_slp', max_slp)
        max_class = str(int(math.ceil(slp_max)))
        slp_class = slp_class_new.replace('max_class', max_class)
252
        # write a copy but still do the pipe
253
254
        with open(recl_slope_file, 'w') as slp_rules:
            slp_rules.write(slp_class)
255
256
257
        pRecode.stdin.write(encode(slp_class))
        pRecode.stdin.close()
        pRecode.wait()
258
    else:
259
        recl_slope_file = os.path.join(configFileDir, 'reclass_rules_slope')
260
261
        if not os.path.isfile(recl_slope_file):
            sys.exit('------------> file reclass_rules_slope is missing')
262
        # improve some subprocess.Popen => run_command with rules=FILE
263
        grass_run_command('r.recode', input=dem_slope, output=slp_rcl, rules=recl_slope_file, stdout=DEVNULL, stderr=DEVNULL)
264
265
        recl_slope_name = os.path.basename(recl_slope_file)
        shutil.copyfile(recl_slope_file, os.path.join(directory_out, recl_slope_name))
266
267

    # Aspect
268
    recl_aspect_file = os.path.join(directory_out, 'reclass_rules_aspect')
269
270
    if (parms.get('reclass_aspect', 'rules_auto_aspect')) == 'yes':
        recl_aspect = os.path.join(DATA_PATH, 'reclass_default_rules_aspect')
271
        shutil.copyfile(recl_aspect, recl_aspect_file)
272
    else:
273
274
275
276
277
        user_aspect = os.path.join(configFileDir, 'reclass_rules_aspect')
        if not os.path.isfile(user_aspect):
            sys.exit('------------> File %s is missing' % user_aspect)
        shutil.copyfile(user_aspect, recl_aspect_file)

278
    asp_rcl = 'asp_rcl'
279
    grass_run_command('r.recode', input=dem_aspect, output=asp_rcl, rules=recl_aspect_file, stdout=DEVNULL, stderr=DEVNULL)
280
281
282
283
284
285
286
287
288

    # Oriented flows generation
    print('---------- HRU-delin Step 1 : Generating Oriented Flows')
    min_size = int(parms.get('basin_min_size', 'size'))
    streams = 'streams_wk'
    subbasins = 'subbasins_wk'
    halfbasins='halfbasins'
    accum = 'accum_wk'
    drain = 'drain_wk'
289
    grass_run_command('r.watershed',
290
291
292
293
294
295
296
297
298
        flags='sb',
        elevation=dem_filled,
        stream=streams,
        basin=subbasins,
        half_basin=halfbasins,
        accumulation=accum,
        drainage=drain,
        threshold=min_size,
        overwrite='True',
299
        stdout=DEVNULL, stderr=DEVNULL
300
301
302
303
304
    )

    # Save rasters
    print('---------- HRU-delin Step 1 : Saving the rasters')
    raster_out = os.path.join(directory_out, 'step1_dem_filled.tif')
305
    grass_run_command('r.out.gdal', input=dem_filled, output=raster_out, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
306
307

    raster_out = os.path.join(directory_out, 'step1_dem_reclass.tif')
308
    grass_run_command('r.out.gdal', input=dem_rcl, output=raster_out, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
309
310

    raster_out = os.path.join(directory_out, 'step1_slope.tif')
311
    grass_run_command('r.out.gdal', input=dem_slope, output=raster_out, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
312
313

    raster_out = os.path.join(directory_out, 'step1_aspect.tif')
314
    grass_run_command('r.out.gdal', input=dem_aspect, output=raster_out, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
315
316

    raster_out = os.path.join(directory_out, 'step1_slope_reclass.tif')
317
    grass_run_command('r.out.gdal', input=slp_rcl, output=raster_out, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
318
319

    raster_out = os.path.join(directory_out, 'step1_aspect_reclass.tif')
320
    grass_run_command('r.out.gdal', input=asp_rcl, output=raster_out, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
321
322

    raster_out = os.path.join(directory_out, 'step1_streams.tif')
323
    grass_run_command('r.out.gdal', input=streams, output=raster_out, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
324
325

    raster_out = os.path.join(directory_out, 'step1_drain.tif')
326
    grass_run_command('r.out.gdal', input=drain, output=raster_out, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
327
328

    raster_out = os.path.join(directory_out, 'step1_accum.tif')
329
    grass_run_command('r.out.gdal', input=accum, output=raster_out, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
330
331

    raster_out = os.path.join(directory_out, 'step1_subbasins.tif')
332
    grass_run_command('r.out.gdal', input=subbasins, output=raster_out, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
333
334

    raster_out = os.path.join(directory_out, 'step1_halfbasins.tif')
335
    grass_run_command('r.out.gdal', input=halfbasins, output=raster_out, overwrite='True', stdout=DEVNULL, stderr=DEVNULL)
336
337
338
    print('---------- HRU-delin Step 1 ended   ---------------------------------------------')

if __name__ == '__main__':
339
340
    from progressColors import *
    from parallel import buildGrassEnv
341
    from ogr2ogr import main as ogr2ogrMain
342
    from utils import grass_run_command, grass_parse_command, grass_feed_command
343
344
345
346
347
348
349
350
351
352
    try:
        from tqdm import tqdm
    except Exception as e:
        print('!!! %s tqdm python3 module not found !!%s\n' % (COLOR_RED, COLOR_RESET))
        print('On Debian/Ubuntu/Linux Mint you can install it with:\n')
        print('%ssudo apt install python3-tqdm%s\n' % (COLOR_GREEN, COLOR_RESET))
        print('or if you don\'t have superuser access:\n')
        print('%spip3 install tqdm%s\n' % (COLOR_GREEN, COLOR_RESET))
        sys.exit(1)

353
354
355
356
357
358
359
360
361
    parms_file = 'hrudelin_config.cfg'
    if len(sys.argv) > 1:
        parms_file = sys.argv[1]

    main(parms_file)

    try:
        os.system('notify-send "hru-delin step 1 complete"')
    except Exception as e:
Julien Veyssier's avatar
Julien Veyssier committed
362
        pass
363
364
365
else:
    from .progressColors import *
    from .parallel import buildGrassEnv
366
    from .ogr2ogr import main as ogr2ogrMain
367
    from .utils import grass_run_command, grass_parse_command, grass_feed_command