Skip to content

Commit

Permalink
ENH: retain index of buildings and plots intersecting sightlines (#667)
Browse files Browse the repository at this point in the history 
* retain building index

* retain plot index

* retain building index on point level

* retain plot index in point data
  • Loading branch information
@martinfleis
martinfleis authored Nov 13, 2024
1 parent 2d9bc9f commit a315499
Show file tree
Hide file tree
Showing 2 changed files with 129 additions and 20 deletions.
101 changes: 91 additions & 10 deletions momepy/streetscape.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,9 @@ class Streetscape:
The function also provides the distribution of measures separately for the two sides
of the street (right and left, assigned arbitrarily). Indicators for the left side
are preceded by ``left_*`` and for the right side by ``right_*``.
are preceded by ``left_*`` and for the right side by ``right_*`` At the same time,
you can retrieve the index of buildings that intersected each sightline via
``*_seq_sb_index``.
Additional street-level measures are: ``street_length`` (Street Length),
``windingness`` (Windingness or Curvature of Streets, equal to 1 -
Expand Down Expand Up @@ -532,6 +534,9 @@ def _compute_sigthlines_indicators(self, street_row, optimize_on=True):
current_street_front_sb = []
current_street_back_sb = []

left_ids_all = []
right_ids_all = []

# [Expanded] each time a sight line or intersight line occured
left_seq_sightlines_end_points = []
right_seq_sightlines_end_points = []
Expand Down Expand Up @@ -563,6 +568,8 @@ def _compute_sigthlines_indicators(self, street_row, optimize_on=True):
current_street_back_sb,
left_seq_sightlines_end_points,
right_seq_sightlines_end_points,
left_ids_all,
right_ids_all,
], None

# ------- SIGHT LINES
Expand All @@ -587,6 +594,9 @@ def _compute_sigthlines_indicators(self, street_row, optimize_on=True):
left_sl_cr_total = 0
right_sl_cr_total = 0

left_ids = []
right_ids = []

# iterate throught each sightline links to the sigh point:
# LEFT(1-*),RIGHT(1-*),FRONT(1), BACK(1)
for row_s in group.itertuples(index=False):
Expand Down Expand Up @@ -672,6 +682,9 @@ def _compute_sigthlines_indicators(self, street_row, optimize_on=True):
current_street_left_seq_sb_categories.append(
match_sl_building_category
)
left_ids.append(match_sl_building_id)
else:
left_ids.append(np.nan)

elif sightline_side == self.SIGHTLINE_RIGHT:
right_sl_count += 1
Expand All @@ -688,12 +701,18 @@ def _compute_sigthlines_indicators(self, street_row, optimize_on=True):
current_street_right_seq_sb_categories.append(
match_sl_building_category
)
right_ids.append(match_sl_building_id)
else:
right_ids.append(np.nan)

elif sightline_side == self.SIGHTLINE_BACK:
back_sl_tan_sb = match_sl_distance
elif sightline_side == self.SIGHTLINE_FRONT:
front_sl_tan_sb = match_sl_distance

left_ids_all.append(left_ids)
right_ids_all.append(right_ids)

# LEFT
left_os_count = left_sl_count
left_os = left_sl_distance_total / left_os_count
Expand Down Expand Up @@ -765,6 +784,8 @@ def _compute_sigthlines_indicators(self, street_row, optimize_on=True):
current_street_back_sb,
left_seq_sightlines_end_points,
right_seq_sightlines_end_points,
left_ids_all,
right_ids_all,
], gdf_sightlines

def _compute_sightline_indicators_full(self):
Expand Down Expand Up @@ -803,6 +824,8 @@ def _compute_sightline_indicators_full(self):
"back_sb",
"left_seq_os_endpoints",
"right_seq_os_endpoints",
"left_ids",
"right_ids",
],
)
df = df.set_index("street_index")
Expand Down Expand Up @@ -843,6 +866,7 @@ def _compute_sigthlines_plot_indicators_one_side(
parcel_seq_sb_ids = []
parcel_seq_sb = []
parcel_seq_sb_depth = []
parcel_ids_all = []

n = len(sightline_points)
if n == 0:
Expand All @@ -852,6 +876,7 @@ def _compute_sigthlines_plot_indicators_one_side(
parcel_seq_sb_ids,
parcel_seq_sb,
parcel_seq_sb_depth,
parcel_ids_all,
]

idx_end_point = 0
Expand Down Expand Up @@ -882,6 +907,7 @@ def _compute_sigthlines_plot_indicators_one_side(
match_distance = _distances.min()
match_geom = gdf_items.geometry[match_id]

parcel_ids = []
# ---------------
# result in intersightline
if match_id is not None:
Expand All @@ -901,13 +927,23 @@ def _compute_sigthlines_plot_indicators_one_side(
):
raise Exception("Not allowed: intersection is not of type Line")
parcel_seq_sb_depth.append(isec.length)
parcel_ids.append(match_id)
else:
parcel_ids.append(np.nan)

# ------- iterate
idx_end_point += 1
parcel_ids_all.append(parcel_ids)

parcel_sb_count.append(n_sightlines_touching)

return [parcel_sb_count, parcel_seq_sb_ids, parcel_seq_sb, parcel_seq_sb_depth]
return [
parcel_sb_count,
parcel_seq_sb_ids,
parcel_seq_sb,
parcel_seq_sb_depth,
parcel_ids_all,
]

def compute_plots(
self, plots: gpd.GeoDataFrame, sightline_plot_depth_extension: float = 300
Expand Down Expand Up @@ -959,10 +995,12 @@ def compute_plots(
"left_parcel_seq_sb_ids",
"left_parcel_seq_sb",
"left_parcel_seq_sb_depth",
"left_parcel_ids",
"right_parcel_sb_count",
"right_parcel_seq_sb_ids",
"right_parcel_seq_sb",
"right_parcel_seq_sb_depth",
"right_parcel_ids",
],
)
df = df.set_index("street_index").join(self._sightline_indicators.street_length)
Expand All @@ -972,7 +1010,9 @@ def compute_plots(
def _aggregate_plots(self):
values = []

for street_uid, row in self._plot_indicators.iterrows():
for street_uid, row in self._plot_indicators.drop(
columns=["left_parcel_ids", "right_parcel_ids"]
).iterrows():
left_parcel_sb_count = row.left_parcel_sb_count
left_parcel_seq_sb_ids = row.left_parcel_seq_sb_ids
left_parcel_seq_sb = row.left_parcel_seq_sb
Expand Down Expand Up @@ -1005,6 +1045,8 @@ def _aggregate_plots(self):
np.nan,
np.nan,
np.nan,
{},
{},
]
)
continue
Expand Down Expand Up @@ -1139,6 +1181,7 @@ def _aggregate_plots(self):
]
+ wd_ratio_list
+ wp_ratio_list
+ [set(left_parcel_seq_sb_ids), set(right_parcel_seq_sb_ids)]
)

columns = [
Expand All @@ -1160,6 +1203,8 @@ def _aggregate_plots(self):
"plot_WP_ratio",
"left_plot_WP_ratio",
"right_plot_WP_ratio",
"left_plot_seq_sb_index",
"right_plot_seq_sb_index",
]

self._aggregate_plot_data = pd.DataFrame(values, columns=columns).set_index(
Expand Down Expand Up @@ -1847,6 +1892,8 @@ def street_level(self) -> gpd.GeoDataFrame:
ind_left_built_coverage,
ind_right_built_coverage,
ind_built_coverage,
set(left_seq_sb_ids),
set(right_seq_sb_ids),
]
)

Expand Down Expand Up @@ -1924,6 +1971,8 @@ def street_level(self) -> gpd.GeoDataFrame:
"left_built_coverage",
"right_built_coverage",
"built_coverage",
"left_seq_sb_index",
"right_seq_sb_index",
],
)
.set_index("street_index")
Expand Down Expand Up @@ -2059,10 +2108,22 @@ def point_level(self) -> gpd.GeoDataFrame:
"right_bc",
"front_sb",
"back_sb",
"left_ids",
"right_ids",
]
]

point_data = point_data.explode(point_data.columns.tolist())
for col in point_data.columns[1:]:
point_data["left_ids"] = point_data["left_ids"].apply(
lambda x: {c for c in x if not pd.isna(c)}
)
point_data["right_ids"] = point_data["right_ids"].apply(
lambda x: {c for c in x if not pd.isna(c)}
)
point_data = point_data.rename(
columns={"left_ids": "left_seq_sb_index", "right_ids": "right_seq_sb_index"}
)
for col in point_data.columns[1:-2]:
point_data[col] = pd.to_numeric(point_data[col])

inds = [
Expand All @@ -2081,6 +2142,8 @@ def point_level(self) -> gpd.GeoDataFrame:
left_parcel_seq_sb_depth = []
right_parcel_seq_sb = []
right_parcel_seq_sb_depth = []
left_ids = []
right_ids = []

# we occasionally have more sightlines per point, so we need to average
# values
Expand Down Expand Up @@ -2116,25 +2179,43 @@ def point_level(self) -> gpd.GeoDataFrame:
for x in np.split(row.right_parcel_seq_sb_depth, right_inds)
]
)
left_ids.append(row.left_parcel_ids)
right_ids.append(row.right_parcel_ids)

point_parcel_data = pd.DataFrame(
{
"left_plot_seq_sb": left_parcel_seq_sb,
"left_plot_seq_sb_depth": left_parcel_seq_sb_depth,
"right_plot_seq_sb": right_parcel_seq_sb,
"right_plot_seq_sb_depth": right_parcel_seq_sb_depth,
"left_plot_seq_sb_index": left_ids,
"right_plot_seq_sb_index": right_ids,
},
index=self._plot_indicators.index,
)
point_parcel_data = point_parcel_data.explode(
point_parcel_data.columns.tolist()
)
point_parcel_data[
point_parcel_data.columns.drop(
["left_plot_seq_sb_index", "right_plot_seq_sb_index"]
)
] = point_parcel_data[
point_parcel_data.columns.drop(
["left_plot_seq_sb_index", "right_plot_seq_sb_index"]
)
].astype(float)
point_data = pd.concat(
[
point_data,
point_parcel_data.explode(
point_parcel_data.columns.tolist()
).astype(float),
],
[point_data, point_parcel_data],
axis=1,
)

point_data["left_plot_seq_sb_index"] = point_data[
"left_plot_seq_sb_index"
].apply(lambda x: {c for c in x if not pd.isna(c)})
point_data["right_plot_seq_sb_index"] = point_data[
"right_plot_seq_sb_index"
].apply(lambda x: {c for c in x if not pd.isna(c)})
inds.extend(
[
"plot_seq_sb",
Expand Down
48 changes: 38 additions & 10 deletions momepy/tests/test_streetscape.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,8 @@ def test_minimal(self):
street_df = sc.street_level()
point_df = sc.point_level()

assert street_df.shape == (35, 75)
assert point_df.shape == (1277, 24)
assert street_df.shape == (35, 77)
assert point_df.shape == (1277, 26)

def test_no_dtm(self):
sc = momepy.Streetscape(self.streets, self.buildings)
Expand All @@ -38,8 +38,8 @@ def test_no_dtm(self):
street_df = sc.street_level()
point_df = sc.point_level()

assert street_df.shape == (35, 92)
assert point_df.shape == (1277, 30)
assert street_df.shape == (35, 96)
assert point_df.shape == (1277, 34)

def test_no_plots(self):
rioxarray = pytest.importorskip("rioxarray")
Expand All @@ -51,8 +51,8 @@ def test_no_plots(self):
street_df = sc.street_level()
point_df = sc.point_level()

assert street_df.shape == (35, 79)
assert point_df.shape == (1277, 24)
assert street_df.shape == (35, 81)
assert point_df.shape == (1277, 26)

def test_all_values(self):
rioxarray = pytest.importorskip("rioxarray")
Expand All @@ -67,11 +67,21 @@ def test_all_values(self):
street_df = sc.street_level()
point_df = sc.point_level()

assert street_df.shape == (35, 102)
assert point_df.shape == (1277, 30)
assert street_df.shape == (35, 106)
assert point_df.shape == (1277, 34)

np.testing.assert_allclose(
street_df.drop(columns="geometry").median().to_numpy(),
street_df.drop(
columns=[
"geometry",
"left_seq_sb_index",
"right_seq_sb_index",
"left_plot_seq_sb_index",
"right_plot_seq_sb_index",
]
)
.median()
.to_numpy(),
[
40.0,
1.0,
Expand Down Expand Up @@ -178,7 +188,17 @@ def test_all_values(self):
)

np.testing.assert_allclose(
point_df.drop(columns="geometry").median().to_numpy(),
point_df.drop(
columns=[
"geometry",
"left_seq_sb_index",
"right_seq_sb_index",
"left_plot_seq_sb_index",
"right_plot_seq_sb_index",
]
)
.median()
.to_numpy(),
[
1.0,
50.0,
Expand Down Expand Up @@ -284,6 +304,8 @@ def test_all_values(self):
"left_built_coverage",
"right_built_coverage",
"built_coverage",
"left_seq_sb_index",
"right_seq_sb_index",
"nodes_degree_1",
"nodes_degree_4",
"nodes_degree_3_5_plus",
Expand Down Expand Up @@ -312,6 +334,8 @@ def test_all_values(self):
"plot_WP_ratio",
"left_plot_WP_ratio",
"right_plot_WP_ratio",
"left_plot_seq_sb_index",
"right_plot_seq_sb_index",
"slope_degree",
"slope_percent",
"n_slopes",
Expand Down Expand Up @@ -340,10 +364,14 @@ def test_all_values(self):
"right_bc",
"front_sb",
"back_sb",
"left_seq_sb_index",
"right_seq_sb_index",
"left_plot_seq_sb",
"left_plot_seq_sb_depth",
"right_plot_seq_sb",
"right_plot_seq_sb_depth",
"left_plot_seq_sb_index",
"right_plot_seq_sb_index",
"os_count",
"os",
"sb_count",
Expand Down

0 comments on commit a315499

Please sign in to comment.