Index: wflow-py/Scripts/wflow_flood.py =================================================================== diff -u -r00f94f519ec3cad07f362bcc5323bf6a324c6c32 -r78869427fc1538fb6f1b316b9d8188e1c66edff2 --- wflow-py/Scripts/wflow_flood.py (.../wflow_flood.py) (revision 00f94f519ec3cad07f362bcc5323bf6a324c6c32) +++ wflow-py/Scripts/wflow_flood.py (.../wflow_flood.py) (revision 78869427fc1538fb6f1b316b9d8188e1c66edff2) @@ -90,6 +90,9 @@ options.ldd_file = inun_lib.configget(config, 'maps', 'ldd_file', True) + options.stream_file = inun_lib.configget(config, 'maps', + 'stream_file', + True) options.riv_length_file = inun_lib.configget(config, 'maps', 'riv_length_file', True) @@ -176,6 +179,7 @@ # Open terrain data for reading ds_dem, rasterband_dem = inun_lib.get_gdal_rasterband(options.dem_file) ds_ldd, rasterband_ldd = inun_lib.get_gdal_rasterband(options.ldd_file) + ds_stream, rasterband_stream = inun_lib.get_gdal_rasterband(options.stream_file) n = 0 for x_loop in range(0, len(x), options.x_tile): x_start = np.maximum(x_loop, 0) @@ -203,28 +207,52 @@ (x_end + x_overlap_max) - (x_start - x_overlap_min), (y_end + y_overlap_max) - (y_start - y_overlap_min) ) + stream = rasterband_stream.ReadAsArray(x_start - x_overlap_min, + y_start - y_overlap_min, + (x_end + x_overlap_max) - (x_start - x_overlap_min), + (y_end + y_overlap_max) - (y_start - y_overlap_min) + ) # write to temporary file terrain_temp_file = os.path.join(options.dest_path, 'terrain_temp.map') drainage_temp_file = os.path.join(options.dest_path, 'drainage_temp.map') - inun_lib.gdal_writemap(terrain_temp_file, 'PCRaster', + stream_temp_file = os.path.join(options.dest_path, 'stream_temp.map') + if rasterband_dem.GetNoDataValue() is not None: + inun_lib.gdal_writemap(terrain_temp_file, 'PCRaster', + np.arange(0, terrain.shape[1]), + np.arange(0, terrain.shape[0]), + terrain, rasterband_dem.GetNoDataValue(), + gdal_type=gdal.GDT_Float32, + logging=logger) + else: + # in case no nodata value is found + logger.warning('No nodata value found in {:s}. assuming -9999'.format(options.dem_file)) + inun_lib.gdal_writemap(terrain_temp_file, 'PCRaster', + np.arange(0, terrain.shape[1]), + np.arange(0, terrain.shape[0]), + terrain, -9999., + gdal_type=gdal.GDT_Float32, + logging=logger) + + inun_lib.gdal_writemap(drainage_temp_file, 'PCRaster', np.arange(0, terrain.shape[1]), np.arange(0, terrain.shape[0]), - terrain, rasterband_dem.GetNoDataValue(), - gdal_type=gdal.GDT_Float32, + drainage, rasterband_ldd.GetNoDataValue(), + gdal_type=gdal.GDT_Int32, logging=logger) - inun_lib.gdal_writemap(drainage_temp_file, 'PCRaster', + inun_lib.gdal_writemap(stream_temp_file, 'PCRaster', np.arange(0, terrain.shape[1]), np.arange(0, terrain.shape[0]), - drainage, rasterband_ldd.GetNoDataValue(), + stream, rasterband_ldd.GetNoDataValue(), gdal_type=gdal.GDT_Int32, logging=logger) # read as pcr objects pcr.setclone(terrain_temp_file) terrain_pcr = pcr.readmap(terrain_temp_file) drainage_pcr = pcr.lddrepair(pcr.ldd(pcr.readmap(drainage_temp_file))) # convert to ldd type map + stream_pcr = pcr.scalar(pcr.readmap(stream_temp_file)) # convert to ldd type map # compute streams - stream_ge, subcatch = inun_lib.subcatch_stream(drainage_pcr, options.hand_strahler) # generate streams + stream_ge, subcatch = inun_lib.subcatch_stream(drainage_pcr, stream_pcr, options.hand_strahler) # generate streams basin = pcr.boolean(subcatch) hand_pcr, dist_pcr = inun_lib.derive_HAND(terrain_pcr, drainage_pcr, 3000, rivers=pcr.boolean(stream_ge), basin=basin) @@ -243,6 +271,7 @@ band_hand.FlushCache() ds_dem = None ds_ldd = None + ds_stream = None band_hand.SetNoDataValue(-9999.) ds_hand = None logger.info('Finalizing {:s}'.format(hand_file)) @@ -271,13 +300,21 @@ inun_file = os.path.join(flood_folder, '{:s}.tif'.format(case_name)) hand_temp_file = os.path.join(flood_folder, 'hand_temp.map') drainage_temp_file = os.path.join(flood_folder, 'drainage_temp.map') + stream_temp_file = os.path.join(flood_folder, 'stream_temp.map') flood_vol_temp_file = os.path.join(flood_folder, 'flood_warp_temp.tif') # load the data with river levels and compute the volumes if options.file_format == 0: + # assume we need the maximum value in a NetCDF time series grid a = nc.Dataset(options.flood_map, 'r') xax = a.variables['x'][:] - yax = np.flipud(a.variables['y'][:]) - flood = np.flipud(a.variables[options.flood_variable][0, :, :]) + yax = a.variables['y'][:] + + flood_series = a.variables[options.flood_variable][:] + flood = flood_series.max(axis=0) + if yax[-1] > yax[0]: + yax = np.flipud(yax) + flood = np.flipud(flood) + a.close() elif options.file_format == 1: xax, yax, flood, flood_fill_value = inun_lib.gdal_readmap(options.flood_map, 'PCRaster') #res_x = x[1]-x[0] @@ -290,8 +327,11 @@ if options.file_format == 0: a = nc.Dataset(options.bankfull_map, 'r') xax = a.variables['x'][:] - yax = np.flipud(a.variables['y'][:]) - bankfull = np.flipud(a.variables[options.flood_variable][0, :, :]) + yax = a.variables['y'][:] + bankfull = a.variables[options.flood_variable][0, :, :] + if yax[-1] > yax[0]: + yax = np.flipud(yax) + bankfull = np.flipud(bankful) a.close() elif options.file_format == 1: xax, yax, bankfull, bankfull_fill_value = inun_lib.gdal_readmap(options.bankfull_map, 'PCRaster') @@ -308,6 +348,7 @@ inun_lib.gdal_writemap(flood_vol_map, 'GTiff', xax, yax, np.maximum(flood_vol_m_data, 0), -999.) ds_hand, rasterband_hand = inun_lib.get_gdal_rasterband(hand_file) ds_ldd, rasterband_ldd = inun_lib.get_gdal_rasterband(options.ldd_file) + ds_stream, rasterband_stream = inun_lib.get_gdal_rasterband(options.stream_file) logger.info('Preparing flood map in {:s} ...please wait...'.format(inun_file)) ds_inun = inun_lib.prepare_gdal(inun_file_tmp, x, y, logging=logger, srs=srs) @@ -344,6 +385,11 @@ (x_end + x_overlap_max) - (x_start - x_overlap_min), (y_end + y_overlap_max) - (y_start - y_overlap_min) ) + stream = rasterband_stream.ReadAsArray(x_start - x_overlap_min, + y_start - y_overlap_min, + (x_end + x_overlap_max) - (x_start - x_overlap_min), + (y_end + y_overlap_max) - (y_start - y_overlap_min) + ) print('len x-ax: {:d} len y-ax {:d} x-shape {:d} y-shape {:d}'.format(len(x_tile_ax), len(y_tile_ax), hand.shape[1], hand.shape[0])) inun_lib.gdal_writemap(hand_temp_file, 'PCRaster', x_tile_ax, @@ -357,13 +403,20 @@ drainage, rasterband_ldd.GetNoDataValue(), gdal_type=gdal.GDT_Int32, logging=logger) + inun_lib.gdal_writemap(stream_temp_file, 'PCRaster', + x_tile_ax, + y_tile_ax, + stream, rasterband_stream.GetNoDataValue(), + gdal_type=gdal.GDT_Int32, + logging=logger) # read as pcr objects pcr.setclone(hand_temp_file) hand_pcr = pcr.readmap(hand_temp_file) drainage_pcr = pcr.lddrepair(pcr.ldd(pcr.readmap(drainage_temp_file))) # convert to ldd type map + stream_pcr = pcr.scalar(pcr.readmap(drainage_temp_file)) # convert to ldd type map # prepare a subcatchment map - stream_ge, subcatch = inun_lib.subcatch_stream(drainage_pcr, options.catchment_strahler) # generate subcatchments + stream_ge, subcatch = inun_lib.subcatch_stream(drainage_pcr, stream_pcr, options.catchment_strahler) # generate subcatchments drainage_surf = pcr.ifthen(stream_ge > 0, pcr.accuflux(drainage_pcr, 1)) # proxy of drainage surface inaccurate at tile edges # compute weights for spreadzone (1/drainage_surf) subcatch = pcr.spreadzone(subcatch, 0, 0) Index: wflow-py/Scripts/wflow_flood_lib.py =================================================================== diff -u -r414daa740d69a0674195d0bc996218a61f76207e -r78869427fc1538fb6f1b316b9d8188e1c66edff2 --- wflow-py/Scripts/wflow_flood_lib.py (.../wflow_flood_lib.py) (revision 414daa740d69a0674195d0bc996218a61f76207e) +++ wflow-py/Scripts/wflow_flood_lib.py (.../wflow_flood_lib.py) (revision 78869427fc1538fb6f1b316b9d8188e1c66edff2) @@ -23,6 +23,7 @@ from osgeo import osr, gdal, gdalconst import pcraster as pcr +import pdb def setlogger(logfilename, logReference, verbose=True): """ @@ -103,7 +104,7 @@ except: Def = True ret = default - configset(config, section, var, str(default), overwrite=False) + #configset(config, section, var, str(default), overwrite=False) default = Def return ret @@ -335,11 +336,12 @@ return hand, dist -def subcatch_stream(ldd, threshold): +def subcatch_stream(ldd, stream, threshold): """ Derive catchments based upon strahler threshold Input: ldd -- pcraster object direction, local drain directions + stream -- pcraster object direction, streamorder threshold -- integer, strahler threshold, subcatchments ge threshold are derived output: @@ -349,7 +351,7 @@ """ # derive stream order - stream = pcr.streamorder(ldd) + # stream = pcr.streamorder(ldd) stream_ge = pcr.ifthen(stream >= threshold, stream) stream_up_sum = pcr.ordinal(pcr.upstream(ldd, pcr.cover(pcr.scalar(stream_ge), 0))) # detect any transfer of strahler order, to a higher strahler order.