import math import numpy as np import pandas as pd import shapely from geopandas import GeoDataFrame, GeoSeries from libpysal.graph import Graph from numpy.typing import NDArray from pandas import DataFrame, Series __all__ = [ "volume", "floor_area", "courtyard_area", "longest_axis_length", "perimeter_wall", "street_profile", "weighted_character", ] try: from numba import njit HAS_NUMBA = True except (ModuleNotFoundError, ImportError): HAS_NUMBA = False from libpysal.common import jit as njit def volume( area: NDArray[np.float64] | Series, height: NDArray[np.float64] | Series, ) -> NDArray[np.float64] | Series: """ Calculates volume of each object in given GeoDataFrame based on its height and area. .. math:: area * height Parameters ---------- area : NDArray[np.float64] | Series array of areas height : NDArray[np.float64] | Series array of heights Returns ------- NDArray[np.float64] | Series array of a type depending on the input Examples -------- >>> import pandas as pd >>> area = pd.Series([100, 30, 40, 75, 230]) >>> height = pd.Series([22, 6.5, 12, 9, 4.5]) >>> momepy.volume(area, height) 0 2200.0 1 195.0 2 480.0 3 675.0 4 1035.0 dtype: float64 """ return area * height def floor_area( area: NDArray[np.float64] | Series, height: NDArray[np.float64] | Series, floor_height: float | NDArray[np.float64] | Series = 3, ) -> NDArray[np.float64] | Series: """Calculates floor area of each object based on height and area. The number of floors is simplified into the formula: ``height // floor_height``. B y default one floor is approximated to 3 metres. .. math:: area * \\frac{height}{floor_height} Parameters ---------- area : NDArray[np.float64] | Series array of areas height : NDArray[np.float64] | Series array of heights floor_height : float | NDArray[np.float64] | Series, optional float denoting the uniform floor height or an aarray reflecting the building height by geometry, by default 3 Returns ------- NDArray[np.float64] | Series array of a type depending on the input Examples -------- >>> import pandas as pd >>> area = pd.Series([100, 30, 40, 75, 230]) >>> height = pd.Series([22, 6.5, 12, 9, 4.5]) >>> momepy.floor_area(area, height) 0 700.0 1 60.0 2 160.0 3 225.0 4 230.0 dtype: float64 If you know average height of floors per each building, you can pass it directly: >>> floor_height = pd.Series([3.2, 3, 4, 3, 4.5]) >>> momepy.floor_area(area, height, floor_height=floor_height) 0 600.0 1 60.0 2 120.0 3 225.0 4 230.0 dtype: float64 """ return area * (height // floor_height) def courtyard_area(geometry: GeoDataFrame | GeoSeries) -> Series: """Calculates area of holes within geometry - area of courtyards. Parameters ---------- geometry : GeoDataFrame | GeoSeries A GeoDataFrame or GeoSeries containing polygons to analyse. Returns ------- Series Examples -------- >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") >>> ca = momepy.courtyard_area(buildings) >>> ca 0 0.0 1 0.0 2 0.0 3 0.0 4 0.0 ... 139 0.0 140 0.0 141 0.0 142 0.0 143 0.0 Name: courtyard_area, Length: 144, dtype: float64 Verify that at least some buildings have courtyards: >>> ca.sum() np.float64(353.33274206543274) """ return Series( shapely.area( shapely.polygons(shapely.get_exterior_ring(geometry.geometry.array)) ) - geometry.area, index=geometry.index, name="courtyard_area", ) def longest_axis_length(geometry: GeoDataFrame | GeoSeries) -> Series: """Calculates the length of the longest axis of object. Axis is defined as a diameter of minimal bounding circle around the geometry. It does not have to be fully inside an object. .. math:: \\max \\left\\{d_{1}, d_{2}, \\ldots, d_{n}\\right\\} Parameters ---------- geometry : GeoDataFrame | GeoSeries A GeoDataFrame or GeoSeries containing polygons to analyse. Returns ------- Series Examples -------- >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") >>> momepy.longest_axis_length(buildings) 0 40.265562 1 191.254382 2 37.247151 3 47.022428 4 37.170142 ... 139 11.587272 140 27.747002 141 52.566435 142 11.091309 143 15.472821 Name: geometry, Length: 144, dtype: float64 """ return shapely.minimum_bounding_radius(geometry.geometry) * 2 def perimeter_wall( geometry: GeoDataFrame | GeoSeries, graph: Graph | None = None, buffer: float = 0.01 ) -> Series: """Calculate the perimeter wall length of the joined structure. Parameters ---------- geometry : GeoDataFrame | GeoSeries A GeoDataFrame or GeoSeries containing polygons to analyse. graph : Graph | None, optional Graph encoding Queen contiguity of ``geometry``. If ``None`` Queen contiguity is built on the fly. buffer: float Buffer value for the geometry. It can be used to account for topological problems. Returns ------- Series Examples -------- >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") >>> momepy.perimeter_wall(buildings) 0 137.186310 1 663.342296 2 663.342296 3 663.342296 4 663.342296 ... 139 42.839590 140 78.562927 141 147.342182 142 118.354123 143 342.909172 Name: perimeter_wall, Length: 144, dtype: float64 By default, ``momepy`` calculates a Queen contiguity graph to determine connected components. Alternatively, you can pass that yourself. This can be useful when the graph is already computed or when you need to use a different method due to topological issues. >>> from libpysal import graph >>> strict_contig = graph.Graph.build_contiguity( ... buildings, rook=False, strict=True, ... ) >>> momepy.perimeter_wall(buildings, graph=strict_contig) 0 137.186310 1 663.342296 2 663.342296 3 663.342296 4 663.342296 ... 139 42.839590 140 78.562927 141 147.342182 142 118.354123 143 342.909172 Name: perimeter_wall, Length: 144, dtype: float64 """ if graph is None: graph = Graph.build_contiguity(geometry, rook=False) isolates = graph.isolates # measure perimeter walls of connected components while ignoring isolates blocks = geometry.drop(isolates) component_perimeter = ( blocks[[blocks.geometry.name]] .set_geometry(blocks.buffer(buffer)) # type: ignore .dissolve(by=graph.component_labels.drop(isolates)) .exterior.length ) # combine components with isolates results = Series(np.nan, index=geometry.index, name="perimeter_wall") results.loc[isolates] = geometry.geometry[isolates].exterior.length results.loc[results.index.drop(isolates)] = component_perimeter.loc[ graph.component_labels.loc[results.index.drop(isolates)] ].values return results def weighted_character( y: NDArray[np.float64] | Series, area: NDArray[np.float64] | Series, graph: Graph ) -> Series: """Calculates the weighted character. Character weighted by the area of the objects within neighbors defined in ``graph``. Results are index based on ``graph``. .. math:: \\frac{\\sum_{i=1}^{n} {character_{i} * area_{i}}}{\\sum_{i=1}^{n} area_{i}} Adapted from :cite:`dibble2017`. Notes ----- The index of ``y`` and ``area`` must match the index along which the ``graph`` is built. Parameters ---------- y : NDArray[np.float64] | Series The character values to be weighted. area : NDArray[np.float64] | Series The area values to be used as weightss graph : libpysal.graph.Graph A spatial weights matrix for values and areas. Returns ------- Series A Series containing the resulting values. Examples -------- Area-weighted elongation within 5 nearest neighbors: >>> from libpysal import graph >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") >>> buildings.head() uID geometry 0 1 POLYGON ((1603599.221 6464369.816, 1603602.984... 1 2 POLYGON ((1603042.88 6464261.498, 1603038.961 ... 2 3 POLYGON ((1603044.65 6464178.035, 1603049.192 ... 3 4 POLYGON ((1603036.557 6464141.467, 1603036.969... 4 5 POLYGON ((1603082.387 6464142.022, 1603081.574... Measure elongation (or anything else): >>> elongation = momepy.elongation(buildings) >>> elongation.head() 0 0.908244 1 0.581318 2 0.726527 3 0.838840 4 0.727294 Name: elongation, dtype: float64 Define spatial graph: >>> knn5 = graph.Graph.build_knn(buildings.centroid, k=5) >>> knn5 Measure the area-weighted character: >>> momepy.weighted_character(elongation, buildings.area, knn5) focal 0 0.808190 1 0.817309 2 0.627589 3 0.794769 4 0.806403 ... 139 0.780744 140 0.875046 141 0.753670 142 0.440000 143 0.901127 Name: sum, Length: 144, dtype: float64 """ stats = graph.describe(y * area, statistics=["sum"])["sum"] agg_area = graph.describe(area, statistics=["sum"])["sum"] return stats / agg_area def street_profile( streets: GeoDataFrame, buildings: GeoDataFrame, distance: float = 10, tick_length: float = 50, height: None | Series = None, ) -> DataFrame: """Calculates the street profile characters. This functions returns a DataFrame with widths, standard deviation of width, openness, heights, standard deviation of height and ratio height/width. The algorithm generates perpendicular lines to the ``streets`` dataframe features every ``distance`` and measures values on intersections with features of ``buildings``. If no feature is reached within ``tick_length`` its value is set as width (being a theoretical maximum). Derived from :cite:`araldi2019`. Parameters ---------- streets : GeoDataFrame A GeoDataFrame containing streets to analyse. buildings : GeoDataFrame A GeoDataFrame containing buildings along the streets. Only Polygon geometries are currently supported. distance : int (default 10) The distance between perpendicular ticks. tick_length : int (default 50) The length of ticks. height: pd.Series (default None) The ``pd.Series`` where building height are stored. If set to ``None``, height and ratio height/width will not be calculated. Returns ------- DataFrame Examples -------- >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") >>> streets = geopandas.read_file(path, layer="streets") >>> streets.head() geometry 0 LINESTRING (1603585.64 6464428.774, 1603413.20... 1 LINESTRING (1603268.502 6464060.781, 1603296.8... 2 LINESTRING (1603607.303 6464181.853, 1603592.8... 3 LINESTRING (1603678.97 6464477.215, 1603675.68... 4 LINESTRING (1603537.194 6464558.112, 1603557.6... >>> result = momepy.street_profile(streets, buildings) >>> result.head() width openness width_deviation 0 47.905964 0.946429 0.020420 1 42.418068 0.615385 2.644521 2 32.131831 0.608696 2.864438 3 50.000000 1.000000 NaN 4 50.000000 1.000000 NaN If you know height of each building, you can pass that along to get back more information: >>> import numpy as np >>> import pandas as pd >>> rng = np.random.default_rng(seed=42) >>> height = pd.Series(rng.integers(low=9, high=30, size=len(buildings))) >>> result = momepy.street_profile(streets, buildings, height=height) >>> result.head() width openness width_deviation height height_deviation hw_ratio 0 47.905964 0.946429 0.020420 12.666667 4.618802 0.264407 1 42.418068 0.615385 2.644521 21.500000 6.467869 0.506859 2 32.131831 0.608696 2.864438 17.555556 4.901647 0.546360 3 50.000000 1.000000 NaN NaN NaN NaN 4 50.000000 1.000000 NaN NaN NaN NaN """ # filter relevant buildings and streets inp, res = shapely.STRtree(streets.geometry).query( buildings.geometry, predicate="dwithin", distance=tick_length // 2 ) buildings_near_streets = np.unique(inp) streets_near_buildings = np.unique(res) relevant_buildings = buildings.iloc[buildings_near_streets].reset_index(drop=True) relevant_streets = streets.iloc[streets_near_buildings].reset_index(drop=True) if height is not None: height = height.iloc[buildings_near_streets].reset_index(drop=True) # calculate street profile on a subset of the data partial_res = _street_profile( relevant_streets, relevant_buildings, distance=distance, tick_length=tick_length, height=height, ) # return full result with defaults final_res = pd.DataFrame(np.nan, index=streets.index, columns=partial_res.columns) final_res.iloc[streets_near_buildings[partial_res.index.values]] = ( partial_res.values ) ## streets with no buildings get the theoretical width and max openness final_res.loc[final_res["width"].isna(), "width"] = tick_length final_res.loc[final_res["openness"].isna(), "openness"] = 1 return final_res def _street_profile( streets: GeoDataFrame, buildings: GeoDataFrame, distance: float = 10, tick_length: float = 50, height: None | Series = None, ) -> DataFrame: """Helper function that actually calculates the street profile characters.""" ## generate points for every street at `distance` intervals segments = streets.segmentize(distance) coords, coord_indxs = shapely.get_coordinates(segments, return_index=True) starts = ~pd.Series(coord_indxs).duplicated(keep="first") ends = ~pd.Series(coord_indxs).duplicated(keep="last") end_markers = starts | ends ## generate tick streings njit_ticks = generate_ticks(coords, end_markers.values, tick_length) ticks = shapely.linestrings(njit_ticks.reshape(-1, 2, 2)) ## find the length of intersection of the nearest building for every tick inp, res = buildings.geometry.sindex.query(ticks, predicate="intersects") intersections = shapely.intersection(ticks[inp], buildings.geometry.array[res]) distances = shapely.distance(intersections, shapely.points(coords[inp // 2])) # streets which intersect buildings have 0 distance to them distances[np.isnan(distances)] = 0 min_distances = pd.Series(distances).groupby(inp).min() dists = np.full((len(ticks),), np.nan) dists[min_distances.index.values] = min_distances.values ## generate tick values and groupby street tick_coords = np.repeat(coord_indxs, 2) ## multiple agg to avoid custom apply left_ticks = ( pd.Series(dists[::2]) .groupby(tick_coords[::2]) .mean() .replace(np.nan, tick_length // 2) ) right_ticks = ( pd.Series(dists[1::2]) .groupby(tick_coords[1::2]) .mean() .replace(np.nan, tick_length // 2) ) w = left_ticks + right_ticks grouper = pd.Series(dists).groupby(tick_coords) openness_agg = grouper.agg(["size", "count"]) # proportion of NAs o = (openness_agg["size"] - openness_agg["count"]) / openness_agg["size"] # needs to be seperate to pass ddof wd = grouper.std(ddof=0) final_result = pd.DataFrame( np.nan, columns=["width", "openness", "width_deviation"], index=streets.index ) final_result["width"] = w final_result["openness"] = o final_result["width_deviation"] = wd ## if heights are available add heights stats to the result if height is not None: min_heights = height.loc[res].groupby(inp).min() tick_heights = np.full((len(ticks),), np.nan) tick_heights[min_heights.index.values] = min_heights.values heights_res = pd.Series(tick_heights).groupby(tick_coords).agg(["mean", "std"]) final_result["height"] = heights_res["mean"] final_result["height_deviation"] = heights_res["std"] final_result["hw_ratio"] = final_result["height"] / final_result[ "width" ].replace(0, np.nan) return final_result # angle between two points @njit def _get_angle_njit(x1, y1, x2, y2): """ pt1, pt2 : tuple """ x_diff = x2 - x1 y_diff = y2 - y1 return math.degrees(math.atan2(y_diff, x_diff)) # get the second end point of a tick # p1 = bearing + 90 @njit def _get_point_njit(x1, y1, bearing, dist): bearing = math.radians(bearing) x = x1 + dist * math.cos(bearing) y = y1 + dist * math.sin(bearing) return np.array((x, y)) @njit def generate_ticks(list_points, end_markers, tick_length): ticks = np.empty((len(list_points) * 2, 4), dtype=float) for i in range(len(list_points)): tick_pos = i * 2 end = end_markers[i] pt = list_points[i] if end: ticks[tick_pos, :] = np.array([pt[0], pt[1], pt[0], pt[1]]) ticks[tick_pos + 1, :] = np.array([pt[0], pt[1], pt[0], pt[1]]) else: next_pt = list_points[i + 1] njit_angle1 = _get_angle_njit(pt[0], pt[1], next_pt[0], next_pt[1]) njit_end_1 = _get_point_njit( pt[0], pt[1], njit_angle1 + 90, tick_length / 2 ) njit_angle2 = _get_angle_njit(njit_end_1[0], njit_end_1[1], pt[0], pt[1]) njit_end_2 = _get_point_njit( njit_end_1[0], njit_end_1[1], njit_angle2, tick_length ) ticks[tick_pos, :] = np.array([njit_end_1[0], njit_end_1[1], pt[0], pt[1]]) ticks[tick_pos + 1, :] = np.array( [njit_end_2[0], njit_end_2[1], pt[0], pt[1]] ) return ticks