added assign-land-use-class method in geospatial process

measprocess/geospatial/process.py

@@ -2,7 +2,7 @@ import geopandas as gpd
 import matplotlib.pyplot as plt
 from shapely.geometry import LineString, Point, Polygon, box
 from shapely.ops import unary_union
-
+from tqdm import tqdm
 
 def project_onto_streets(
     point_series: gpd.GeoSeries,
@@ -203,3 +203,76 @@ def delete_indoor_meas(
     outdoor_meas_gdf = meas_gdf.drop(indoor_meas_gdf.index).reset_index()
 
     return outdoor_meas_gdf.geometry
+
+def assign_landuse_classes(
+    point_geodataframe: gpd.GeoDataFrame,
+    polygon_geodataframe: gpd.GeoDataFrame,
+    cut_polygons_by_buffer: bool = False,
+    buffer_size: integer = 0,
+    epsg: str = "EPSG:31287",
+) -> gpd.GeoDataFrame:
+    '''
+    Returns a DataFrame with the landuse class distribution in the area defined by the buffer_size around the passed locations
+    Landuseclasses are defined according to the Corine Land Cover Project:
+    https://land.copernicus.eu/user-corner/technical-library/corine-land-cover-nomenclature-guidelines/html
+
+    For Austria the land_use_classes can be obtained via:
+    https://www.data.gv.at/katalog/dataset/76617316-b9e6-4bcd-ba09-e328b578fed2
+
+    For Germany the source is:
+    https://gis.uba.de/catalog/OpenSearch.do?type=/Query/OpenSearch.do&search=CORINE
+
+    :param point_geodataframe: A Geodataframe with geometry of type point for which the landuse classes will be evaluated.
+    :param polygon_geodataframe: A GeoDataFrame with the landuse class polygons. Can be retrieved from
+    :param cut_polygons_by_buffer: If True the polygon_geodataframe will be cut by the total bounds of the point series. Can result in a speedup.
+    :param buffer_size: The buffer_size in meters used to construct the area around points for which the distribution of classes will be computed.
+    :param epsg: The projection used for the computation of the area. Defaults to EPSG:31287 for Austria.
+
+    :return: A GeoDataFrame with columns describing the percentage of the area assigned to the respective land use class.
+    '''
+
+    if point_geodataframe.crs is None or polygon_geodataframe.crs is None:
+        raise ValueError(
+            "CRS needs to be set explicitly for passed Geoseries: Use .set_crs()."
+        )
+
+    point_geodataframe['first5codes'] = None
+    point_geodataframe['first5percent'] = None
+
+    if buffer_size < 0:
+        raise ValueError(
+            "Buffer size has to be a positive integer. Default value is zero."
+        )
+    if buffer_size != 0:
+        point_geodataframe.geometry = point_geodataframe.geometry.buffer(buffer_size)
+
+    if cut_polygons_by_buffer == True:
+        xmin, ymin, xmax, ymax = point_geodataframe.to_crs(epsg).total_bounds
+        polygon_geodataframe = polygon_geodataframe.cx[xmin-buffer_size:xmax+buffer_size, ymin-buffer_size:ymax+buffer_size]
+
+    point_geodataframe = point_geodataframe.to_crs(epsg)
+    polygon_geodataframe = polygon_geodataframe.to_crs(epsg)
+
+    total = len(point_geodataframe)
+    results = {}
+
+    for i, element in tqdm(point_geodataframe.iterrows(), total=total):
+        point_buffer = element.geometry
+
+        percentages = []
+        codes = []
+        for j, poly in enumerate(polygon_geodataframe.geometry):
+            if point_buffer.intersects(poly):
+                percentages.append(point_buffer.intersection(poly).area/point_buffer.area*100)
+                codes.append(polygon_geodataframe.iloc[j]['CODE_18'])
+        code_percent_dict = {}
+        for key, value in zip(codes, percentages):
+            code_percent_dict[key] = code_percent_dict.get(key, 0) + value
+        code_percent_dict = dict(sorted(code_percent_dict.items(), key=lambda x: x[1], reverse =True))
+        first5pairs = {k: dict(code_percent_dict)[k] for k in list(dict(code_percent_dict))[:5]}
+
+        results[i] = first5pairs
+
+    result_df = pd.DataFrame(results).transpose()
+    result_df.index = point_geodataframe.index
+    return result_df