330 lines
15 KiB
Python
330 lines
15 KiB
Python
from __future__ import annotations
|
|
|
|
import json
|
|
import uuid
|
|
from pathlib import Path
|
|
from typing import Any
|
|
|
|
from geoalchemy2.shape import to_shape
|
|
from shapely.geometry import GeometryCollection, MultiPolygon, shape
|
|
from shapely.geometry.base import BaseGeometry
|
|
from shapely.geometry import mapping
|
|
from shapely.ops import unary_union
|
|
from shapely.validation import make_valid
|
|
from sqlalchemy.orm import Session
|
|
|
|
from app.core.errors import AppError
|
|
from app.models import Area, Dataset
|
|
from app.schemas.operations import VectorOperationResult
|
|
from app.services.geojson_service import parse_geojson_payload
|
|
from app.services.storage_service import StorageService
|
|
|
|
|
|
class VectorOperationsService:
|
|
@staticmethod
|
|
def _require_vector_dataset(dataset: Dataset) -> None:
|
|
if dataset.dataset_type not in {"vector", "geojson"}:
|
|
raise AppError(code="INVALID_DATASET_TYPE", message="Dataset is not a vector dataset", status_code=400)
|
|
|
|
@staticmethod
|
|
def _load_dataset_payload(dataset: Dataset) -> tuple[dict[str, Any], list[dict[str, Any]]]:
|
|
if not dataset.storage_path:
|
|
raise AppError(code="DATASET_FILE_MISSING", message="Stored file missing", status_code=404)
|
|
path = Path(dataset.storage_path)
|
|
if not path.exists():
|
|
raise AppError(code="DATASET_FILE_MISSING", message="Stored file missing", status_code=404)
|
|
try:
|
|
payload = json.loads(path.read_text(encoding="utf-8"))
|
|
except Exception as exc:
|
|
raise AppError(code="INVALID_GEOJSON", message="Stored dataset is not valid JSON", status_code=400) from exc
|
|
|
|
if not isinstance(payload, dict) or payload.get("type") != "FeatureCollection":
|
|
raise AppError(code="INVALID_GEOJSON", message="Dataset payload is not a FeatureCollection", status_code=400)
|
|
|
|
features = payload.get("features")
|
|
if not isinstance(features, list):
|
|
raise AppError(code="INVALID_GEOJSON", message="Dataset payload is missing features", status_code=400)
|
|
return payload, [feature for feature in features if isinstance(feature, dict)]
|
|
|
|
@staticmethod
|
|
def _extract_geometries(features: list[dict[str, Any]]) -> list[tuple[dict[str, Any], BaseGeometry]]:
|
|
geometries: list[tuple[dict[str, Any], BaseGeometry]] = []
|
|
for feature in features:
|
|
if not isinstance(feature, dict):
|
|
continue
|
|
geometry = feature.get("geometry")
|
|
if not geometry:
|
|
continue
|
|
try:
|
|
shapely_geom = shape(geometry)
|
|
except Exception as exc:
|
|
raise AppError(code="INVALID_GEOMETRY", message="Feature geometry invalid", status_code=400) from exc
|
|
if not shapely_geom.is_valid:
|
|
shapely_geom = make_valid(shapely_geom)
|
|
if not shapely_geom.is_valid:
|
|
raise AppError(code="INVALID_GEOMETRY", message="Feature geometry cannot be repaired", status_code=400)
|
|
|
|
geometries.append((feature, shapely_geom))
|
|
|
|
if not geometries:
|
|
raise AppError(code="EMPTY_VECTOR_DATASET", message="Dataset has no parseable geometries", status_code=422)
|
|
return geometries
|
|
|
|
@staticmethod
|
|
def inspect(db: Session, dataset_id: uuid.UUID) -> VectorOperationResult:
|
|
dataset = db.get(Dataset, dataset_id)
|
|
if not dataset:
|
|
raise AppError(code="DATASET_NOT_FOUND", message="Dataset not found", status_code=404)
|
|
VectorOperationsService._require_vector_dataset(dataset)
|
|
|
|
payload, features = VectorOperationsService._load_dataset_payload(dataset)
|
|
geometries = VectorOperationsService._extract_geometries(features)
|
|
|
|
geometry_type_summary: dict[str, int] = {}
|
|
for _, geometry in geometries:
|
|
geometry_type_summary[geometry.geom_type] = geometry_type_summary.get(geometry.geom_type, 0) + 1
|
|
|
|
unioned = unary_union([geometry for _, geometry in geometries])
|
|
bounds = unioned.bounds
|
|
return VectorOperationResult(
|
|
source_dataset_id=str(dataset_id),
|
|
feature_count=len(geometries),
|
|
geometry_type_summary=geometry_type_summary,
|
|
bounds_json={"min_x": float(bounds[0]), "min_y": float(bounds[1]), "max_x": float(bounds[2]), "max_y": float(bounds[3])},
|
|
crs=payload.get("crs") if isinstance(payload.get("crs"), str) else dataset.crs,
|
|
)
|
|
|
|
@staticmethod
|
|
def bbox(db: Session, dataset_id: uuid.UUID) -> dict[str, Any]:
|
|
summary = VectorOperationsService.inspect(db, dataset_id)
|
|
return {
|
|
"dataset_id": str(dataset_id),
|
|
"bounds_json": summary.bounds_json,
|
|
"feature_count": summary.feature_count,
|
|
"crs": summary.crs,
|
|
}
|
|
|
|
@staticmethod
|
|
def stats(db: Session, dataset_id: uuid.UUID) -> dict[str, Any]:
|
|
summary = VectorOperationsService.inspect(db, dataset_id)
|
|
return {
|
|
"dataset_id": str(dataset_id),
|
|
"feature_count": summary.feature_count,
|
|
"geometry_type_summary": summary.geometry_type_summary,
|
|
"bounds_json": summary.bounds_json,
|
|
"crs": summary.crs,
|
|
}
|
|
|
|
@staticmethod
|
|
def clip_by_area(db: Session, dataset_id: uuid.UUID, area_id: uuid.UUID, output_name: str | None) -> uuid.UUID:
|
|
source_dataset = db.get(Dataset, dataset_id)
|
|
if not source_dataset:
|
|
raise AppError(code="DATASET_NOT_FOUND", message="Dataset not found", status_code=404)
|
|
VectorOperationsService._require_vector_dataset(source_dataset)
|
|
|
|
area = db.get(Area, area_id)
|
|
if not area:
|
|
raise AppError(code="AREA_NOT_FOUND", message="Area not found", status_code=404)
|
|
if area.project_id != source_dataset.project_id:
|
|
raise AppError(code="INVALID_DATASET_SCOPE", message="Area does not belong to dataset project", status_code=400)
|
|
|
|
payload, features = VectorOperationsService._load_dataset_payload(source_dataset)
|
|
geometries = VectorOperationsService._extract_geometries(features)
|
|
area_geom = to_shape(area.geometry)
|
|
if area_geom.is_empty:
|
|
raise AppError(code="INVALID_AREA_GEOMETRY", message="Area geometry is empty", status_code=400)
|
|
|
|
if isinstance(area_geom, GeometryCollection):
|
|
area_geom = unary_union(area_geom.geoms)
|
|
if area_geom.geom_type == "MultiPolygon":
|
|
area_geom = MultiPolygon(area_geom.geoms)
|
|
|
|
if not area_geom.is_valid:
|
|
area_geom = make_valid(area_geom)
|
|
if not area_geom.is_valid:
|
|
raise AppError(code="INVALID_AREA_GEOMETRY", message="Area geometry cannot be repaired", status_code=400)
|
|
|
|
output_features: list[dict[str, Any]] = []
|
|
for feature, source_geom in geometries:
|
|
clipped = source_geom.intersection(area_geom)
|
|
if clipped.is_empty:
|
|
continue
|
|
if not clipped.is_valid:
|
|
clipped = make_valid(clipped)
|
|
if not clipped.is_valid:
|
|
raise AppError(code="INVALID_GEOMETRY", message="Clipped geometry became invalid", status_code=400)
|
|
output_features.append({
|
|
"type": "Feature",
|
|
"geometry": mapping(clipped),
|
|
"properties": feature.get("properties", {}) or {},
|
|
})
|
|
|
|
if not output_features:
|
|
raise AppError(code="VECTOR_OPERATION_EMPTY_RESULT", message="Clip operation produced no output features", status_code=422)
|
|
|
|
return VectorOperationsService._persist_derived_dataset(
|
|
db=db,
|
|
source_dataset=source_dataset,
|
|
source_id=dataset_id,
|
|
operation="clip",
|
|
feature_collection={"type": "FeatureCollection", "features": output_features},
|
|
output_name=output_name,
|
|
default_name="vector_clipped",
|
|
)
|
|
|
|
@staticmethod
|
|
def buffer(db: Session, dataset_id: uuid.UUID, distance_m: float, dissolve: bool, output_name: str | None) -> uuid.UUID:
|
|
source_dataset = db.get(Dataset, dataset_id)
|
|
if not source_dataset:
|
|
raise AppError(code="DATASET_NOT_FOUND", message="Dataset not found", status_code=404)
|
|
VectorOperationsService._require_vector_dataset(source_dataset)
|
|
if distance_m <= 0:
|
|
raise AppError(code="INVALID_PARAMETERS", message="distance_m must be greater than 0", status_code=400)
|
|
|
|
_, features = VectorOperationsService._load_dataset_payload(source_dataset)
|
|
geometries = VectorOperationsService._extract_geometries(features)
|
|
buffered_features = [(feature, geometry.buffer(distance_m)) for feature, geometry in geometries]
|
|
|
|
output_features: list[dict[str, Any]] = []
|
|
for feature, geometry in buffered_features:
|
|
if geometry.is_empty:
|
|
continue
|
|
if not geometry.is_valid:
|
|
geometry = make_valid(geometry)
|
|
if not geometry.is_valid:
|
|
raise AppError(code="INVALID_GEOMETRY", message="Buffer geometry became invalid", status_code=400)
|
|
output_features.append({
|
|
"type": "Feature",
|
|
"geometry": mapping(geometry),
|
|
"properties": feature.get("properties", {}) or {},
|
|
})
|
|
|
|
if dissolve:
|
|
dissolved = unary_union([shape(feature["geometry"]) for feature in output_features])
|
|
output_features = [{
|
|
"type": "Feature",
|
|
"geometry": mapping(dissolved),
|
|
"properties": {"operation": "vector_buffer", "distance_m": distance_m, "dissolve": True},
|
|
}]
|
|
|
|
if not output_features:
|
|
raise AppError(code="VECTOR_OPERATION_EMPTY_RESULT", message="Buffer operation produced no output features", status_code=422)
|
|
|
|
return VectorOperationsService._persist_derived_dataset(
|
|
db=db,
|
|
source_dataset=source_dataset,
|
|
source_id=dataset_id,
|
|
operation="buffer",
|
|
feature_collection={"type": "FeatureCollection", "features": output_features},
|
|
output_name=output_name,
|
|
default_name="vector_buffered",
|
|
)
|
|
|
|
@staticmethod
|
|
def intersect(
|
|
db: Session,
|
|
source_dataset_id: uuid.UUID,
|
|
target_dataset_id: uuid.UUID,
|
|
output_name: str | None,
|
|
) -> uuid.UUID:
|
|
if source_dataset_id == target_dataset_id:
|
|
raise AppError(code="INVALID_PARAMETERS", message="other_dataset_id must be different from source dataset", status_code=400)
|
|
|
|
source_dataset = db.get(Dataset, source_dataset_id)
|
|
if not source_dataset:
|
|
raise AppError(code="DATASET_NOT_FOUND", message="Source dataset not found", status_code=404)
|
|
VectorOperationsService._require_vector_dataset(source_dataset)
|
|
|
|
target_dataset = db.get(Dataset, target_dataset_id)
|
|
if not target_dataset:
|
|
raise AppError(code="DATASET_NOT_FOUND", message="Target dataset not found", status_code=404)
|
|
VectorOperationsService._require_vector_dataset(target_dataset)
|
|
if target_dataset.project_id != source_dataset.project_id:
|
|
raise AppError(code="INVALID_DATASET_SCOPE", message="Datasets must belong to same project", status_code=400)
|
|
|
|
source_payload, source_features = VectorOperationsService._load_dataset_payload(source_dataset)
|
|
target_payload, _ = VectorOperationsService._load_dataset_payload(target_dataset)
|
|
source_geometries = VectorOperationsService._extract_geometries(source_features)
|
|
target_geometries = VectorOperationsService._extract_geometries(target_payload.get("features", []))
|
|
target_union = unary_union([geometry for _, geometry in target_geometries])
|
|
|
|
output_features: list[dict[str, Any]] = []
|
|
for source_feature, source_geometry in source_geometries:
|
|
intersection = source_geometry.intersection(target_union)
|
|
if intersection.is_empty:
|
|
continue
|
|
if not intersection.is_valid:
|
|
intersection = make_valid(intersection)
|
|
if not intersection.is_valid:
|
|
raise AppError(code="INVALID_GEOMETRY", message="Intersection geometry became invalid", status_code=400)
|
|
output_features.append({
|
|
"type": "Feature",
|
|
"geometry": mapping(intersection),
|
|
"properties": source_feature.get("properties", {}) or {},
|
|
})
|
|
|
|
if not output_features:
|
|
raise AppError(code="VECTOR_OPERATION_EMPTY_RESULT", message="Intersection operation produced no output features", status_code=422)
|
|
|
|
return VectorOperationsService._persist_derived_dataset(
|
|
db=db,
|
|
source_dataset=source_dataset,
|
|
source_id=source_dataset_id,
|
|
operation="intersect",
|
|
feature_collection={"type": "FeatureCollection", "features": output_features},
|
|
output_name=output_name,
|
|
default_name="vector_intersect",
|
|
)
|
|
|
|
@staticmethod
|
|
def _persist_derived_dataset(
|
|
db: Session,
|
|
source_dataset: Dataset,
|
|
source_id: uuid.UUID,
|
|
operation: str,
|
|
feature_collection: dict[str, Any],
|
|
output_name: str | None,
|
|
default_name: str,
|
|
) -> uuid.UUID:
|
|
derived_id = uuid.uuid4()
|
|
output_name_value = f"{(output_name or default_name)}.geojson"
|
|
if not output_name_value.strip():
|
|
output_name_value = f"{default_name}.geojson"
|
|
|
|
stored = json.dumps(feature_collection, ensure_ascii=False, separators=(",", ":")).encode("utf-8")
|
|
storage_info = StorageService.persist_dataset_file(
|
|
project_id=str(source_dataset.project_id),
|
|
dataset_id=str(derived_id),
|
|
dataset_type="vector",
|
|
original_filename=output_name_value,
|
|
content=stored,
|
|
content_type="application/geo+json",
|
|
)
|
|
|
|
metadata = parse_geojson_payload(json.dumps(feature_collection, ensure_ascii=False, separators=(",", ":")))
|
|
derived_dataset = Dataset(
|
|
id=derived_id,
|
|
project_id=source_dataset.project_id,
|
|
area_id=source_dataset.area_id,
|
|
name=output_name_value,
|
|
dataset_type="vector",
|
|
source=f"operation:{operation}",
|
|
storage_path=storage_info["storage_path"],
|
|
original_filename=storage_info["original_filename"],
|
|
stored_filename=storage_info["stored_filename"],
|
|
content_type=storage_info["content_type"],
|
|
size_bytes=storage_info["size_bytes"],
|
|
checksum_sha256=storage_info["checksum_sha256"],
|
|
derived_from_dataset_id=source_id,
|
|
crs=metadata.get("crs"),
|
|
bounds_json=metadata.get("bounds_json"),
|
|
resolution_json=metadata.get("resolution_json"),
|
|
bands_json=metadata.get("bands_json"),
|
|
metadata_json=metadata,
|
|
status="ready",
|
|
)
|
|
db.add(derived_dataset)
|
|
db.commit()
|
|
db.refresh(derived_dataset)
|
|
return derived_id
|