Index: wflow-py/Scripts/wflow_flood.py =================================================================== diff -u -r77a2ccbb602adf1c3698bfa8c50b79320c8362cd -rf6d346c8a0883f44e74b99d2dce672a3c881b30e --- wflow-py/Scripts/wflow_flood.py (.../wflow_flood.py) (revision 77a2ccbb602adf1c3698bfa8c50b79320c8362cd) +++ wflow-py/Scripts/wflow_flood.py (.../wflow_flood.py) (revision f6d346c8a0883f44e74b99d2dce672a3c881b30e) @@ -438,9 +438,7 @@ 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, stream_pcr, options.hand_strahler) # generate streams - pcr.report(stream_ge, 'stream_ge.map') - pcr.report(pcr.scalar(subcatch), 'subcatch.map') + stream_ge, subcatch = inun_lib.subcatch_stream(drainage_pcr, options.hand_strahler, stream=stream_pcr) # 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) @@ -473,9 +471,9 @@ # HAND file has now been prepared, moving to flood mapping part # ##################################################################################### # load the staticmaps needed to estimate volumes across all - xax, yax, riv_length, fill_value = inun_lib.gdal_readmap(options.riv_length_file, 'GTiff') + xax, yax, riv_length, fill_value = inun_lib.gdal_readmap(options.riv_length_file, 'GTiff', logging=logger) riv_length = np.ma.masked_where(riv_length==fill_value, riv_length) - xax, yax, riv_width, fill_value = inun_lib.gdal_readmap(options.riv_width_file, 'GTiff') + xax, yax, riv_width, fill_value = inun_lib.gdal_readmap(options.riv_width_file, 'GTiff', logging=logger) riv_width[riv_width == fill_value] = 0 x_res = np.abs((xax[-1]-xax[0])/(len(xax)-1)) @@ -510,7 +508,7 @@ a.close() elif options.file_format == 1: logger.info('Reading flood from {:s} PCRaster file'.format(options.flood_map)) - xax, yax, flood, flood_fill_value = inun_lib.gdal_readmap(options.flood_map, 'PCRaster') + xax, yax, flood, flood_fill_value = inun_lib.gdal_readmap(options.flood_map, 'PCRaster', logging=logger) flood[flood==flood_fill_value] = 0. #res_x = x[1]-x[0] #res_y = y[1]-y[0] @@ -535,7 +533,7 @@ a.close() elif options.file_format == 1: logger.info('Reading bankfull from {:s} PCRaster file'.format(options.bankfull_map)) - xax, yax, bankfull, bankfull_fill_value = inun_lib.gdal_readmap(options.bankfull_map, 'PCRaster') + xax, yax, bankfull, bankfull_fill_value = inun_lib.gdal_readmap(options.bankfull_map, 'PCRaster', logging=logger) # flood = bankfull*2 # res_x = 2000 # res_y = 2000 @@ -544,9 +542,9 @@ flood_vol_m = riv_length*riv_width*flood_vol/(x_res * y_res) # volume expressed in meters water disc (1e6 is the surface area of one wflow grid cell) flood_vol_m_data = flood_vol_m.data flood_vol_m_data[flood_vol_m.mask] = -999. - print('Saving water layer map to {:s}'.format(flood_vol_map)) + logger.info('Saving water layer map to {:s}'.format(flood_vol_map)) # write to a tiff file - inun_lib.gdal_writemap(flood_vol_map, 'GTiff', xax, yax, np.maximum(flood_vol_m_data, 0), -999.) + inun_lib.gdal_writemap(flood_vol_map, 'GTiff', xax, yax, np.maximum(flood_vol_m_data, 0), -999., logging=logger) 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) @@ -614,35 +612,34 @@ 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 + stream_pcr = pcr.scalar(pcr.readmap(stream_temp_file)) # convert to ldd type map # prepare a subcatchment map # TODO check weighting scheme, perhaps not necessary # weight = 1./pcr.scalar(pcr.spreadzone(pcr.cover(pcr.ordinal(drainage_surf), 0), 0, 0)) # subcatch_fill = pcr.scalar(pcr.spreadzone(subcatch, 0, weight)) # # cover subcatch with subcatch_fill - # pcr.report(weight, 'weight_{:02d}.map'.format(n)) - # pcr.report(subcatch, 'subcatch_{:02d}.map'.format(n)) - # pcr.report(pcr.nominal(subcatch_fill), 'subcatch_fill_{:02d}.map'.format(n)) inun_lib.gdal_warp(flood_vol_map, hand_temp_file, flood_vol_temp_file, gdal_interp=gdalconst.GRA_NearestNeighbour) # , - x_tile_ax, y_tile_ax, flood_meter, fill_value = inun_lib.gdal_readmap(flood_vol_temp_file, 'GTiff') + x_tile_ax, y_tile_ax, flood_meter, fill_value = inun_lib.gdal_readmap(flood_vol_temp_file, 'GTiff', logging=logger) # convert meter depth to volume [m3] flood_vol = pcr.numpy2pcr(pcr.Scalar, flood_meter, fill_value)*((x_tile_ax[1] - x_tile_ax[0]) * (y_tile_ax[0] - y_tile_ax[1])) # resolution of SRTM *1166400000. + # first prepare a basin map, belonging to the HAND order we are looking at + # TODO: replace options.hand_strahler by current HAND strahler order + stream_ge_hand, subcatch_hand = inun_lib.subcatch_stream(drainage_pcr, options.hand_strahler, stream=stream_pcr) # TODO: loop over river order options.catchment_strahler to maximum order found in tile - stream_ge, subcatch = inun_lib.subcatch_stream(drainage_pcr, stream_pcr, options.catchment_strahler) # generate subcatchments + stream_ge, subcatch = inun_lib.subcatch_stream(drainage_pcr, options.catchment_strahler, stream=stream_pcr, basin=pcr.boolean(pcr.cover(subcatch_hand, 0)), assign_existing=True) # generate subcatchments, only within basin for HAND # TODO: mask subcatchments with river order lower than thres or higher than thres 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) + # subcatch = pcr.spreadzone(subcatch, 0, 0) # TODO: put areas outside subcatch to zero ## now we have some nice volume. Now we need to redistribute! # TODO: insert selected HAND map inundation_pcr = inun_lib.volume_spread(drainage_pcr, hand_pcr, subcatch, flood_vol, volume_thres=0., iterations=options.iterations, - area_multiplier=options.area_multiplier) # 1166400000. + area_multiplier=options.area_multiplier, logging=logger) # 1166400000. inundation = pcr.pcr2numpy(inundation_pcr, -9999.) # cut relevant part if y_overlap_max == 0: Index: wflow-py/Scripts/wflow_flood_lib.py =================================================================== diff -u -r77a2ccbb602adf1c3698bfa8c50b79320c8362cd -rf6d346c8a0883f44e74b99d2dce672a3c881b30e --- wflow-py/Scripts/wflow_flood_lib.py (.../wflow_flood_lib.py) (revision 77a2ccbb602adf1c3698bfa8c50b79320c8362cd) +++ wflow-py/Scripts/wflow_flood_lib.py (.../wflow_flood_lib.py) (revision f6d346c8a0883f44e74b99d2dce672a3c881b30e) @@ -340,19 +340,23 @@ up_elevation = pcr.ifthenelse(basin, pcr.scalar(pcr.subcatchment(ldd, height_river)), pcr.scalar(pcr.spreadzone(height_river, 0, friction))) # replace areas outside of basin by a spread zone calculation. hand = pcr.max(pcr.scalar(pcr.ordinal(dem*100))-up_elevation, 0)/100 # convert back to float in DEM units + # hand = (pcr.scalar(pcr.ordinal(dem*100))-up_elevation)/100 # convert back to float in DEM units dist = pcr.ldddist(ldd, stream, 1) # compute horizontal distance estimate return hand, dist -def subcatch_stream(ldd, stream, threshold, min_strahler=-999, max_strahler=999, assign_edge=False, assign_existing=False, up_area=None): +def subcatch_stream(ldd, threshold, stream=None, min_strahler=-999, max_strahler=999, assign_edge=False, assign_existing=False, up_area=None, basin=None): """ Derive catchments based upon strahler threshold Input: ldd -- pcraster object direction, local drain directions threshold -- integer, strahler threshold, subcatchments ge threshold are derived - min_strahler -- integer, minimum strahler threshold of river catchments + stream=None -- pcraster object ordinal, stream order map (made with pcr.streamorder), if provided, stream order + map is not generated on the fly but used from this map. Useful when a subdomain within a catchment is + provided, which would cause edge effects in the stream order map + min_strahler=-999 -- integer, minimum strahler threshold of river catchments to return - max_strahler -- integer, maximum strahler threshold of river catchments + max_strahler=999 -- integer, maximum strahler threshold of river catchments to return assign_unique=False -- if set to True, unassigned connected areas at the edges of the domain are assigned a unique id as well. If set @@ -367,7 +371,9 @@ """ # derive stream order - # stream = pcr.streamorder(ldd) + if stream is None: + 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. @@ -379,6 +385,9 @@ # derive upstream catchment areas (write to file) subcatch = pcr.nominal(pcr.subcatchment(ldd, transition_unique)) + # mask out areas outside basin + if basin is not None: + subcatch = pcr.ifthen(basin, subcatch) if assign_edge: # fill unclassified areas (in pcraster equal to zero) with a unique id, above the maximum id assigned so far @@ -402,7 +411,7 @@ return stream_ge, pcr.ordinal(subcatch) -def volume_spread(ldd, hand, subcatch, volume, volume_thres=0., area_multiplier=1., iterations=15): +def volume_spread(ldd, hand, subcatch, volume, volume_thres=0., area_multiplier=1., iterations=15, logging=logging): """ Estimate 2D flooding from a 1D simulation per subcatchment reach Input: @@ -419,33 +428,25 @@ #initial values pcr.setglobaloption("unittrue") dem_min = pcr.areaminimum(hand, subcatch) # minimum elevation in subcatchments - # pcr.report(dem_min, 'dem_min.map') dem_norm = hand - dem_min - # pcr.report(dem_norm, 'dem_norm.map') # surface of each subcatchment surface = pcr.areaarea(subcatch)*area_multiplier - pcr.report(surface, 'surface.map') error_abs = pcr.scalar(1e10) # initial error (very high) volume_catch = pcr.areatotal(volume, subcatch) - # pcr.report(volume_catch, 'volume_catch.map') depth_catch = volume_catch/surface # meters water disc averaged over subcatchment - pcr.report(depth_catch, 'depth_catch.map') dem_max = pcr.ifthenelse(volume_catch > volume_thres, pcr.scalar(32.), pcr.scalar(0)) # bizarre high inundation depth dem_min = pcr.scalar(0.) for n in range(iterations): - print('Iteration: {:02d}'.format(n + 1)) + logging.debug('Iteration: {:02d}'.format(n + 1)) #####while np.logical_and(error_abs > error_thres, dem_min < dem_max): dem_av = (dem_min + dem_max)/2 - # pcr.report(dem_av, 'dem_av00.{:03d}'.format(n + 1)) # compute value at dem_av average_depth_catch = pcr.areaaverage(pcr.max(dem_av - dem_norm, 0), subcatch) - # pcr.report(average_depth_catch, 'depth_c0.{:03d}'.format(n + 1)) error = pcr.cover((depth_catch-average_depth_catch)/depth_catch, depth_catch*0) - # pcr.report(error, 'error000.{:03d}'.format(n + 1)) dem_min = pcr.ifthenelse(error > 0, dem_av, dem_min) dem_max = pcr.ifthenelse(error <= 0, dem_av, dem_max) # error_abs = np.abs(error) # TODO: not needed probably, remove @@ -523,7 +524,7 @@ TempDataset = None os.remove(temp_file_name) -def gdal_readmap(file_name, file_format, give_geotrans=False): +def gdal_readmap(file_name, file_format, give_geotrans=False, logging=logging): """ Read geographical file into memory Dependencies are osgeo.gdal and numpy Input: