hoàn thành chức năng change detection
This commit is contained in:
Victor Phan
+773
-9
@@ -3,7 +3,7 @@ API Server for Land Classification Model Training
|
||||
Cho phép chọn dữ liệu và cấu hình training qua giao diện web
|
||||
"""
|
||||
|
||||
from fastapi import FastAPI, BackgroundTasks, HTTPException
|
||||
from fastapi import FastAPI, BackgroundTasks, HTTPException, UploadFile, File
|
||||
from fastapi.middleware.cors import CORSMiddleware
|
||||
from fastapi.staticfiles import StaticFiles
|
||||
from fastapi.responses import HTMLResponse, FileResponse
|
||||
@@ -15,10 +15,27 @@ import json
|
||||
from datetime import datetime
|
||||
from pathlib import Path
|
||||
import sys
|
||||
import numpy as np
|
||||
import xarray as xr
|
||||
import rasterio
|
||||
from rasterio.transform import from_bounds
|
||||
import asyncio
|
||||
import hashlib
|
||||
import traceback
|
||||
|
||||
# Import report generator
|
||||
from report_generator import generate_training_report, generate_prediction_report
|
||||
|
||||
# Import planetary computer libraries (conditional)
|
||||
try:
|
||||
from pystac_client import Client
|
||||
import planetary_computer
|
||||
import odc.stac
|
||||
except ImportError:
|
||||
Client = None
|
||||
planetary_computer = None
|
||||
odc = None
|
||||
|
||||
app = FastAPI(title="Land Classification Training API", version="1.0.0")
|
||||
|
||||
# Enable CORS
|
||||
@@ -131,6 +148,28 @@ class NDVIConfig(BaseModel):
|
||||
resolution: int = 20
|
||||
|
||||
|
||||
class ChangeDetectionWorkflowRequest(BaseModel):
|
||||
"""Request for change detection workflow"""
|
||||
prediction_result: dict
|
||||
bbox: List[float]
|
||||
|
||||
|
||||
class ComparePeriodsPredictionConfig(BaseModel):
|
||||
"""Compare predictions between two time periods"""
|
||||
model_filename: str
|
||||
min_lon: float
|
||||
min_lat: float
|
||||
max_lon: float
|
||||
max_lat: float
|
||||
current_period: dict # {start_date, end_date}
|
||||
prediction_period: dict # {start_date, end_date}
|
||||
max_scenes: int = 12
|
||||
cloud_cover: int = 30
|
||||
resolution: int = 20
|
||||
export_ndvi: bool = True
|
||||
export_classification: bool = True
|
||||
|
||||
|
||||
class PredictionWithNDVIConfig(BaseModel):
|
||||
"""Cấu hình predict kết hợp land classification và NDVI"""
|
||||
model_filename: str
|
||||
@@ -147,6 +186,16 @@ class PredictionWithNDVIConfig(BaseModel):
|
||||
export_classification: bool = True # Export classification raster
|
||||
|
||||
|
||||
# Serve change detection interface page (moved here after app is defined)
|
||||
@app.get("/change-detection", response_class=HTMLResponse)
|
||||
async def change_detection_page():
|
||||
html_file = Path(__file__).parent / "change_detection_interface.html"
|
||||
if html_file.exists():
|
||||
return FileResponse(html_file)
|
||||
else:
|
||||
return HTMLResponse("<h2>Change Detection Interface not found.</h2>")
|
||||
|
||||
|
||||
@app.get("/", response_class=HTMLResponse)
|
||||
async def root():
|
||||
"""Serve main index page with tabs"""
|
||||
@@ -1846,6 +1895,683 @@ def run_batch_prediction(job: dict, config: PredictionConfig):
|
||||
print(f"[BATCH ERROR] Job {job['job_id']}: {traceback.format_exc()}")
|
||||
|
||||
|
||||
# ============ CHANGE DETECTION API ============
|
||||
|
||||
def rasterize_ground_truth(shapefile_path, out_shape, bbox, class_column="class"):
|
||||
"""Rasterize ground truth shapefile to match prediction raster shape."""
|
||||
try:
|
||||
import geopandas as gpd
|
||||
from rasterio import features as rio_features
|
||||
|
||||
gdf = gpd.read_file(shapefile_path)
|
||||
minx, miny, maxx, maxy = bbox
|
||||
|
||||
# Crop to bbox
|
||||
gdf = gdf.cx[minx:maxx, miny:maxy]
|
||||
|
||||
if class_column not in gdf.columns:
|
||||
raise ValueError(f"Shapefile missing '{class_column}' column. Available: {list(gdf.columns)}")
|
||||
|
||||
# Create transform for rasterization
|
||||
transform = from_bounds(minx, miny, maxx, maxy, out_shape[1], out_shape[0])
|
||||
|
||||
# Prepare geometries and values for rasterization
|
||||
shapes = zip(gdf.geometry, gdf[class_column])
|
||||
|
||||
# Rasterize
|
||||
gt_raster = rio_features.rasterize(
|
||||
shapes,
|
||||
out_shape=out_shape,
|
||||
fill=-1,
|
||||
transform=transform,
|
||||
dtype="int16"
|
||||
)
|
||||
|
||||
return gt_raster
|
||||
except Exception as e:
|
||||
print(f"[RASTERIZE ERROR] {e}")
|
||||
raise
|
||||
|
||||
|
||||
@app.post("/api/change-detection/predict")
|
||||
async def change_detection_predict_workflow(
|
||||
model_filename: str,
|
||||
min_lon: float,
|
||||
min_lat: float,
|
||||
max_lon: float,
|
||||
max_lat: float,
|
||||
start_date: str,
|
||||
end_date: str,
|
||||
max_scenes: int = 12,
|
||||
cloud_cover: int = 30,
|
||||
resolution: int = 20
|
||||
):
|
||||
"""
|
||||
Complete workflow: Predict + Compare with Ground Truth
|
||||
1. Load Sentinel-2 data for bbox and date range
|
||||
2. Run prediction using trained model
|
||||
3. Rasterize ground truth from training shapefile
|
||||
4. Compare and generate change detection results
|
||||
"""
|
||||
try:
|
||||
# --- STEP 1: LOAD MODEL ---
|
||||
model_path = Path("model_train") / model_filename
|
||||
if not model_path.exists():
|
||||
raise HTTPException(status_code=404, detail=f"Model not found: {model_filename}")
|
||||
|
||||
model_data = joblib.load(model_path)
|
||||
|
||||
if isinstance(model_data, dict):
|
||||
model = model_data.get('model')
|
||||
label_encoder = model_data.get('label_encoder')
|
||||
else:
|
||||
model = model_data
|
||||
label_encoder = None
|
||||
|
||||
print(f"[CHANGE DETECTION] Loaded model: {model_filename}")
|
||||
|
||||
# --- STEP 2: LOAD SENTINEL-2 DATA ---
|
||||
bbox = [min_lon, min_lat, max_lon, max_lat]
|
||||
time_range = f"{start_date}/{end_date}"
|
||||
|
||||
catalog = Client.open(
|
||||
"https://planetarycomputer.microsoft.com/api/stac/v1",
|
||||
modifier=planetary_computer.sign_inplace
|
||||
)
|
||||
|
||||
search = catalog.search(
|
||||
collections=["sentinel-2-l2a"],
|
||||
bbox=bbox,
|
||||
datetime=time_range,
|
||||
query={"eo:cloud_cover": {"lt": cloud_cover}}
|
||||
)
|
||||
|
||||
items = list(search.items())[:max_scenes]
|
||||
print(f"[CHANGE DETECTION] Found {len(items)} Sentinel-2 scenes")
|
||||
|
||||
if len(items) == 0:
|
||||
raise HTTPException(status_code=404, detail="No Sentinel-2 data found for the given area and date range")
|
||||
|
||||
# Load data
|
||||
signed_items = [planetary_computer.sign(item) for item in items]
|
||||
data = odc.stac.load(
|
||||
signed_items,
|
||||
bbox=bbox,
|
||||
bands=["B02", "B03", "B04", "B08"],
|
||||
resolution=resolution,
|
||||
chunks={"x": 2048, "y": 2048}
|
||||
).compute()
|
||||
|
||||
# --- STEP 3: CALCULATE NDVI ---
|
||||
print("[CHANGE DETECTION] Calculating NDVI...")
|
||||
nir = data["B08"].astype('float32')
|
||||
red = data["B04"].astype('float32')
|
||||
ndvi = (nir - red) / (nir + red + 1e-8)
|
||||
|
||||
# Handle clouds if SCL available
|
||||
if "SCL" in data:
|
||||
scl = data["SCL"]
|
||||
cloud_mask = (scl == 3) | (scl == 8) | (scl == 9) | (scl == 10)
|
||||
ndvi = ndvi.where(~cloud_mask)
|
||||
|
||||
# --- STEP 4: PREPARE FEATURES ---
|
||||
ndvi_filled = ndvi.ffill(dim='time').bfill(dim='time')
|
||||
ndvi_mean = np.nanmean(ndvi_filled.values, axis=0)
|
||||
|
||||
height, width = ndvi_mean.shape
|
||||
n_pixels = height * width
|
||||
|
||||
# Prepare features
|
||||
features = ndvi_mean.flatten().reshape(-1, 1)
|
||||
valid_mask = ~np.isnan(features[:, 0])
|
||||
features_clean = features[valid_mask]
|
||||
|
||||
# --- STEP 5: PREDICT ---
|
||||
print("[CHANGE DETECTION] Running prediction...")
|
||||
predictions = model.predict(features_clean)
|
||||
|
||||
# Decode labels if needed
|
||||
if label_encoder is not None:
|
||||
try:
|
||||
predictions = label_encoder.inverse_transform(predictions)
|
||||
except:
|
||||
pass
|
||||
|
||||
# Reshape to raster
|
||||
prediction_raster = np.full(n_pixels, -1, dtype=np.int16)
|
||||
prediction_raster[valid_mask] = predictions.astype(np.int16)
|
||||
prediction_raster = prediction_raster.reshape(height, width)
|
||||
|
||||
# --- STEP 6: COMPARE WITH GROUND TRUTH ---
|
||||
print("[CHANGE DETECTION] Comparing with ground truth...")
|
||||
gt_shapefile = "train/ST_training data_updated_1130points_new.shp"
|
||||
gt_raster = rasterize_ground_truth(gt_shapefile, (height, width), bbox, class_column="class")
|
||||
|
||||
# Calculate changes
|
||||
mask_valid = (gt_raster >= 0) & (prediction_raster >= 0)
|
||||
changes = gt_raster[mask_valid] != prediction_raster[mask_valid]
|
||||
n_total = np.count_nonzero(mask_valid)
|
||||
n_changed = np.count_nonzero(changes)
|
||||
|
||||
# Create change matrix
|
||||
from collections import Counter
|
||||
change_pairs = list(zip(gt_raster[mask_valid][changes], prediction_raster[mask_valid][changes]))
|
||||
change_counter = Counter(change_pairs)
|
||||
change_matrix = {f"{int(gt)}->{int(pred)}": int(cnt) for (gt, pred), cnt in change_counter.items()}
|
||||
|
||||
# Create change map
|
||||
change_map = np.full((height, width), -1, dtype=np.int8)
|
||||
change_map[mask_valid] = changes.astype(np.int8)
|
||||
|
||||
# --- STEP 7: SAVE RESULTS ---
|
||||
output_dir = Path("predictions")
|
||||
output_dir.mkdir(exist_ok=True)
|
||||
|
||||
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
|
||||
change_file = output_dir / f"change_map_{timestamp}.tif"
|
||||
|
||||
transform = from_bounds(min_lon, min_lat, max_lon, max_lat, width, height)
|
||||
|
||||
with rasterio.open(
|
||||
change_file, 'w',
|
||||
driver='GTiff',
|
||||
height=height,
|
||||
width=width,
|
||||
count=1,
|
||||
dtype=change_map.dtype,
|
||||
crs='EPSG:4326',
|
||||
transform=transform
|
||||
) as dst:
|
||||
dst.write(change_map, 1)
|
||||
|
||||
print(f"[CHANGE DETECTION] Saved change map to {change_file}")
|
||||
|
||||
# --- RETURN RESULTS ---
|
||||
return {
|
||||
"success": True,
|
||||
"n_scenes": len(items),
|
||||
"ndvi_stats": {
|
||||
"mean": float(np.nanmean(ndvi_mean)),
|
||||
"min": float(np.nanmin(ndvi_mean)),
|
||||
"max": float(np.nanmax(ndvi_mean)),
|
||||
"std": float(np.nanstd(ndvi_mean))
|
||||
},
|
||||
"class_distribution": {
|
||||
int(cls): int(count)
|
||||
for cls, count in zip(*np.unique(predictions, return_counts=True))
|
||||
},
|
||||
"change_detection": {
|
||||
"n_total_pixels": int(n_total),
|
||||
"n_changed_pixels": int(n_changed),
|
||||
"change_rate": float(n_changed) / n_total if n_total > 0 else 0.0,
|
||||
"change_matrix": change_matrix,
|
||||
"message": f"Detected {n_changed} changes out of {n_total} valid pixels ({(n_changed/n_total*100):.1f}%)" if n_total > 0 else "No valid pixels for comparison"
|
||||
},
|
||||
"change_map_file": str(change_file),
|
||||
"timestamp": timestamp
|
||||
}
|
||||
|
||||
except HTTPException:
|
||||
raise
|
||||
except Exception as e:
|
||||
print(f"[CHANGE DETECTION ERROR] {e}")
|
||||
import traceback
|
||||
traceback.print_exc()
|
||||
raise HTTPException(status_code=500, detail=f"Change detection workflow failed: {str(e)}")
|
||||
|
||||
|
||||
@app.post("/api/change-detection/workflow")
|
||||
async def change_detection_workflow(request: ChangeDetectionWorkflowRequest):
|
||||
"""Workflow: compare prediction with ground truth training data."""
|
||||
from collections import Counter
|
||||
|
||||
try:
|
||||
prediction_result = request.prediction_result
|
||||
bbox = request.bbox
|
||||
|
||||
if not prediction_result or "output_files" not in prediction_result:
|
||||
raise ValueError("Invalid prediction result")
|
||||
|
||||
# Get classification raster from prediction
|
||||
class_file = None
|
||||
for f in prediction_result.get("output_files", []):
|
||||
if f.get("type") == "classification":
|
||||
class_file = f.get("path")
|
||||
break
|
||||
|
||||
if not class_file:
|
||||
raise ValueError("No classification raster in prediction result")
|
||||
|
||||
# Load prediction raster
|
||||
with rasterio.open(class_file) as pred_ds:
|
||||
pred_arr = pred_ds.read(1)
|
||||
pred_crs = pred_ds.crs
|
||||
pred_transform = pred_ds.transform
|
||||
|
||||
# Rasterize ground truth training data
|
||||
gt_shapefile = "train/ST_training data_updated_1130points_new.shp"
|
||||
gt_raster = rasterize_ground_truth(gt_shapefile, pred_arr.shape, bbox, class_column="class")
|
||||
|
||||
# Calculate change detection
|
||||
mask_valid = (gt_raster >= 0) & (pred_arr >= 0) & ~np.isnan(gt_raster) & ~np.isnan(pred_arr)
|
||||
changes = gt_raster[mask_valid] != pred_arr[mask_valid]
|
||||
n_total = np.count_nonzero(mask_valid)
|
||||
n_changed = np.count_nonzero(changes)
|
||||
|
||||
# Create change pairs matrix
|
||||
change_pairs = list(zip(gt_raster[mask_valid][changes], pred_arr[mask_valid][changes]))
|
||||
change_counter = Counter(change_pairs)
|
||||
change_matrix = {f"{int(gt)}->{int(pred)}": int(cnt) for (gt, pred), cnt in change_counter.items()}
|
||||
|
||||
# Create change map
|
||||
change_map = np.full(pred_arr.shape, -1, dtype=np.int8)
|
||||
change_map[mask_valid] = changes.astype(np.int8)
|
||||
|
||||
# Change rate
|
||||
change_rate = float(n_changed) / n_total if n_total > 0 else 0.0
|
||||
|
||||
# Save change map
|
||||
change_dir = Path("predictions")
|
||||
change_dir.mkdir(exist_ok=True)
|
||||
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
|
||||
change_file = change_dir / f"change_map_{timestamp}.tif"
|
||||
|
||||
with rasterio.open(
|
||||
change_file, 'w',
|
||||
driver='GTiff',
|
||||
height=change_map.shape[0],
|
||||
width=change_map.shape[1],
|
||||
count=1,
|
||||
dtype=change_map.dtype,
|
||||
crs=pred_crs,
|
||||
transform=pred_transform
|
||||
) as dst:
|
||||
dst.write(change_map, 1)
|
||||
|
||||
return {
|
||||
"success": True,
|
||||
"change_detection": {
|
||||
"n_total_pixels": int(n_total),
|
||||
"n_changed_pixels": int(n_changed),
|
||||
"change_rate": change_rate,
|
||||
"change_matrix": change_matrix,
|
||||
"message": f"Detected {n_changed} changes out of {n_total} valid pixels ({change_rate*100:.1f}%)"
|
||||
},
|
||||
"change_map_file": str(change_file),
|
||||
"timestamp": timestamp
|
||||
}
|
||||
|
||||
except Exception as e:
|
||||
print(f"[CHANGE DETECTION WORKFLOW ERROR] {str(e)}")
|
||||
import traceback
|
||||
traceback.print_exc()
|
||||
raise HTTPException(status_code=500, detail=f"Change detection workflow failed: {str(e)}")
|
||||
|
||||
|
||||
@app.post("/api/change-detection/compare-periods")
|
||||
async def compare_periods(request: ComparePeriodsPredictionConfig, background_tasks: BackgroundTasks):
|
||||
"""Compare land use classification between two time periods."""
|
||||
from collections import Counter
|
||||
|
||||
try:
|
||||
# Extract parameters
|
||||
model_filename = request.model_filename
|
||||
bbox = [request.min_lon, request.min_lat, request.max_lon, request.max_lat]
|
||||
|
||||
current_start = request.current_period["start_date"]
|
||||
current_end = request.current_period["end_date"]
|
||||
pred_start = request.prediction_period["start_date"]
|
||||
pred_end = request.prediction_period["end_date"]
|
||||
|
||||
max_scenes = request.max_scenes
|
||||
cloud_cover = request.cloud_cover
|
||||
resolution = request.resolution
|
||||
|
||||
print(f"[COMPARE PERIODS] Current: {current_start} to {current_end} | Prediction: {pred_start} to {pred_end}")
|
||||
|
||||
# Step 1: Predict on current period
|
||||
print(f"[COMPARE PERIODS] Step 1: Predicting current period...")
|
||||
current_result = await predict_with_ndvi(PredictionWithNDVIConfig(
|
||||
model_filename=model_filename,
|
||||
min_lon=request.min_lon,
|
||||
min_lat=request.min_lat,
|
||||
max_lon=request.max_lon,
|
||||
max_lat=request.max_lat,
|
||||
start_date=current_start,
|
||||
end_date=current_end,
|
||||
max_scenes=max_scenes,
|
||||
cloud_cover=cloud_cover,
|
||||
resolution=resolution,
|
||||
export_classification=True,
|
||||
export_ndvi=False
|
||||
), background_tasks)
|
||||
|
||||
# Get current classification raster
|
||||
current_class_file = None
|
||||
for f in current_result.get("output_files", []):
|
||||
if f.get("type") == "classification":
|
||||
current_class_file = f.get("path")
|
||||
break
|
||||
|
||||
if not current_class_file:
|
||||
raise ValueError("No classification raster for current period")
|
||||
|
||||
# Step 2: Predict on prediction period
|
||||
print(f"[COMPARE PERIODS] Step 2: Predicting future period...")
|
||||
pred_result = await predict_with_ndvi(PredictionWithNDVIConfig(
|
||||
model_filename=model_filename,
|
||||
min_lon=request.min_lon,
|
||||
min_lat=request.min_lat,
|
||||
max_lon=request.max_lon,
|
||||
max_lat=request.max_lat,
|
||||
start_date=pred_start,
|
||||
end_date=pred_end,
|
||||
max_scenes=max_scenes,
|
||||
cloud_cover=cloud_cover,
|
||||
resolution=resolution,
|
||||
export_classification=True,
|
||||
export_ndvi=False
|
||||
), background_tasks)
|
||||
|
||||
# Get prediction classification raster
|
||||
pred_class_file = None
|
||||
for f in pred_result.get("output_files", []):
|
||||
if f.get("type") == "classification":
|
||||
pred_class_file = f.get("path")
|
||||
break
|
||||
|
||||
if not pred_class_file:
|
||||
raise ValueError("No classification raster for prediction period")
|
||||
|
||||
# Step 3: Load both rasters
|
||||
print(f"[COMPARE PERIODS] Step 3: Comparing classifications...")
|
||||
with rasterio.open(current_class_file) as src:
|
||||
current_arr = src.read(1)
|
||||
crs = src.crs
|
||||
transform = src.transform
|
||||
|
||||
with rasterio.open(pred_class_file) as src:
|
||||
pred_arr = src.read(1)
|
||||
|
||||
# Ensure same shape
|
||||
if current_arr.shape != pred_arr.shape:
|
||||
raise ValueError(f"Shape mismatch: current {current_arr.shape} vs prediction {pred_arr.shape}")
|
||||
|
||||
# Calculate changes
|
||||
mask_valid = ~np.isnan(current_arr) & ~np.isnan(pred_arr)
|
||||
changes = current_arr[mask_valid] != pred_arr[mask_valid]
|
||||
n_total = np.count_nonzero(mask_valid)
|
||||
n_changed = np.count_nonzero(changes)
|
||||
|
||||
# Create change pairs
|
||||
change_pairs = list(zip(current_arr[mask_valid][changes], pred_arr[mask_valid][changes]))
|
||||
change_counter = Counter(change_pairs)
|
||||
change_matrix = {f"{int(curr)}->{int(pred)}": int(cnt)
|
||||
for (curr, pred), cnt in change_counter.items()}
|
||||
|
||||
# Change rate
|
||||
change_rate = float(n_changed) / n_total if n_total > 0 else 0.0
|
||||
|
||||
# Save change map
|
||||
change_dir = Path("predictions")
|
||||
change_dir.mkdir(exist_ok=True)
|
||||
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
|
||||
change_file = change_dir / f"change_map_{timestamp}.tif"
|
||||
|
||||
change_map = np.zeros(current_arr.shape, dtype=np.int8)
|
||||
change_map[mask_valid] = changes.astype(np.int8)
|
||||
|
||||
with rasterio.open(
|
||||
change_file, 'w',
|
||||
driver='GTiff',
|
||||
height=change_map.shape[0],
|
||||
width=change_map.shape[1],
|
||||
count=1,
|
||||
dtype=change_map.dtype,
|
||||
crs=crs,
|
||||
transform=transform
|
||||
) as dst:
|
||||
dst.write(change_map, 1)
|
||||
|
||||
# Extract class distributions
|
||||
current_classes = np.unique(current_arr[~np.isnan(current_arr)]).astype(int)
|
||||
current_dist = {int(c): int(np.count_nonzero(current_arr == c)) for c in current_classes}
|
||||
|
||||
pred_classes = np.unique(pred_arr[~np.isnan(pred_arr)]).astype(int)
|
||||
pred_dist = {int(c): int(np.count_nonzero(pred_arr == c)) for c in pred_classes}
|
||||
|
||||
return {
|
||||
"success": True,
|
||||
"current_classification": {
|
||||
"n_scenes": current_result.get("n_scenes"),
|
||||
"resolution": current_result.get("resolution"),
|
||||
"class_distribution": current_dist
|
||||
},
|
||||
"prediction_classification": {
|
||||
"n_scenes": pred_result.get("n_scenes"),
|
||||
"resolution": pred_result.get("resolution"),
|
||||
"class_distribution": pred_dist
|
||||
},
|
||||
"change_detection": {
|
||||
"n_total_pixels": int(n_total),
|
||||
"n_changed_pixels": int(n_changed),
|
||||
"change_rate": change_rate,
|
||||
"change_matrix": change_matrix,
|
||||
"message": f"Detected {n_changed} changes out of {n_total} valid pixels ({change_rate*100:.1f}%)"
|
||||
},
|
||||
"change_map_file": str(change_file),
|
||||
"timestamp": timestamp
|
||||
}
|
||||
|
||||
except Exception as e:
|
||||
print(f"[COMPARE PERIODS ERROR] {str(e)}")
|
||||
import traceback
|
||||
traceback.print_exc()
|
||||
raise HTTPException(status_code=500, detail=f"Period comparison failed: {str(e)}")
|
||||
|
||||
|
||||
@app.post("/api/change-detection")
|
||||
async def change_detection_api(
|
||||
prediction_file: UploadFile = File(...),
|
||||
gt_file: UploadFile = File(...)
|
||||
):
|
||||
"""Detect changes between two raster files (prediction and ground truth)."""
|
||||
import tempfile
|
||||
from collections import Counter
|
||||
|
||||
try:
|
||||
# Save uploaded files temporarily
|
||||
pred_tmp = tempfile.NamedTemporaryFile(suffix='.tif', delete=False)
|
||||
gt_tmp = tempfile.NamedTemporaryFile(suffix='.tif', delete=False)
|
||||
|
||||
try:
|
||||
# Write uploaded files to temp
|
||||
pred_content = await prediction_file.read()
|
||||
gt_content = await gt_file.read()
|
||||
|
||||
pred_tmp.write(pred_content)
|
||||
gt_tmp.write(gt_content)
|
||||
pred_tmp.close()
|
||||
gt_tmp.close()
|
||||
|
||||
# Read prediction raster
|
||||
with rasterio.open(pred_tmp.name) as pred_ds:
|
||||
pred_arr = pred_ds.read(1)
|
||||
pred_crs = pred_ds.crs
|
||||
pred_transform = pred_ds.transform
|
||||
|
||||
# Read ground truth raster
|
||||
with rasterio.open(gt_tmp.name) as gt_ds:
|
||||
gt_arr = gt_ds.read(1)
|
||||
|
||||
# Ensure same shape
|
||||
if pred_arr.shape != gt_arr.shape:
|
||||
raise ValueError(f"Raster shapes don't match: prediction {pred_arr.shape} vs ground truth {gt_arr.shape}")
|
||||
|
||||
# Calculate change detection
|
||||
mask_valid = (gt_arr >= 0) & (pred_arr >= 0) & ~np.isnan(gt_arr) & ~np.isnan(pred_arr)
|
||||
changes = gt_arr[mask_valid] != pred_arr[mask_valid]
|
||||
n_total = np.count_nonzero(mask_valid)
|
||||
n_changed = np.count_nonzero(changes)
|
||||
|
||||
# Create change pairs matrix
|
||||
change_pairs = list(zip(gt_arr[mask_valid][changes], pred_arr[mask_valid][changes]))
|
||||
change_counter = Counter(change_pairs)
|
||||
change_matrix = {f"{int(gt)}->{int(pred)}": int(cnt) for (gt, pred), cnt in change_counter.items()}
|
||||
|
||||
# Create change map (0=same, 1=changed, -1=invalid)
|
||||
change_map = np.full(pred_arr.shape, -1, dtype=np.int8)
|
||||
change_map[mask_valid] = changes.astype(np.int8)
|
||||
|
||||
# Change rate
|
||||
change_rate = float(n_changed) / n_total if n_total > 0 else 0.0
|
||||
|
||||
# Save change map as GeoTIFF
|
||||
change_dir = Path("predictions")
|
||||
change_dir.mkdir(exist_ok=True)
|
||||
|
||||
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
|
||||
change_file = change_dir / f"change_map_{timestamp}.tif"
|
||||
|
||||
with rasterio.open(
|
||||
change_file, 'w',
|
||||
driver='GTiff',
|
||||
height=change_map.shape[0],
|
||||
width=change_map.shape[1],
|
||||
count=1,
|
||||
dtype=change_map.dtype,
|
||||
crs=pred_crs,
|
||||
transform=pred_transform
|
||||
) as dst:
|
||||
dst.write(change_map, 1)
|
||||
|
||||
# Return results
|
||||
return {
|
||||
"success": True,
|
||||
"change_detection": {
|
||||
"n_total_pixels": int(n_total),
|
||||
"n_changed_pixels": int(n_changed),
|
||||
"change_rate": change_rate,
|
||||
"change_matrix": change_matrix,
|
||||
"message": f"Detected {n_changed} changes out of {n_total} valid pixels ({change_rate*100:.1f}%)"
|
||||
},
|
||||
"change_map_file": str(change_file),
|
||||
"timestamp": timestamp
|
||||
}
|
||||
|
||||
finally:
|
||||
# Cleanup temp files
|
||||
try:
|
||||
Path(pred_tmp.name).unlink()
|
||||
Path(gt_tmp.name).unlink()
|
||||
except:
|
||||
pass
|
||||
|
||||
except Exception as e:
|
||||
print(f"[CHANGE DETECTION ERROR] {str(e)}")
|
||||
raise HTTPException(status_code=500, detail=f"Change detection failed: {str(e)}")
|
||||
async def change_detection_api(
|
||||
prediction_file: UploadFile = File(...),
|
||||
gt_file: UploadFile = File(...)
|
||||
):
|
||||
"""Detect changes between prediction raster and ground truth raster."""
|
||||
import tempfile
|
||||
from collections import Counter
|
||||
|
||||
try:
|
||||
# Save uploaded files temporarily
|
||||
pred_tmp = tempfile.NamedTemporaryFile(suffix='.tif', delete=False)
|
||||
gt_tmp = tempfile.NamedTemporaryFile(suffix='.tif', delete=False)
|
||||
|
||||
try:
|
||||
# Write uploaded files to temp
|
||||
pred_content = await prediction_file.read()
|
||||
gt_content = await gt_file.read()
|
||||
|
||||
pred_tmp.write(pred_content)
|
||||
gt_tmp.write(gt_content)
|
||||
pred_tmp.close()
|
||||
gt_tmp.close()
|
||||
|
||||
# Read prediction raster
|
||||
import rasterio
|
||||
with rasterio.open(pred_tmp.name) as pred_ds:
|
||||
pred_arr = pred_ds.read(1)
|
||||
pred_crs = pred_ds.crs
|
||||
pred_transform = pred_ds.transform
|
||||
|
||||
# Read ground truth raster
|
||||
with rasterio.open(gt_tmp.name) as gt_ds:
|
||||
gt_arr = gt_ds.read(1)
|
||||
|
||||
# Ensure same shape
|
||||
if pred_arr.shape != gt_arr.shape:
|
||||
raise ValueError(f"Raster shapes don't match: prediction {pred_arr.shape} vs ground truth {gt_arr.shape}")
|
||||
|
||||
# Calculate change detection
|
||||
mask_valid = (gt_arr >= 0) & (pred_arr >= 0) & ~np.isnan(gt_arr) & ~np.isnan(pred_arr)
|
||||
changes = gt_arr[mask_valid] != pred_arr[mask_valid]
|
||||
n_total = np.count_nonzero(mask_valid)
|
||||
n_changed = np.count_nonzero(changes)
|
||||
|
||||
# Create change pairs matrix
|
||||
change_pairs = list(zip(gt_arr[mask_valid][changes], pred_arr[mask_valid][changes]))
|
||||
change_counter = Counter(change_pairs)
|
||||
change_matrix = {f"{int(gt)}->{int(pred)}": int(cnt) for (gt, pred), cnt in change_counter.items()}
|
||||
|
||||
# Create change map (0=same, 1=changed, -1=invalid)
|
||||
change_map = np.full(pred_arr.shape, -1, dtype=np.int8)
|
||||
change_map[mask_valid] = changes.astype(np.int8)
|
||||
|
||||
# Change rate
|
||||
change_rate = float(n_changed) / n_total if n_total > 0 else 0.0
|
||||
|
||||
# Save change map as GeoTIFF
|
||||
change_dir = Path("predictions")
|
||||
change_dir.mkdir(exist_ok=True)
|
||||
|
||||
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
|
||||
change_file = change_dir / f"change_map_{timestamp}.tif"
|
||||
|
||||
with rasterio.open(
|
||||
change_file, 'w',
|
||||
driver='GTiff',
|
||||
height=change_map.shape[0],
|
||||
width=change_map.shape[1],
|
||||
count=1,
|
||||
dtype=change_map.dtype,
|
||||
crs=pred_crs,
|
||||
transform=pred_transform
|
||||
) as dst:
|
||||
dst.write(change_map, 1)
|
||||
|
||||
# Return results
|
||||
return {
|
||||
"success": True,
|
||||
"change_detection": {
|
||||
"n_total_pixels": int(n_total),
|
||||
"n_changed_pixels": int(n_changed),
|
||||
"change_rate": change_rate,
|
||||
"change_matrix": change_matrix,
|
||||
"message": f"Detected {n_changed} changes out of {n_total} valid pixels ({change_rate*100:.1f}%)"
|
||||
},
|
||||
"change_map_file": str(change_file),
|
||||
"timestamp": timestamp
|
||||
}
|
||||
|
||||
finally:
|
||||
# Cleanup temp files
|
||||
try:
|
||||
Path(pred_tmp.name).unlink()
|
||||
Path(gt_tmp.name).unlink()
|
||||
except:
|
||||
pass
|
||||
|
||||
except Exception as e:
|
||||
print(f"[CHANGE DETECTION ERROR] {str(e)}")
|
||||
raise HTTPException(status_code=500, detail=f"Change detection failed: {str(e)}")
|
||||
|
||||
|
||||
# ============ PREDICTION WITH NDVI API ============
|
||||
|
||||
@app.post("/api/predict/with-ndvi")
|
||||
@@ -1853,13 +2579,6 @@ async def predict_with_ndvi(config: PredictionWithNDVIConfig, background_tasks:
|
||||
"""Predict land classification và NDVI cho một khu vực"""
|
||||
try:
|
||||
import numpy as np
|
||||
import xarray as xr
|
||||
from pystac_client import Client
|
||||
import planetary_computer
|
||||
import odc.stac
|
||||
import rasterio
|
||||
from rasterio.transform import from_bounds
|
||||
import hashlib, os
|
||||
|
||||
# Load model
|
||||
model_path = Path(f"model_train/{config.model_filename}")
|
||||
@@ -1985,6 +2704,50 @@ async def predict_with_ndvi(config: PredictionWithNDVIConfig, background_tasks:
|
||||
prediction_raster = np.full(n_pixels, -1, dtype=np.int16)
|
||||
prediction_raster[valid_mask] = predictions
|
||||
prediction_raster = prediction_raster.reshape(height, width)
|
||||
|
||||
# --- CHANGE DETECTION ---
|
||||
change_summary = None
|
||||
change_map = None
|
||||
try:
|
||||
# Use training shapefile as ground truth
|
||||
gt_shapefile = "train/ST_training data_updated_1130points_new.shp"
|
||||
gt_raster = rasterize_ground_truth(gt_shapefile, (height, width), bbox, class_column="class")
|
||||
# Compare prediction and ground truth
|
||||
mask_valid = (gt_raster >= 0) & (prediction_raster >= 0)
|
||||
changes = gt_raster[mask_valid] != prediction_raster[mask_valid]
|
||||
n_total = np.count_nonzero(mask_valid)
|
||||
n_changed = np.count_nonzero(changes)
|
||||
# Per-class change matrix
|
||||
from collections import Counter
|
||||
change_pairs = list(zip(gt_raster[mask_valid][changes], prediction_raster[mask_valid][changes]))
|
||||
change_counter = Counter(change_pairs)
|
||||
change_matrix = {f"{int(gt)}->{int(pred)}": int(cnt) for (gt, pred), cnt in change_counter.items()}
|
||||
change_summary = {
|
||||
"n_total": int(n_total),
|
||||
"n_changed": int(n_changed),
|
||||
"change_rate": float(n_changed) / n_total if n_total > 0 else 0.0,
|
||||
"change_matrix": change_matrix
|
||||
}
|
||||
# Optionally, create a change map (1=changed, 0=same, -1=invalid)
|
||||
change_map = np.full((height, width), -1, dtype=np.int8)
|
||||
change_map[mask_valid] = changes.astype(np.int8)
|
||||
# Save change map as GeoTIFF
|
||||
change_file = output_dir / f"change_map_{timestamp}.tif"
|
||||
with rasterio.open(
|
||||
change_file, 'w',
|
||||
driver='GTiff',
|
||||
height=height,
|
||||
width=width,
|
||||
count=1,
|
||||
dtype=change_map.dtype,
|
||||
crs='EPSG:4326',
|
||||
transform=from_bounds(bbox[0], bbox[1], bbox[2], bbox[3], width, height)
|
||||
) as dst:
|
||||
dst.write(change_map, 1)
|
||||
output_files.append({"type": "change_map", "path": str(change_file)})
|
||||
print(f"[CHANGE DETECTION] Saved change map to {change_file}")
|
||||
except Exception as change_exc:
|
||||
print(f"[CHANGE DETECTION] Warning: {change_exc}")
|
||||
|
||||
# Prepare outputs
|
||||
timestamp = datetime.now().strftime('%Y%m%d_%H%M%S')
|
||||
@@ -2055,7 +2818,8 @@ async def predict_with_ndvi(config: PredictionWithNDVIConfig, background_tasks:
|
||||
"class_distribution": class_distribution,
|
||||
"n_scenes": len(items),
|
||||
"resolution": config.resolution,
|
||||
"bbox": bbox
|
||||
"bbox": bbox,
|
||||
"change_detection": change_summary
|
||||
}
|
||||
|
||||
except Exception as e:
|
||||
|
||||
Reference in New Issue
Block a user