From 0932a076cb1923698c5ecd45deca8c3dd7d1e15d Mon Sep 17 00:00:00 2001 From: Joris Nettelstroth Date: Fri, 17 Jan 2025 15:00:53 +0100 Subject: [PATCH 1/7] Replace print statements with logger info and warning --- dhnx/optimization/optimization_models.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/dhnx/optimization/optimization_models.py b/dhnx/optimization/optimization_models.py index 5cb242e9..c85ab4ea 100644 --- a/dhnx/optimization/optimization_models.py +++ b/dhnx/optimization/optimization_models.py @@ -438,7 +438,7 @@ def get_invest_val(lab): if lab == str(x[0].label)] if len(outflow) > 1: - print('Multiple IDs!') + logger.info('Multiple IDs!') try: invest = res[outflow[0]]['scalars']['invest'] @@ -513,7 +513,7 @@ def get_hp_results(p): for r, c in df.iterrows(): if df.at[r, hp_lab + '.' + 'status-1'] + \ df.at[r, hp_lab + '.' + 'status-2'] > 1: - print( + logger.warning( "Investment status of pipe id {} is 1 for both dircetions!" " This is not allowed!".format(r) ) @@ -522,7 +522,7 @@ def get_hp_results(p): or\ (df.at[r, hp_lab + '.' + 'status-2'] == 1 and df.at[ r, hp_lab + '.' + 'size-2'] == 0): - print( + logger.warning( "Investment status of pipe id {} is 1, and capacity is 0!" "What happend?!".format(r) ) @@ -766,7 +766,7 @@ def solve_optimisation_investment(model): if model.settings['dump_path'] is not None: my_es = model.es my_es.dump(dpath=model.settings['dump_path'], filename=model.settings['dump_name']) - print('oemof Energysystem stored in "{}"'.format(model.settings['dump_path'])) + logger.info('oemof Energysystem stored in "{}"'.format(model.settings['dump_path'])) edges_results = model.get_results_edges() From 34d96645d362764805a7a690d42bbeaa2b18f19d Mon Sep 17 00:00:00 2001 From: Joris Nettelstroth Date: Fri, 17 Jan 2025 15:01:50 +0100 Subject: [PATCH 2/7] Fix deprecations in pandas>=2.0 --- dhnx/optimization/optimization_models.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/dhnx/optimization/optimization_models.py b/dhnx/optimization/optimization_models.py index c85ab4ea..b91dafca 100644 --- a/dhnx/optimization/optimization_models.py +++ b/dhnx/optimization/optimization_models.py @@ -446,10 +446,10 @@ def get_invest_val(lab): try: # that's in case of a one timestep optimisation due to # an oemof bug in outputlib - invest = res[outflow[0]]['sequences']['invest'][0] + invest = res[outflow[0]]['sequences']['invest'].iloc[0] except (KeyError, IndexError): - # this is in case there is no bi-directional heatpipe, e.g. at - # forks-consumers, producers-forks + # this is in case there is no bi-directional heatpipe, + # e.g. at forks-consumers, producers-forks invest = 0 # the rounding is performed due to numerical issues @@ -469,10 +469,10 @@ def get_invest_status(lab): try: # that's in case of a one timestep optimisation due to # an oemof bug in outputlib - invest_status = res[outflow[0]]['sequences']['invest_status'][0] + invest_status = res[outflow[0]]['sequences']['invest_status'].iloc[0] except (KeyError, IndexError): - # this is in case there is no bi-directional heatpipe, e.g. at - # forks-consumers, producers-forks + # this is in case there is no bi-directional heatpipe, + # e.g. at forks-consumers, producers-forks invest_status = 0 return invest_status @@ -667,7 +667,7 @@ def optimize_operation(thermal_network): def setup_optimise_investment( thermal_network, invest_options, heat_demand='scalar', num_ts=1, - time_res=1, start_date='1/1/2018', frequence='H', solver='cbc', + time_res=1, start_date='1/1/2018', frequence='h', solver='cbc', solve_kw=None, solver_cmdline_options=None, simultaneity=1, bidirectional_pipes=False, dump_path=None, dump_name='dump.oemof', print_logging_info=False, write_lp_file=False): From d9a05a5e90e0c5be96b63548f2aec55fdbc415c7 Mon Sep 17 00:00:00 2001 From: Joris Nettelstroth Date: Fri, 17 Jan 2025 15:13:32 +0100 Subject: [PATCH 3/7] Enforce reduction of 3D geometry to 2D Most underlying functions (e.g. distance measurement in shapely) do not support 3D geometries. Having z-coordinates in the data can cause issues in several places in the code, especially when e.g. buildings contain z coordinates, but streets do not. The savest option currently is to force 2d on all input geometries, which is now the default in check_crs() --- dhnx/gistools/geometry_operations.py | 32 +++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/dhnx/gistools/geometry_operations.py b/dhnx/gistools/geometry_operations.py index 38e21a44..e7950a2a 100644 --- a/dhnx/gistools/geometry_operations.py +++ b/dhnx/gistools/geometry_operations.py @@ -20,6 +20,7 @@ print("Need to install geopandas to process geometry data.") try: + import shapely from shapely import wkt from shapely.geometry import LineString from shapely.geometry import MultiLineString @@ -510,16 +511,45 @@ def drop_parallel_lines(gdf): return gdf -def check_crs(gdf, crs=4647): +def check_crs(gdf, crs=4647, force_2d=True): """Convert CRS to EPSG:4647 - ETRS89 / UTM zone 32N (zE-N). This is the (only?) Coordinate Reference System that gives the correct results for distance calculations. + Note about enforcing 2D geometries: + + Most underlying functions (e.g. distance measurement in shapely) do + not support 3D geometries. Having z-coordinates in the + data can cause issues in several places in the code, especially when + e.g. buildings contain z coordinates, but streets do not. + + The savest option currently is to force 2d on all input geometries. + + Parameters + ---------- + gdf : GeoDataFrame + The GeoDataFrame to update. + crs : int (optional) + EPSG code for a coordinate reference system to convert to. + Default is 4647. + force_2d : boolean (optional) + If True, enforce reduction of 3D geometry to 2D. Default is True. + + Returns + ------- + gdf : GeoDataFrame + Updated GeoDataFrame. + """ if gdf.crs.to_epsg() != crs: gdf.to_crs(epsg=crs, inplace=True) logger.info('CRS of GeoDataFrame converted to EPSG:{0}'.format(crs)) + if force_2d and gdf.has_z.any(): + logger.debug("Reducing 3D geometry to 2D for compatibility") + # Alternatively, use "gdf.force_2d()" with geopandas>0.14.3 + gdf.geometry = shapely.force_2d(gdf.geometry) # requires shapely>=2.0 + return gdf From 546e8ec04cac36d5b1bf23d2cb6abd06d6ece771 Mon Sep 17 00:00:00 2001 From: Joris Nettelstroth Date: Fri, 17 Jan 2025 15:19:51 +0100 Subject: [PATCH 4/7] Fix some edge cases in weld_segments() --- dhnx/gistools/geometry_operations.py | 55 +++++++++++++++++++++++----- 1 file changed, 46 insertions(+), 9 deletions(-) diff --git a/dhnx/gistools/geometry_operations.py b/dhnx/gistools/geometry_operations.py index e7950a2a..2eafcbab 100644 --- a/dhnx/gistools/geometry_operations.py +++ b/dhnx/gistools/geometry_operations.py @@ -184,7 +184,7 @@ def check_double_points(gdf, radius=0.001, id_column=None): count += 1 if count > 0: - logger.info('Number of duplicated points: ', count) + logger.info('Number of duplicated points: %s', count) else: logger.info( 'Check passed: No points with a distance closer than {}'.format(radius)) @@ -279,7 +279,7 @@ def weld_segments(gdf_line_net, gdf_line_gen, gdf_line_houses, """Weld continuous line segments together and cut loose ends. This is a public function that recursively calls the internal function - weld_line_segments_(), until the problem cannot be simplified further. + _weld_segments(), until the problem cannot be simplified further. Find all lines that only connect to one other line and connect those to a single MultiLine object. Points that connect to Generators and @@ -366,11 +366,26 @@ def debug_plot(neighbours, color='red'): # Drop this object, because it is contained within a merged object continue # Continue with the next line segment - # Find all neighbours of the current segment - neighbours = gdf_line_net[gdf_line_net.geometry.touches(geom)] - # If all of the neighbours intersect with each other, it is the - # last segement before an intersection, which can be removed - if all([all(neighbours.geometry.intersects(neighbour)) + if geom.is_ring: + # Drop this object, because rings are not valid options for street + # segments, since they have no start or end + continue # Continue with the next line segment + + # Find all neighbours of the current segment. + # Neighbours are geometries whose boundary "touches" the current + # geometry. We need to use the boundary to allow roads that cross + # themselves, e.g. for spiral ramps. Also they must not be equal + # to the current segment. + neighbours = gdf_line_net[ + (gdf_line_net.geometry.boundary.touches(geom) + & ~gdf_line_net.geometry.geom_equals(geom) + )] + # If all of the neighbours touch each other, it is the + # last segment before an intersection, which can be removed. + # The tests needs to be "touches" OR "equals", since per definition + # a line geometry cannot "touch" itself + if all([all(neighbours.geometry.touches(neighbour) + | neighbours.geometry.geom_equals(neighbour)) for neighbour in neighbours.geometry]): # Treat as if there was only one neighbour (like end segment) neighbours = neighbours.head(1) @@ -385,12 +400,34 @@ def debug_plot(neighbours, color='red'): p2 = geom.boundary.geoms[-1] p1_neighbours = neighbours.geometry.intersects(p1).to_list() p2_neighbours = neighbours.geometry.intersects(p2).to_list() + if (any(gdf_line_ext.geometry.touches(p1)) and p2_neighbours.count(True) > 0): unused = False elif (any(gdf_line_ext.geometry.touches(p2)) and p1_neighbours.count(True) > 0): unused = False + elif (any(gdf_line_ext.geometry.touches(geom)) + and not any(gdf_line_net.geometry.touches(geom))): + # The current segment is touched by an external line, but not + # by any other network segment. Select connected external line + geom_ext = unary_union(gdf_line_ext[ + gdf_line_ext.geometry.touches(geom)].geometry) + # Test if the network line touched by the external line + # is equal to the current line segement, i.e. the external + # line is not connected to any other network line + if gdf_line_net[gdf_line_net.geometry.touches(geom_ext) + ].geometry.geom_equals(geom).all(): + logger.warning( + "Welding is about to remove a street network line " + "segment that is connected to nothing but a building " + "connection line. This would leave the building " + "unconnected, so the segment is not removed. " + "This indicates a bug in the welding logic or an " + "issue in the input data, e.g. an initial street " + "network where not all lines are connected. If the " + "remaining process fails, this may be the cause.") + unused = False # Keep line, despite not being useful if unused: # If truly unused, we can discard it to simplify the network @@ -427,7 +464,7 @@ def debug_plot(neighbours, color='red'): # There are excactly two separate neighbours that can be merged pass # Run the rest of the loop - # Before merging, we need to futher clean up the list of neighbours + # Before merging, we need to further clean up the list of neighbours neighbours_list = [] for neighbour in neighbours.geometry: if any(gdf_deleted.geometry.geom_equals(neighbour)): @@ -491,7 +528,7 @@ def drop_parallel_lines(gdf): These can be two actually distinct paths between two points, or two identical lines on top of each other. When downloading streets with osmnx, - this can intruduce such duplicates where the two have the attributes + this can introduce such duplicates where the two have the attributes 'reversed=True' and 'reversed=False' This function modifies the GeoDataFrame in place and resets the index. From 153cd11fc8c4ad42803274bef0f96474a82cc4c3 Mon Sep 17 00:00:00 2001 From: Joris Nettelstroth Date: Fri, 17 Jan 2025 15:41:26 +0100 Subject: [PATCH 5/7] Simplify and improve run_point_method_boundary() Previously, run_point_method_boundary() required consumer and producer data as input and repeated all steps for both data sets. To avoid unnecessary code repetition, the duplicates were removed and the function now has to be called once with consumer data and once with producer data as input. Additionally, incorrect behaviour for some edge cases was fixed and support for the multiple connection lines implemented in create_object_connections() was added. Usage of the function weld_segments(), which welds continuous line segments together and cuts loose ends has been made optional via the parameter 'welding' (default is True) The function process_geometry() had to be adapted to support the changes in and run_point_method_boundary() --- dhnx/gistools/connect_points.py | 159 +++++++++++++++++++------------- 1 file changed, 93 insertions(+), 66 deletions(-) diff --git a/dhnx/gistools/connect_points.py b/dhnx/gistools/connect_points.py index 2e4186f4..7e409e0d 100644 --- a/dhnx/gistools/connect_points.py +++ b/dhnx/gistools/connect_points.py @@ -322,11 +322,12 @@ def create_points_from_polygons(gdf, method='midpoint'): ) -def run_point_method_boundary( - consumers_poly, consumers, producers_poly, producers, - lines_consumers, lines_producers): +def run_point_method_boundary(consumers_poly, consumers, lines_consumers): """Run 'boundary' method for finding the building connection point. + This is meant to be called once with the consumers as input and + once with the producers (producers_poly, producers, lines_producers). + The 'midpoint' method (using the centroid) must already have been run, generating the default connection lines from street to centroid. @@ -352,29 +353,34 @@ def run_point_method_boundary( but they are not changed. consumers : geopandas.GeoDataFrame Points of the consumer buildings (as returned by 'midpoint' method). - producers_poly : geopandas.GeoDataFrame - Polygons of the producer buildings. Point geometries are also allowed, - but they are not changed. - producers : geopandas.GeoDataFrame - Points of the producer buildings (as returned by 'midpoint' method). lines_consumers : geopandas.GeoDataFrame Connection lines from street to each consumer point. - lines_producers : geopandas.GeoDataFrame - Connection lines from street to each producer point. Returns ------- consumers : geopandas.GeoDataFrame Updated points of the consumer buildings. - producers : geopandas.GeoDataFrame - Updated points of the producer buildings. lines_consumers : geopandas.GeoDataFrame Updated connection lines from street to each consumer point. - lines_producers : geopandas.GeoDataFrame - Updated connection lines from street to each producer point. """ logger.info('Run "boundary" method for finding the building connections') + # lines_consumers may represent multiple lines per consumer + # Duplicate geometries in consumers_poly have to be created accordingly + consumers_poly['id_full'] = consumers['id_full'] + consumers_poly = pd.merge( + left=consumers_poly, + right=lines_consumers.drop(columns=[lines_consumers.geometry.name]), + how='right', on='id_full') + + # When using the original consumer points as a fallback later, we + # require it to have the same index as lines_consumers. Therefore + # we create the 'duplicate' consumers object + consumers_d = pd.merge( + left=consumers, + right=lines_consumers.drop(columns=[lines_consumers.geometry.name]), + how='right', on='id_full') + # Cut the part off of each "line_consumer" that is within the building # polygon. As a result, the heating grid will only reach to the wall of # the building. @@ -387,13 +393,10 @@ def run_point_method_boundary( lines_consumers_n.loc[lines_consumers_n.type == "MultiLineString"] = \ gpd.GeoDataFrame(geometry=lines_consumers.difference( consumers_poly.convex_hull, align=False)) - - # Repeat for the producer lines - lines_producers_n = gpd.GeoDataFrame(geometry=lines_producers.difference( - producers_poly, align=False)) - lines_producers_n.loc[lines_producers_n.type == "MultiLineString"] = \ - gpd.GeoDataFrame(geometry=lines_producers.difference( - producers_poly.convex_hull, align=False)) + # Only keep the new consumer lines if they have a useful minimum length. + # There was an edgecase where a street 'almost' touched a building, + # and the cut consumer line had a length of 1e-9 m + lines_consumers_n[lines_consumers_n.length < 1e-3] = LineString() # Now the "consumers" (point objects for each building) need to be # updated to touch the end of the consumer_lines @@ -403,13 +406,6 @@ def run_point_method_boundary( gpd.GeoDataFrame(geometry=lines_consumers.intersection( consumers_poly.convex_hull.boundary, align=False)) - # Repeat for the producers - producers_n = gpd.GeoDataFrame(geometry=lines_producers.intersection( - producers_poly.boundary, align=False)) - producers_n.loc[producers_n.type == "MultiPoint"] = \ - gpd.GeoDataFrame(geometry=lines_producers.intersection( - producers_poly.convex_hull.boundary, align=False)) - # Sometimes the centroid does not lie within a building and there may be # no intersetions, i.e. the new point is an 'empty' geometry. This can # happen if buildings are multipolygons, which is not forbidden. @@ -426,25 +422,48 @@ def run_point_method_boundary( # equals the starting point of the connection line. mask2 = consumers_n.geom_equals( lines_consumers.geometry.apply(lambda line: line.boundary.geoms[0])) - mask = mask1 | mask2 - consumers_n.loc[mask] = consumers.loc[mask].geometry - lines_consumers_n.loc[mask] = lines_consumers.loc[mask].geometry - # Repeat for the producers - mask1 = (producers_n.is_empty | lines_producers_n.is_empty) - mask2 = producers_n.geom_equals( - lines_producers.geometry.apply(lambda line: line.boundary.geoms[0])) - mask = mask1 | mask2 - producers_n.loc[mask] = producers.loc[mask].geometry - lines_producers_n.loc[mask] = lines_producers.loc[mask].geometry + # If for whatever reason the street-side "start" of the new connection + # line is not the same point as the original connection line start, use + # the original line. This may happen for complex geometries, where the + # street line lies within the building geometry + lines_consumers_n_start = lines_consumers_n.copy() + lines_consumers_n_start.geometry = ( + lines_consumers_n_start[~mask1 & ~mask2].boundary.apply( + lambda g: g.geoms[0])) + mask3 = lines_consumers_n_start.geom_equals( + lines_consumers.geometry.apply(lambda line: line.boundary.geoms[0])) + + # Now apply all the filters above to reset the geometries + mask = mask1 | mask2 | ~mask3 + consumers_n.loc[mask] = consumers_d.loc[mask].geometry + lines_consumers_n.loc[mask] = lines_consumers.loc[mask].geometry # Now update all the original variables with the new data + lines_consumers_n['id_full'] = lines_consumers['id_full'] + consumers_n['id_full'] = lines_consumers_n['id_full'] + + # If multiple building connection lines existed before, we now also + # have created multiple building points for each building. + # We need to group those points into one MultiPoint per initial unique + # building, to keep the original index structure intact. + consumers_n = (pd.DataFrame(consumers_n) # convert gdf to df + .groupby('id_full', sort=False, as_index=False) + .agg(lambda x: MultiPoint(x.values)) # returns df + .set_geometry(consumers_n.geometry.name, # convert to gdf + crs=consumers_n.crs) + ) + + # For each new consumer point(s), test if they actually touch + # the new conumser line(s) that have the same 'id_full' assigned + for id_full, points in zip(consumers_n['id_full'], consumers_n.geometry): + if not points.touches(lines_consumers_n[ + lines_consumers_n['id_full'] == id_full].geometry).all(): + raise ValueError(f"Points from {id_full} have no matching lines") + consumers.geometry = consumers_n.geometry - producers.geometry = producers_n.geometry lines_consumers = lines_consumers_n - lines_producers = lines_producers_n - - return consumers, producers, lines_consumers, lines_producers + return consumers, lines_consumers def check_duplicate_geometries(gdf): @@ -467,7 +486,8 @@ def check_duplicate_geometries(gdf): def process_geometry(lines, consumers, producers, method='midpoint', projected_crs=4647, - tol_distance=2, reset_index=True): + tol_distance=2, reset_index=True, + welding=True): """ This function connects the consumers and producers to the line network, and prepares the attributes of the geopandas.GeoDataFrames for importing as dhnx.ThermalNetwork. @@ -501,6 +521,10 @@ def process_geometry(lines, consumers, producers, If True, reset the index and ignore the existing index. If False, use the existing index for consumer and producer identificators. Default: True + welding : bool, optional + Weld continuous line segments together and cut loose ends. This + can improve the performance of the optimization, as it decreases + the total number of line elements. Default is True. Returns ------- @@ -509,14 +533,17 @@ def process_geometry(lines, consumers, producers, 'producers', 'pipes'. """ - if method == 'boundary': - # copies of the original polygons are needed for method 'boundary' - consumers_poly = go.check_crs(consumers, crs=projected_crs).copy() - producers_poly = go.check_crs(producers, crs=projected_crs).copy() + if not reset_index: + raise ValueError("Keeping the orginal index is currently not " + "supported. Use 'reset_index=True'.") + + # Copies of the original polygons are needed for method 'boundary' + consumers_poly = go.check_crs(consumers, crs=projected_crs).copy() + producers_poly = go.check_crs(producers, crs=projected_crs).copy() # check whether the expected geometry is used for geo dataframes check_geometry_type(lines, types=['LineString', 'MultiLineString']) - for gdf in [producers, consumers]: + for gdf in [producers, consumers, producers_poly, consumers_poly]: check_geometry_type(gdf, types=['Polygon', 'Point', 'MultiPolygon']) check_duplicate_geometries(gdf) @@ -529,17 +556,18 @@ def process_geometry(lines, consumers, producers, for layer in [producers, consumers]: layer = go.check_crs(layer, crs=projected_crs) layer = create_points_from_polygons(layer, method=method) + layer['lat'] = layer['geometry'].apply(lambda x: x.y) + layer['lon'] = layer['geometry'].apply(lambda x: x.x) + + for layer in [producers, consumers, producers_poly, consumers_poly]: if reset_index: layer.reset_index(inplace=True, drop=True) layer.index.name = 'id' - if 'id' in layer.columns: - layer.drop(['id'], axis=1, inplace=True) + layer.drop(columns=['id'], inplace=True, errors='ignore') else: if layer.index.has_duplicates: raise ValueError("The index of input data has duplicate " "values, which is not allowed") - layer['lat'] = layer['geometry'].apply(lambda x: x.y) - layer['lon'] = layer['geometry'].apply(lambda x: x.x) producers['id_full'] = 'producers-' + producers.index.astype('str') producers['type'] = 'G' @@ -558,16 +586,18 @@ def process_geometry(lines, consumers, producers, if method == 'boundary': # Can only be performed after 'midpoint' method - consumers, producers, lines_consumers, lines_producers = ( - run_point_method_boundary( - consumers_poly, consumers, producers_poly, producers, - lines_consumers, lines_producers)) - - # Weld continuous line segments together and cut loose ends - lines = go.weld_segments( - lines, lines_producers, lines_consumers, - # debug_plotting=True, - ) + consumers, lines_consumers = run_point_method_boundary( + consumers_poly, consumers, lines_consumers) + producers, lines_producers = run_point_method_boundary( + producers_poly, producers, lines_producers) + + if welding: + # Weld continuous line segments together and cut loose ends + lines = go.weld_segments( + lines, lines_producers, lines_consumers, + # debug_plotting=True, + ) + # Keep only the shortest of all lines connecting the same two points lines = go.drop_parallel_lines(lines) @@ -584,8 +614,7 @@ def process_geometry(lines, consumers, producers, lines_all.reset_index(inplace=True, drop=True) if reset_index: lines_all.index.name = 'id' - if 'id' in lines_all.columns: - lines_all.drop(['id'], axis=1, inplace=True) + lines_all.drop(columns=['id'], inplace=True, errors='ignore') # concat point layer points_all = pd.concat([ @@ -594,8 +623,6 @@ def process_geometry(lines, consumers, producers, forks[['id_full', 'geometry']]], sort=False ) - points_all['geo_wkt'] = points_all['geometry'].apply(lambda x: x.wkt) - points_all.set_index('geo_wkt', drop=True, inplace=True) # add from_node, to_node to lines layer lines_all = go.insert_node_ids(lines_all, points_all) From 7da0502f30bf9345490ea7faf37a3e9ae7490492 Mon Sep 17 00:00:00 2001 From: Joris Nettelstroth Date: Fri, 17 Jan 2025 15:48:33 +0100 Subject: [PATCH 6/7] Allow multiple potential building connections DHNx requires connecting each building to the line network, e.g. the streets. To do this, create_object_connections() normally finds the perpendicular line which represents the closest connection between building and street. This commit updates create_object_connections() to allow more than 1 (default value) number of connections 'n_conn' to be created. During the optimization phase, the solver now has multiple options to choose from when finding the best district heating network. This may be desired when the nearest street segment is not actually the best way to connect a building. In cases where a building has a street or path in the front and back, both may be valid options. For each line segment that a possible connection has been attached to, its direct connections are ignored when searching for the next possible line segment. This mostly avoids obviously inferior alternative connections and helps find more relevant alternatives. To allow this functionality, create_object_connections() had to be completely rewritten. But the default settings should retain the old behaviour. This feature increases the optimization time. --- dhnx/gistools/connect_points.py | 241 ++++++++++++----------- dhnx/gistools/geometry_operations.py | 45 +++-- dhnx/optimization/optimization_models.py | 9 +- 3 files changed, 167 insertions(+), 128 deletions(-) diff --git a/dhnx/gistools/connect_points.py b/dhnx/gistools/connect_points.py index 7e409e0d..0bd1c26b 100644 --- a/dhnx/gistools/connect_points.py +++ b/dhnx/gistools/connect_points.py @@ -116,8 +116,9 @@ def calc_lot_foot(line, point): return lot_foot_point -def create_object_connections(points, lines, tol_distance=1): - """Connects points to a line network. +def create_object_connections(points, lines, tol_distance=1, n_conn=1, + drop_neighbours=True): + """Connect points to a line network. Generally, the nearest point of the next line is used as connection the point. Depending on the geometry, there are 3 options, the connection is created: @@ -147,6 +148,17 @@ def create_object_connections(points, lines, tol_distance=1): which contain more than 2 points are not allowed. tol_distance : float Tolerance distance for choosing the end of the line instead of the nearest point. + n_conn : int, optional + Number of connection lines created from each consumer/producer to + the nearest line segments in the street network. This allows the + placement of the connection lines to be part of the optimization + process. The default is 1. + drop_neighbours : bool, optional + (Only relevant for n_conn>1). When searching the next connection + line, ignore the neighbour segements of the last connection segment. + Chances are that they do not provide an advantage + over the nearest segment. This allows the next connection + line to find a more relevant alternative. Default is True. Returns ------- @@ -155,120 +167,115 @@ def create_object_connections(points, lines, tol_distance=1): All lines should only touch at the line endings. """ - # check linestrings - for _, c in lines.iterrows(): - if len(c['geometry'].coords) > 2: - raise ValueError("The Linestrings must consists of simple lines," - " with only two coordinates!") - - # Pepare merging all the street lines - all_lines = lines['geometry'] - - # There seems to be a conflict between shapely and pygeos, - # which both use 'geos' internally, if both are installed. - # This can cause - # 'OSError exception: access violation reading 0xFFFFFFFFFFFFFFFF'. - # - # With shapely 1.8.0 and pygeos 0.12.0 it was observed that - # this sometimes even fails without error. In such a case - # mergedlines might only contain a single LineString (one street - # segment) instead of a MultiLineString (the combined network - # of all street segments). This completely messes up the - # following nearest_points(). - # - # Wrapping the argument in 'list()' seems to be a valid workaround. - # It may come with a performance cost, as noted here: - # https://github.com/geopandas/geopandas/issues/1820 - # https://github.com/geopandas/geopandas/issues/2171 - # This issue may disappear when shapely 2.0 is released (then pygeos - # is merged with shapely). - mergedlines = unary_union(list(all_lines)) - # mergedlines = unary_union(all_lines) # TODO Try this with shapely 2.0 - - # empty geopandas dataframe for house connections - conn_lines = gpd.GeoDataFrame(geometry=[], crs=lines.crs) - - # iterate over all houses - for index, row in points.iterrows(): - - house_geo = row['geometry'] - - # new nearest point method ############ ######### - n_p = nearest_points(mergedlines, house_geo)[0] - - # get index of line which is closest to the house - line_index = line_of_point(n_p, lines) - - # get geometry of supply line - supply_line = lines.loc[line_index, 'geometry'] - - # get end points of line - supply_line_p0 = Point(list(supply_line.coords)[0]) - supply_line_p1 = Point(list(supply_line.coords)[1]) - supply_line_points = [supply_line_p0, supply_line_p1] - supply_line_mulitpoints = MultiPoint(supply_line_points) - - if n_p in supply_line_points: - # case that nearest point is a line ending - - logger.debug( - 'Connect buildings... id {}: ' - 'Connected to supply line ending (nearest point)'.format(index) - ) - - con_line = LineString([n_p, house_geo]) - - conn_lines = pd.concat( - [gpd.GeoDataFrame(conn_lines, crs=lines.crs), - gpd.GeoDataFrame(geometry=[con_line], crs=lines.crs)], - ignore_index=True) - + logger.debug("Create connections from street to buildings") + + def create_object_connection(point_geom, lines, tol_distance=tol_distance): + # Find the nearest line and its nearest point + lines_merged = unary_union(lines.geometry) + nearest_line_point = nearest_points(point_geom, lines_merged)[1] + nearest_line_idx = lines[nearest_line_point.distance(lines.geometry) + < 1e-8].index + nearest_line = lines.geometry[nearest_line_idx[0]] + + # Check if the nearest point is an end point of the line + line_start, line_end = nearest_line.boundary.geoms + if (nearest_line_point.equals(line_start) + or nearest_line_point.equals(line_end)): + connection_point = nearest_line_point else: - - dist_to_endings = [x.distance(n_p) for x in supply_line_points] - - if min(dist_to_endings) >= tol_distance: - # line is split, no line ending is close to the nearest point - # this also means the original supply line needs to be deleted - - logger.debug( - 'Connect buildings... id {}: Supply line split'.format(index)) - - con_line = LineString([n_p, house_geo]) - - conn_lines = pd.concat( - [gpd.GeoDataFrame(conn_lines, crs=lines.crs), - gpd.GeoDataFrame(geometry=[con_line], crs=lines.crs)], - ignore_index=True) - - lines.drop([line_index], inplace=True) - + # Check if the distance of nearest_point_on_line is close + # to an existing point on the line + points_on_line = nearest_line.boundary + closest_existing_point = nearest_points( + nearest_line_point, points_on_line)[1] + dist_on_line = nearest_line_point.distance(closest_existing_point) + if dist_on_line <= tol_distance: + connection_point = closest_existing_point + # If connection_point changes, the nearest lines have to be + # updated. There are probably two nearest lines instead of one + nearest_line_idx = lines[ + connection_point.distance(lines.geometry) < 1e-8].index + else: + # Split the line and use the nearest point as connection point + connection_point = nearest_line_point + + # Create connection line (Direction: From street to building) + conn_line = LineString([connection_point, point_geom]) + return conn_line, nearest_line_idx + + conn_lines_list = [] + for id_full, point_geom in zip(points['id_full'], points.geometry): + # For each point (building), find the n closest connection lines + # to the street lines, by dropping the previous closest street + # sections before searching the next connection line + lines_drop = [] + conn_lines_id = [] + for i in range(n_conn): + conn_line, nearest_line_idx = create_object_connection( + point_geom, + lines.drop(lines_drop), + tol_distance=tol_distance) + + lines_drop.extend(nearest_line_idx) + if drop_neighbours and i+1 < n_conn: + # Find all neighbours of the nearest segment(s) and delete them + # as well. Chances are that they do not provide an advantage + # over the nearest segment. This allows the next connection + # line to find a more relevant alternative + neighbours = lines[lines.touches( + unary_union(lines.geometry[nearest_line_idx]))] + lines_drop.extend(neighbours.index) + neighbours2 = lines[lines.touches( + unary_union(neighbours.geometry))] + lines_drop.extend(neighbours2.index) + + conn_lines_id.append(conn_line) + + # Create a GeoDataFrame with the connection lines of the current id + gdf_conn_lines_id = gpd.GeoDataFrame( + data={'id_full': [id_full]*len(conn_lines_id)}, + geometry=conn_lines_id, crs=lines.crs) + + # Drop duplicate geometries in connection lines + gdf_conn_lines_id = gdf_conn_lines_id.drop_duplicates( + subset=gdf_conn_lines_id.geometry.name).reset_index(drop=True) + # Make sure that the street network is split into new sections + # where the connection lines meet the street lines + for conn_line in gdf_conn_lines_id.geometry: + conn_point = conn_line.boundary.geoms[0] + + if conn_point.within(lines.geometry.boundary).any(): + # Connection point is one of the existing street line points + continue + else: + nearest_line_idx = lines[conn_point.distance(lines.geometry) + < 1e-8].index[0] + nearest_line = lines.geometry[nearest_line_idx] + # Identify start and end points for new lines + line_start, line_end = nearest_line.boundary.geoms + # Store any existing attributes from original data + attributes = lines.loc[[nearest_line_idx]].drop( + columns=[lines.geometry.name]) + # Remove the line that is to be replaced + lines.drop([nearest_line_idx], inplace=True) + # Combine remaining lines with two new replacement lines lines = pd.concat( [gpd.GeoDataFrame(lines, crs=lines.crs), - gpd.GeoDataFrame(geometry=[ - LineString([supply_line_p0, n_p]), - LineString([n_p, supply_line_p1])], crs=lines.crs)], + gpd.GeoDataFrame( + geometry=[ + LineString([line_start, conn_point]), + LineString([conn_point, line_end])], + crs=lines.crs, + data=pd.concat([attributes, attributes]), # keep data + ), + ], ignore_index=True) - else: - # case that one or both line endings are closer than tolerance - # thus, the next line ending is chosen - logger.debug( - 'Connect buildings... id {}: Connected to Supply line ' - 'ending due to tolerance'.format(index)) - - conn_point = nearest_points(supply_line_mulitpoints, n_p)[0] - - con_line = LineString([conn_point, house_geo]) - - conn_lines = pd.concat( - [gpd.GeoDataFrame(conn_lines, crs=lines.crs), - gpd.GeoDataFrame(geometry=[con_line], crs=lines.crs)], - ignore_index=True) + conn_lines_list.append(gdf_conn_lines_id) - logger.info('Connection of buildings completed.') + gdf_conn_lines = pd.concat(conn_lines_list, ignore_index=True) - return conn_lines, lines + return gdf_conn_lines, lines def check_geometry_type(gdf, types): @@ -486,7 +493,7 @@ def check_duplicate_geometries(gdf): def process_geometry(lines, consumers, producers, method='midpoint', projected_crs=4647, - tol_distance=2, reset_index=True, + tol_distance=2, reset_index=True, n_conn=1, n_conn_prod=1, welding=True): """ This function connects the consumers and producers to the line network, and prepares the @@ -521,6 +528,16 @@ def process_geometry(lines, consumers, producers, If True, reset the index and ignore the existing index. If False, use the existing index for consumer and producer identificators. Default: True + n_conn : int, optional + Number of connection lines created from each consumer to + the nearest line segments in the street network. This allows the + placement of the connection lines to be part of the optimization + process. The default is 1. + n_conn_prod : int, optional + Number of connection lines created from each producer to + the nearest line segments in the street network. This allows the + placement of the connection lines to be part of the optimization + process. The default is 1. welding : bool, optional Weld continuous line segments together and cut loose ends. This can improve the performance of the optimization, as it decreases @@ -576,9 +593,9 @@ def process_geometry(lines, consumers, producers, # Add lines to consumers and producers lines_consumers, lines = create_object_connections( - consumers, lines, tol_distance=tol_distance) + consumers, lines, tol_distance=tol_distance, n_conn=n_conn) lines_producers, lines = create_object_connections( - producers, lines, tol_distance=tol_distance) + producers, lines, tol_distance=tol_distance, n_conn=n_conn_prod) if not reset_index: # Connection lines are ordered the same as points. Match their indexes lines_consumers.index = consumers.index diff --git a/dhnx/gistools/geometry_operations.py b/dhnx/gistools/geometry_operations.py index 2eafcbab..e1c90080 100644 --- a/dhnx/gistools/geometry_operations.py +++ b/dhnx/gistools/geometry_operations.py @@ -103,30 +103,53 @@ def insert_node_ids(lines, nodes): ------- geopandas.GeoDataFrame """ - # add id to gdf_lines for starting and ending node - # point as wkt - lines['b0_wkt'] = lines["geometry"].apply( - lambda geom: geom.boundary.geoms[0].wkt) - lines['b1_wkt'] = lines["geometry"].apply( - lambda geom: geom.boundary.geoms[-1].wkt) + nodes['geo_wkt'] = nodes.geometry.apply( + lambda x: wkt.dumps(x, output_dimension=2)) + nodes.set_index('geo_wkt', drop=True, inplace=True) + + # add id to gdf_lines for starting and ending node point as wkt + lines['b0_wkt'] = lines.geometry.apply( + lambda geom: wkt.dumps(geom.boundary.geoms[0], output_dimension=2)) + lines['b1_wkt'] = lines.geometry.apply( + lambda geom: wkt.dumps(geom.boundary.geoms[-1], output_dimension=2)) + + def match_multipoint(point_wkt): + """Return id_full from matching 'nodes' for each point in point_wkt. + + This is necessary if point_wkt contains MultiPoint objects, which + result from multiple connection lines, and not only single Point + objects. + """ + point_line = wkt.loads(point_wkt) + for point_node in nodes.index: + if point_line.within(wkt.loads(point_node)): + return nodes.loc[point_node, 'id_full'] + logger.error("Point not found: %s", point_wkt) + return False try: lines['from_node'] = lines['b0_wkt'].apply( lambda x: nodes.at[x, 'id_full']) - lines['to_node'] = lines['b1_wkt'].apply( - lambda x: nodes.at[x, 'id_full']) + lines['to_node'] = lines['b1_wkt'].apply(lambda x: match_multipoint(x)) except KeyError as e: errors = ([wkt.loads(x) for x in lines['b0_wkt'] if x not in nodes['id_full']]) errors.extend([wkt.loads(x) for x in lines['b1_wkt'] - if x not in nodes['id_full']]) + if not match_multipoint(x)]) gdf_errors = gpd.GeoDataFrame(geometry=errors, crs=lines.crs) ax = lines.plot() gdf_errors.plot(ax=ax, color='red', label='Point(s) causing error') plt.legend() plt.show() - # gdf_errors.to_file('debug_points.geojson') - # lines.to_file('debug_lines.geojson') + try: + nodes.to_file('debug_nodes.geojson') + gdf_errors.to_file('debug_error_points.geojson') + # lines[[lines.geometry.name, 'type', 'b0_wkt', 'b1_wkt'] + # ].to_file('debug_lines.geojson') + except Exception as e: + breakpoint() + logger.error(e) + raise KeyError("This error indicates specific problems with the data. " "A plot of the problematic point(s) is shown.") from e diff --git a/dhnx/optimization/optimization_models.py b/dhnx/optimization/optimization_models.py index b91dafca..19f9e2c1 100644 --- a/dhnx/optimization/optimization_models.py +++ b/dhnx/optimization/optimization_models.py @@ -184,11 +184,10 @@ def check_input(self): "The consumer {} of pipe id {} does not exist!" .format(cons_id, p)) - pipe_to_cons_ids = list(self.thermal_network.components['pipes']['to_node'].values) - pipe_to_cons_ids = [x.split('-', 1)[1] for x in pipe_to_cons_ids - if x.split('-', 1)[0] == 'consumers'] - - for id in list(self.thermal_network.components['consumers'].index): + pipe_to_cons_ids = list( + self.thermal_network.components['pipes']['to_node'].values) + pipe_to_cons_ids = [x for x in pipe_to_cons_ids if 'consumers' in x] + for id in list(self.thermal_network.components['consumers']['id_full']): if id not in pipe_to_cons_ids: raise ValueError( "The consumer id {} has no connection the the grid!".format(id)) From 44bf390b8839835e0067b9c7a292d6d7215e751b Mon Sep 17 00:00:00 2001 From: Joris Nettelstroth Date: Fri, 17 Jan 2025 16:17:26 +0100 Subject: [PATCH 7/7] Fix flake8 issues --- dhnx/gistools/connect_points.py | 21 +++++++++------------ 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/dhnx/gistools/connect_points.py b/dhnx/gistools/connect_points.py index 0bd1c26b..1215e5c4 100644 --- a/dhnx/gistools/connect_points.py +++ b/dhnx/gistools/connect_points.py @@ -217,7 +217,7 @@ def create_object_connection(point_geom, lines, tol_distance=tol_distance): tol_distance=tol_distance) lines_drop.extend(nearest_line_idx) - if drop_neighbours and i+1 < n_conn: + if drop_neighbours and i + 1 < n_conn: # Find all neighbours of the nearest segment(s) and delete them # as well. Chances are that they do not provide an advantage # over the nearest segment. This allows the next connection @@ -233,7 +233,7 @@ def create_object_connection(point_geom, lines, tol_distance=tol_distance): # Create a GeoDataFrame with the connection lines of the current id gdf_conn_lines_id = gpd.GeoDataFrame( - data={'id_full': [id_full]*len(conn_lines_id)}, + data={'id_full': [id_full] * len(conn_lines_id)}, geometry=conn_lines_id, crs=lines.crs) # Drop duplicate geometries in connection lines @@ -259,16 +259,13 @@ def create_object_connection(point_geom, lines, tol_distance=tol_distance): # Remove the line that is to be replaced lines.drop([nearest_line_idx], inplace=True) # Combine remaining lines with two new replacement lines - lines = pd.concat( - [gpd.GeoDataFrame(lines, crs=lines.crs), - gpd.GeoDataFrame( - geometry=[ - LineString([line_start, conn_point]), - LineString([conn_point, line_end])], - crs=lines.crs, - data=pd.concat([attributes, attributes]), # keep data - ), - ], + lines = pd.concat([ + gpd.GeoDataFrame(lines, crs=lines.crs), + gpd.GeoDataFrame( + geometry=[LineString([line_start, conn_point]), + LineString([conn_point, line_end])], + crs=lines.crs, + data=pd.concat([attributes, attributes]))], ignore_index=True) conn_lines_list.append(gdf_conn_lines_id)