from __future__ import annotations import uuid from datetime import datetime, timezone from hashlib import sha256 from pathlib import Path from typing import Any from geoalchemy2.shape import to_shape from shapely.geometry import mapping from shapely.ops import transform as shapely_transform from shapely.validation import make_valid from app.core.errors import AppError from app.models import Area, Dataset, DatasetVersion from app.services.derived_dataset_governance_service import DerivedDatasetGovernanceService from app.services.raster_service import extract_raster_metadata from app.services.storage_service import StorageService from app.services.tile_manifest_service import TileManifestService, canonical_manifest_json def _import_rasterio(): import importlib rasterio = importlib.import_module("rasterio") errors = importlib.import_module("rasterio.errors") return rasterio, errors def _import_numpy(): import importlib return importlib.import_module("numpy") def _import_pillow(): import importlib return importlib.import_module("PIL") class RasterOperationsService: RASTER_UNAVAILABLE_MESSAGE = ( "Raster processing unavailable. Install rasterio and GDAL-compatible drivers to enable raster processing operations." ) RASTER_STATS_UNAVAILABLE_MESSAGE = ( "Raster statistics unavailable. Install rasterio and numpy to enable raster band statistics." ) RASTER_INDEX_UNAVAILABLE_MESSAGE = ( "Raster processing unavailable. Install rasterio and numpy to enable raster index operations." ) PREVIEW_UNAVAILABLE_MESSAGE = "Raster preview unavailable. Install rasterio, numpy and pillow to enable preview generation." DEFAULT_REPROJECT_CRS = "EPSG:31370" DEFAULT_STATS_HISTOGRAM_BINS = 16 @staticmethod def _require_raster_dataset(dataset: Dataset) -> None: if dataset.dataset_type != "raster": raise AppError(code="INVALID_DATASET_TYPE", message="Dataset is not a raster dataset", status_code=400) if not dataset.storage_path: raise AppError(code="DATASET_FILE_MISSING", message="Stored raster file is missing", status_code=404) @staticmethod def _load_dataset(db, dataset_id: uuid.UUID) -> Dataset: dataset = db.get(Dataset, dataset_id) if not dataset: raise AppError(code="DATASET_NOT_FOUND", message="Dataset not found", status_code=404) RasterOperationsService._require_raster_dataset(dataset) source_path = Path(dataset.storage_path) if not source_path.exists(): raise AppError(code="DATASET_FILE_MISSING", message="Stored raster file missing", status_code=404) return dataset @staticmethod def _validate_storage_checksum(dataset: Dataset) -> None: """Refuse tiling pixels that no longer match the governed Dataset row.""" expected_checksum = str(dataset.checksum_sha256 or "").strip().lower() governed_artifact = bool(getattr(dataset, "data_contract_key", None)) valid_checksum = len(expected_checksum) == 64 and all(character in "0123456789abcdef" for character in expected_checksum) if not valid_checksum: if expected_checksum or governed_artifact: raise AppError( code="DATASET_STORAGE_CHECKSUM_UNVERIFIABLE", message="Governed raster storage requires a valid SHA-256 checksum before tiling.", details={"checksum_sha256": dataset.checksum_sha256}, status_code=409, ) return digest = sha256() with Path(str(dataset.storage_path)).open("rb") as stream: for chunk in iter(lambda: stream.read(8 * 1024 * 1024), b""): digest.update(chunk) actual_checksum = digest.hexdigest() if actual_checksum != expected_checksum: raise AppError( code="DATASET_STORAGE_CHECKSUM_MISMATCH", message="Raster dataset storage no longer matches its validated checksum.", details={ "expected_checksum_sha256": expected_checksum, "actual_checksum_sha256": actual_checksum, }, status_code=409, ) @staticmethod def _raster_dependencies() -> tuple[Any, Any]: try: return _import_rasterio() except Exception as exc: raise AppError( code="RASTER_PROCESSING_UNAVAILABLE", message=RasterOperationsService.RASTER_UNAVAILABLE_MESSAGE, status_code=503, ) from exc @staticmethod def _stats_dependencies() -> tuple[Any, Any]: rasterio, _ = RasterOperationsService._raster_dependencies() try: numpy = _import_numpy() except Exception as exc: raise AppError( code="RASTER_PROCESSING_UNAVAILABLE", message=RasterOperationsService.RASTER_STATS_UNAVAILABLE_MESSAGE, status_code=503, ) from exc return rasterio, numpy @staticmethod def _index_dependencies() -> tuple[Any, Any]: rasterio, _ = RasterOperationsService._raster_dependencies() try: numpy = _import_numpy() except Exception as exc: raise AppError( code="RASTER_PROCESSING_UNAVAILABLE", message=RasterOperationsService.RASTER_INDEX_UNAVAILABLE_MESSAGE, status_code=503, ) from exc return rasterio, numpy @staticmethod def _validate_positive_band_index(value: int, label: str) -> int: if not isinstance(value, int): raise AppError(code="INVALID_PARAMETERS", message=f"{label} must be a positive integer", status_code=400) if value <= 0: raise AppError(code="INVALID_PARAMETERS", message=f"{label} must be greater than 0", status_code=400) return value @staticmethod def _normalize_nodata(value: Any) -> float | int | None: if value is None: return None if isinstance(value, (list, tuple)): if not value: return None value = value[0] if value == "": return None try: return float(value) except (TypeError, ValueError): return None @staticmethod def _normalize_nodata_for_band(nodata: Any, band_index: int) -> float | int | None: if isinstance(nodata, (list, tuple)): if band_index <= 0 or band_index > len(nodata): return None return RasterOperationsService._normalize_nodata(nodata[band_index - 1]) return RasterOperationsService._normalize_nodata(nodata) @staticmethod def _validate_tile_request(tile_size: int, overlap: int) -> None: if tile_size <= 0: raise AppError(code="INVALID_PARAMETERS", message="tile_size must be greater than 0", status_code=400) if overlap < 0: raise AppError(code="INVALID_PARAMETERS", message="overlap must be greater or equal to 0", status_code=400) if overlap >= tile_size: raise AppError(code="INVALID_PARAMETERS", message="overlap must be smaller than tile_size", status_code=400) @staticmethod def _dataset_metadata(dataset_id: uuid.UUID, storage: dict[str, Any], extra: dict[str, Any] | None = None) -> dict[str, Any]: metadata = { "dataset_id": str(dataset_id), "size_bytes": storage.get("size_bytes"), "checksum_sha256": storage.get("checksum_sha256"), "path": storage.get("storage_path"), } if extra: metadata.update(extra) return metadata @staticmethod def _validate_band_mapping(dataset: Dataset, source_band_count: int, mapping: dict[str, int]) -> dict[str, int]: if source_band_count <= 0: raise AppError(code="INVALID_DATASET", message="Source raster has no bands", status_code=400) validated: dict[str, int] = {} for name, value in mapping.items(): band_index = RasterOperationsService._validate_positive_band_index(value, name) if band_index > source_band_count: raise AppError( code="INVALID_PARAMETERS", message=f"{name} exceeds available band count ({band_index} > {source_band_count})", status_code=400, ) validated[name] = band_index return validated @staticmethod def _coerce_rasterio_crs(rasterio: Any, value: str | None) -> Any: if not value: raise ValueError("CRS value is missing") crs_namespace = getattr(rasterio, "crs", rasterio) crs_class = getattr(crs_namespace, "CRS", crs_namespace) if hasattr(crs_class, "from_user_input"): return crs_class.from_user_input(value) raise AttributeError("rasterio CRS converter unavailable") @staticmethod def _preview_dimensions(source_width: int, source_height: int, max_dimension: int = 2048) -> tuple[int, int]: width = max(1, int(source_width)) height = max(1, int(source_height)) preview_width = min(width, max_dimension) preview_height = int(height * (preview_width / width)) if preview_height <= 0: preview_height = 1 if preview_height > max_dimension: preview_height = max_dimension preview_width = int(width * (preview_height / height)) if preview_width <= 0: preview_width = 1 return preview_width, preview_height @staticmethod def _window_bounds_to_list(bounds: Any) -> list[float]: if isinstance(bounds, (list, tuple)) and len(bounds) == 4: left, bottom, right, top = bounds return [float(left), float(bottom), float(right), float(top)] return [float(bounds.left), float(bounds.bottom), float(bounds.right), float(bounds.top)] @staticmethod def _normalize_preview_data(data: Any) -> Any: try: if isinstance(data, (list, tuple)) and len(data) > 0: return data[0] except Exception: pass return data @staticmethod def _write_preview_image( data: Any, preview_path: Path, preview_width: int | None = None, preview_height: int | None = None, ) -> tuple[int, int]: _import_pillow() numpy = _import_numpy() image_data = numpy.asarray(RasterOperationsService._normalize_preview_data(data)) if image_data.size == 0: raise AppError(code="RASTER_PREVIEW_ERROR", message="Cannot generate preview for empty raster", status_code=422) if image_data.ndim > 2: image_data = image_data[0] if image_data.ndim != 2: raise AppError(code="RASTER_PREVIEW_ERROR", message="Cannot generate preview for raster shape", status_code=500) valid = numpy.isfinite(image_data) if valid.any(): valid_values = image_data.astype("float64")[valid] minimum = float(valid_values.min()) maximum = float(valid_values.max()) scale = maximum - minimum if scale == 0: scale = 1.0 normalized = ((image_data.astype("float64") - minimum) / scale * 255).clip(0, 255) normalized = normalized.astype("uint8") else: normalized = numpy.zeros(image_data.shape, dtype="uint8") from PIL import Image image = Image.fromarray(normalized, mode="L") if preview_width is not None and preview_height is not None and ( preview_width != image.width or preview_height != image.height ): image = image.resize( (int(preview_width), int(preview_height)), resample=getattr(Image.Resampling, "LANCZOS", Image.BICUBIC), ) preview_path.parent.mkdir(parents=True, exist_ok=True) image.save(preview_path) return int(image.width), int(image.height) @staticmethod def _transform_area_area_geometry(area_geom, area: Area, source_crs_str: str) -> Any: if not area.original_crs: raise AppError( code="INVALID_CRS", message="Area CRS is required to align clipping geometry with raster CRS.", status_code=400, ) if area.original_crs == source_crs_str: return area_geom try: import pyproj except Exception as exc: raise AppError(code="INVALID_CRS", message="pyproj is required to reproject clip area", status_code=400) from exc try: transformer = pyproj.Transformer.from_crs(area.original_crs, source_crs_str, always_xy=True) return shapely_transform(transformer.transform, area_geom) except Exception as exc: raise AppError(code="INVALID_CRS", message="Unable to align area CRS to raster CRS", status_code=400) from exc @staticmethod def _persist_derived_dataset( db, source_dataset: Dataset, source_dataset_id: uuid.UUID, operation: str, output_path: str, output_name: str, metadata: dict[str, Any], ) -> uuid.UUID: derived_id = uuid.uuid4() output_file = Path(output_path) if output_file.suffix.lower() not in {".tif", ".tiff", ".geotiff"}: output_file = output_file.with_suffix(".tif") storage_metadata: dict[str, Any] = {} if output_file.exists(): storage_metadata = { "size_bytes": output_file.stat().st_size, "checksum_sha256": StorageService.calculate_checksum_sha256(output_file.read_bytes()), } storage_metadata.update( { "original_filename": output_file.name, "stored_filename": output_file.name, "content_type": "image/tiff", "storage_path": str(output_file), }, ) metadata_payload = dict(metadata or {}) operation_name = operation if operation.startswith("raster.") else f"raster.{operation}" provenance = { "operation": operation_name, "source_dataset_id": str(source_dataset_id), "input_dataset_id": str(source_dataset_id), "operation_parameters": metadata_payload.get("operation_parameters", {}), } metadata_payload.setdefault("operation", operation_name) metadata_payload.update(provenance) metadata_payload.setdefault("output_dataset_id", str(derived_id)) derived_dataset = Dataset( id=derived_id, project_id=source_dataset.project_id, area_id=source_dataset.area_id, name=output_name, dataset_type="raster", source=f"operation:{operation_name}", dataset_role="derived", # The source identity of this artifact is the server-owned # derived-operation registry entry. The parent remains explicit # in provenance and the lineage edge below. source_name="derived", source_metadata=source_dataset.source_metadata, provenance_metadata=provenance, imported_at=datetime.now(timezone.utc), temporal_series_key=( f"{source_dataset.temporal_series_key}:{operation_name}" if source_dataset.temporal_series_key else None ), observed_at=source_dataset.observed_at, valid_from=source_dataset.valid_from, valid_to=source_dataset.valid_to, temporal_granularity=source_dataset.temporal_granularity, source_version=source_dataset.source_version, storage_path=str(output_file), original_filename=storage_metadata["original_filename"], stored_filename=storage_metadata["stored_filename"], content_type=storage_metadata["content_type"], size_bytes=storage_metadata.get("size_bytes"), checksum_sha256=storage_metadata.get("checksum_sha256"), derived_from_dataset_id=source_dataset_id, crs=metadata_payload.get("crs"), bounds_json=metadata_payload.get("bounds"), resolution_json=metadata_payload.get("resolution"), bands_json={"dtype": metadata_payload.get("dtype")} if metadata_payload.get("dtype") is not None else None, metadata_json=metadata_payload, status="validating", ) db.add(derived_dataset) dataset_version = DatasetVersion( dataset_id=derived_dataset.id, version=1, storage_path=derived_dataset.storage_path, source_version=derived_dataset.source_version, observed_at=derived_dataset.observed_at, valid_from=derived_dataset.valid_from, valid_to=derived_dataset.valid_to, checksum_sha256=derived_dataset.checksum_sha256, source_metadata=derived_dataset.source_metadata, provenance_metadata=derived_dataset.provenance_metadata, ) db.add(dataset_version) DerivedDatasetGovernanceService.govern_raster( db, dataset=derived_dataset, dataset_version=dataset_version, raster_metadata=metadata_payload, source_key="derived", operation=operation_name, parent_dataset=source_dataset, operation_parameters=metadata_payload.get("operation_parameters"), ) db.commit() db.refresh(derived_dataset) return derived_id @staticmethod def metadata(db, dataset_id: uuid.UUID) -> dict[str, Any]: dataset = RasterOperationsService._load_dataset(db, dataset_id) metadata = extract_raster_metadata(dataset.storage_path) metadata["dataset_id"] = str(dataset.id) metadata["size_bytes"] = dataset.size_bytes metadata["checksum_sha256"] = dataset.checksum_sha256 metadata["path"] = dataset.storage_path return metadata @staticmethod def inspect(db, dataset_id: uuid.UUID) -> dict[str, Any]: dataset = RasterOperationsService._load_dataset(db, dataset_id) profile = RasterOperationsService.metadata(db, dataset_id) return { "dataset_id": str(dataset.id), "ready": True, "metadata": profile, "operation": "raster.inspect", "output_dataset_id": None, "source_dataset_id": None, } @staticmethod def preview(db, dataset_id: uuid.UUID) -> dict[str, Any]: dataset = RasterOperationsService._load_dataset(db, dataset_id) rasterio, _ = RasterOperationsService._raster_dependencies() preview_dir = StorageService.preview_root(str(dataset.project_id), str(dataset.id)) preview_dir.mkdir(parents=True, exist_ok=True) preview_path = preview_dir / "preview.png" try: with rasterio.open(dataset.storage_path) as source: width = int(source.width) height = int(source.height) preview_width, preview_height = RasterOperationsService._preview_dimensions(width, height) if not preview_path.exists(): data = source.read(1) try: preview_width, preview_height = RasterOperationsService._write_preview_image( data=data, preview_path=preview_path, preview_width=preview_width, preview_height=preview_height, ) except TypeError: preview_width, preview_height = RasterOperationsService._write_preview_image(data, preview_path) else: try: from PIL import Image with Image.open(preview_path) as cached: preview_width = int(cached.width) preview_height = int(cached.height) except Exception: # best effort fallback; keep computed dimensions. pass except AppError: raise except Exception as exc: # pragma: no cover if isinstance(exc, AppError): raise raise AppError(code="RASTER_PREVIEW_ERROR", message="Unable to generate raster preview", status_code=500) from exc metadata = RasterOperationsService._dataset_metadata( dataset.id, { "storage_path": dataset.storage_path, "size_bytes": dataset.size_bytes, "checksum_sha256": dataset.checksum_sha256, }, extra=extract_raster_metadata(dataset.storage_path), ) return { "dataset_id": str(dataset.id), "ready": True, "preview": { "path": str(preview_path), "format": "PNG", "width": preview_width, "height": preview_height, }, "metadata": metadata, "operation": "raster.preview", "source_dataset_id": str(dataset.id), } @staticmethod def _compute_spectral_index( db, dataset_id: uuid.UUID, mapping: dict[str, int], operation_name: str, formula: str, output_name: str, subtraction_order: str = "second_minus_first", ) -> uuid.UUID: dataset = RasterOperationsService._load_dataset(db, dataset_id) rasterio, numpy = RasterOperationsService._index_dependencies() output_id = uuid.uuid4() output_filename = f"{output_name or operation_name}.tif" output_path = Path(StorageService.derived_raster_root(str(dataset.project_id), str(output_id)) / output_filename) output_path.parent.mkdir(parents=True, exist_ok=True) with rasterio.open(dataset.storage_path) as source: source_band_count = int(source.count) validated_mapping = RasterOperationsService._validate_band_mapping( dataset=dataset, source_band_count=source_band_count, mapping={str(key): int(value) for key, value in mapping.items()}, ) first_key = [key for key in ("red_band", "green_band", "swir_band") if key in validated_mapping][0] first_band = validated_mapping[first_key] second_band = validated_mapping["nir_band"] source_profile = source.profile.copy() source_profile.update( { "count": 1, "dtype": "float32", "nodata": float("nan"), }, ) block_size = max(1, min(1024, int(source.width), int(source.height))) with rasterio.open(output_path, "w", **source_profile) as destination: for yoff in range(0, int(source.height), block_size): row_count = min(block_size, int(source.height) - yoff) for xoff in range(0, int(source.width), block_size): column_count = min(block_size, int(source.width) - xoff) window = rasterio.windows.Window(xoff, yoff, column_count, row_count) first_data = numpy.asarray( source.read(first_band, window=window, out_dtype="float32"), ).astype("float32") second_data = numpy.asarray( source.read(second_band, window=window, out_dtype="float32"), ).astype("float32") nodata = source.nodata first_nodata = RasterOperationsService._normalize_nodata_for_band(nodata, first_band) second_nodata = RasterOperationsService._normalize_nodata_for_band(nodata, second_band) valid = numpy.isfinite(first_data) & numpy.isfinite(second_data) if first_nodata is not None: valid &= first_data != first_nodata if second_nodata is not None: valid &= second_data != second_nodata denominator = first_data + second_data computed = numpy.full_like(first_data, float("nan"), dtype="float32") if not numpy.all(~valid): np_valid = valid.astype(bool) if np_valid.any(): safe_denominator = denominator.copy() safe_denominator[~np_valid] = 1.0 with numpy.errstate(divide="ignore", invalid="ignore", over="ignore", under="ignore"): if subtraction_order == "first_minus_second": difference = first_data - second_data else: difference = second_data - first_data computed_values = difference / safe_denominator computed[~np_valid] = float("nan") computed[np_valid] = numpy.where( (first_data[np_valid] + second_data[np_valid]) == 0.0, float("nan"), computed_values[np_valid], ) destination.write(computed, indexes=1, window=window) output_metadata = extract_raster_metadata(str(output_path)) output_metadata["operation"] = f"raster.{operation_name}" output_metadata["source_dataset_id"] = str(dataset.id) output_metadata["operation_parameters"] = { **validated_mapping, "formula": formula, "nodata_strategy": "nan", "source_band_count": source_band_count, } output_metadata["band_mapping"] = validated_mapping output_metadata["formula"] = formula output_metadata["output_dtype"] = "float32" output_metadata["nodata_strategy"] = { "mode": "nan", "value_range_note": "Expected index range is approximately [-1, 1] before optional clipping.", } output_metadata["created_at"] = datetime.now(timezone.utc).isoformat() output_metadata["path"] = str(output_path) output_metadata["output_dataset_id"] = str(output_id) derived_id = RasterOperationsService._persist_derived_dataset( db=db, source_dataset=dataset, source_dataset_id=dataset.id, operation=f"raster.{operation_name}", output_path=str(output_path), output_name=output_filename, metadata=output_metadata, ) return derived_id @staticmethod def ndvi(db, dataset_id: uuid.UUID, nir_band: int, red_band: int, output_name: str | None = None) -> uuid.UUID: return RasterOperationsService._compute_spectral_index( db=db, dataset_id=dataset_id, mapping={"nir_band": nir_band, "red_band": red_band}, operation_name="ndvi", formula="(nir - red) / (nir + red)", output_name=(output_name or "ndvi"), ) @staticmethod def ndwi(db, dataset_id: uuid.UUID, green_band: int, nir_band: int, output_name: str | None = None) -> uuid.UUID: return RasterOperationsService._compute_spectral_index( db=db, dataset_id=dataset_id, mapping={"green_band": green_band, "nir_band": nir_band}, operation_name="ndwi", formula="(nir - green) / (nir + green)", output_name=(output_name or "ndwi"), ) @staticmethod def ndbi(db, dataset_id: uuid.UUID, swir_band: int, nir_band: int, output_name: str | None = None) -> uuid.UUID: return RasterOperationsService._compute_spectral_index( db=db, dataset_id=dataset_id, mapping={"swir_band": swir_band, "nir_band": nir_band}, operation_name="ndbi", formula="(swir - nir) / (swir + nir)", output_name=(output_name or "ndbi"), subtraction_order="first_minus_second", ) @staticmethod def stats(db, dataset_id: uuid.UUID) -> dict[str, Any]: dataset = RasterOperationsService._load_dataset(db, dataset_id) rasterio, numpy = RasterOperationsService._stats_dependencies() with rasterio.open(dataset.storage_path) as source: height = int(source.height) width = int(source.width) count = int(source.count) dataset_profile = extract_raster_metadata(dataset.storage_path) metadata = { "dataset_id": str(dataset.id), "source_dataset_id": str(dataset.id), "size_bytes": dataset.size_bytes, "checksum_sha256": dataset.checksum_sha256, "profile": dataset_profile, } bands = [] chunk_rows = max(1, min(2048, height)) for band_index in range(1, count + 1): nodata = RasterOperationsService._normalize_nodata_for_band(source.nodata, band_index) dtype = str(source.dtypes[band_index - 1]) if source.dtypes else None band_min = None band_max = None valid_count = 0 total_sum = 0.0 total_sq = 0.0 nodata_count = 0 hist = None hist_bins = None for row_offset in range(0, height, chunk_rows): row_count = min(chunk_rows, height - row_offset) data = source.read(band_index, window=rasterio.windows.Window(0, row_offset, width, row_count)) values = numpy.asarray(data) if values.size == 0: continue finite = numpy.isfinite(values) if nodata is not None: valid = finite & (values != nodata) nodata_count += int(values.size - valid.sum()) else: valid = finite band_values = values[valid].astype("float64") if band_values.size == 0: continue current_min = float(band_values.min()) current_max = float(band_values.max()) if band_min is None or current_min < band_min: band_min = current_min if band_max is None or current_max > band_max: band_max = current_max valid_count += int(band_values.size) total_sum += float(band_values.sum()) total_sq += float((band_values**2).sum()) if hist is None: hist, hist_bins = numpy.histogram(band_values, bins=RasterOperationsService.DEFAULT_STATS_HISTOGRAM_BINS) else: additional, _ = numpy.histogram(band_values, bins=hist_bins) hist = hist + additional total_pixels = width * height if valid_count == 0: bands.append( { "band_index": band_index, "dtype": dtype, "min": None, "max": None, "mean": None, "std": None, "nodata_count": int(nodata_count), "nodata_ratio": 1.0 if total_pixels else 0.0, "valid_pixel_count": 0, "histogram": None, "histogram_bins": None, }, ) continue mean = total_sum / valid_count variance = max(0.0, (total_sq / valid_count) - (mean * mean)) std = float(numpy.sqrt(variance)) bands.append( { "band_index": band_index, "dtype": dtype, "min": float(band_min) if band_min is not None else None, "max": float(band_max) if band_max is not None else None, "mean": float(mean), "std": float(std), "nodata_count": int(nodata_count), "nodata_ratio": float(nodata_count) / max(1, total_pixels), "valid_pixel_count": int(valid_count), "histogram": hist.astype(int).tolist() if hist is not None else None, "histogram_bins": [float(item) for item in hist_bins] if hist_bins is not None else None, }, ) return { "dataset_id": str(dataset.id), "source_dataset_id": str(dataset.id), "bands": bands, "generated_at": datetime.now(timezone.utc).isoformat(), "metadata": metadata, } @staticmethod def reproject( db, dataset_id: uuid.UUID, target_crs: str | None, output_name: str | None, resampling: str = "nearest", ) -> uuid.UUID: dataset = RasterOperationsService._load_dataset(db, dataset_id) target_crs = target_crs or RasterOperationsService.DEFAULT_REPROJECT_CRS rasterio, _ = RasterOperationsService._raster_dependencies() try: target = RasterOperationsService._coerce_rasterio_crs(rasterio, target_crs) except Exception as exc: raise AppError(code="INVALID_PARAMETERS", message="Invalid target CRS", status_code=400) from exc if not hasattr(target, "to_string"): raise AppError(code="INVALID_CRS", message="Invalid target CRS", status_code=400) resampling_map = { "nearest": getattr(rasterio.enums.Resampling, "nearest", None), "bilinear": getattr(rasterio.enums.Resampling, "bilinear", None), "cubic": getattr(rasterio.enums.Resampling, "cubic", None), } selected_resampling = resampling_map.get(resampling or "nearest") if selected_resampling is None: raise AppError(code="INVALID_PARAMETERS", message="Unsupported resampling method", status_code=400) output_name = (output_name or "raster_reprojected").strip() or "raster_reprojected" output_id = uuid.uuid4() output_filename = f"{output_name}.tif" output_path = Path(StorageService.derived_raster_root(str(dataset.project_id), str(output_id)) / output_filename) output_path.parent.mkdir(parents=True, exist_ok=True) with rasterio.open(dataset.storage_path) as source: if not source.crs: raise AppError(code="INVALID_DATASET_CRS", message="Source raster CRS is missing", status_code=400) source_transform = source.transform source_crs = source.crs output_kwargs = source.meta.copy() source_bounds = getattr(source, "bounds", None) try: if source_bounds is not None: source_bounds_tuple = ( source_bounds.left, source_bounds.bottom, source_bounds.right, source_bounds.top, ) else: raise AttributeError except Exception: source_bounds_tuple = ( 0.0, 0.0, float(source.width), float(source.height), ) transform, width, height = rasterio.warp.calculate_default_transform( source_crs, target, source.width, source.height, *source_bounds_tuple, ) output_kwargs.update( { "crs": target, "transform": transform, "width": int(width), "height": int(height), "count": source.count, }, ) with rasterio.open(output_path, "w", **output_kwargs) as destination: for band_index in range(1, source.count + 1): source_band_reader = rasterio.band destination_band_reader = rasterio.band if hasattr(source_band_reader, "__self__"): source_band_reader = getattr(rasterio.__class__, "band", source_band_reader) if hasattr(destination_band_reader, "__self__"): destination_band_reader = getattr(rasterio.__class__, "band", destination_band_reader) source_band = source_band_reader(source, band_index) destination_band = destination_band_reader(destination, band_index) rasterio.warp.reproject( source=source_band, destination=destination_band, src_transform=source_transform, src_crs=source_crs, dst_transform=transform, dst_crs=target, resampling=selected_resampling, ) output_metadata = extract_raster_metadata(str(output_path)) output_metadata["operation"] = "raster.reproject" output_metadata["source_dataset_id"] = str(dataset.id) output_metadata["operation_parameters"] = { "target_crs": target_crs, "resampling": resampling, } output_metadata["target_crs"] = target_crs output_metadata["output_dataset_id"] = str(output_id) derived_id = RasterOperationsService._persist_derived_dataset( db=db, source_dataset=dataset, source_dataset_id=dataset.id, operation="raster.reproject", output_path=str(output_path), output_name=output_filename, metadata=output_metadata, ) return derived_id @staticmethod def clip(db, dataset_id: uuid.UUID, area_id: uuid.UUID, output_name: str | None) -> uuid.UUID: dataset = RasterOperationsService._load_dataset(db, dataset_id) rasterio, _ = RasterOperationsService._raster_dependencies() 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 != dataset.project_id: raise AppError(code="INVALID_DATASET_SCOPE", message="Area does not belong to dataset project", status_code=400) if not area.geometry: raise AppError(code="INVALID_AREA_GEOMETRY", message="Area geometry is missing", status_code=400) area_geom = to_shape(area.geometry) if not area_geom.is_valid: area_geom = make_valid(area_geom) if not area_geom.is_valid: raise AppError(code="INVALID_GEOMETRY", message="Area geometry cannot be repaired", status_code=400) if area_geom.is_empty: raise AppError(code="INVALID_AREA_GEOMETRY", message="Area geometry is empty", status_code=400) with rasterio.open(dataset.storage_path) as source: raw_source_crs = source.crs source_crs = raw_source_crs.to_string() if hasattr(raw_source_crs, "to_string") else ( str(raw_source_crs) if raw_source_crs else None ) if not source_crs: raise AppError(code="INVALID_DATASET_CRS", message="Source raster CRS is missing", status_code=400) transformed_area = RasterOperationsService._transform_area_area_geometry(area_geom, area, source_crs) if not transformed_area.is_valid: transformed_area = make_valid(transformed_area) if not transformed_area.is_valid: raise AppError(code="INVALID_AREA_GEOMETRY", message="Area geometry cannot be repaired", status_code=400) mask_input = [mapping(transformed_area)] try: clipped_data, clipped_transform = rasterio.mask.mask(source, mask_input, crop=True, nodata=source.nodata, filled=True) except Exception as exc: raise AppError(code="RASTER_OPERATION_ERROR", message="Raster clipping failed", status_code=500) from exc clipped_has_data = True try: import numpy clipped_array = numpy.asarray(clipped_data) clipped_has_data = bool(clipped_array.size and numpy.isfinite(clipped_array).any()) except Exception: clipped_has_data = clipped_data.size > 0 if not clipped_has_data: raise AppError(code="RASTER_OPERATION_EMPTY_RESULT", message="Raster clip produced no output data", status_code=422) source_count = getattr(source, "count", None) if not source_count: if hasattr(clipped_data, "shape") and len(clipped_data.shape) >= 1: source_count = int(clipped_data.shape[0]) else: source_count = 1 source_count = int(source_count) profile = source.profile.copy() profile.update( { "count": source_count, "height": int(clipped_data.shape[1]), "width": int(clipped_data.shape[2]), "transform": clipped_transform, }, ) output_id = uuid.uuid4() output_filename = f"{output_name or 'raster_clipped'}.tif" output_path = Path(StorageService.derived_raster_root(str(dataset.project_id), str(output_id)) / output_filename) output_path.parent.mkdir(parents=True, exist_ok=True) with rasterio.open(output_path, "w", **profile) as destination: destination.write(clipped_data) if not output_path.exists(): raise AppError(code="RASTER_OPERATION_ERROR", message="Failed to write clip output", status_code=500) derived_metadata = extract_raster_metadata(str(output_path)) derived_metadata["operation"] = "raster.clip" derived_metadata["source_dataset_id"] = str(dataset.id) derived_metadata["operation_parameters"] = { "area_id": str(area_id), "source_crs": source_crs, "area_crs": area.original_crs, } derived_metadata["output_dataset_id"] = str(output_id) derived_id = RasterOperationsService._persist_derived_dataset( db=db, source_dataset=dataset, source_dataset_id=dataset.id, operation="raster.clip", output_path=str(output_path), output_name=output_filename, metadata=derived_metadata, ) return derived_id @staticmethod def tile( db, dataset_id: uuid.UUID, tile_size: int = 512, overlap: int = 64, output_name: str | None = None, max_tiles: int | None = None, ) -> dict[str, Any]: RasterOperationsService._validate_tile_request(tile_size=tile_size, overlap=overlap) if max_tiles is not None and max_tiles <= 0: raise AppError(code="INVALID_PARAMETERS", message="max_tiles must be positive", status_code=400) dataset = RasterOperationsService._load_dataset(db, dataset_id) RasterOperationsService._validate_storage_checksum(dataset) rasterio, _ = RasterOperationsService._raster_dependencies() tile_set_id = str(uuid.uuid4()) tile_root = StorageService.raster_tiles_root(str(dataset.project_id), str(dataset.id), tile_set_id) manifest_tiles: list[dict[str, Any]] = [] tile_paths: list[str] = [] source_crs: str | None = None with rasterio.open(dataset.storage_path) as source: raw_source_crs = getattr(source, "crs", None) source_crs = raw_source_crs.to_string() if hasattr(raw_source_crs, "to_string") else ( str(raw_source_crs) if raw_source_crs else dataset.crs ) source_count = getattr(source, "count", 0) if not source_count: source_count = 1 if source_count == 0: raise AppError(code="INVALID_DATASET", message="Dataset has no raster bands", status_code=400) source_width = int(source.width) source_height = int(source.height) step = max(1, tile_size - overlap) x_offsets = RasterOperationsService._tile_offsets(source_width, tile_size, step) y_offsets = RasterOperationsService._tile_offsets(source_height, tile_size, step) expected_tile_count = len(x_offsets) * len(y_offsets) if max_tiles is not None and expected_tile_count > max_tiles: raise AppError( code="RASTER_TILE_LIMIT_EXCEEDED", message="Raster tile generation exceeds the guest analysis limit.", details={ "expected_tile_count": expected_tile_count, "max_tiles": max_tiles, }, status_code=422, ) tile_root.mkdir(parents=True, exist_ok=True) tile_index = 0 for yoff in y_offsets: for xoff in x_offsets: tile_width = min(tile_size, source_width - xoff) tile_height = min(tile_size, source_height - yoff) if tile_width <= 0 or tile_height <= 0: continue window = rasterio.windows.Window(xoff, yoff, tile_width, tile_height) tile_data = source.read(window=window) if tile_data.size == 0: continue bounds = rasterio.windows.bounds(window, source.transform) transform = rasterio.windows.transform(window, source.transform) tile_path = tile_root / f"tile_{tile_index:04d}.tif" profile = source.profile.copy() profile.update(width=int(tile_width), height=int(tile_height), transform=transform) profile.pop("transform", None) profile["transform"] = transform with rasterio.open(tile_path, "w", **profile) as tile_dest: tile_dest.write(tile_data) tile_paths.append(str(tile_path)) tile_integrity = TileManifestService.tile_integrity(tile_path) manifest_tiles.append( { "path": str(tile_path), "pixel_window": [int(xoff), int(yoff), int(tile_width), int(tile_height)], "bounds": RasterOperationsService._window_bounds_to_list(bounds), "transform": [float(item) for item in transform.to_gdal()], "crs": source_crs, "index": tile_index, **tile_integrity, }, ) tile_index += 1 if not manifest_tiles: raise AppError(code="RASTER_OPERATION_EMPTY_RESULT", message="Raster tile generation produced no tiles", status_code=422) try: source_metadata = extract_raster_metadata(dataset.storage_path) except AppError: source_metadata = {"bounds": [0.0, 0.0, 0.0, 0.0]} bounds = source_metadata.get("bounds", [0.0, 0.0, 0.0, 0.0]) manifest_crs = source_crs or source_metadata.get("crs") or dataset.crs manifest_payload = { **TileManifestService.dataset_binding(db, dataset), "tile_set_id": tile_set_id, "crs": manifest_crs, "source_crs": manifest_crs, "dataset_crs": dataset.crs, "bounds": [float(value) for value in bounds], "tile_size": int(tile_size), "overlap": int(overlap), "parameters": { "tile_size": int(tile_size), "overlap": int(overlap), "output_name": output_name, }, "created_at": datetime.now(timezone.utc).isoformat(), "tile_paths": tile_paths, "count": len(manifest_tiles), "tiles": manifest_tiles, "ai_inference": False, "tile_server": None, } manifest_path = tile_root / "manifest.json" manifest_path.write_text(canonical_manifest_json(manifest_payload), encoding="utf-8") return { "dataset_id": str(dataset.id), "ready": True, "operation": "raster.tile", "tile_set_id": tile_set_id, "tile_size": tile_size, "overlap": overlap, "manifest_path": str(manifest_path), "count": len(manifest_tiles), "manifest": manifest_payload, } @staticmethod def _tile_offsets(dimension: int, tile_size: int, step: int) -> list[int]: """Return full-tile starts plus one unique edge-aligned final start.""" if dimension <= 0 or tile_size <= 0 or step <= 0: raise AppError(code="INVALID_PARAMETERS", message="Raster tile dimensions must be positive", status_code=400) if dimension <= tile_size: return [0] final_start = dimension - tile_size offsets = list(range(0, final_start + 1, step)) if offsets[-1] != final_start: offsets.append(final_start) return offsets