import geopandas as gpd import numpy as np import shapely from geopandas import GeoDataFrame, GeoSeries from libpysal.graph import Graph from numpy.typing import NDArray from packaging.version import Version from pandas import Series from scipy import sparse __all__ = [ "orientation", "shared_walls", "alignment", "neighbor_distance", "mean_interbuilding_distance", "building_adjacency", "neighbors", "street_alignment", "cell_alignment", ] GPD_GE_013 = Version(gpd.__version__) >= Version("0.13.0") def orientation(geometry: GeoDataFrame | GeoSeries) -> Series: """Calculate the orientation of objects. The 'orientation' is defined as the deviation of orientation of the bounding rectangle from cardinal directions. As such it is within a range 0 - 45. The orientation of LineStrings is represented by the orientation of the line connecting the first and the last point of the segment. 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.orientation(buildings) 0 41.051468 1 20.829200 2 20.969610 3 20.811525 4 20.412895 ... 139 21.107057 140 13.669062 141 13.671217 142 21.463914 143 13.754569 Name: orientation, Length: 144, dtype: float64 """ geom_types = geometry.geom_type poly_mask = geom_types.str.contains("Polygon") line_mask = geom_types.str.contains("Line") result = np.full(len(geometry), np.nan, dtype=float) if poly_mask.any(): bboxes = shapely.minimum_rotated_rectangle(geometry.geometry.loc[poly_mask]) coords = shapely.get_coordinates(bboxes) pt0 = coords[::5] pt1 = coords[1::5] angle = np.arctan2(pt1[:, 0] - pt0[:, 0], pt1[:, 1] - pt0[:, 1]) result[poly_mask] = np.degrees(angle) if line_mask.any(): first = shapely.get_point(geometry.geometry, 0) last = shapely.get_point(geometry.geometry, -1) pt0 = shapely.get_coordinates(first) pt1 = shapely.get_coordinates(last) angle = np.arctan2(pt1[:, 0] - pt0[:, 0], pt1[:, 1] - pt0[:, 1]) result[line_mask] = np.degrees(angle) return Series( np.abs((result + 45) % 90 - 45), index=geometry.index, dtype=float, name="orientation", ) def shared_walls( geometry: GeoDataFrame | GeoSeries, strict: bool = True, tolerance: float = 0.01 ) -> Series: """Calculate the length of shared walls of adjacent elements (typically buildings). Note that data needs to be topologically correct. Overlapping polygons will lead to incorrect results. Adapted from :cite:`hamaina2012a`. Parameters ---------- geometry : GeoDataFrame | GeoSeries A GeoDataFrame or GeoSeries containing polygons to analyse. strict : bool Perform calculations based on strict contiguity. If set to `False`, consider overlapping or nearly overlapping polygons as touching. tolerance: float Tolerance for non-strict calculations, if strict is True, tolerance has no effect on the results. Returns ------- Series Examples -------- >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") >>> momepy.shared_walls(buildings) 0 0.000000 1 33.805154 2 39.766297 3 40.604643 4 40.455735 ... 139 9.886397 140 0.000000 141 0.000000 142 19.332735 143 10.876113 Name: shared_walls, Length: 144, dtype: float64 """ predicate = "touches" if not strict: orig_lengths = geometry.length geometry = geometry.buffer(tolerance) predicate = "intersects" if GPD_GE_013: inp, res = geometry.sindex.query(geometry.geometry, predicate=predicate) else: inp, res = geometry.sindex.query_bulk(geometry.geometry, predicate=predicate) mask = inp != res inp, res = inp[mask], res[mask] left = geometry.geometry.take(inp).reset_index(drop=True) right = geometry.geometry.take(res).reset_index(drop=True) intersections = left.intersection(right).length walls = intersections.groupby(inp).sum() walls.index = geometry.index.take(walls.index) results = Series(0.0, index=geometry.index, name="shared_walls") results.loc[walls.index] = walls if not strict: results = (results / 2) - (2 * tolerance) results = results.clip(0, orig_lengths) return results def alignment(orientation: Series, graph: Graph) -> Series: """Calculate the mean deviation of orientation adjacent elements .. math:: \\frac{1}{n}\\sum_{i=1}^n dev_i=\\frac{dev_1+dev_2+\\cdots+dev_n}{n} Takes orientation of adjacent elements defined in ``graph`` and calculates the mean deviation. Notes ----- The index of ``orientation`` must match the index along which the ``graph`` is built. Parameters ---------- orientation : Series A series containing orientation (e.g. measured by the :func:`orientation` function) indexed using the same index that has been used to build the graph. graph : libpysal.graph.Graph Graph representing spatial relationships between elements. Returns ------- Series Examples -------- >>> from libpysal import graph >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") >>> orientation = momepy.orientation(buildings) Define a spatial graph that includes observations within its own neighborhood: >>> delaunay = graph.Graph.build_triangulation( ... buildings.centroid ... ).assign_self_weight() >>> delaunay Alignment of orienation within triangulated neighbors: >>> momepy.alignment(orientation, delaunay) 0 14.639585 1 0.217417 2 0.205626 3 0.151730 4 0.352692 ... 139 1.970642 140 0.127322 141 0.161906 142 4.350890 143 0.084884 Name: alignment, Length: 144, dtype: float64 """ orientation_expanded = orientation.loc[graph._adjacency.index.get_level_values(1)] orientation_expanded.index = graph._adjacency.index.get_level_values(0) r = (orientation_expanded - orientation).abs().groupby(level=0).mean() r.name = "alignment" return r def neighbor_distance(geometry: GeoDataFrame | GeoSeries, graph: Graph) -> Series: """Calculate the mean distance to adjacent elements. Takes geometry of adjacent elements defined in ``graph`` and calculates the mean distance. Notes ----- The index of ``geometry`` must match the index along which the ``graph`` is built. Parameters ---------- geometry : GeoDataFrame | GeoSeries A GeoDataFrame or GeoSeries containing geometries to analyse. graph : libpysal.graph.Graph Graph representing spatial relationships between elements. Returns ------- Series Examples -------- >>> from libpysal import graph >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") Define a spatial graph: >>> delaunay = graph.Graph.build_triangulation(buildings.centroid) >>> delaunay Mean distance to adjacent buildings within triangulated neighbors: >>> momepy.neighbor_distance(buildings, delaunay) 0 29.185890 1 30.244905 2 47.052305 3 22.831824 4 16.183615 ... 139 39.698734 140 20.634252 141 38.208668 142 6.304569 143 14.551355 Length: 144, dtype: float64 """ geoms = geometry.geometry.loc[graph._adjacency.index.get_level_values(1)] geoms.index = graph._adjacency.index.get_level_values(0) mean_distance = ( (geoms.distance(geometry.geometry, align=True)).groupby(level=0).mean() ) mean_distance.loc[graph.isolates] = np.nan return mean_distance def mean_interbuilding_distance( geometry: GeoDataFrame | GeoSeries, adjacency_graph: Graph, neighborhood_graph: Graph, ) -> Series: """Calculate the mean distance between adjacent geometries within a set neighborhood For each building, this function takes a neighborhood based on the neighbors within a ``neighborhood_graph`` and calculates the mean distance between adjacent buildings within this neighborhood where adjacency is captured by ``adjacency_graph``. Notes ----- The index of ``geometry`` must match the index along which both of the graphs are built. Parameters ---------- geometry : GeoDataFrame | GeoSeries A GeoDataFrame or GeoSeries containing geometries to analyse. adjacency_graph : libpysal.graph.Graph Graph representing the adjacency of geometries. Typically, this is a contiguity graph derived from tessellation cells linked to buildings. neighborhood_graph : libpysal.graph.Graph Graph representing the extent around each geometry within which to calculate the mean interbuilding distance. This can be a distance based graph, KNN graph, higher order contiguity, etc. Returns ------- Series Examples -------- >>> from libpysal import graph >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") Define a spatial graph denoting building adjacency: >>> delaunay = graph.Graph.build_triangulation(buildings.centroid) >>> delaunay Define a spatial graph denoting the neighborhood: >>> knn15 = graph.Graph.build_knn(buildings.centroid, k=15) >>> knn15 Measure mean interbuilding distance: >>> momepy.mean_interbuilding_distance(buildings, delaunay, knn15) 0 29.516506 1 18.673132 2 23.277728 3 25.409034 4 18.454463 ... 139 21.642580 140 16.427126 141 17.792155 142 10.844367 143 14.896066 Name: mean_interbuilding_distance, Length: 144, dtype: float64 """ distance = Series( shapely.distance( geometry.geometry.loc[ adjacency_graph._adjacency.index.get_level_values(0) ].values, geometry.geometry.loc[ adjacency_graph._adjacency.index.get_level_values(1) ].values, ), index=adjacency_graph._adjacency.index, name="distance", ) distance_matrix = ( distance.astype("Sparse[float]").sparse.to_coo(sort_labels=True)[0].tocsr() ) neighborhood_matrix = sparse.coo_matrix(neighborhood_graph.sparse).tocsr() mean_distances = np.zeros(distance_matrix.shape[0], dtype=float) for i in range(distance_matrix.shape[0]): neighborhood_indices = np.append(neighborhood_matrix[i].indices, i) sub_matrix = distance_matrix[neighborhood_indices][:, neighborhood_indices] mean_distances[i] = sub_matrix.sum() / sub_matrix.nnz return Series( mean_distances, index=geometry.index, name="mean_interbuilding_distance" ) def building_adjacency( contiguity_graph: Graph, neighborhood_graph: Graph, ) -> Series: """Calculate the level of building adjacency. Building adjacency reflects how much buildings tend to join together into larger structures. It is calculated as a ratio of joined built-up structures captured by ``contiguity_graph`` and buildings within the neighborhood defined in ``neighborhood_graph``. Adapted from :cite:`vanderhaegen2017`. Notes ----- Both graphs must be built on the same index. Parameters ---------- contiguity_graph : libpysal.graph.Graph Graph representing contiguity between geometries, typically a rook contiguity graph derived from buildings. neighborhood_graph : libpysal.graph.Graph Graph representing the extent around each geometry within which to calculate the level of building adjacency. This can be a distance based graph, KNN graph, higher order contiguity, etc. Returns ------- Series Examples -------- >>> from libpysal import graph >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") Define a spatial graph denoting building contiguity: >>> contig = graph.Graph.build_contiguity(buildings) >>> contig Define a spatial graph denoting the neighborhood: >>> knn15 = graph.Graph.build_knn(buildings.centroid, k=15) >>> knn15 Measure mean interbuilding distance: >>> momepy.building_adjacency(contig, knn15) 0 0.6875 1 0.1875 2 0.1875 3 0.2500 4 0.1875 ... 139 0.4375 140 0.1875 141 0.1875 142 0.1875 143 0.2500 Name: building_adjacency, Length: 144, dtype: float64 """ components = contiguity_graph.component_labels # check if self-weights are present, otherwise assign them to treat self as part of # the neighborhood has_self_weights = ( neighborhood_graph._adjacency.index.get_level_values("focal") == neighborhood_graph._adjacency.index.get_level_values("neighbor") ).sum() == neighborhood_graph.n if not has_self_weights: neighborhood_graph = neighborhood_graph.assign_self_weight() grouper = components.loc[ neighborhood_graph._adjacency.index.get_level_values(1) ].groupby(neighborhood_graph._adjacency.index.get_level_values(0)) result = grouper.agg("nunique") / grouper.agg("count") result.name = "building_adjacency" result.index.name = None return result def neighbors( geometry: GeoDataFrame | GeoSeries, graph: Graph, weighted=False ) -> Series: """Calculate the number of neighbours captured by ``graph``. If ``weighted=True``, the number of neighbours will be divided by the perimeter of the object to return a relative value (neighbors per meter). Adapted from :cite:`hermosilla2012`. Notes ----- The index of ``geometry`` must match the index along which the ``graph`` is built. Parameters ---------- gdf : GeoDataFrame | GeoSeries GeoDataFrame containing geometries to analyse. graph : libpysal.graph.Graph Graph representing spatial relationships between elements. weighted : bool If True, the number of neighbours will be divided by the perimeter of the object to return a relative value (neighbors per meter). Returns ------- Series Examples -------- >>> from libpysal import graph >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") >>> tessellation = momepy.morphological_tessellation(buildings) Define a spatial graph denoting adjacency: >>> contig = graph.Graph.build_contiguity(tessellation) >>> contig Number of neighbors of each tessellation cell: >>> momepy.neighbors(tessellation, contig) focal 0 4 1 9 2 3 3 3 4 7 .. 139 3 140 6 141 12 142 2 143 5 Name: neighbors, Length: 144, dtype: int64 Weighted by the tessellation area: >>> momepy.neighbors(tessellation, contig, weighted=True) focal 0 0.012732 1 0.010116 2 0.013350 3 0.010172 4 0.038916 ... 139 0.020037 140 0.036766 141 0.045287 142 0.044147 143 0.051799 Name: neighbors, Length: 144, dtype: float64 """ if weighted: r = graph.cardinalities / geometry.length else: r = graph.cardinalities r.name = "neighbors" return r def street_alignment( building_orientation: Series, street_orientation: Series, street_index: Series, ) -> Series: """Calulate the deviation of the building orientation from the street orientation. Parameters ---------- building_orientation : Series Series with the orientation of buildings. Can be measured using :func:`orientation`. street_orientation : Series Series with the orientation of streets. Can be measured using :func:`orientation`. street_index : Series Series with the index of the street to which the building belongs. Can be retrieved using :func:`momepy.get_nearest_street`. Returns ------- Series Examples -------- >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") >>> streets = geopandas.read_file(path, layer="streets") Get street index. >>> buildings["street_index"] = momepy.get_nearest_street(buildings, streets) >>> buildings.head() uID geometry street_index 0 1 POLYGON ((1603599.221 6464369.816, 1603602.984... 0.0 1 2 POLYGON ((1603042.88 6464261.498, 1603038.961 ... 33.0 2 3 POLYGON ((1603044.65 6464178.035, 1603049.192 ... 10.0 3 4 POLYGON ((1603036.557 6464141.467, 1603036.969... 8.0 4 5 POLYGON ((1603082.387 6464142.022, 1603081.574... 8.0 Compute orientations. >>> blg_orient = momepy.orientation(buildings) >>> str_orient = momepy.orientation(streets) Compute street alignment. >>> momepy.street_alignment(blg_orient, str_orient, buildings["street_index"]) 0 0.290739 1 4.542071 2 0.105745 3 0.424903 4 0.823533 ... 139 1.779876 140 0.109254 141 0.466453 142 1.223387 143 0.455081 Name: street_alignment, Length: 144, dtype: float64 """ r = (building_orientation - street_orientation.loc[street_index].values).abs() r.name = "street_alignment" return r def cell_alignment( left_orientation: NDArray[np.float64] | Series, right_orientation: NDArray[np.float64] | Series, ) -> Series: """ Calculate the difference between cell orientation and the orientation of object. .. math:: \\left|{\\textit{building orientation} - \\textit{cell orientation}}\\right| Notes ----- ``left_orientation`` and ``right_orientation`` must be aligned or have an index. Parameters ---------- left_orientation : np.array, pd.Series The ``np.array``, or `pd.Series`` with orientation of cells. This can be calculated using :func:`orientation`. right_orientation : np.array, pd.Series The ``np.array`` or ``pd.Series`` with orientation of objects. This can be calculated using :func:`orientation`. Returns ------- Series Examples -------- >>> from libpysal import graph >>> path = momepy.datasets.get_path("bubenec") >>> buildings = geopandas.read_file(path, layer="buildings") >>> tessellation = momepy.morphological_tessellation(buildings) Measure orientations: >>> blg_orient = momepy.orientation(buildings) >>> tess_orient = momepy.orientation(tessellation) Compute alignment: >>> momepy.cell_alignment(blg_orient, tess_orient) 0 0.854788 1 20.829200 2 5.552120 3 4.052674 4 0.159289 ... 139 0.189750 140 17.920139 141 0.393708 142 0.024618 143 0.252122 Name: orientation, Length: 144, dtype: float64 """ if not isinstance(left_orientation, Series): left_orientation = Series(left_orientation) if not isinstance(right_orientation, Series): right_orientation = Series(right_orientation) return (left_orientation - right_orientation).abs()