diff --git a/measprocess/geospatial/process.py b/measprocess/geospatial/process.py index 77452a5ede18f4e97e1da0f0c169ab2d75b54be2..73c78646882722385522002b6f86097764f22ae2 100644 --- a/measprocess/geospatial/process.py +++ b/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