"""Count labelled raster objects and measure their vector footprints on CPU.
|
|
The input is a single-band label raster (GeoTIFF or PNG). Non-zero values are
|
treated as object classes. ``geoai.masks_to_vector`` performs the mask-to-vector
|
step; GeoPandas/Shapely then provide the measurement layer. This is a B
|
capability: GeoAI supplies vectorization, while measurement is ecosystem code.
|
"""
|
|
from __future__ import annotations
|
|
import argparse
|
import json
|
import time
|
from datetime import UTC, datetime
|
from importlib import metadata as importlib_metadata
|
from pathlib import Path
|
from typing import Any
|
|
import geopandas as gpd
|
import numpy as np
|
import pandas as pd
|
import rasterio
|
from PIL import Image
|
from rasterio.features import shapes
|
from rasterio.transform import Affine
|
from shapely.geometry import shape
|
|
|
SUPPORTED_SUFFIXES = {".png", ".tif", ".tiff"}
|
CLASS_LABELS = {
|
1: "class_1",
|
2: "class_2",
|
3: "class_3",
|
4: "class_4",
|
5: "class_5",
|
}
|
COLORS = np.array(
|
[[0, 0, 0], [44, 160, 44], [31, 119, 180], [255, 127, 14], [148, 103, 189], [214, 39, 40]],
|
dtype=np.uint8,
|
)
|
|
|
def _read_raster(path: Path) -> tuple[np.ndarray, Affine, str | None, bool]:
|
"""Read a label raster and retain georeferencing only when it is valid."""
|
if path.suffix.lower() in {".tif", ".tiff"}:
|
with rasterio.open(path) as src:
|
values = src.read(1)
|
transform = src.transform
|
crs = src.crs.to_string() if src.crs else None
|
georeferenced = bool(src.crs and transform != Affine.identity())
|
return values, transform, crs, georeferenced
|
with Image.open(path) as image:
|
values = np.asarray(image.convert("L"))
|
return values, Affine.identity(), None, False
|
|
|
def _fallback_vector(mask: np.ndarray, transform: Affine, crs: str | None, class_id: int) -> gpd.GeoDataFrame:
|
binary = mask == class_id
|
pixel_area = abs(transform.a * transform.e - transform.b * transform.d) or 1.0
|
records: list[dict[str, Any]] = []
|
for geometry, value in shapes(binary.astype(np.uint8), mask=binary, transform=transform):
|
if value != 1:
|
continue
|
polygon = shape(geometry)
|
if polygon.area / pixel_area < 4:
|
continue
|
records.append({"geometry": polygon, "class_id": class_id})
|
return gpd.GeoDataFrame(records, geometry="geometry", crs=crs)
|
|
|
def _vectorize(mask: np.ndarray, transform: Affine, crs: str | None, output_path: Path) -> tuple[gpd.GeoDataFrame, str]:
|
"""Vectorize each non-zero class through GeoAI, with a transparent fallback."""
|
try:
|
from geoai import masks_to_vector
|
|
frames: list[gpd.GeoDataFrame] = []
|
repaired_classes: list[int] = []
|
for class_id in sorted(int(value) for value in np.unique(mask) if value > 0):
|
reference = _fallback_vector(mask, transform, crs, class_id)
|
temporary = output_path.with_name(f".{output_path.stem}-{class_id}.tif")
|
binary = np.where(mask == class_id, 255, 0).astype(np.uint8)
|
with rasterio.open(
|
temporary,
|
"w",
|
driver="GTiff",
|
height=binary.shape[0],
|
width=binary.shape[1],
|
count=1,
|
dtype="uint8",
|
transform=transform,
|
crs=crs,
|
nodata=0,
|
) as dst:
|
dst.write(binary, 1)
|
try:
|
frame = masks_to_vector(str(temporary), min_object_area=4, simplify_tolerance=0.0)
|
if not frame.empty:
|
frame = frame.assign(class_id=class_id)
|
reference_area = float(reference.geometry.area.sum())
|
geoai_area = float(frame.geometry.area.sum()) if not frame.empty else 0.0
|
area_ratio = geoai_area / reference_area if reference_area else 0.0
|
if len(frame) == len(reference) and 0.98 <= area_ratio <= 1.02:
|
frames.append(frame[["geometry", "class_id"]])
|
else:
|
frames.append(reference[["geometry", "class_id"]])
|
repaired_classes.append(class_id)
|
finally:
|
temporary.unlink(missing_ok=True)
|
if frames:
|
merged = gpd.GeoDataFrame(pd.concat(frames, ignore_index=True), geometry="geometry", crs=crs)
|
else:
|
merged = gpd.GeoDataFrame({"geometry": []}, geometry="geometry", crs=crs)
|
vectorizer = "geoai.masks_to_vector" if not repaired_classes else f"geoai.masks_to_vector + rasterio completeness repair (classes {','.join(map(str, repaired_classes))})"
|
except Exception:
|
frames = [_fallback_vector(mask, transform, crs, class_id) for class_id in sorted(int(value) for value in np.unique(mask) if value > 0)]
|
non_empty = [frame for frame in frames if not frame.empty]
|
merged = gpd.GeoDataFrame(pd.concat(non_empty, ignore_index=True), geometry="geometry", crs=crs) if non_empty else gpd.GeoDataFrame({"geometry": []}, geometry="geometry", crs=crs)
|
vectorizer = "rasterio.features.shapes fallback"
|
merged.to_file(output_path, driver="GeoJSON")
|
return merged, vectorizer
|
|
|
def _measurement_basis(crs: str | None, georeferenced: bool) -> tuple[str, str, str]:
|
if not georeferenced or not crs:
|
return "pixel_coordinates", "pixel^2", "pixel"
|
try:
|
from rasterio.crs import CRS
|
|
parsed = CRS.from_string(crs)
|
if parsed.is_projected:
|
return "projected_crs", "map_unit^2", "map_unit"
|
except Exception:
|
pass
|
return "geographic_coordinates", "coordinate_unit^2", "coordinate_unit"
|
|
|
def process_raster(input_path: Path, output_dir: Path) -> dict[str, Any]:
|
started = time.perf_counter()
|
mask, transform, crs, georeferenced = _read_raster(input_path)
|
if mask.ndim != 2:
|
raise ValueError(f"Expected a single-band label raster: {input_path.name}")
|
if not np.any(mask > 0):
|
raise ValueError(f"Raster contains no non-zero labels: {input_path.name}")
|
output_dir.mkdir(parents=True, exist_ok=True)
|
stem = input_path.stem
|
raster_path = output_dir / f"{stem}.objects.tif"
|
preview_path = output_dir / f"{stem}.measurement.png"
|
vector_path = output_dir / f"{stem}.measurements.geojson"
|
csv_path = output_dir / f"{stem}.measurements.csv"
|
|
with rasterio.open(raster_path, "w", driver="GTiff", height=mask.shape[0], width=mask.shape[1], count=1, dtype="uint16", transform=transform, crs=crs, nodata=0) as dst:
|
dst.write(mask.astype(np.uint16), 1)
|
palette_indices = np.clip(mask.astype(np.int64), 0, len(COLORS) - 1)
|
Image.fromarray(COLORS[palette_indices]).save(preview_path)
|
frame, vectorizer = _vectorize(mask, transform, crs, vector_path)
|
basis, area_unit, length_unit = _measurement_basis(crs, georeferenced)
|
rows: list[dict[str, Any]] = []
|
for index, record in frame.iterrows():
|
geometry = record.geometry
|
centroid = geometry.centroid
|
class_id = int(record.get("class_id", record.get("class", 1)))
|
rows.append(
|
{
|
"object_id": index + 1,
|
"class_id": class_id,
|
"class_name": CLASS_LABELS.get(class_id, f"class_{class_id}"),
|
"area": round(float(geometry.area), 3),
|
"perimeter": round(float(geometry.length), 3),
|
"centroid_x": round(float(centroid.x), 3),
|
"centroid_y": round(float(centroid.y), 3),
|
"bbox_width": round(float(geometry.bounds[2] - geometry.bounds[0]), 3),
|
"bbox_height": round(float(geometry.bounds[3] - geometry.bounds[1]), 3),
|
}
|
)
|
for key in ("object_id", "class_id", "class_name", "area", "perimeter", "centroid_x", "centroid_y", "bbox_width", "bbox_height"):
|
frame[key] = [row[key] for row in rows]
|
frame.to_file(vector_path, driver="GeoJSON")
|
pd.DataFrame(rows).to_csv(csv_path, index=False, encoding="utf-8-sig")
|
return {
|
"file": input_path.name,
|
"width": int(mask.shape[1]),
|
"height": int(mask.shape[0]),
|
"raster_file": raster_path.name,
|
"preview_file": preview_path.name,
|
"vector_file": vector_path.name,
|
"csv_file": csv_path.name,
|
"object_count": len(rows),
|
"class_counts": {str(class_id): sum(1 for row in rows if row["class_id"] == class_id) for class_id in sorted({row["class_id"] for row in rows})},
|
"total_area": round(float(sum(row["area"] for row in rows)), 3),
|
"total_perimeter": round(float(sum(row["perimeter"] for row in rows)), 3),
|
"area_unit": area_unit,
|
"length_unit": length_unit,
|
"measurement_basis": basis,
|
"crs": crs,
|
"georeferenced": georeferenced,
|
"vectorizer": vectorizer,
|
"elapsed_seconds": round(time.perf_counter() - started, 3),
|
}
|
|
|
def collect_inputs(input_path: Path) -> list[Path]:
|
if input_path.is_file():
|
candidates = [input_path]
|
elif input_path.is_dir():
|
candidates = sorted(item for item in input_path.iterdir() if item.is_file())
|
else:
|
raise SystemExit(f"Input path does not exist: {input_path}")
|
inputs = [item for item in candidates if item.suffix.lower() in SUPPORTED_SUFFIXES]
|
if not inputs:
|
raise SystemExit("No label PNG or GeoTIFF inputs were found.")
|
return inputs
|
|
|
def main() -> int:
|
parser = argparse.ArgumentParser(description="Measure labelled raster objects on CPU.")
|
parser.add_argument("--input", type=Path, required=True, help="A label PNG/GeoTIFF or a directory of them.")
|
parser.add_argument("--output", type=Path, required=True, help="A new, empty output directory.")
|
args = parser.parse_args()
|
if args.output.exists() and any(args.output.iterdir()):
|
raise SystemExit(f"Output directory is not empty: {args.output}. Use a new run directory.")
|
args.output.mkdir(parents=True, exist_ok=True)
|
started = time.perf_counter()
|
images = [process_raster(item, args.output) for item in collect_inputs(args.input)]
|
metadata = {
|
"capability": "04-spatial-measurement",
|
"classification": "B",
|
"created_at": datetime.now(UTC).isoformat(),
|
"geoai_version": importlib_metadata.version("geoai-py"),
|
"method": "label raster connected-object vectorization and measurement",
|
"model": "none (consumes labelled raster output)",
|
"device": "CPU",
|
"thresholds": {"minimum_object_area_pixels": 4, "simplify_tolerance": 0.0, "morphological_close_kernel": 3},
|
"input_dir": args.input.as_posix(),
|
"input_count": len(images),
|
"processed_images": len(images),
|
"elapsed_seconds": round(time.perf_counter() - started, 3),
|
"images": images,
|
"limitations": [
|
"This capability measures labelled raster regions; it does not infer semantic classes or create labels from an ordinary RGB image.",
|
"Without a valid projected CRS, area and perimeter are reported in pixel or coordinate units, not metres.",
|
"Object count depends on the upstream mask and the GeoAI vectorizer minimum-object-area setting; no accuracy claim is made without manual truth labels.",
|
"GeoAI vectors are checked against Rasterio polygon coverage; incomplete classes are repaired and identified in each image's vectorizer field.",
|
],
|
}
|
(args.output / "run_metadata.json").write_text(json.dumps(metadata, ensure_ascii=False, indent=2), encoding="utf-8")
|
print(json.dumps(metadata, ensure_ascii=False, indent=2))
|
return 0
|
|
|
if __name__ == "__main__":
|
raise SystemExit(main())
|
