Index: wflow-py/Scripts/wflow_flood.py =================================================================== diff -u -r8b0b000b126b3d3dd5b3145ba9e19175d10c7866 -r1d05f7d4cda92d45dc5f6400d5ab88020766d883 --- wflow-py/Scripts/wflow_flood.py (.../wflow_flood.py) (revision 8b0b000b126b3d3dd5b3145ba9e19175d10c7866) +++ wflow-py/Scripts/wflow_flood.py (.../wflow_flood.py) (revision 1d05f7d4cda92d45dc5f6400d5ab88020766d883) @@ -228,6 +228,7 @@ import wflow_flood_lib as inun_lib import wflow.pcrut as pcrut import wflow.wflow_lib as wflow_lib +import datetime as dt import pdb @@ -255,6 +256,9 @@ parser.add_option('-c', '--catchment', dest='catchment_strahler', default=7, type='int', help='Strahler order threshold >= are selected as catchment boundaries') + parser.add_option('-t', '--time', + dest='time', default='', + help='time in YYYYMMDDHHMMSS, overrides time in NetCDF input if set') # parser.add_option('-s', '--hand_strahler', # dest='hand_strahler', default=7, type='int', # help='Strahler order threshold >= selected as riverine') @@ -313,6 +317,8 @@ True) options.file_format = inun_lib.configget(config, 'file_settings', 'file_format', 0, datatype='int') + options.out_format = inun_lib.configget(config, 'file_settings', + 'out_format', 0, datatype='int') options.latlon = inun_lib.configget(config, 'file_settings', 'latlon', 0, datatype='int') options.x_tile = inun_lib.configget(config, 'tiling', @@ -327,6 +333,9 @@ 'iterations', 20, datatype='int') options.initial_level = inun_lib.configget(config, 'inundation', 'initial_level', 32., datatype='float') + options.flood_volume_type = inun_lib.configget(config, 'inundation', + 'flood_volume_type', 0, datatype='int') + # options.area_multiplier = inun_lib.configget(config, 'inundation', # 'area_multiplier', 1., datatype='float') logger.info('DEM file: {:s}'.format(options.dem_file)) @@ -396,8 +405,8 @@ # hand file does not exist yet! Generate it, otherwise skip! logger.info('HAND file {:s} not found, start setting up...please wait...'.format(hand_file)) hand_file_tmp = os.path.join(options.dest_path, '{:s}_hand_strahler_{:02d}.tif.tmp'.format(dem_name, hand_strahler)) - ds_hand = inun_lib.prepare_gdal(hand_file_tmp, x, y, logging=logger, srs=srs) - band_hand = ds_hand.GetRasterBand(1) + ds_hand, band_hand = inun_lib.prepare_gdal(hand_file_tmp, x, y, logging=logger, srs=srs) + # band_hand = ds_hand.GetRasterBand(1) # Open terrain data for reading ds_dem, rasterband_dem = inun_lib.get_gdal_rasterband(options.dem_file) @@ -476,6 +485,8 @@ # compute streams stream_ge, subcatch = inun_lib.subcatch_stream(drainage_pcr, hand_strahler, stream=stream_pcr) # generate streams + # compute basins + stream_ge_dummy, subcatch = inun_lib.subcatch_stream(drainage_pcr, options.catchment_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) @@ -507,13 +518,8 @@ ##################################################################################### # 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', 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', logging=logger) - riv_width[riv_width == fill_value] = 0 + # set the clone pcr.setclone(options.ldd_wflow) - # read wflow ldd as pcraster object ldd_pcr = pcr.readmap(options.ldd_wflow) xax, yax, riv_width, fill_value = inun_lib.gdal_readmap(options.riv_width_file, 'GTiff', logging=logger) @@ -522,19 +528,34 @@ x_res, y_res, reallength_wflow = pcrut.detRealCellLength(pcr.scalar(ldd_pcr), not(bool(options.latlon))) cell_surface_wflow = pcr.pcr2numpy(x_res * y_res, 0) - xax, yax, riv_length_fact, fill_value = inun_lib.gdal_readmap(options.riv_length_fact_file, 'GTiff', logging=logger) - riv_length_fact = np.ma.masked_where(riv_length_fact==fill_value, riv_length_fact) - drain_length = wflow_lib.detdrainlength(ldd_pcr, x_res, y_res) + if options.flood_volume_type == 0: + # load the staticmaps needed to estimate volumes across all + # 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', logging=logger) + riv_width[riv_width == fill_value] = 0 - riv_length = pcr.pcr2numpy(drain_length, 0) * riv_length_fact - # riv_length_pcr = pcr.numpy2pcr(pcr.Scalar, riv_length, 0) + # read river length factor file (multiplier) + xax, yax, riv_length_fact, fill_value = inun_lib.gdal_readmap(options.riv_length_fact_file, 'GTiff', logging=logger) + riv_length_fact = np.ma.masked_where(riv_length_fact==fill_value, riv_length_fact) + drain_length = wflow_lib.detdrainlength(ldd_pcr, x_res, y_res) + # compute river length in each cell + riv_length = pcr.pcr2numpy(drain_length, 0) * riv_length_fact + # riv_length_pcr = pcr.numpy2pcr(pcr.Scalar, riv_length, 0) + + flood_folder = os.path.join(options.dest_path, case_name) flood_vol_map = os.path.join(flood_folder, '{:s}_vol.tif'.format(os.path.split(options.flood_map)[1].split('.')[0])) if not(os.path.isdir(flood_folder)): os.makedirs(flood_folder) - inun_file_tmp = os.path.join(flood_folder, '{:s}.tif.tmp'.format(case_name)) - inun_file = os.path.join(flood_folder, '{:s}.tif'.format(case_name)) + if options.out_format == 0: + inun_file_tmp = os.path.join(flood_folder, '{:s}.tif.tmp'.format(case_name)) + inun_file = os.path.join(flood_folder, '{:s}.tif'.format(case_name)) + else: + inun_file_tmp = os.path.join(flood_folder, '{:s}.nc.tmp'.format(case_name)) + inun_file = os.path.join(flood_folder, '{:s}.nc'.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') @@ -546,6 +567,11 @@ a = nc.Dataset(options.flood_map, 'r') xax = a.variables['x'][:] yax = a.variables['y'][:] + if options.time == '': + time_list = nc.num2date(a.variables['time'][:], units = a.variables['time'].units, calendar=a.variables['time'].calendar) + time = [time_list[len(time_list)/2]] + else: + time = [dt.datetime.strptime(options.time, '%Y%m%d%H%M%S')] flood_series = a.variables[options.flood_variable][:] flood_data = flood_series.max(axis=0) @@ -559,6 +585,11 @@ 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', logging=logger) + flood = np.ma.masked_equal(flood, flood_fill_value) + if options.time == '': + options.time = '20000101000000' + time = [dt.datetime.strptime(options.time, '%Y%m%d%H%M%S')] + flood[flood==flood_fill_value] = 0. # load the bankfull depths if options.bankfull_map == '': @@ -569,6 +600,7 @@ a = nc.Dataset(options.bankfull_map, 'r') xax = a.variables['x'][:] yax = a.variables['y'][:] + bankfull_series = a.variables[options.flood_variable][:] bankfull_data = bankfull_series.max(axis=0) if np.ma.is_masked(bankfull_data): @@ -581,12 +613,16 @@ 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', logging=logger) + bankfull = np.ma.masked_equal(bankfull, bankfull_fill_value) # flood = bankfull*2 # res_x = 2000 # res_y = 2000 # subtract the bankfull water level to get flood levels (above bankfull) flood_vol = np.maximum(flood-bankfull, 0) - flood_vol_m = riv_length*riv_width*flood_vol/cell_surface_wflow # volume expressed in meters water disc + if options.flood_volume_type == 0: + flood_vol_m = riv_length*riv_width*flood_vol/cell_surface_wflow # volume expressed in meters water disc + else: + flood_vol_m = flood_vol/cell_surface_wflow flood_vol_m_data = flood_vol_m.data flood_vol_m_data[flood_vol_m.mask] = -999. logger.info('Saving water layer map to {:s}'.format(flood_vol_map)) @@ -598,9 +634,12 @@ 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) - band_inun = ds_inun.GetRasterBand(1) - + if options.out_format == 0: + ds_inun, band_inun = inun_lib.prepare_gdal(inun_file_tmp, x, y, logging=logger, srs=srs) + # band_inun = ds_inun.GetRasterBand(1) + else: + ds_inun, band_inun = inun_lib.prepare_nc(inun_file_tmp, time, x, np.flipud(y), metadata=metadata_global, + metadata_var=metadata_var, logging=logger) # loop over all the tiles n = 0 for x_loop in range(0, len(x), options.x_tile): @@ -654,12 +693,8 @@ pcr.setclone(stream_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 - # 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 + # warp of flood volume to inundation resolution inun_lib.gdal_warp(flood_vol_map, stream_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', logging=logger) x_res_tile, y_res_tile, reallength = pcrut.detRealCellLength(pcr.scalar(stream_pcr), not(bool(options.latlon))) @@ -668,7 +703,7 @@ # convert meter depth to volume [m3] flood_vol = pcr.numpy2pcr(pcr.Scalar, flood_meter*cell_surface_tile, fill_value) - # first prepare a basin map, belonging to the HAND order we are looking at + # first prepare a basin map, belonging to the lowest order we are looking at inundation_pcr = pcr.scalar(stream_pcr) * 0 for hand_strahler in range(options.catchment_strahler, stream_max + 1, 1): # hand_temp_file = os.path.join(flood_folder, 'hand_temp.map') @@ -693,7 +728,9 @@ hand_pcr = pcr.readmap(hand_temp_file) - stream_ge_hand, subcatch_hand = inun_lib.subcatch_stream(drainage_pcr, hand_strahler, stream=stream_pcr) + stream_ge_hand, subcatch_hand = inun_lib.subcatch_stream(drainage_pcr, options.catchment_strahler, stream=stream_pcr) + # stream_ge_hand, subcatch_hand = inun_lib.subcatch_stream(drainage_pcr, hand_strahler, stream=stream_pcr) + # pcr.report(subcatch_hand, 'subcatch_hand.map') stream_ge, subcatch = inun_lib.subcatch_stream(drainage_pcr, options.catchment_strahler, stream=stream_pcr, @@ -702,13 +739,20 @@ min_strahler=hand_strahler, max_strahler=hand_strahler) # generate subcatchments, only within basin for HAND flood_vol_strahler = pcr.ifthenelse(pcr.boolean(pcr.cover(subcatch, 0)), flood_vol, 0) # mask the flood volume map with the created subcatch map for strahler order = hand_strahler + # pdb.set_trace() + # pcr.report(pcr.scalar(subcatch), 'subcatch_{:02d}.map'.format(hand_strahler)) + # pcr.report(pcr.scalar(pcr.subcatchment(drainage_pcr, subcatch)), 'subcatch_upstream_{:02d}.map'.format(hand_strahler)) + + inundation_pcr_step = inun_lib.volume_spread(drainage_pcr, hand_pcr, pcr.subcatchment(drainage_pcr, subcatch), # to make sure backwater effects can occur from higher order rivers to lower order rivers flood_vol_strahler, volume_thres=0., iterations=options.iterations, cell_surface=pcr.numpy2pcr(pcr.Scalar, cell_surface_tile, -9999), - logging=logger) # 1166400000. + logging=logger, + order=hand_strahler) # 1166400000. + # pcr.report(inundation_pcr_step, 'inun_step_{:02d}.map'.format(hand_strahler)) # use maximum value of inundation_pcr_step and new inundation for higher strahler order inundation_pcr = pcr.max(inundation_pcr, inundation_pcr_step) inundation = pcr.pcr2numpy(inundation_pcr, -9999.) @@ -719,8 +763,12 @@ x_overlap_max = -inundation.shape[1] inundation_cut = inundation[0+y_overlap_min:-y_overlap_max, 0+x_overlap_min:-x_overlap_max] # inundation_cut - band_inun.WriteArray(inundation_cut, x_start, y_start) - band_inun.FlushCache() + if options.out_format == 0: + band_inun.WriteArray(inundation_cut, x_start, y_start) + band_inun.FlushCache() + else: + # with netCDF, data is up-side-down. + inun_lib.write_tile_nc(band_inun, inundation_cut, x_start, y_start) # clean up os.unlink(flood_vol_temp_file) os.unlink(drainage_temp_file) @@ -731,15 +779,19 @@ # band_inun.SetNoDataValue(-9999.) # ds_inun = None # sys.exit(0) - os.unlink(flood_vol_map) + # os.unlink(flood_vol_map) logger.info('Finalizing {:s}'.format(inun_file)) # add the metadata to the file and band - band_inun.SetNoDataValue(-9999.) - ds_inun.SetMetadata(metadata_global) - band_inun.SetMetadata(metadata_var) - ds_inun = None - ds_hand = None + # band_inun.SetNoDataValue(-9999.) + # ds_inun.SetMetadata(metadata_global) + # band_inun.SetMetadata(metadata_var) + if options.out_format == 0: + ds_inun = None + ds_hand = None + else: + ds_inun.close() + ds_ldd = None # rename temporary file to final hand file if os.path.isfile(inun_file): Index: wflow-py/Scripts/wflow_flood_lib.py =================================================================== diff -u -rd0835ddf132994a8a772d23ec5f5936f22aafea8 -r1d05f7d4cda92d45dc5f6400d5ab88020766d883 --- wflow-py/Scripts/wflow_flood_lib.py (.../wflow_flood_lib.py) (revision d0835ddf132994a8a772d23ec5f5936f22aafea8) +++ wflow-py/Scripts/wflow_flood_lib.py (.../wflow_flood_lib.py) (revision 1d05f7d4cda92d45dc5f6400d5ab88020766d883) @@ -23,6 +23,8 @@ from osgeo import osr, gdal, gdalconst import pcraster as pcr +import netCDF4 as nc +import datetime as dt import pdb def setlogger(logfilename, logReference, verbose=True): @@ -191,7 +193,48 @@ # Retrieve geoTransform info return ds, ds.GetRasterBand(band) # there's only 1 band, starting from 1 +def prepare_nc(trg_file, times, x, y, metadata={}, logging=logging, units='Days since 1900-01-01 00:00:00', calendar='gregorian'): + """ + This function prepares a NetCDF file with given metadata, for a certain year, daily basis data + The function assumes a gregorian calendar and a time unit 'Days since 1900-01-01 00:00:00' + """ + logger.info('Setting up "' + trg_file + '"') + times_list = nc.date2num(times, units=units, calendar=calendar) + nc_trg = nc.Dataset(trg_file, 'w') + logger.info('Setting up dimensions and attributes') + nc_trg.createDimension('time', 0) #NrOfDays*8 + nc_trg.createDimension('lat', len(y)) + nc_trg.createDimension('lon', len(x)) + times_nc = nc_trg.createVariable('time', 'f8', ('time',)) + times_nc.units = units + times_nc.calendar = calendar + times_nc.standard_name = 'time' + times_nc.long_name = 'time' + times_nc[:] = times_list + y_var = nc_trg.createVariable('lat', 'f4', ('lat',)) + y_var.standard_name = 'latitude' + y_var.long_name = 'latitude' + y_var.units = 'degrees_north' + x_var = nc_trg.createVariable('lon', 'f4', ('lon',)) + x_var.standard_name = 'longitude' + x_var.long_name = 'longitude' + x_var.units = 'degrees_east' + y_var[:] = y + x_var[:] = x + projection= nc_trg.createVariable('projection', 'c') + projection.long_name = 'wgs84' + projection.EPSG_code = 'EPSG:4326' + projection.proj4_params = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs' + projection.grid_mapping_name = 'latitude_longitude' + + # now add all attributes from user-defined metadata + for attr in metadata: + nc_trg.setncattr(attr, metadata[attr]) + nc_trg.sync() + return nc_trg + def prepare_gdal(filename, x, y, format='GTiff', logging=logging, + metadata={}, metadata_var={}, gdal_type=gdal.GDT_Float32, zlib=True, srs=None): # prepare geotrans xul = x[0] - (x[1] - x[0]) / 2 @@ -215,10 +258,37 @@ if srs: ds.SetProjection(srs.ExportToWkt()) # get rasterband entry + band = ds.GetRasterBand(1) + band.SetNoDataValue(-9999.) + ds.SetMetadata(metadata) + band.SetMetadata(metadata_var) + logging.info('Prepared {:s}'.format(filename)) - return ds + return ds, band +def write_tile_nc(var, data, x_start, y_start, flipud=False): + """ + + :param var: netcdf variable to write into + :param data: + :param x_start: column to start writing + :param y_start: row to start writing (is flipped up-side-down) + :return: var + """ + # determine x end and y end + x_end = x_start + data.shape[1] + y_end = y_start + data.shape[0] + if flipud: + if y_start == 0: + var[-y_end:, x_start:x_end] = data + else: + var[-y_end:-y_start, x_start:x_end] = data + else: + var[y_start:y_end] = data + return var + + def gdal_warp(src_filename, clone_filename, dst_filename, gdal_type=gdalconst.GDT_Float32, gdal_interp=gdalconst.GRA_Bilinear, format='GTiff', ds_in=None, override_src_proj=None): """ @@ -411,7 +481,7 @@ return stream_ge, pcr.ordinal(subcatch) -def volume_spread(ldd, hand, subcatch, volume, volume_thres=0., cell_surface=1., iterations=15, logging=logging): +def volume_spread(ldd, hand, subcatch, volume, volume_thres=0., cell_surface=1., iterations=15, logging=logging, order=0): """ Estimate 2D flooding from a 1D simulation per subcatchment reach Input: @@ -433,9 +503,9 @@ surface = pcr.areaarea(subcatch)*pcr.areaaverage(cell_surface, subcatch) # area_multiplier error_abs = pcr.scalar(1e10) # initial error (very high) volume_catch = pcr.areatotal(volume, subcatch) - depth_catch = volume_catch/surface # meters water disc averaged over subcatchment - + # pcr.report(depth_catch, 'depth_catch_{:02d}.map'.format(order)) + # pcr.report(volume, 'volume_{:02d}.map'.format(order)) dem_max = pcr.ifthenelse(volume_catch > volume_thres, pcr.scalar(32.), pcr.scalar(0)) # bizarre high inundation depth dem_min = pcr.scalar(0.)