Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature: Allow multiple potential building connections #147

Merged
merged 8 commits into from
Feb 17, 2025
399 changes: 220 additions & 179 deletions dhnx/gistools/connect_points.py
Original file line number Diff line number Diff line change
@@ -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,112 @@ 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)

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)],
ignore_index=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:
# 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)],
# 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([line_start, conn_point]),
LineString([conn_point, line_end])],
crs=lines.crs,
data=pd.concat([attributes, attributes]))],
ignore_index=True)

logger.info('Connection of buildings completed.')
conn_lines_list.append(gdf_conn_lines_id)

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):
@@ -322,11 +326,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 +357,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 +397,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 +410,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 +426,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 +490,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, n_conn=1, n_conn_prod=1,
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 +525,20 @@ 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
the total number of line elements. Default is True.
Returns
-------
@@ -509,14 +547,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 +570,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'
@@ -548,26 +590,28 @@ 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
lines_producers.index = producers.index

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 +628,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 +637,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)
132 changes: 111 additions & 21 deletions dhnx/gistools/geometry_operations.py
Original file line number Diff line number Diff line change
@@ -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
@@ -102,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

@@ -183,7 +207,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))
@@ -278,7 +302,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
@@ -365,11 +389,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)
@@ -384,12 +423,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
@@ -426,7 +487,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)):
@@ -490,7 +551,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.
@@ -510,16 +571,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
31 changes: 15 additions & 16 deletions dhnx/optimization/optimization_models.py
Original file line number Diff line number Diff line change
@@ -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))
@@ -438,18 +437,18 @@ 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']
except (KeyError, IndexError):
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 +468,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
@@ -513,7 +512,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 +521,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)
)
@@ -667,7 +666,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):
@@ -766,7 +765,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()