Skip to content

Commit

Permalink
fix intersection cleaning
Browse files Browse the repository at this point in the history
  • Loading branch information
@Mary-h86
Mary-h86 committed Apr 17, 2024
1 parent f77a674 commit 96cab4c
Showing 1 changed file with 21 additions and 67 deletions.
88 changes: 21 additions & 67 deletions src/tile2net/raster/pednet.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
import datetime
import pandas as pd
import os

os.environ['USE_PYGEOS'] = '0'
import geopandas as gpd

Expand Down Expand Up @@ -137,37 +138,6 @@ def make_longer(self, line, thr):
merged = shapely.ops.linemerge(new_l)
return merged

def update_sw_intersecting(self, cw_pol):
"""
Cleans up sidewalk lines extending over crosswalks polygons and cut them
Parameters
----------
cw_pol : :class:`gpd.GeoDataFrame`
:class:`GeoDataFrame` of crosswalks
"""
# finding sidewalk lines that extends to the crosswalks polygons and cut them
swl = self.sidewalk.copy()
swl = swl.reset_index(drop=True).explode().reset_index(drop=True)
# print(swl.head(10))
cwrect = [i.minimum_rotated_rectangle for i in cw_pol.geometry]
cwrectgdf = geo2geodf(cwrect)

cwin, swin = swl.geometry.values.sindex.query_bulk(cwrectgdf.geometry.values,
predicate="intersects")
indsect = list(zip(cwin, swin))
inds = []
updated_geom = []
for v, k in indsect:
if swl.iloc[k, 1].intersects(cwrectgdf.iloc[v, 0]):
inds.append(k)
updated_geom.append(
swl.iloc[k, swl.columns.get_loc(swl.geometry.name)].difference(
cwrectgdf.iloc[v, 0]))
swl.iloc[inds, swl.columns.get_loc(swl.geometry.name)] = updated_geom
smoothed = wrinkle_remover(swl, 1)
self.sidewalk = smoothed

def create_crosswalk(self):
"""
Create crosswalks centerlines from polygons
Expand All @@ -184,27 +154,24 @@ def create_crosswalk(self):
cw_explode = cw_union.explode().reset_index(drop=True)
cw_explode = cw_explode[cw_explode.geometry.notna()].reset_index(drop=True)
cw_explode_ = morpho_atts(cw_explode)
polak = []
for c, geom in enumerate(cw_explode_.geometry):
if geom.area < 5:
continue
# logging.info(c, 'continue')
else:
# if crosswalks are attached to each other and form a T or U
if cw_explode_.iloc[
c, cw_explode_.columns.get_loc(cw_explode_.convexity.name)] < 0.8:
if cw_explode_.iloc[c, cw_explode_.columns.get_loc(cw_explode_.convexity.name)] < 0.8:
av_width = 4 * geom.area / geom.length
geom_er = geom.buffer(-av_width / 4)
polak.append(geom_er)
if geom_er.type == "MultiPolygon":
for g in list(geom_er.geoms): # shapely 2
if g.area > 2:
cnl = to_cline(g, 0.3, 1)
tr_line_ = trim_checkempty(cnl, 4.5, 2)
if tr_line_.length < 6:
extended = self.make_longer(tr_line_, 3/4)
if tr_line_.length < 8:
extended = self.make_longer(tr_line_, 0.8)
extended_line = extend_lines(geo2geodf([extended]),
tolerance=5,
tolerance=8,
target=geo2geodf([geom.boundary]), extension=0)
for gi in extended_line.geometry:
cw_lin_geom.append(gi)
Expand All @@ -214,13 +181,12 @@ def create_crosswalk(self):
continue
elif geom_er.type == "Polygon":
if geom_er.area > 2:
# print(geom_er.type, geom_er.area)
cnl = to_cline(geom_er, 0.2, 1)
tr_line_ = trim_checkempty(cnl, 4.5, 2)
if tr_line_.length < 6:
extended = self.make_longer(tr_line_, 3 / 4)
if tr_line_.length < 8:
extended = self.make_longer(tr_line_, 0.8)
extended_line = extend_lines(geo2geodf([extended]),
target=geo2geodf([geom.boundary]), tolerance=5,
target=geo2geodf([geom.boundary]), tolerance=8,
extension=0)
for g in extended_line.geometry:
cw_lin_geom.append(g)
Expand All @@ -231,26 +197,20 @@ def create_crosswalk(self):
else:
continue
else:
polak.append(geom)
line = get_crosswalk_cnl(geom)
if line.length < 6:
extended = self.make_longer(line, 3 / 4)
if line.length < 8:
extended = self.make_longer(line, 0.8)
extended_line = extend_lines(geo2geodf([extended]),
target=geo2geodf([geom.boundary]), tolerance=5, extension=0)
target=geo2geodf([geom.boundary]), tolerance=8, extension=0)
for g in extended_line.geometry:
cw_lin_geom.append(g)
else:
cw_lin_geom.append(line)

cwpol = gpd.GeoDataFrame(geometry=polak)
cwpol.geometry = cwpol.geometry.set_crs(3857)

self.update_sw_intersecting(cwpol)

cw_ntw = geo2geodf(cw_lin_geom)
cw_ntw['f_type'] = 'crosswalk'
cw_ntw.geometry = cw_ntw.geometry.set_crs(3857)
smoothed = wrinkle_remover(cw_ntw, 1.2)
smoothed = wrinkle_remover(cw_ntw, 1.3)
self.crosswalk = smoothed

def create_lines(self, gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
Expand Down Expand Up @@ -330,20 +290,14 @@ def create_sidewalks(self):
None
"""
sw_all = self.prepare_class_gdf('sidewalk')
cwp = self.prepare_class_gdf('crosswalk')

if len(sw_all) > 0:
# sw_filtered = self.find_medianisland(sw_all, cwp)
# if not isinstance(sw_filtered, int):
# swntw = self.create_lines(sw_filtered) # union_df
# else:
# swntw = self.create_lines(sw_all) # swuinon_df
swntw = self.create_lines(sw_all)
logging.info('..... creating the processed sidewalk network')
#
swntw.geometry = swntw.simplify(0.6)
sw_modif_uni = gpd.GeoDataFrame(
geometry=gpd.GeoSeries([geom for geom in swntw.unary_union.geoms]))
geometry=gpd.GeoSeries([geom for geom in swntw.unary_union.geoms]))
sw_modif_uni_met = set_gdf_crs(sw_modif_uni, 3857)
sw_uni_lines = sw_modif_uni_met.explode()
sw_uni_lines.reset_index(drop=True, inplace=True)
Expand Down Expand Up @@ -380,7 +334,7 @@ def convert_whole_poly2line(self):

# query LineString geometry to identify points intersecting 2 geometries
inp, res = self.crosswalk.sindex.query(geo2geodf(points).geometry,
predicate="intersects")
predicate="intersects")
unique, counts = np.unique(inp, return_counts=True)
ends = np.unique(res[np.isin(inp, unique[counts == 1])])

Expand Down Expand Up @@ -413,7 +367,7 @@ def convert_whole_poly2line(self):
# create a dataframe of points
if len(new_geoms_s) > 0:
ps = [g[1] for g in new_geoms_s]
ls = [g[0] for g in new_geoms_s]
# ls = [g[0] for g in new_geoms_s]
pdfs = gpd.GeoDataFrame(geometry=ps)
pdfs.set_crs(3857, inplace=True)

Expand All @@ -422,7 +376,7 @@ def convert_whole_poly2line(self):

if len(new_geoms_e) > 0:
pe = [g[1] for g in new_geoms_e]
le = [g[0] for g in new_geoms_e]
# le = [g[0] for g in new_geoms_e]
pdfe = gpd.GeoDataFrame(geometry=pe)
pdfe.set_crs(3857, inplace=True)

Expand All @@ -431,14 +385,13 @@ def convert_whole_poly2line(self):

if len(new_geoms_both) > 0:
pb = [g[1] for g in new_geoms_both]
lb = [g[0] for g in new_geoms_both] # crosswalk lines where both ends do not intersect
# lb = [g[0] for g in new_geoms_both] # crosswalk lines where both ends do not intersect
pdfb = gpd.GeoDataFrame(geometry=pb)
pdfb.set_crs(3857, inplace=True)

connect_b = get_shortest(self.sidewalk, pdfb, f_type='sidewalk_connection')
all_connections.append(connect_b)
if len(all_connections) > 1:

connect = pd.concat(all_connections)
elif len(all_connections) == 1:
connect = all_connections[0]
Expand All @@ -456,8 +409,8 @@ def convert_whole_poly2line(self):

for k, v in indcwnear:
island_lines.append(
shapely.shortest_line(self.island.geometry.values[v],
pdfb.geometry.values[k]))
shapely.shortest_line(self.island.geometry.values[v],
pdfb.geometry.values[k]))

island = gpd.GeoDataFrame(geometry=island_lines)

Expand All @@ -474,12 +427,13 @@ def convert_whole_poly2line(self):
combined.geometry = combined.geometry.set_crs(3857)
combined.geometry = combined.geometry.to_crs(4326)
combined = combined[~combined.geometry.isna()]
combined.drop_duplicates(subset='geometry', inplace=True)
combined.reset_index(drop=True, inplace=True)
path = self.project.network.path

path.mkdir(parents=True, exist_ok=True)
path = path.joinpath(
f'{self.project.name}-Network-{datetime.datetime.now().strftime("%d-%m-%Y_%H")}'
f'{self.project.name}-Network-{datetime.datetime.now().strftime("%d-%m-%Y_%H")}'
)
combined.to_file(path)

Expand Down

0 comments on commit 96cab4c

Please sign in to comment.