1069 lines
46 KiB
Python
1069 lines
46 KiB
Python
from __future__ import annotations
|
|
|
|
import json
|
|
import uuid
|
|
from datetime import datetime, timezone
|
|
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.raster_service import extract_raster_metadata
|
|
from app.services.storage_service import StorageService
|
|
|
|
|
|
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 _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",
|
|
source_name=source_dataset.source_name,
|
|
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="ready",
|
|
)
|
|
db.add(derived_dataset)
|
|
db.add(
|
|
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.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,
|
|
) -> dict[str, Any]:
|
|
RasterOperationsService._validate_tile_request(tile_size=tile_size, overlap=overlap)
|
|
dataset = RasterOperationsService._load_dataset(db, dataset_id)
|
|
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)
|
|
tile_root.mkdir(parents=True, exist_ok=True)
|
|
|
|
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)
|
|
tile_index = 0
|
|
for yoff in range(0, source_height, step):
|
|
for xoff in range(0, source_width, step):
|
|
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))
|
|
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_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 = {
|
|
"tile_set_id": tile_set_id,
|
|
"source_dataset_id": str(dataset.id),
|
|
"source_raster_id": str(dataset.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(json.dumps(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,
|
|
}
|