import geopandas as gpd import numpy as np import shapely from geopandas import GeoDataFrame, GeoSeries from numpy.typing import NDArray from packaging.version import Version from pandas import DataFrame, MultiIndex, Series from momepy.functional import _dimension __all__ = [ "form_factor", "fractal_dimension", "facade_ratio", "circular_compactness", "square_compactness", "convexity", "courtyard_index", "rectangularity", "shape_index", "corners", "squareness", "equivalent_rectangular_index", "elongation", "centroid_corner_distance", "linearity", "compactness_weighted_axis", ] GPD_013 = Version(gpd.__version__) >= Version("0.13") def form_factor( geometry: GeoDataFrame | GeoSeries, height: NDArray[np.float64] | Series, ) -> Series: """Calculates the form factor of each object given its geometry and height. .. math:: surface \\over {volume^{2 \\over 3}} where .. math:: surface = (perimeter * height) + area Adapted from :cite:`bourdic2012`. Parameters ---------- geometry : GeoDataFrame | GeoSeries A GeoDataFrame or GeoSeries containing polygons to analyse. height : NDArray[np.float64] | Series array of heights Returns ------- Series Examples -------- >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") Synthesize some building height information. >>> import numpy as np >>> rng = np.random.default_rng(seed=42) >>> height = rng.integers(low=9, high=30, size=len(buildings)) >>> momepy.form_factor(buildings, height) 0 5.588952 1 8.403204 2 5.149302 3 5.381587 4 5.030861 ... 139 6.039771 140 5.904980 141 5.508910 142 5.869335 143 5.378663 Name: form_factor, Length: 144, dtype: float64 """ area = geometry.area volume = area * height surface = (geometry.length * height) + geometry.area zeros = volume == 0 res = np.empty(len(geometry)) res[zeros] = np.nan res[~zeros] = surface[~zeros] / (volume[~zeros] ** (2 / 3)) return Series(res, index=geometry.index, name="form_factor") def fractal_dimension( geometry: GeoDataFrame | GeoSeries, ) -> Series: """Calculates fractal dimension based on area and perimeter. .. math:: {2log({{perimeter} \\over {4}})} \\over log(area) Based on :cite:`mcgarigal1995fragstats`. 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.fractal_dimension(buildings) 0 1.072678 1 1.182350 2 1.018422 3 1.048314 4 1.017328 ... 139 1.014975 140 1.033581 141 1.064103 142 1.022617 143 1.000008 Length: 144, dtype: float64 """ return (2 * np.log(geometry.length / 4)) / np.log(geometry.area) def facade_ratio( geometry: GeoDataFrame | GeoSeries, ) -> Series: """ Calculates the facade ratio of each object given its geometry. .. math:: area \\over perimeter Adapted from :cite:`schirmer2015`. 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.facade_ratio(buildings) 0 5.310716 1 11.314008 2 5.963959 3 6.376086 4 5.987687 ... 139 1.868981 140 4.046407 141 5.963454 142 1.711740 143 2.734681 Length: 144, dtype: float64 """ return geometry.area / geometry.length def circular_compactness(geometry: GeoDataFrame | GeoSeries) -> Series: """Calculates the circular compactness of each object given its geometry. .. math:: area \\over \\textit{area of enclosing circle} Adapted from :cite:`dibble2017`. 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.circular_compactness(buildings) 0 0.572145 1 0.390417 2 0.588332 3 0.520411 4 0.591297 ... 139 0.563433 140 0.525733 141 0.404872 142 0.530405 143 0.636387 Length: 144, dtype: float64 """ return geometry.area / ( np.pi * shapely.minimum_bounding_radius(geometry.geometry.array) ** 2 ) def square_compactness(geometry: GeoDataFrame | GeoSeries) -> Series: """Calculates the square compactness of each object given its geometry. .. math:: \\begin{equation*} \\left(\\frac{4 \\sqrt{area}}{perimeter}\\right) ^ 2 \\end{equation*} Adapted from :cite:`feliciotti2018`. 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.square_compactness(buildings) 0 0.619387 1 0.182604 2 0.887750 3 0.719750 4 0.894034 ... 139 0.940666 140 0.824085 141 0.647576 142 0.914813 143 0.999961 Length: 144, dtype: float64 """ return ((np.sqrt(geometry.area) * 4) / geometry.length) ** 2 def convexity(geometry: GeoDataFrame | GeoSeries) -> Series: """Calculates the convexity of each object given its geometry. .. math:: \\frac{\\textit{area}}{\\textit{area of convex hull}} Adapted from :cite:`dibble2017`. 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.convexity(buildings) 0 0.815196 1 0.703008 2 0.953398 3 0.890489 4 0.957323 ... 139 1.000000 140 0.904176 141 0.795774 142 1.000000 143 1.000000 Length: 144, dtype: float64 """ return geometry.area / geometry.geometry.convex_hull.area def courtyard_index( geometry: GeoDataFrame | GeoSeries, courtyard_area: NDArray[np.float64] | Series | None = None, ) -> Series: """Calculates the courtyard index of each object given its geometry. .. math:: \\textit{area of courtyards} \\over \\textit{total area} Adapted from :cite:`schirmer2015`. Parameters ---------- geometry : GeoDataFrame | GeoSeries A GeoDataFrame or GeoSeries containing polygons to analyse. courtyard_area : NDArray[np.float64] | Series | None, optional array of courtyard areas. If None, it will be calculated, by default None Returns ------- Series Examples -------- >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") >>> ci = momepy.courtyard_index(buildings) >>> ci 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 Length: 144, dtype: float64 >>> ci.max() np.float64(0.16605915738643523) If you know the courtyard area, you can pass it to skip the computation step. >>> courtyard_area = momepy.courtyard_area(buildings) >>> momepy.courtyard_index(buildings, courtyard_area=courtyard_area) 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 Length: 144, dtype: float64 """ if courtyard_area is None: courtyard_area = _dimension.courtyard_area(geometry) return courtyard_area / geometry.area def rectangularity(geometry: GeoDataFrame | GeoSeries) -> Series: """Calculates the rectangularity of each object given its geometry. .. math:: \\frac{\\textit{area}}{\\textit{area of minimum bounding rectangle}} Adapted from :cite:`dibble2017`. 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.rectangularity(buildings) 0 0.694268 1 0.702242 2 0.901582 3 0.821797 4 0.912858 ... 139 0.996876 140 0.820865 141 0.659281 142 0.971600 143 0.999400 Length: 144, dtype: float64 """ return geometry.area / shapely.area( shapely.minimum_rotated_rectangle(geometry.geometry.array) ) def shape_index( geometry: GeoDataFrame | GeoSeries, longest_axis_length: NDArray[np.float64] | Series | None = None, ) -> Series: """Calculates the shape index of each object given its geometry. .. math:: {\\sqrt{{area} \\over {\\pi}}} \\over {0.5 * \\textit{longest axis}} Parameters ---------- geometry : GeoDataFrame | GeoSeries A GeoDataFrame or GeoSeries containing polygons to analyse. longest_axis_length : NDArray[np.float64] | Series | None, optional array of longest axis lengths. If None, it will be calculated, by default None Returns ------- Series Examples -------- >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") >>> momepy.shape_index(buildings) 0 0.756403 1 0.624834 2 0.767028 3 0.721395 4 0.768958 ... 139 0.750622 140 0.725074 141 0.636296 142 0.728289 143 0.797739 Length: 144, dtype: float64 If you know the longest axis length, you can pass it to skip the computation step. >>> lal = momepy.longest_axis_length(buildings) >>> momepy.shape_index(buildings, longest_axis_length=lal) 0 0.756403 1 0.624834 2 0.767028 3 0.721395 4 0.768958 ... 139 0.750622 140 0.725074 141 0.636296 142 0.728289 143 0.797739 Length: 144, dtype: float64 """ if longest_axis_length is None: longest_axis_length = _dimension.longest_axis_length(geometry) return np.sqrt(geometry.area / np.pi) / (0.5 * longest_axis_length) def corners( geometry: GeoDataFrame | GeoSeries, eps: float = 10, include_interiors: bool = False, ) -> Series: """Calculates the number of corners of each object given its geometry. As a corner is considered a point where the angle between two consecutive segments deviates from 180 degrees by more than ``eps``. Parameters ---------- geometry : GeoDataFrame | GeoSeries A GeoDataFrame or GeoSeries containing polygons to analyse. eps : float, optional Deviation from 180 degrees to consider a corner, by default 10 include_interiors : bool, optional If True, polygon interiors are included in the calculation. If False, only exterior is considered, by default False Returns ------- Series Examples -------- >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") >>> momepy.corners(buildings) 0 24 1 43 2 8 3 16 4 8 .. 139 4 140 6 141 6 142 4 143 4 Length: 144, dtype: int64 """ if not GPD_013: raise ImportError("momepy.corners requires geopandas 0.13 or later. ") def _count_corners(points: DataFrame, eps: float) -> int: pts = points.values[:-1] true_angles = _true_angles_mask(pts, eps=eps) corners = np.count_nonzero(true_angles) return corners if include_interiors: coords = geometry.reset_index(drop=True).get_coordinates(index_parts=False) else: coords = geometry.reset_index(drop=True).exterior.get_coordinates( index_parts=False ) cc = coords.groupby(level=0).apply(_count_corners, eps=eps) cc.index = geometry.index return cc def squareness( geometry: GeoDataFrame | GeoSeries, eps: float = 10, include_interiors: bool = False, ) -> Series: """Calculates the squareness of each object given its geometry. Squareness is a mean deviation of angles at corners from 90 degrees. As a corner is considered a point where the angle between two consecutive segments deviates from 180 degrees by more than ``eps``. Parameters ---------- geometry : GeoDataFrame | GeoSeries A GeoDataFrame or GeoSeries containing polygons to analyse. eps : float, optional Deviation from 180 degrees to consider a corner, by default 10 include_interiors : bool, optional If True, polygon interiors are included in the calculation. If False, only exterior is considered, by default False Returns ------- Series Examples -------- >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") >>> momepy.squareness(buildings) 0 3.707582 1 2.990318 2 0.437987 3 4.573564 4 0.382146 ... 139 0.214733 140 0.307916 141 0.343259 142 0.899731 143 0.028140 Length: 144, dtype: float64 """ if not GPD_013: raise ImportError("momepy.squareness requires geopandas 0.13 or later. ") def _squareness(points: DataFrame, eps: float): pts = points.values[:-1] true_angles, degrees = _true_angles_mask(pts, eps=eps, return_degrees=True) return np.nanmean(np.abs(90 - degrees[true_angles])) if include_interiors: coords = geometry.reset_index(drop=True).get_coordinates(index_parts=False) else: coords = geometry.reset_index(drop=True).exterior.get_coordinates( index_parts=False ) sq = coords.groupby(level=0).apply(_squareness, eps=eps) sq.index = geometry.index return sq def equivalent_rectangular_index(geometry: GeoDataFrame | GeoSeries) -> Series: """Calculates the equivalent rectangular index of each object given its geometry. .. math:: \\sqrt{{area} \\over \\textit{area of bounding rectangle}} * {\\textit{perimeter of bounding rectangle} \\over {perimeter}} Based on :cite:`basaraner2017`. 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.equivalent_rectangular_index(buildings) 0 0.787923 1 0.443137 2 0.954252 3 0.851658 4 0.957543 ... 139 1.000050 140 0.907837 141 0.813269 142 0.995926 143 0.999999 Length: 144, dtype: float64 """ bbox = shapely.minimum_rotated_rectangle(geometry.geometry.array) return np.sqrt(geometry.area / shapely.area(bbox)) * ( shapely.length(bbox) / geometry.length ) def elongation(geometry: GeoDataFrame | GeoSeries) -> Series: """Calculates the elongation of each object given its geometry. The elongation is defined as the elongation of the minimum bounding rectangle. .. math:: {{p - \\sqrt{p^2 - 16a}} \\over {4}} \\over {{{p} \\over {2}} - {{p - \\sqrt{p^2 - 16a}} \\over {4}}} where `a` is the area of the object and `p` its perimeter. Based on :cite:`gil2012`. 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.elongation(buildings) 0 0.908244 1 0.581318 2 0.726527 3 0.838840 4 0.727294 ... 139 0.608004 140 0.979998 141 0.747326 142 0.564060 143 0.987953 Name: elongation, Length: 144, dtype: float64 """ bbox = shapely.minimum_rotated_rectangle(geometry.geometry.array) a = shapely.area(bbox) p = shapely.length(bbox) sqrt = np.maximum(p**2 - 16 * a, 0) elo1 = ((p - np.sqrt(sqrt)) / 4) / ((p / 2) - ((p - np.sqrt(sqrt)) / 4)) elo2 = ((p + np.sqrt(sqrt)) / 4) / ((p / 2) - ((p + np.sqrt(sqrt)) / 4)) res = np.where(elo1 <= elo2, elo1, elo2) return Series(res, index=geometry.index, name="elongation") def centroid_corner_distance( geometry: GeoDataFrame | GeoSeries, eps: float = 10, include_interiors: bool = False, ) -> DataFrame: """Calculates the centroid-corner distance of each object given its geometry. As a corner is considered a point where the angle between two consecutive segments deviates from 180 degrees by more than ``eps``. Parameters ---------- geometry : GeoDataFrame | GeoSeries A GeoDataFrame or GeoSeries containing polygons to analyse. eps : float, optional Deviation from 180 degrees to consider a corner, by default 10 include_interiors : bool, optional If True, polygon interiors are included in the calculation. If False, only exterior is considered, by default False Returns ------- DataFrame DataFrame with columns 'mean' and 'std' Examples -------- >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") >>> momepy.centroid_corner_distance(buildings).head() mean std 0 15.961532 3.081063 1 58.763388 22.922368 2 14.988106 3.648731 3 15.000439 4.999226 4 14.965557 3.660826 """ if not GPD_013: raise ImportError( "momepy.centroid_corner_distance requires geopandas 0.13 or later. " ) result_index = geometry.index if isinstance(geometry.index, MultiIndex): geometry = geometry.reset_index(drop=True) def _ccd(points: DataFrame, eps: float) -> Series: centroid = points.values[0, 2:] pts = points.values[:-1, :2] true_angles = _true_angles_mask(pts, eps=eps) dists = np.linalg.norm(pts[true_angles] - centroid, axis=1) return Series({"mean": np.nanmean(dists), "std": np.nanstd(dists)}) if include_interiors: coords = geometry.get_coordinates(index_parts=False) else: coords = geometry.exterior.get_coordinates(index_parts=False) coords[["cent_x", "cent_y"]] = geometry.centroid.get_coordinates(index_parts=False) ccd = coords.groupby(level=0).apply(_ccd, eps=eps) ccd.index = result_index return ccd def linearity(geometry: GeoDataFrame | GeoSeries) -> Series: """Calculates the linearity of each LineString The linearity is defined as the ratio of the length of the segment between the first and last point to the length of the LineString. While other geometry types are accepted, the result is not well defined. .. math:: \\frac{l_{euclidean}}{l_{segment}} where `l` is the length of the LineString. Adapted from :cite:`araldi2019`. Parameters ---------- geometry : GeoDataFrame | GeoSeries A GeoDataFrame or GeoSeries containing lines to analyse. Returns ------- Series Examples -------- >>> path = momepy.datasets.get_path("bubenec") >>> streets = geopandas.read_file(path, layer="streets") >>> momepy.linearity(streets).head() 0 1.000000 1 0.995987 2 0.999653 3 0.999997 4 1.000000 dtype: float64 """ return ( shapely.distance( shapely.get_point(geometry.geometry.array, 0), shapely.get_point(geometry.geometry.array, -1), ) / geometry.length ) def compactness_weighted_axis( geometry: GeoDataFrame | GeoSeries, longest_axis_length: NDArray[np.float64] | Series | None = None, ) -> Series: """Calculates the compactness-weighted axis of each object in a given GeoDataFrame. .. math:: d_{i} \\times\\left(\\frac{4}{\\pi}-\\frac{16 (area_{i})} {perimeter_{i}^{2}}\\right) Parameters ---------- geometry : GeoDataFrame | GeoSeries A GeoDataFrame or GeoSeries containing polygons to analyse. longest_axis_length : NDArray[np.float64] | Series | None, optional array of longest axis lengths. If None, it will be calculated, by default None Returns ------- Series Examples -------- >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") >>> momepy.compactness_weighted_axis(buildings) 0 26.327730 1 208.588747 2 14.358370 3 26.026429 4 14.095118 ... 139 3.853623 140 12.462700 141 32.888900 142 3.975422 143 4.228395 Length: 144, dtype: float64 If you know the longest axis length, you can pass it to skip the computation step. >>> lal = momepy.longest_axis_length(buildings) >>> momepy.compactness_weighted_axis(buildings, longest_axis_length=lal) 0 26.327730 1 208.588747 2 14.358370 3 26.026429 4 14.095118 ... 139 3.853623 140 12.462700 141 32.888900 142 3.975422 143 4.228395 Length: 144, dtype: float64 """ if longest_axis_length is None: longest_axis_length = _dimension.longest_axis_length(geometry) return longest_axis_length * ( (4 / np.pi) - (16 * geometry.area) / (geometry.length**2) ) # helper functions def _true_angles_mask( points: NDArray[np.float64], eps: float, return_degrees: bool = False ) -> NDArray[np.bool_] | tuple[NDArray[np.bool_], NDArray[np.float64]]: """Calculates the mask of true angles. Parameters ---------- points : NDArray[np.float64] array of points eps : float Deviation from 180 degrees to consider a corner return_degrees : bool, optional If True, returns also degrees, by default False Returns ------- NDArray[np.bool_] | tuple[NDArray[np.bool_], NDArray[np.float64]] boolean array or a tuple of boolean array and float array of degrees """ a = np.roll(points, 1, axis=0) b = points c = np.roll(points, -1, axis=0) ba = a - b bc = c - b cosine_angle = np.sum(ba * bc, axis=1) / ( np.linalg.norm(ba, axis=1) * np.linalg.norm(bc, axis=1) ) angles = np.arccos(cosine_angle) degrees = np.degrees(angles) if return_degrees: return np.logical_or(degrees <= 180 - eps, degrees >= 180 + eps), degrees return np.logical_or(degrees <= 180 - eps, degrees >= 180 + eps)