Merge pull request #83 from hotosm/feat/clipping-no-fly-zones

Feat/clipping no fly zones

src/backend/app/db/db_models.py

@@ -89,7 +89,7 @@ class DbProject(Base):
     # GEOMETRY
     outline = cast(WKBElement, Column(Geometry("POLYGON", srid=4326)))
     centroid = cast(WKBElement, Column(Geometry("POINT", srid=4326)))
-    no_fly_zones = cast(WKBElement, Column(Geometry("POLYGON", srid=4326)))
+    no_fly_zones = cast(WKBElement, Column(Geometry("MULTIPOLYGON", srid=4326)))
 
     organisation_id = cast(
         int,

src/backend/app/migrations/versions/9d01411fd221_.py

@@ -0,0 +1,64 @@
+"""
+
+Revision ID: 9d01411fd221
+Revises: acee47666167
+Create Date: 2024-07-20 13:35:29.634665
+
+"""
+from typing import Sequence, Union
+
+from alembic import op
+import geoalchemy2
+
+# revision identifiers, used by Alembic.
+revision: str = "9d01411fd221"
+down_revision: Union[str, None] = "acee47666167"
+branch_labels: Union[str, Sequence[str], None] = None
+depends_on: Union[str, Sequence[str], None] = None
+
+
+def upgrade() -> None:
+    # ### commands auto generated by Alembic - please adjust! ###
+
+    op.alter_column(
+        "projects",
+        "no_fly_zones",
+        existing_type=geoalchemy2.types.Geometry(
+            geometry_type="POLYGON",
+            srid=4326,
+            from_text="ST_GeomFromEWKT",
+            name="geometry",
+            _spatial_index_reflected=True,
+        ),
+        type_=geoalchemy2.types.Geometry(
+            geometry_type="MULTIPOLYGON",
+            srid=4326,
+            from_text="ST_GeomFromEWKT",
+            name="geometry",
+        ),
+        existing_nullable=True,
+    )
+    # ### end Alembic commands ###
+
+
+def downgrade() -> None:
+    # ### commands auto generated by Alembic - please adjust! ###
+    op.alter_column(
+        "projects",
+        "no_fly_zones",
+        existing_type=geoalchemy2.types.Geometry(
+            geometry_type="MULTIPOLYGON",
+            srid=4326,
+            from_text="ST_GeomFromEWKT",
+            name="geometry",
+        ),
+        type_=geoalchemy2.types.Geometry(
+            geometry_type="POLYGON",
+            srid=4326,
+            from_text="ST_GeomFromEWKT",
+            name="geometry",
+            _spatial_index_reflected=True,
+        ),
+        existing_nullable=True,
+    )
+    # ### end Alembic commands ###

src/backend/app/projects/project_routes.py

@@ -5,11 +5,9 @@
 from app.users.user_schemas import AuthUser
 import geojson
 from datetime import timedelta
-
 from fastapi import APIRouter, HTTPException, Depends, UploadFile, File, Form
 from sqlalchemy.orm import Session
 from loguru import logger as log
-
 from app.projects import project_schemas, project_crud
 from app.db import database
 from app.models.enums import HTTPStatus
@@ -18,6 +16,8 @@
 from app.config import settings
 from databases import Database
 from app.db import db_models
+from shapely.geometry import shape, mapping
+from shapely.ops import unary_union
 
 router = APIRouter(
     prefix=f"{settings.API_PREFIX}/projects",
@@ -125,6 +125,7 @@ async def upload_project_task_boundaries(
 @router.post("/preview-split-by-square/", tags=["Projects"])
 async def preview_split_by_square(
     project_geojson: UploadFile = File(...),
+    no_fly_zones: UploadFile = File(default=None),
     dimension: int = Form(100),
     user: AuthUser = Depends(login_required),
 ):
@@ -140,8 +141,25 @@ async def preview_split_by_square(
     # read entire file
     content = await project_geojson.read()
     boundary = geojson.loads(content)
+    project_shape = shape(boundary["features"][0]["geometry"])
+
+    # If no_fly_zones is provided, read and parse it
+    if no_fly_zones:
+        no_fly_content = await no_fly_zones.read()
+        no_fly_zones_geojson = geojson.loads(no_fly_content)
+        no_fly_shapes = [
+            shape(feature["geometry"]) for feature in no_fly_zones_geojson["features"]
+        ]
+        no_fly_union = unary_union(no_fly_shapes)
+
+        # Calculate the difference between the project shape and no-fly zones
+        new_outline = project_shape.difference(no_fly_union)
+    else:
+        new_outline = project_shape
+    result_geojson = geojson.Feature(geometry=mapping(new_outline))
+
+    result = await project_crud.preview_split_by_square(result_geojson, dimension)
 
-    result = await project_crud.preview_split_by_square(boundary, dimension)
     return result
 
 

src/backend/app/utils.py

@@ -8,7 +8,7 @@
 from geojson_pydantic import FeatureCollection as FeatCol
 from geoalchemy2 import WKBElement
 from geoalchemy2.shape import from_shape, to_shape
-from shapely.geometry import mapping, shape
+from shapely.geometry import mapping, shape, MultiPolygon as ShapelyMultiPolygon
 from shapely.ops import unary_union
 from fastapi import HTTPException
 from shapely import wkb
@@ -75,13 +75,15 @@ def geojson_to_geometry(
     features = parsed_geojson.get("features", [])
 
     if len(features) > 1:
-        # TODO code to merge all geoms into multipolygon
-        # TODO do not use convex hull
-        pass
-
-    geometry = features[0].get("geometry")
+        geometries = [shape(feature.get("geometry")) for feature in features]
+        merged_geometry = unary_union(geometries)
+        if not isinstance(merged_geometry, ShapelyMultiPolygon):
+            merged_geometry = ShapelyMultiPolygon([merged_geometry])
+        shapely_geom = merged_geometry
+    else:
+        geometry = features[0].get("geometry")
 
-    shapely_geom = shape(geometry)
+        shapely_geom = shape(geometry)
 
     return from_shape(shapely_geom)