From a315499e995da942d14cadb7c1bac1bf9a6ddc66 Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Wed, 13 Nov 2024 16:22:42 +0100 Subject: [PATCH] ENH: retain index of buildings and plots intersecting sightlines (#667) * retain building index * retain plot index * retain building index on point level * retain plot index in point data --- momepy/streetscape.py | 101 ++++++++++++++++++++++++++++--- momepy/tests/test_streetscape.py | 48 ++++++++++++--- 2 files changed, 129 insertions(+), 20 deletions(-) diff --git a/momepy/streetscape.py b/momepy/streetscape.py index bed68af0..e4d9f6a1 100644 --- a/momepy/streetscape.py +++ b/momepy/streetscape.py @@ -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 - @@ -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 = [] @@ -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 @@ -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): @@ -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 @@ -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 @@ -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): @@ -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") @@ -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: @@ -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 @@ -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: @@ -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 @@ -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) @@ -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 @@ -1005,6 +1045,8 @@ def _aggregate_plots(self): np.nan, np.nan, np.nan, + {}, + {}, ] ) continue @@ -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 = [ @@ -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( @@ -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), ] ) @@ -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") @@ -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 = [ @@ -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 @@ -2116,6 +2179,8 @@ 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( { @@ -2123,18 +2188,34 @@ def point_level(self) -> gpd.GeoDataFrame: "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", diff --git a/momepy/tests/test_streetscape.py b/momepy/tests/test_streetscape.py index bc4857f4..86cdee49 100644 --- a/momepy/tests/test_streetscape.py +++ b/momepy/tests/test_streetscape.py @@ -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) @@ -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") @@ -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") @@ -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, @@ -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, @@ -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", @@ -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", @@ -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",