From 4f56cf9cf6438e9dcaf4658175f4f6a4a73c75ad Mon Sep 17 00:00:00 2001 From: Krasen Samardzhiev Date: Fri, 3 May 2024 16:23:44 +0200 Subject: [PATCH 1/6] functional area_ratio --- momepy/functional/_intensity.py | 47 ++++++++- momepy/functional/tests/test_intensity.py | 122 ++++++++++++++++++++++ 2 files changed, 168 insertions(+), 1 deletion(-) diff --git a/momepy/functional/_intensity.py b/momepy/functional/_intensity.py index 0b4fe6bb..6e1cc0a1 100644 --- a/momepy/functional/_intensity.py +++ b/momepy/functional/_intensity.py @@ -1,9 +1,11 @@ +import numpy as np +import pandas as pd import shapely from geopandas import GeoDataFrame, GeoSeries from libpysal.graph import Graph from pandas import Series -__all__ = ["courtyards"] +__all__ = ["courtyards", "area_ratio"] def courtyards(geometry: GeoDataFrame | GeoSeries, graph: Graph) -> Series: @@ -48,3 +50,46 @@ def _calculate_courtyards(group): ) return result + + +def area_ratio( + left: Series, right: Series, right_group_key: Series | np.ndarray +) -> pd.Series: + """ + Calculate covered area ratio or floor area ratio of objects. + .. math:: + \\textit{covering object area} \\over \\textit{covered object area} + + Adapted from :cite:`schirmer2015`. + + Parameters + ---------- + left : Series + A GeoDataFrame with the areas of the objects being covered (e.g. land unit). + right : Series + A GeoDataFrame with the areas of the covering objects (e.g. building). + right_group_key: np.array | pd.Series + The group key that assigns objects from ``right`` to ``left``. + + Returns + ------- + Series + + + Examples + -------- + >>> tessellation_df['CAR'] = mm.area_ratio(tessellation_df['area'], + ... buildings_df['area'], + ... buildings_df['uID']) + """ + + if isinstance(right_group_key, np.ndarray): + right_group_key = pd.Series(right_group_key, index=right.index) + + results = pd.Series(np.nan, left.index) + stats = ( + right.loc[right_group_key.index.values].groupby(right_group_key.values).sum() + ) + results.loc[stats.index.values] = stats.values + + return results / left.values diff --git a/momepy/functional/tests/test_intensity.py b/momepy/functional/tests/test_intensity.py index 18e25122..f36a7284 100644 --- a/momepy/functional/tests/test_intensity.py +++ b/momepy/functional/tests/test_intensity.py @@ -1,5 +1,6 @@ import geopandas as gpd import numpy as np +import pandas as pd from libpysal.graph import Graph from pandas.testing import assert_series_equal @@ -38,6 +39,65 @@ def test_courtyards(self): expected = {"mean": 0.6805555555555556, "sum": 98, "min": 0, "max": 1} assert_result(courtyards, expected, self.df_buildings) + def test_area_ratio(self): + car_block = mm.area_ratio( + self.blocks.geometry.area, + self.df_buildings["area"], + self.df_buildings["bID"], + ) + car_block_expected = { + "mean": 0.27619743196980123, + "max": 0.35699143584461146, + "min": 0.12975039475826336, + "count": 8, + } + assert_result(car_block, car_block_expected, self.blocks) + + car = mm.area_ratio( + self.df_tessellation.geometry.area, + self.df_buildings["area"], + self.df_buildings["uID"] - 1, + ) + car2 = mm.area_ratio( + self.df_tessellation.set_index("uID").area, + self.df_buildings.set_index("uID").area, + self.df_buildings["uID"].values, + ) + car_expected = { + "mean": 0.3206556897709747, + "max": 0.8754071653707558, + "min": 0.029097983413141276, + "count": 144, + } + assert_result(car, car_expected, self.df_tessellation) + assert_result(car2, car_expected, self.df_tessellation.set_index("uID")) + + car_sel = mm.area_ratio( + self.df_tessellation.iloc[10:20]["area"], + self.df_buildings["area"], + self.df_tessellation.iloc[10:20]["uID"] - 1, + ) + car_sel_expected = { + "mean": 0.3892868062654601, + "max": 0.5428192449477212, + "min": 0.22057633949526625, + "count": 10, + } + assert_result(car_sel, car_sel_expected, self.df_tessellation.iloc[10:20]) + + far = mm.area_ratio( + self.df_tessellation.geometry.area, + self.df_buildings["fl_area"], + self.df_buildings["uID"] - 1, + ) + far_expected = { + "mean": 1.910949846262234, + "max": 7.003257322966046, + "min": 0.26188185071827147, + "count": 144, + } + assert_result(far, far_expected, self.df_tessellation) + class TestIntensityEquality: def setup_method(self): @@ -70,3 +130,65 @@ def test_courtyards(self): assert_series_equal( new_courtyards, old_courtyards, check_names=False, check_dtype=False ) + + def test_area_ratio(self): + self.blocks["area"] = self.blocks.geometry.area + car_block_new = mm.area_ratio( + self.blocks.geometry.area, + self.df_buildings["area"], + self.df_buildings["bID"], + ) + car_block_old = mm.AreaRatio( + self.blocks, self.df_buildings, "area", "area", "bID" + ).series + assert_series_equal( + car_block_new, car_block_old, check_dtype=False, check_names=False + ) + + car_new = mm.area_ratio( + self.df_tessellation.geometry.area, + self.df_buildings["area"], + self.df_buildings["uID"] - 1, + ) + car2_new = mm.area_ratio( + self.df_tessellation.set_index("uID").area, + self.df_buildings.set_index("uID").area, + self.df_buildings["uID"].values, + ) + car_old = mm.AreaRatio( + self.df_tessellation, self.df_buildings, "area", "area", "uID" + ).series + assert_series_equal(car_new, car_old, check_dtype=False, check_names=False) + assert_series_equal( + pd.Series(car_old.values, car2_new.index), + car2_new, + check_dtype=False, + check_names=False, + ) + + car_sel = mm.AreaRatio( + self.df_tessellation.iloc[10:20], self.df_buildings, "area", "area", "uID" + ).series + car_sel_new = mm.area_ratio( + self.df_tessellation.iloc[10:20]["area"], + self.df_buildings["area"], + self.df_tessellation.iloc[10:20]["uID"] - 1, + ) + + assert_series_equal(car_sel_new, car_sel, check_dtype=False, check_names=False) + + far_new = mm.area_ratio( + self.df_tessellation.geometry.area, + self.df_buildings["fl_area"], + self.df_buildings["uID"] - 1, + ) + + far_old = mm.AreaRatio( + self.df_tessellation, + self.df_buildings, + self.df_tessellation.area, + self.df_buildings.fl_area, + "uID", + ).series + + assert_series_equal(far_new, far_old, check_dtype=False, check_names=False) From b0065ea3c55a1b36593119d0ac6e9bbb4a60814f Mon Sep 17 00:00:00 2001 From: Krasen Samardzhiev Date: Fri, 31 May 2024 11:30:30 +0200 Subject: [PATCH 2/6] Update momepy/functional/_intensity.py Co-authored-by: Martin Fleischmann --- momepy/functional/_intensity.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/momepy/functional/_intensity.py b/momepy/functional/_intensity.py index 97be0a88..58b2be59 100644 --- a/momepy/functional/_intensity.py +++ b/momepy/functional/_intensity.py @@ -116,7 +116,7 @@ def _calc_nodedensity(group, edges): def area_ratio( - left: Series, right: Series, right_group_key: Series | np.ndarray + covered: Series, covering: Series, aggregation_key: Series | np.ndarray ) -> pd.Series: """Calculate covered area ratio or floor area ratio of objects. .. math:: From b530ade08ca7845cf731cdd89b28e4b21112a8f8 Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Mon, 20 May 2024 13:10:57 +0200 Subject: [PATCH 3/6] BUG/ENH: support custom enclosure index is in tessellation and GeoDataFrame as enclosure input (#582) * BUG: ensure custom enclosure index is supported in tessellation * support both series and datarame as enclosures input * tests * typing --- momepy/functional/_elements.py | 17 +++++++++++------ momepy/functional/tests/test_elements.py | 17 +++++++++++++++++ 2 files changed, 28 insertions(+), 6 deletions(-) diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py index 7fe1c411..9c466410 100644 --- a/momepy/functional/_elements.py +++ b/momepy/functional/_elements.py @@ -91,7 +91,7 @@ def morphological_tessellation( def enclosed_tessellation( geometry: GeoSeries | GeoDataFrame, - enclosures: GeoSeries, + enclosures: GeoSeries | GeoDataFrame, shrink: float = 0.4, segment: float = 0.5, threshold: float = 0.05, @@ -129,7 +129,7 @@ def enclosed_tessellation( ---------- geometry : GeoSeries | GeoDataFrame A GeoDataFrame or GeoSeries containing buildings to tessellate the space around. - enclosures : GeoSeries + enclosures : GeoSeries | GeoDataFrame The enclosures geometry, which can be generated using :func:`momepy.enclosures`. shrink : float, optional The distance for negative buffer to generate space between adjacent polygons). @@ -161,7 +161,7 @@ def enclosed_tessellation( """ # convert to GeoDataFrame and add position (we will need it later) - enclosures = enclosures.to_frame() + enclosures = enclosures.geometry.to_frame() enclosures["position"] = range(len(enclosures)) # preserve index name if exists @@ -206,17 +206,22 @@ def enclosed_tessellation( unchanged_in_new = new_df.loc[[-1]] new_df = new_df.drop(-1) clean_blocks = pd.concat( - [enclosures.drop(altered).drop(columns="position"), unchanged_in_new] + [ + enclosures.drop(enclosures.index[altered]).drop(columns="position"), + unchanged_in_new, + ] ) else: - clean_blocks = enclosures.drop(altered).drop(columns="position") + clean_blocks = enclosures.drop(enclosures.index[altered]).drop( + columns="position" + ) # assign negative index to enclosures with no buildings clean_blocks.index = range(-len(clean_blocks), 0, 1) # get building index for enclosures with single building singles = enclosures.iloc[single] - singles.index = singles.position.loc[single].apply( + singles.index = singles.position.loc[singles.index].apply( lambda ix: geometry.iloc[res[inp == ix]].index[0] ) # combine results diff --git a/momepy/functional/tests/test_elements.py b/momepy/functional/tests/test_elements.py index 4480f1a0..e4f11b77 100644 --- a/momepy/functional/tests/test_elements.py +++ b/momepy/functional/tests/test_elements.py @@ -142,6 +142,23 @@ def test_enclosed_tessellation(self): threshold_elimination.sort_values("geometry").reset_index(drop=True), ) + tessellation_df = mm.enclosed_tessellation( + self.df_buildings, + self.enclosures, + ) + assert_geodataframe_equal(tessellation, tessellation_df) + + custom_index = self.enclosures + custom_index.index = (custom_index.index + 100).astype(str) + tessellation_custom_index = mm.enclosed_tessellation( + self.df_buildings, + custom_index, + ) + assert (tessellation_custom_index.geom_type == "Polygon").all() + assert tessellation_custom_index.crs == self.df_buildings.crs + assert (self.df_buildings.index.isin(tessellation_custom_index.index)).all() + assert tessellation_custom_index.enclosure_index.isin(custom_index.index).all() + def test_verify_tessellation(self): df = self.df_buildings b = df.total_bounds From 12b85e35147826f93671e8807ef9ce69baa60a17 Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Mon, 20 May 2024 20:02:15 +0200 Subject: [PATCH 4/6] ENH: adaptive buffer as a tessellation limit (#590) * ENH: adaptive buffer as a tessellation limit * pass kwargs * version checks and tests * fix expected --- momepy/elements.py | 30 ----------- momepy/functional/_elements.py | 67 ++++++++++++++++++++++++ momepy/functional/tests/test_elements.py | 36 +++++++++++++ 3 files changed, 103 insertions(+), 30 deletions(-) diff --git a/momepy/elements.py b/momepy/elements.py index 9698e6e5..798050ce 100644 --- a/momepy/elements.py +++ b/momepy/elements.py @@ -16,7 +16,6 @@ from tqdm.auto import tqdm __all__ = [ - "buffered_limit", "Tessellation", "Blocks", "get_network_id", @@ -29,35 +28,6 @@ GPD_GE_10 = Version(gpd.__version__) >= Version("1.0dev") -def buffered_limit(gdf, buffer=100): - """ - Define limit for :class:`momepy.Tessellation` as a buffer around buildings. - - See :cite:`fleischmann2020` for details. - - Parameters - ---------- - gdf : GeoDataFrame - A GeoDataFrame containing building footprints. - buffer : float - A buffer around buildings limiting the extend of tessellation. - - Returns - ------- - MultiPolygon - A MultiPolygon or Polygon defining the study area. - - Examples - -------- - >>> limit = mm.buffered_limit(buildings_df) - >>> type(limit) - shapely.geometry.polygon.Polygon - """ - return ( - gdf.buffer(buffer).union_all() if GPD_GE_10 else gdf.buffer(buffer).unary_union - ) - - class Tessellation: """ Generates tessellation. Three versions of tessellation can be created: diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py index 9c466410..97cfe585 100644 --- a/momepy/functional/_elements.py +++ b/momepy/functional/_elements.py @@ -1,22 +1,26 @@ import warnings import geopandas as gpd +import libpysal import numpy as np import pandas as pd import shapely from geopandas import GeoDataFrame, GeoSeries from joblib import Parallel, delayed from libpysal.cg import voronoi_frames +from libpysal.graph import Graph from packaging.version import Version GPD_GE_013 = Version(gpd.__version__) >= Version("0.13.0") GPD_GE_10 = Version(gpd.__version__) >= Version("1.0dev") +LPS_GE_411 = Version(libpysal.__version__) >= Version("4.11.dev") __all__ = [ "morphological_tessellation", "enclosed_tessellation", "verify_tessellation", "get_nearest_street", + "buffered_limit", ] @@ -348,3 +352,66 @@ def get_nearest_street( ids[blg_idx] = streets.index[str_idx] return ids + + +def buffered_limit( + gdf: GeoDataFrame | GeoSeries, + buffer: float | str = 100, + min_buffer: float = 0, + max_buffer: float = 100, + **kwargs, +) -> shapely.Geometry: + """ + Define limit for tessellation as a buffer around buildings. + + The function calculates a buffer around buildings and returns a MultiPolygon or + Polygon defining the study area. The buffer can be either a fixed number or + "adaptive" which calculates the buffer based on Gabriel graph. + + See :cite:`fleischmann2020` for details. + + Parameters + ---------- + gdf : GeoDataFrame | GeoSeries + A GeoDataFrame containing building footprints. + buffer : float | str, optional + A buffer around buildings limiting the extend of tessellation. If "adaptive", + the buffer is calculated based on Gabriel graph as the half of the maximum + distance between neighbors (represented as centroids) of each node + 10% of + such the maximum distance. The lower and upper bounds can be furhter specified + by ``min_buffer`` and ``max_buffer``. By default 100. + min_buffer : float, optional + The minimum adaptive buffer distance. By default 0. + max_buffer : float, optional + The maximum adaptive buffer distance. By default 100. + **kwargs + Keyword arguments passed to :meth:`geopandas.GeoSeries.buffer`. + + Returns + ------- + MultiPolygon + A MultiPolygon or Polygon defining the study area. + + Examples + -------- + >>> limit = mm.buffered_limit(buildings_df) + >>> type(limit) + shapely.geometry.polygon.Polygon + """ + if buffer == "adaptive": + if not LPS_GE_411: + raise ImportError( + "Adaptive buffer requires libpysal 4.11 or higher." + ) # because https://github.com/pysal/libpysal/pull/709 + gabriel = Graph.build_triangulation(gdf.centroid, "gabriel", kernel="identity") + max_dist = gabriel.aggregate("max") + buffer = np.clip(max_dist / 2 + max_dist * 0.1, min_buffer, max_buffer).values + + elif not isinstance(buffer, int | float): + raise ValueError("`buffer` must be either 'adaptive' or a number.") + + return ( + gdf.buffer(buffer, **kwargs).union_all() + if GPD_GE_10 + else gdf.buffer(buffer, **kwargs).unary_union + ) diff --git a/momepy/functional/tests/test_elements.py b/momepy/functional/tests/test_elements.py index e4f11b77..a5440be0 100644 --- a/momepy/functional/tests/test_elements.py +++ b/momepy/functional/tests/test_elements.py @@ -1,4 +1,5 @@ import geopandas as gpd +import libpysal import numpy as np import pandas as pd import pytest @@ -11,6 +12,7 @@ import momepy as mm GPD_GE_013 = Version(gpd.__version__) >= Version("0.13.0") +LPS_GE_411 = Version(libpysal.__version__) >= Version("4.11.dev") class TestElements: @@ -220,3 +222,37 @@ def test_get_nearest_street(self): streets.index = streets.index.astype(str) nearest = mm.get_nearest_street(self.df_buildings, streets, 10) assert (nearest == None).sum() == 137 # noqa: E711 + + def test_buffered_limit(self): + limit = mm.buffered_limit(self.df_buildings, 50) + assert limit.geom_type == "Polygon" + assert pytest.approx(limit.area) == 366525.967849688 + + @pytest.mark.skipif(not LPS_GE_411, reason="libpysal>=4.11 required") + def test_buffered_limit_adaptive(self): + limit = mm.buffered_limit(self.df_buildings, "adaptive") + assert limit.geom_type == "Polygon" + assert pytest.approx(limit.area) == 355819.18954170 + + limit = mm.buffered_limit(self.df_buildings, "adaptive", max_buffer=30) + assert limit.geom_type == "Polygon" + assert pytest.approx(limit.area) == 304200.301833294 + + limit = mm.buffered_limit( + self.df_buildings, "adaptive", min_buffer=30, max_buffer=300 + ) + assert limit.geom_type == "Polygon" + assert pytest.approx(limit.area) == 357671.831894244 + + @pytest.mark.skipif(LPS_GE_411, reason="libpysal>=4.11 required") + def test_buffered_limit_adaptive_error(self): + with pytest.raises( + ImportError, match="Adaptive buffer requires libpysal 4.11 or higher." + ): + mm.buffered_limit(self.df_buildings, "adaptive") + + def test_buffered_limit_error(self): + with pytest.raises( + ValueError, match="`buffer` must be either 'adaptive' or a number." + ): + mm.buffered_limit(self.df_buildings, "invalid") From 9da632211c15be12bb0e8fc140be9144af85c151 Mon Sep 17 00:00:00 2001 From: Krasen Samardzhiev Date: Fri, 31 May 2024 11:59:16 +0200 Subject: [PATCH 5/6] pr changes --- momepy/functional/_intensity.py | 21 ++++++++++++--------- momepy/functional/tests/test_intensity.py | 5 +++++ 2 files changed, 17 insertions(+), 9 deletions(-) diff --git a/momepy/functional/_intensity.py b/momepy/functional/_intensity.py index 58b2be59..db0943fe 100644 --- a/momepy/functional/_intensity.py +++ b/momepy/functional/_intensity.py @@ -116,7 +116,7 @@ def _calc_nodedensity(group, edges): def area_ratio( - covered: Series, covering: Series, aggregation_key: Series | np.ndarray + area: Series, overlay: Series, aggregation_key: Series | np.ndarray = None ) -> pd.Series: """Calculate covered area ratio or floor area ratio of objects. .. math:: @@ -126,12 +126,12 @@ def area_ratio( Parameters ---------- - left : Series + area : Series A GeoDataFrame with the areas of the objects being covered (e.g. land unit). - right : Series + overlay : Series A GeoDataFrame with the areas of the covering objects (e.g. building). - right_group_key: np.array | pd.Series - The group key that assigns objects from ``right`` to ``left``. + aggregation_key: np.array | pd.Series + The group key that assigns objects from ``overlay`` to ``area``. Returns ------- @@ -145,8 +145,11 @@ def area_ratio( ... buildings_df['uID']) """ - # if isinstance(right_group_key, np.ndarray): - # right_group_key = pd.Series(right_group_key, index=right.index) + if aggregation_key is None: + return overlay / area - results = describe_reached(right, right_group_key, result_index=left.index) - return results["sum"] / left.values + if isinstance(aggregation_key, np.ndarray): + aggregation_key = pd.Series(aggregation_key, index=overlay.index) + + results = describe_reached(overlay, aggregation_key, result_index=area.index) + return results["sum"] / area diff --git a/momepy/functional/tests/test_intensity.py b/momepy/functional/tests/test_intensity.py index cc470f9e..b8fb4b9d 100644 --- a/momepy/functional/tests/test_intensity.py +++ b/momepy/functional/tests/test_intensity.py @@ -115,6 +115,10 @@ def test_area_ratio(self): self.df_buildings.set_index("uID").area, self.df_tessellation.set_index("uID").index, ) + car3 = mm.area_ratio( + self.df_tessellation.geometry.area, + self.df_buildings.geometry.area, + ) car_expected = { "mean": 0.3206556897709747, "max": 0.8754071653707558, @@ -123,6 +127,7 @@ def test_area_ratio(self): } assert_result(car, car_expected, self.df_tessellation) assert_result(car2, car_expected, self.df_tessellation.set_index("uID")) + assert_result(car3, car_expected, self.df_tessellation) car_sel = mm.area_ratio( self.df_tessellation.iloc[10:20]["area"], From 54b6165286d3fcc776016b4c8abf2d74f8f7a00d Mon Sep 17 00:00:00 2001 From: Krasen Samardzhiev Date: Fri, 31 May 2024 15:12:27 +0200 Subject: [PATCH 6/6] test versioning --- momepy/functional/tests/test_intensity.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/momepy/functional/tests/test_intensity.py b/momepy/functional/tests/test_intensity.py index b8fb4b9d..b12fe738 100644 --- a/momepy/functional/tests/test_intensity.py +++ b/momepy/functional/tests/test_intensity.py @@ -1,13 +1,17 @@ import geopandas as gpd import numpy as np +import pandas as pd import pytest from libpysal.graph import Graph +from packaging.version import Version from pandas.testing import assert_series_equal import momepy as mm from .conftest import assert_result +PD_210 = Version(pd.__version__) >= Version("2.1.0") + class TestIntensity: def setup_method(self): @@ -90,6 +94,9 @@ def test_node_density(self): ): mm.node_density(edges, edges, g, weighted=True) + @pytest.mark.skipif( + not PD_210, reason="aggregation is different in previous pandas versions" + ) def test_area_ratio(self): car_block = mm.area_ratio( self.blocks.geometry.area, @@ -188,6 +195,9 @@ def test_courtyards(self): new_courtyards, old_courtyards, check_names=False, check_dtype=False ) + @pytest.mark.skipif( + not PD_210, reason="aggregation is different in previous pandas versions" + ) def test_area_ratio(self): self.blocks["area"] = self.blocks.geometry.area car_block_new = mm.area_ratio(