cập nhật chức năng predict các loại đất
This commit is contained in:
Victor Phan
+153
-36
@@ -330,7 +330,7 @@ def update_prediction_progress(message: str):
|
||||
|
||||
|
||||
async def run_prediction(config: PredictionConfig):
|
||||
"""Chạy prediction process"""
|
||||
"""Chạy prediction process - Áp dụng phương pháp từ 02.predict_ODC.ipynb"""
|
||||
global prediction_status
|
||||
|
||||
try:
|
||||
@@ -341,6 +341,7 @@ async def run_prediction(config: PredictionConfig):
|
||||
import numpy as np
|
||||
from datetime import datetime as dt
|
||||
import rioxarray
|
||||
import dask.array as da
|
||||
|
||||
prediction_status["progress"] = "Đang load model..."
|
||||
|
||||
@@ -364,6 +365,7 @@ async def run_prediction(config: PredictionConfig):
|
||||
# Import and use Microsoft Planetary Computer STAC API
|
||||
import pystac_client
|
||||
import planetary_computer
|
||||
from odc.stac import load
|
||||
|
||||
catalog = pystac_client.Client.open(
|
||||
"https://planetarycomputer.microsoft.com/api/stac/v1",
|
||||
@@ -373,65 +375,168 @@ async def run_prediction(config: PredictionConfig):
|
||||
bbox = [config.min_lon, config.min_lat, config.max_lon, config.max_lat]
|
||||
time_range = f"{config.start_date}/{config.end_date}"
|
||||
|
||||
# ============ BƯỚC 1: TẢI DỮ LIỆU SENTINEL-2 ============
|
||||
prediction_status["progress"] = "Đang tải dữ liệu Sentinel-2..."
|
||||
|
||||
# Search Sentinel-2 data
|
||||
search = catalog.search(
|
||||
s2_search = catalog.search(
|
||||
collections=["sentinel-2-l2a"],
|
||||
bbox=bbox,
|
||||
datetime=time_range,
|
||||
query={"eo:cloud_cover": {"lt": config.cloud_cover}}
|
||||
)
|
||||
|
||||
items = list(search.items()) # Changed from items_as_dicts() to items()
|
||||
if not items:
|
||||
s2_items = list(s2_search.items())
|
||||
if not s2_items:
|
||||
raise ValueError("Không tìm thấy dữ liệu Sentinel-2 cho khu vực và thời gian này")
|
||||
|
||||
items = items[:config.max_scenes]
|
||||
|
||||
prediction_status["progress"] = f"Đang xử lý {len(items)} scenes Sentinel-2..."
|
||||
|
||||
# Load and process Sentinel-2 data (simplified)
|
||||
# Note: This is a simplified version. Full implementation would need more processing
|
||||
from odc.stac import load
|
||||
s2_items = s2_items[:config.max_scenes]
|
||||
prediction_status["progress"] = f"Đang xử lý {len(s2_items)} scenes Sentinel-2..."
|
||||
|
||||
# Load Sentinel-2 data
|
||||
s2_data = load(
|
||||
items,
|
||||
s2_items,
|
||||
bbox=bbox,
|
||||
chunks={"time": 1, "x": 2048, "y": 2048},
|
||||
groupby="solar_day",
|
||||
resolution=config.resolution
|
||||
)
|
||||
|
||||
prediction_status["progress"] = "Đang tính toán các chỉ số..."
|
||||
# ============ BƯỚC 2: TÍNH NDVI VÀ XỬ LÝ MÂY ============
|
||||
prediction_status["progress"] = "Đang tính toán NDVI và xử lý mây..."
|
||||
|
||||
# Calculate NDVI using Sentinel-2 band names
|
||||
# B08 = NIR, B04 = Red
|
||||
nir = s2_data["B08"] # NIR band
|
||||
red = s2_data["B04"] # Red band
|
||||
ndvi = (nir - red) / (nir + red + 1e-8) # Add small value to avoid division by zero
|
||||
# Calculate NDVI using Sentinel-2 band names (B08 = NIR, B04 = Red)
|
||||
nir = s2_data["B08"].astype('float32')
|
||||
red = s2_data["B04"].astype('float32')
|
||||
ndvi = (nir - red) / (nir + red + 1e-8)
|
||||
|
||||
# Resample to monthly
|
||||
ndvi_monthly = ndvi.resample(time="1M").mean()
|
||||
# Mask clouds using SCL band if available
|
||||
if "SCL" in s2_data:
|
||||
scl = s2_data["SCL"]
|
||||
# SCL values: 4=vegetation, 5=bare soil, 6=water - these are clear
|
||||
# 3=cloud shadow, 8=cloud medium, 9=cloud high, 10=cirrus - mask these
|
||||
cloud_mask = (scl == 3) | (scl == 8) | (scl == 9) | (scl == 10)
|
||||
ndvi = ndvi.where(~cloud_mask)
|
||||
|
||||
prediction_status["progress"] = "Đang dự đoán..."
|
||||
# ============ BƯỚC 3: ĐIỀN GIÁ TRỊ NAN (FILL NAN) ============
|
||||
prediction_status["progress"] = "Đang điền giá trị bị che mây..."
|
||||
|
||||
# Prepare features for prediction
|
||||
features_list = []
|
||||
for t in range(len(ndvi_monthly.time)):
|
||||
ndvi_t = ndvi_monthly.isel(time=t).values
|
||||
features_list.append(ndvi_t.flatten())
|
||||
# Fill NaN using forward fill and backward fill
|
||||
ndvi_filled = ndvi.ffill(dim='time').bfill(dim='time')
|
||||
|
||||
# Stack features
|
||||
features = np.column_stack(features_list)
|
||||
# Resample to monthly average
|
||||
prediction_status["progress"] = "Đang tính trung bình NDVI theo tháng..."
|
||||
ndvi_monthly = ndvi_filled.resample(time="1ME").mean()
|
||||
|
||||
# Compute NDVI (convert from dask to numpy)
|
||||
ndvi_monthly = ndvi_monthly.compute()
|
||||
|
||||
# ============ BƯỚC 4: TẢI DỮ LIỆU SENTINEL-1 (VH, VV) ============
|
||||
prediction_status["progress"] = "Đang tải dữ liệu Sentinel-1 (Radar)..."
|
||||
|
||||
# Search Sentinel-1 data
|
||||
s1_search = catalog.search(
|
||||
collections=["sentinel-1-rtc"],
|
||||
bbox=bbox,
|
||||
datetime=time_range,
|
||||
)
|
||||
|
||||
s1_items = list(s1_search.items())
|
||||
|
||||
if s1_items:
|
||||
s1_items = s1_items[:config.max_scenes]
|
||||
prediction_status["progress"] = f"Đang xử lý {len(s1_items)} scenes Sentinel-1..."
|
||||
|
||||
# Load Sentinel-1 data
|
||||
s1_data = load(
|
||||
s1_items,
|
||||
bbox=bbox,
|
||||
chunks={"time": 1, "x": 2048, "y": 2048},
|
||||
groupby="sat:absolute_orbit",
|
||||
resolution=config.resolution,
|
||||
like=ndvi_monthly # Align with NDVI grid
|
||||
)
|
||||
|
||||
# Extract VH and VV bands
|
||||
if "vh" in s1_data and "vv" in s1_data:
|
||||
vh = s1_data["vh"].astype('float32')
|
||||
vv = s1_data["vv"].astype('float32')
|
||||
|
||||
# Resample to monthly average
|
||||
prediction_status["progress"] = "Đang tính trung bình VH/VV theo tháng..."
|
||||
vh_monthly = vh.resample(time="1ME").mean().compute()
|
||||
vv_monthly = vv.resample(time="1ME").mean().compute()
|
||||
|
||||
use_radar = True
|
||||
else:
|
||||
prediction_status["progress"] = "Không tìm thấy bands VH/VV, tiếp tục với NDVI..."
|
||||
use_radar = False
|
||||
else:
|
||||
prediction_status["progress"] = "Không có dữ liệu Sentinel-1, tiếp tục với NDVI..."
|
||||
use_radar = False
|
||||
|
||||
# ============ BƯỚC 5: CHUẨN BỊ FEATURES CHO DỰ ĐOÁN ============
|
||||
prediction_status["progress"] = "Đang chuẩn bị features cho dự đoán..."
|
||||
|
||||
# Get shape information
|
||||
n_times_ndvi = len(ndvi_monthly.time)
|
||||
y_size = len(ndvi_monthly.y)
|
||||
x_size = len(ndvi_monthly.x)
|
||||
n_pixels = y_size * x_size
|
||||
|
||||
# Prepare NDVI features (flatten each time step)
|
||||
ndvi_features = []
|
||||
for t in range(n_times_ndvi):
|
||||
ndvi_t = ndvi_monthly.isel(time=t).values.flatten()
|
||||
ndvi_features.append(ndvi_t)
|
||||
|
||||
# Stack NDVI features
|
||||
features = np.column_stack(ndvi_features)
|
||||
|
||||
# Add radar features if available
|
||||
if use_radar:
|
||||
n_times_vh = len(vh_monthly.time)
|
||||
n_times_vv = len(vv_monthly.time)
|
||||
|
||||
# Add VH features
|
||||
for t in range(min(n_times_vh, n_times_ndvi)):
|
||||
vh_t = vh_monthly.isel(time=t).values.flatten()
|
||||
# Resize if needed
|
||||
if len(vh_t) != n_pixels:
|
||||
vh_t = np.resize(vh_t, n_pixels)
|
||||
features = np.column_stack([features, vh_t])
|
||||
|
||||
# Add VV features
|
||||
for t in range(min(n_times_vv, n_times_ndvi)):
|
||||
vv_t = vv_monthly.isel(time=t).values.flatten()
|
||||
# Resize if needed
|
||||
if len(vv_t) != n_pixels:
|
||||
vv_t = np.resize(vv_t, n_pixels)
|
||||
features = np.column_stack([features, vv_t])
|
||||
|
||||
# Handle NaN values in features
|
||||
features = np.nan_to_num(features, nan=0.0)
|
||||
|
||||
# ============ BƯỚC 6: DỰ ĐOÁN ============
|
||||
prediction_status["progress"] = f"Đang dự đoán với {features.shape[1]} features..."
|
||||
|
||||
# Make prediction
|
||||
predictions = model.predict(features)
|
||||
|
||||
# Decode labels if label_encoder exists
|
||||
if label_encoder is not None:
|
||||
try:
|
||||
predictions = label_encoder.inverse_transform(predictions)
|
||||
except:
|
||||
pass # Keep numeric predictions if inverse_transform fails
|
||||
|
||||
# Reshape to original shape
|
||||
pred_shape = ndvi_monthly.isel(time=0).shape
|
||||
pred_shape = (y_size, x_size)
|
||||
predictions_2d = predictions.reshape(pred_shape)
|
||||
|
||||
# ============ BƯỚC 7: TẠO OUTPUT VÀ LƯU KẾT QUẢ ============
|
||||
prediction_status["progress"] = "Đang tạo bản đồ phân loại..."
|
||||
|
||||
# Create output xarray
|
||||
prediction_da = xr.DataArray(
|
||||
predictions_2d,
|
||||
@@ -450,21 +555,33 @@ async def run_prediction(config: PredictionConfig):
|
||||
timestamp = dt.now().strftime("%Y%m%d_%H%M%S")
|
||||
output_file = output_dir / f"prediction_{timestamp}.tif"
|
||||
|
||||
prediction_status["progress"] = "Đang lưu kết quả..."
|
||||
prediction_status["progress"] = "Đang lưu kết quả GeoTIFF..."
|
||||
|
||||
# Save as GeoTIFF
|
||||
prediction_da.rio.write_crs(s2_data.rio.crs, inplace=True)
|
||||
prediction_da.rio.to_raster(output_file, driver="GTiff")
|
||||
# Set CRS and save as GeoTIFF
|
||||
if hasattr(s2_data, 'rio') and s2_data.rio.crs is not None:
|
||||
prediction_da.rio.write_crs(s2_data.rio.crs, inplace=True)
|
||||
else:
|
||||
prediction_da.rio.write_crs("EPSG:4326", inplace=True)
|
||||
|
||||
prediction_da.rio.to_raster(str(output_file), driver="GTiff")
|
||||
|
||||
# Get unique classes for result
|
||||
unique_classes = np.unique(predictions_2d)
|
||||
unique_classes = unique_classes[~np.isnan(unique_classes)].tolist()
|
||||
|
||||
prediction_status["is_predicting"] = False
|
||||
prediction_status["progress"] = "Hoàn thành!"
|
||||
prediction_status["output_file"] = str(output_file)
|
||||
prediction_status["result"] = {
|
||||
"output_file": str(output_file),
|
||||
"shape": pred_shape,
|
||||
"unique_classes": np.unique(predictions).tolist(),
|
||||
"shape": list(pred_shape),
|
||||
"unique_classes": unique_classes,
|
||||
"bbox": bbox,
|
||||
"time_range": time_range
|
||||
"time_range": time_range,
|
||||
"n_features": features.shape[1],
|
||||
"n_times_ndvi": n_times_ndvi,
|
||||
"used_radar": use_radar,
|
||||
"model_used": config.model_filename
|
||||
}
|
||||
prediction_status["end_time"] = dt.now().isoformat()
|
||||
|
||||
|
||||
Reference in New Issue
Block a user