diff --git a/delft3dfmpy/converters/osm_to_dflowfm.py b/delft3dfmpy/converters/osm_to_dflowfm.py new file mode 100644 index 0000000..94fe429 --- /dev/null +++ b/delft3dfmpy/converters/osm_to_dflowfm.py @@ -0,0 +1,111 @@ +import logging + +import geopandas as gpd +import numpy as np +import pandas as pd +from shapely.geometry import LineString + +from delft3dfmpy.core import checks, geometry +from delft3dfmpy.datamodels.common import ExtendedDataFrame, ExtendedGeoDataFrame + +def generate_culverts(culverts, id_col='id', roughness_values=None, logger=logging): + """ + Generate culverts from OSM to D-Flow FM + """ + + culverts_dfm = culverts.copy() + culverts_dfm['crosssection'] = [{} for _ in range(len(culverts_dfm))] + + # Set culvert id + culverts_dfm[id_col] = 'C_' + culverts_dfm[id_col] + + # Rename to id + culverts_dfm.set_index(id_col, inplace=True, drop=False) + + # Set cross-section culvert + for culvert in culverts_dfm.itertuples(): + # Generate cross section definition name + if culvert.profile == 'round': + crosssection = {'shape': 'circle', 'diameter': culvert.diameter} + + elif culvert.profile == 'boxed_rectangular': + crosssection = {'shape': 'rectangle', 'height': culvert.depth, 'width': culvert.width, + 'closed': 1} + else: + logger.debug(f'Culvert {culvert.Index} on branch {culvert.branch_id} has an unknown shape: {culvert.profile}. Applying a default profile (round - 50cm)') + crosssection = {'shape': 'circle', 'diameter': 0.50} + + # Set culvert parameters + culverts_dfm.at[culvert.Index, 'allowedflowdir'] = 'both' + culverts_dfm.at[culvert.Index, 'valveonoff'] = int(0) + culverts_dfm.at[culvert.Index, 'numlosscoeff'] = int(0) + culverts_dfm.at[culvert.Index, 'valveopeningheight'] = 0 + culverts_dfm.at[culvert.Index, 'relopening'] = 0 + culverts_dfm.at[culvert.Index, 'losscoeff'] = 0 + culverts_dfm.at[culvert.Index, 'crosssection'] = crosssection + + # Set roughness values of culvert + if culvert.material in roughness_values.keys(): + culverts_dfm.at[culvert.Index,'friction_value'] = roughness_values[culvert.material] + elif culvert.material=='concret': + culverts_dfm.at[culvert.Index, 'friction_value'] = roughness_values['concrete'] + else: + culverts_dfm.at[culvert.Index, 'friction_value'] = roughness_values['default'] + logger.debug(f'Material is not known for {culvert.id}, therefore the default friction values is used') + return culverts_dfm + + +def move_structure(struc, struc_dict, branch, offset): + """ + Function the move a structure if needed for a compound event. + + Parameters + ---------- + struc : string + current sub-structure id + struc_dict : dict + dict with all structures of a certain type + branch : string + branch id of the first structure in the compound + offset : float + chainage of the first structure in the compound + + Returns + ------- + Dict with shifted coordinates. + + """ + branch2 = struc_dict[struc]['branchid'] + if branch2!=branch: + logger.warning(f'Structures of not on the same branche. Moving structure {struc} to branch {branch}.') + struc_dict[struc]['branchid'] = branch + struc_dict[struc]['chainage'] = offset + return struc_dict + +def generate_compounds(idlist, structurelist, structures, id_col): + # probably the coordinates should all be set to those of the first structure (still to do) + compounds_dfm = ExtendedDataFrame(required_columns=[id_col,'structurelist']) + compounds_dfm.set_data(pd.DataFrame(np.zeros((len(idlist),3)), columns=[id_col,'numstructures','structurelist'], dtype='str'),index_col=id_col) + compounds_dfm.index = idlist + for ii,compound in enumerate(compounds_dfm.itertuples()): + compounds_dfm.at[compound.Index, id_col] = idlist[ii] + compounds_dfm.at[compound.Index, 'numstructures'] = len(structurelist[ii]) + + # check the substructure coordinates. If they do not coincide, move subsequent structures to the coordinates of the first + for s_i, struc in enumerate(structurelist[ii]): + if s_i == 0: + # find out what type the first structure it is and get its coordinates + if struc in structures.culverts.keys(): + branch = structures.culverts[struc]['branchid'] + offset = structures.culverts[struc]['chainage'] + else: + raise IndexError('Structure id not found. Make sure all other structures have been added to the model.') + else: + # move a subsequent structure to the location of the first + if struc in structures.culverts.keys(): + structures.culverts = move_structure(struc, structures.culverts, branch, offset) + + compounds_dfm.at[compound.Index, 'structurelist'] = ';'.join([f'{s}' for s in structurelist[ii]]) + + return compounds_dfm + diff --git a/delft3dfmpy/io/dfmreader.py b/delft3dfmpy/io/dfmreader.py index 4d7f709..a6863f7 100644 --- a/delft3dfmpy/io/dfmreader.py +++ b/delft3dfmpy/io/dfmreader.py @@ -5,7 +5,7 @@ import pandas as pd from scipy.spatial import KDTree -from delft3dfmpy.converters import hydamo_to_dflowfm +from delft3dfmpy.converters import hydamo_to_dflowfm, osm_to_dflowfm from delft3dfmpy.core import checks from delft3dfmpy.datamodels.common import ExtendedGeoDataFrame from shapely.geometry import LineString, MultiPolygon, Point, Polygon @@ -169,6 +169,38 @@ def culverts_from_hydamo(self, culverts, afsluitmiddel=None): frictionvalue=culvert.ruwheidswaarde ) + def culverts_from_osm(self, culverts, id_col='id', roughness_type = None, roughness_values=None,logger=logging): + """ + Method to convert osm culverts to dflowfm culverts. + """ + + # Convert to dflowfm input + logger.info(f'Culverts are generated from OSM data') + generated_culverts = osm_to_dflowfm.generate_culverts(culverts,id_col, roughness_values, logger=logger) + + # Add to dict + logger.info(f'Add culverts to D-Flow FM structure') + for culvert in generated_culverts.itertuples(): + self.structures.add_culvert( + id=culvert.Index, + branchid=culvert.branch_id, + chainage=culvert.branch_offset, + leftlevel=culvert.leftlevel, + rightlevel=culvert.rightlevel, + crosssection=culvert.crosssection, + length=culvert.geometry.length, + inletlosscoeff=culvert.inletlosscoeff, + outletlosscoeff=culvert.outletlosscoeff, + allowedflowdir=culvert.allowedflowdir, + valveonoff=culvert.valveonoff, + numlosscoeff=culvert.numlosscoeff, + valveopeningheight=culvert.valveopeningheight, + relopening=culvert.relopening, + losscoeff=culvert.losscoeff, + frictiontype=roughness_type, + frictionvalue=culvert.friction_value + ) + def compound_structures(self, idlist, structurelist): """ Method to add compound structures to the model. @@ -185,7 +217,6 @@ def compound_structures(self, idlist, structurelist): ) - class CrossSectionsIO: def __init__(self, crosssections): diff --git a/examples/test_osm.py b/examples/test_osm.py index b89abfe..687035e 100644 --- a/examples/test_osm.py +++ b/examples/test_osm.py @@ -50,7 +50,7 @@ # Friction type and values friction_type = parameters['frictiontype'] -friction_values = dict(zip(parameters['frictionmaterials'].split(','), list(map(float,parameters['frictionvalues'].split(','))))) +friction_values = dict(zip(parameters['frictionmaterials'].replace(' ','').split(','), list(map(float,parameters['frictionvalues'].split(','))))) # Read branches and store in OSM data model logger.info(f'Read branches') @@ -92,6 +92,7 @@ # Snap culvert to branches and determine centroid. osm.culverts.snap_to_branch(osm.branches, snap_method='ends') +# TODO: review adjusted sampling method # Determine Elevation at profiles and culvert ends with rasterio.open(fn_dem) as ds: # Elevation at profiles @@ -133,19 +134,20 @@ dfmmodel = DFlowFMModel() # Collect structures -#dfmmodel.structures.io.culverts_from_osm(osm.culverts, friction_type, friction_values, logger=logger) +dfmmodel.structures.io.culverts_from_osm(osm.culverts,id_col=id, roughness_type=friction_type, roughness_values=friction_values, logger=logger) + +# TODO: stop review here + +# Add cross-sections from OSM +#dfmmodel.crosssections.io. # TODO: CROSS SECTION DEFINITION - create circular profiles --> prof_idofbranch # TODO: CROSS SECTION DEFINITION - create rectangular profiles --> prof_idofbranch # TODO: CROSS SECTION DEFINITION - create trapezoid profiles --> prof_idofbranch +# TODO: CROSS SECTION DEFINITION - create trapezoid profiles --> prof_idofbranch # TODO: CROSS SECTION LOCATION - select rows with drain type other than culvert -# TODO: STRUCTURE - determine length -# TODO: STRUCTURE - snapping of open drains over culverts. May not be needed as wel only use parameterized profiles -# TODO: STRUCTURE - where cross sections are meeting a culvert, it may make sense to straighten the elevation value up and downstream of culvert # TODO: plot branches + cross sections locations + structures -# TODO: create DFM parsers for drains, culverts and cross sections based on dfmmodel.structures.io.xxxx_from_hydamo as example - # TODO: create 1D network logger.info(f'Create 1D network') dfmmodel.network.set_branches(osm.branches, id_col=id)