From 740504cc17411bcefe3574e736e0aec2696e8309 Mon Sep 17 00:00:00 2001 From: Sujan Adhikari <109404840+Sujanadh@users.noreply.github.com> Date: Mon, 9 Dec 2024 18:23:55 +0545 Subject: [PATCH] fix(backend): polygon geometries; remove holes, fix right hand rule (#1961) * feat: clean polygon geometries; remove holes, fix right hand rule, create single aoi * refactor: change bad request http status to unprocessable entity --- src/backend/app/db/postgis_utils.py | 87 ++++++++++++++++----- src/backend/app/projects/project_routes.py | 5 +- src/backend/app/projects/project_schemas.py | 4 +- 3 files changed, 72 insertions(+), 24 deletions(-) diff --git a/src/backend/app/db/postgis_utils.py b/src/backend/app/db/postgis_utils.py index a84221758d..2dac3b9f02 100644 --- a/src/backend/app/db/postgis_utils.py +++ b/src/backend/app/db/postgis_utils.py @@ -32,7 +32,7 @@ from osm_rawdata.postgres import PostgresClient from psycopg import Connection, ProgrammingError from psycopg.rows import class_row -from shapely.geometry import mapping, shape +from shapely.geometry import MultiPolygon, Polygon, mapping, shape from shapely.geometry.base import BaseGeometry from shapely.ops import unary_union @@ -698,7 +698,7 @@ async def javarosa_to_geojson_geom(javarosa_geom_string: str, geom_type: str) -> def multigeom_to_singlegeom( featcol: geojson.FeatureCollection, ) -> geojson.FeatureCollection: - """Converts any Multi(xxx) geometry types to individual geometries. + """Converts any Multi(xxx) geometry types to list of individual geometries. Args: featcol : A GeoJSON FeatureCollection of geometries. @@ -723,8 +723,6 @@ def split_multigeom(geom, properties): geom = shape(feature["geometry"]) except ValueError: log.warning(f"Geometry is not valid, so was skipped: {feature['geometry']}") - continue - if geom.geom_type.startswith("Multi"): # Handle all MultiXXX types final_features.extend(split_multigeom(geom, properties)) @@ -740,12 +738,55 @@ def split_multigeom(geom, properties): return geojson.FeatureCollection(final_features) +def remove_holes(polygon: Polygon): + """Detect and remove holes within a polygon.""" + if polygon.interiors: + return Polygon(polygon.exterior) # Keep only the exterior ring + return polygon + + +def create_single_polygon(multipolygon: MultiPolygon, dissolve_polygon: bool): + """If a MultiPolygon can create a common exterior ring, return a single AOI Polygon. + + Otherwise, dissolve the polygons with convex hull. + """ + unified = [Polygon(poly.exterior) for poly in multipolygon.geoms] + merged_polygon = unary_union(unified) + + if merged_polygon.geom_type == "MultiPolygon": + polygons = [ + Polygon(poly.exterior) for poly in merged_polygon.geoms if poly.is_valid + ] + union_poly = unary_union(polygons) + + if union_poly.geom_type == "Polygon": + return union_poly + + if union_poly.geom_type == "MultiPolygon" and dissolve_polygon: + # disjoint polygons + return union_poly.convex_hull + + return merged_polygon + + +def ensure_right_hand_rule(polygon: Polygon): + """Check if a polygon follows the right-hand rule, fix it if not.""" + if polygon.exterior.is_ccw: # If counter-clockwise, reverse it + return Polygon( + polygon.exterior.coords[::-1], + [interior.coords[::-1] for interior in polygon.interiors], + ) + return polygon + + def merge_polygons( featcol: geojson.FeatureCollection, - dissolve_polygon: bool = True, + dissolve_polygon: bool = False, ) -> geojson.FeatureCollection: """Merge multiple Polygons or MultiPolygons into a single Polygon. + It is used to create a single polygon boundary. + Args: featcol: a FeatureCollection containing geometries. dissolve_polygon: True to dissolve polygons to single polygon. @@ -754,29 +795,35 @@ def merge_polygons( geojson.FeatureCollection: a FeatureCollection of a single Polygon. """ geom_list = [] + properties = {} try: features = featcol.get("features", []) for feature in features: + properties = feature["properties"] polygon = shape(feature["geometry"]) - geom_list.append(polygon) - - merged_polygon = unary_union(geom_list) - merged_geojson = mapping(merged_polygon) - - # MultiPolygons are stripped out earlier - if dissolve_polygon: - merged_polygon = merged_polygon.convex_hull - merged_geojson = mapping(merged_polygon) - log.warning( - "Resulted GeoJSON contains disjoint Polygons. " - "Adjacent polygons are preferred." - ) - return geojson.FeatureCollection([geojson.Feature(geometry=merged_geojson)]) + if isinstance(polygon, MultiPolygon): + # Remove holes in each polygon + polygons_without_holes = [remove_holes(poly) for poly in polygon.geoms] + valid_polygons = [ + ensure_right_hand_rule(poly) for poly in polygons_without_holes + ] + else: + polygon_without_holes = remove_holes(polygon) + valid_polygons = [ensure_right_hand_rule(polygon_without_holes)] + geom_list.extend(valid_polygons) + + merged_geom = create_single_polygon(MultiPolygon(geom_list), dissolve_polygon) + merged_geojson = mapping(merged_geom) + + # Create FeatureCollection + return geojson.FeatureCollection( + [geojson.Feature(geometry=merged_geojson, properties=properties)] + ) except Exception as e: raise HTTPException( - HTTPStatus.BAD_REQUEST, + HTTPStatus.UNPROCESSABLE_ENTITY, detail=f"Couldn't merge the multipolygon to polygon: {str(e)}", ) from e diff --git a/src/backend/app/projects/project_routes.py b/src/backend/app/projects/project_routes.py index d4b7ef7243..ecfe348778 100644 --- a/src/backend/app/projects/project_routes.py +++ b/src/backend/app/projects/project_routes.py @@ -614,7 +614,7 @@ async def preview_split_by_square( else: log.warning("Parsed geojson file contained no geometries") - if len(boundary_featcol["features"]) == 0: + if len(boundary_featcol["features"]) > 0: boundary_featcol = merge_polygons(boundary_featcol) return split_by_square( @@ -639,6 +639,7 @@ async def get_data_extract( TODO alternatively, direct to raw-data-api to generate first, then upload """ boundary_geojson = json.loads(await geojson_file.read()) + clean_boundary_geojson = merge_polygons(boundary_geojson) # Get extract config file from existing data_models if form_category: @@ -650,7 +651,7 @@ async def get_data_extract( extract_config = None fgb_url = await project_crud.generate_data_extract( - boundary_geojson, + clean_boundary_geojson, extract_config, ) diff --git a/src/backend/app/projects/project_schemas.py b/src/backend/app/projects/project_schemas.py index 9eba15b46e..4b565330d7 100644 --- a/src/backend/app/projects/project_schemas.py +++ b/src/backend/app/projects/project_schemas.py @@ -136,8 +136,8 @@ def parse_input_geojson( # geojson_pydantic.GeometryCollection # FIXME update this to remove the Featcol parsing at some point featcol = geojson_to_featcol(value) - merged = merge_polygons(featcol) - return merged.get("features")[0].get("geometry") + merged_geojson = merge_polygons(featcol, True) + return merged_geojson.get("features")[0].get("geometry") @model_validator(mode="after") def append_fmtm_hashtag_and_slug(self) -> Self: