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

ENH: retain index of buildings and plots intersecting sightlines #667

Merged
merged 4 commits into from
Nov 13, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
101 changes: 91 additions & 10 deletions momepy/streetscape.py
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
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
Loading