basketballcantho/remote-sensing
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

bổ sung thêm hàm tự resign token SAS

Browse Source
This commit is contained in:
Victor Phan
2026-01-19 12:31:40 +07:00
parent 3a7d5bb21b
commit 61646de647
1 changed files with 565 additions and 192 deletions
Download Patch File Download Diff File Expand all files Collapse all files
api_server.py
+565 -192
View File
@@ -22,6 +22,8 @@ from rasterio.transform import from_bounds
import asyncio
import hashlib
import traceback
import socket
import urllib.request
# Import report generator
from report_generator import generate_training_report, generate_prediction_report
@@ -107,6 +109,58 @@ DEFAULT_LABEL_NAMES = {
}
# ============ NETWORK CONNECTIVITY HELPERS ============
def check_internet_connectivity(timeout=5):
"""
Check if we have internet connectivity to common hosts
Returns: (is_connected: bool, error_message: str)
"""
test_hosts = [
("8.8.8.8", 53), # Google DNS
("1.1.1.1", 53), # Cloudflare DNS
]
for host, port in test_hosts:
try:
socket.setdefaulttimeout(timeout)
socket.socket(socket.AF_INET, socket.SOCK_STREAM).connect((host, port))
return True, None
except (socket.timeout, socket.error):
continue
return False, "No internet connectivity detected. Please check your network connection."
def check_planetary_computer_access(timeout=10):
"""
Check if we can access Microsoft Planetary Computer
Returns: (is_accessible: bool, error_message: str)
"""
test_urls = [
"https://planetarycomputer.microsoft.com/api/stac/v1",
]
for url in test_urls:
try:
req = urllib.request.Request(url)
response = urllib.request.urlopen(req, timeout=timeout)
if response.status == 200:
return True, None
except urllib.error.HTTPError as e:
# 405 Method Not Allowed is OK - means the service is reachable
if e.code in [200, 405]:
return True, None
return False, f"HTTP Error {e.code}: {e.reason}"
except urllib.error.URLError as e:
if "Could not resolve host" in str(e):
return False, f"Cannot resolve Planetary Computer hosts. DNS issue detected: {str(e)}"
return False, f"URL Error: {str(e)}"
except Exception as e:
return False, f"Unexpected error: {str(e)}"
return False, "Cannot access Microsoft Planetary Computer. Service may be down or blocked by firewall."
class TrainingConfig(BaseModel):
"""Cấu hình training - Tất cả bắt buộc nhập từ giao diện"""
# Khu vực (bbox)
@@ -233,6 +287,44 @@ async def change_detection_page():
return HTMLResponse("<h2>Change Detection Interface not found.</h2>")
@app.get("/api/network/check")
async def check_network():
"""
Check network connectivity and access to Planetary Computer
Useful for diagnosing connection issues before attempting data downloads
"""
result = {
"timestamp": datetime.now().isoformat(),
"internet_connected": False,
"planetary_computer_accessible": False,
"internet_error": None,
"pc_error": None,
"status": "unknown"
}
# Check basic internet connectivity
is_connected, conn_error = check_internet_connectivity(timeout=5)
result["internet_connected"] = is_connected
result["internet_error"] = conn_error
# Check Planetary Computer access
if is_connected:
is_accessible, access_error = check_planetary_computer_access(timeout=10)
result["planetary_computer_accessible"] = is_accessible
result["pc_error"] = access_error
if is_accessible:
result["status"] = "ready"
else:
result["status"] = "pc_unavailable"
else:
result["status"] = "no_internet"
result["planetary_computer_accessible"] = False
result["pc_error"] = "Cannot check - no internet connection"
return result
@app.get("/", response_class=HTMLResponse)
async def root():
"""Serve main index page with tabs"""
@@ -1234,14 +1326,51 @@ async def run_prediction(config: PredictionConfig):
else: # temporal or extended
bands_to_load = ["B02", "B03", "B04", "B08", "B11", "SCL"]
s2_data = load(
s2_items,
bbox=bbox,
bands=bands_to_load,
chunks={"time": 1, "x": 2048, "y": 2048},
groupby="solar_day",
resolution=config.resolution
).compute()
# Retry logic for data loading
max_retries = 3
s2_data = None
for attempt in range(max_retries):
try:
s2_data = load(
s2_items,
bbox=bbox,
bands=bands_to_load,
chunks={"time": 1, "x": 2048, "y": 2048},
groupby="solar_day",
resolution=config.resolution
).compute()
break
except Exception as e:
error_msg = str(e)
is_auth_error = any(err in error_msg for err in [
"AuthenticationFailed", "Signature fields not well formed",
"Server failed to authenticate", "403", "401"
])
if is_auth_error:
print(f"⚠️ Authentication error. Re-signing items...")
try:
import planetary_computer
s2_items = [planetary_computer.sign(item) for item in s2_items]
if attempt < max_retries - 1:
import time
time.sleep(1)
continue
except Exception:
pass
if "does not exist" in error_msg or "RasterioIOError" in error_msg:
if len(s2_items) > 3 and attempt < max_retries - 1:
remove_count = max(1, len(s2_items) // 5)
s2_items = s2_items[:-remove_count]
print(f"⚠️ Access error. Reduced to {len(s2_items)} scenes, retrying...")
continue
if attempt == max_retries - 1:
raise ValueError(f"Failed to load Sentinel-2 data: {error_msg[:200]}")
if s2_data is None:
raise ValueError("Failed to load Sentinel-2 data after retries")
prediction_status["progress"] = "Đã load Sentinel-2 data"
@@ -1261,20 +1390,25 @@ async def run_prediction(config: PredictionConfig):
if s1_items:
s1_items = s1_items[:config.max_scenes]
s1_data = load(
s1_items,
bbox=bbox,
bands=["vh", "vv"],
chunks={"time": 1, "x": 2048, "y": 2048},
groupby="solar_day",
resolution=config.resolution
).compute()
# Convert to dB
vh_data = 10 * np.log10(s1_data['vh'].where(s1_data['vh'] > 0))
vv_data = 10 * np.log10(s1_data['vv'].where(s1_data['vv'] > 0))
use_radar = True
prediction_status["progress"] = f"Đã load Sentinel-1 data ({len(s1_items)} scenes)"
# Try to load Sentinel-1 with error handling
try:
s1_data = load(
s1_items,
bbox=bbox,
bands=["vh", "vv"],
chunks={"time": 1, "x": 2048, "y": 2048},
groupby="solar_day",
resolution=config.resolution
).compute()
# Convert to dB
vh_data = 10 * np.log10(s1_data['vh'].where(s1_data['vh'] > 0))
vv_data = 10 * np.log10(s1_data['vv'].where(s1_data['vv'] > 0))
use_radar = True
prediction_status["progress"] = f"Đã load Sentinel-1 data ({len(s1_items)} scenes)"
except Exception as s1_error:
print(f"⚠️ Failed to load Sentinel-1 data: {str(s1_error)[:100]}")
prediction_status["progress"] = "Lỗi Sentinel-1, bỏ qua radar features"
else:
prediction_status["progress"] = "Không có dữ liệu Sentinel-1, bỏ qua radar features"
except Exception as e:
@@ -2647,15 +2781,55 @@ async def change_detection_predict_workflow(
if len(items) == 0:
raise HTTPException(status_code=404, detail="No Sentinel-2 data found for the given area and date range")
# Load data
# Load data with error handling
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()
max_retries = 3
data = None
for attempt in range(max_retries):
try:
data = odc.stac.load(
signed_items,
bbox=bbox,
bands=["B02", "B03", "B04", "B08"],
resolution=resolution,
chunks={"x": 2048, "y": 2048}
).compute()
break
except Exception as e:
error_msg = str(e)
is_auth_error = any(err in error_msg for err in [
"AuthenticationFailed", "Signature fields not well formed",
"Server failed to authenticate", "403", "401"
])
if is_auth_error:
print(f"⚠️ Authentication error. Re-signing items...")
try:
signed_items = [planetary_computer.sign(item) for item in items]
if attempt < max_retries - 1:
import time
time.sleep(1)
continue
except Exception:
pass
if ("does not exist" in error_msg or "RasterioIOError" in error_msg) and len(signed_items) > 3:
remove_count = max(1, len(signed_items) // 5)
signed_items = signed_items[:-remove_count]
items = items[:-remove_count]
print(f"⚠️ Access error. Reduced to {len(signed_items)} scenes, retrying...")
if attempt < max_retries - 1:
continue
if attempt == max_retries - 1:
raise HTTPException(
status_code=500,
detail=f"Failed to load Sentinel-2 data: {error_msg[:200]}"
)
if data is None:
raise HTTPException(status_code=500, detail="Failed to load Sentinel-2 data")
# --- STEP 3: CALCULATE NDVI ---
print("[CHANGE DETECTION] Calculating NDVI...")
@@ -3931,6 +4105,24 @@ async def ndvi_predict_timeseries(config: NDVIPredictionConfig):
bbox = [min_lon, min_lat, max_lon, max_lat]
time_range = f"{config.start_date}/{config.end_date}"
# Check network connectivity first
print(f"\n🌐 Checking network connectivity...")
is_connected, conn_error = check_internet_connectivity()
if not is_connected:
raise HTTPException(
status_code=503,
detail=f"❌ {conn_error}\n\n"
f"Please check:\n"
f"1. Your internet connection is active\n"
f"2. DNS servers are configured correctly\n"
f"3. Firewall/proxy settings allow outbound connections"
)
is_accessible, access_error = check_planetary_computer_access()
if not is_accessible:
print(f"⚠️ {access_error}")
print(f" Will attempt to load data anyway (may fail)...")
# Load Sentinel-2 time series using Planetary Computer
print(f"\n📡 Loading Sentinel-2 data from Planetary Computer...")
@@ -3984,8 +4176,85 @@ async def ndvi_predict_timeseries(config: NDVIPredictionConfig):
except Exception as e:
error_msg = str(e)
if "Could not resolve host" in error_msg or "CURL error" in error_msg:
# Check for various error types that might indicate problematic scenes
is_network_error = any(err in error_msg for err in [
"Could not resolve host", "CURL error", "Connection", "Timeout"
])
is_access_error = any(err in error_msg for err in [
"does not exist in the file system", "RasterioIOError",
"not recognized as a supported dataset", "Aborting load"
])
is_auth_error = any(err in error_msg for err in [
"AuthenticationFailed", "Signature fields not well formed",
"Server failed to authenticate", "403", "401"
])
if is_auth_error:
print(f"⚠️ Authentication error detected. Re-signing items and retrying...")
# Re-sign all items with fresh tokens
try:
import planetary_computer
s2_items = [planetary_computer.sign(item) for item in s2_items]
print(f" Re-signed {len(s2_items)} items")
if attempt < max_retries - 1:
import time
time.sleep(1) # Brief pause before retry
continue
else:
raise HTTPException(
status_code=500,
detail=f"❌ Authentication failed after {max_retries} attempts.\n\n"
f"Azure Blob Storage authentication is failing. This may be due to:\n"
f"1. Expired or invalid SAS tokens from Planetary Computer\n"
f"2. Planetary Computer API service issues\n"
f"3. System clock synchronization issues (check date/time)\n\n"
f"Please try again in a few minutes. If the problem persists,\n"
f"Microsoft Planetary Computer may be experiencing issues.\n\n"
f"Error: {error_msg[:300]}"
)
except Exception as resign_error:
print(f" Re-signing failed: {resign_error}")
if attempt == max_retries - 1:
raise HTTPException(
status_code=500,
detail=f"Failed to re-sign authentication tokens: {str(resign_error)}"
)
elif is_access_error and len(s2_items) > 5:
# If we have access errors and multiple scenes, try with fewer scenes
print(f"⚠️ Data access error detected. Reducing scene count and retrying...")
# Remove last 20% of scenes and try again
remove_count = max(1, len(s2_items) // 5)
s2_items = s2_items[:-remove_count]
print(f" Reduced to {len(s2_items)} scenes")
if attempt < max_retries - 1:
continue
else:
# Final attempt with reduced scenes
print(f" Final attempt with {len(s2_items)} scenes...")
elif is_network_error:
print(f"⚠️ Network error on attempt {attempt + 1}: {error_msg[:100]}")
# Re-check connectivity on network errors
if "Could not resolve host" in error_msg:
is_conn, conn_msg = check_internet_connectivity()
if not is_conn:
raise HTTPException(
status_code=503,
detail=f"❌ Network connectivity lost!\n\n"
f"{conn_msg}\n\n"
f"Please:\n"
f"1. Check your internet connection\n"
f"2. Verify DNS settings (try 8.8.8.8 or 1.1.1.1)\n"
f"3. Check firewall/proxy settings\n"
f"4. Try again once connection is restored\n\n"
f"Original error: {error_msg[:150]}"
)
if attempt < max_retries - 1:
import time
print(f" Retrying in {retry_delay} seconds...")
@@ -3994,13 +4263,30 @@ async def ndvi_predict_timeseries(config: NDVIPredictionConfig):
else:
raise HTTPException(
status_code=503,
detail=f"Network error: Unable to download Sentinel-2 data after {max_retries} attempts. "
f"Please check your internet connection or try again later. "
detail=f"❌ Network error persists after {max_retries} attempts.\n\n"
f"Unable to download Sentinel-2 data from Microsoft Planetary Computer.\n\n"
f"Possible causes:\n"
f"1. Unstable internet connection\n"
f"2. DNS resolution issues\n"
f"3. Planetary Computer service temporarily unavailable\n"
f"4. Firewall blocking access to sentinel2l2a01.blob.core.windows.net\n\n"
f"Please check your network and try again.\n\n"
f"Error: {error_msg[:200]}"
)
else:
# Non-network error, raise immediately
raise
# If still failing after retries, try with minimal data
if attempt == max_retries - 1 and len(s2_items) > 3:
print(f"⚠️ Multiple failures. Attempting with minimal scenes (3)...")
s2_items = s2_items[:3]
continue
else:
# Last resort: provide meaningful error
raise HTTPException(
status_code=500,
detail=f"Failed to load Sentinel-2 data. This may be due to expired access tokens "
f"or unavailable data. Please try a different time range or area. "
f"Error: {error_msg[:300]}"
)
if s2_data is None:
raise HTTPException(status_code=500, detail="Failed to load Sentinel-2 data")
@@ -4365,17 +4651,104 @@ async def ndvi_forecast(config: NDVIForecastConfig):
print(f"✅ Found {len(s2_items)} historical scenes")
# Load data
s2_data = load(
s2_items,
bbox=bbox,
bands=["B02", "B03", "B04", "B05", "B08", "B11", "SCL"],
chunks={"time": 1, "x": 2048, "y": 2048},
groupby="solar_day",
resolution=config.resolution
).compute()
# Load data with retry logic for network and access errors
max_retries = 3
retry_delay = 2
s2_data = None
print(f"✅ Loaded {len(s2_data.time)} time steps")
for attempt in range(max_retries):
try:
print(f"📥 Loading Sentinel-2 data (attempt {attempt + 1}/{max_retries})...")
s2_data = load(
s2_items,
bbox=bbox,
bands=["B02", "B03", "B04", "B05", "B08", "B11", "SCL"],
chunks={"time": 1, "x": 2048, "y": 2048},
groupby="solar_day",
resolution=config.resolution
).compute()
print(f"✅ Loaded {len(s2_data.time)} time steps")
break
except Exception as e:
error_msg = str(e)
# Check for various error types
is_network_error = any(err in error_msg for err in [
"Could not resolve host", "CURL error", "Connection", "Timeout"
])
is_access_error = any(err in error_msg for err in [
"does not exist in the file system", "RasterioIOError",
"not recognized as a supported dataset", "Aborting load"
])
is_auth_error = any(err in error_msg for err in [
"AuthenticationFailed", "Signature fields not well formed",
"Server failed to authenticate", "403", "401"
])
if is_auth_error:
print(f"⚠️ Authentication error. Re-signing items...")
try:
import planetary_computer
s2_items = [planetary_computer.sign(item) for item in s2_items]
if attempt < max_retries - 1:
import time
time.sleep(1)
continue
else:
raise HTTPException(
status_code=500,
detail=f"Authentication failed: {error_msg[:200]}"
)
except Exception as resign_error:
if attempt == max_retries - 1:
raise HTTPException(
status_code=500,
detail=f"Failed to re-sign tokens: {str(resign_error)}"
)
elif is_access_error and len(s2_items) > 5:
# If we have access errors and multiple scenes, try with fewer scenes
print(f"⚠️ Data access error detected. Reducing scene count and retrying...")
remove_count = max(1, len(s2_items) // 5)
s2_items = s2_items[:-remove_count]
print(f" Reduced to {len(s2_items)} scenes")
if attempt < max_retries - 1:
continue
else:
print(f" Final attempt with {len(s2_items)} scenes...")
elif is_network_error:
print(f"⚠️ Network error on attempt {attempt + 1}: {error_msg[:100]}")
if attempt < max_retries - 1:
import time
print(f" Retrying in {retry_delay} seconds...")
time.sleep(retry_delay)
retry_delay *= 2
else:
raise HTTPException(
status_code=503,
detail=f"Network error: Unable to download Sentinel-2 data after {max_retries} attempts. "
f"Error: {error_msg[:200]}"
)
else:
# Try with minimal scenes on final attempt
if attempt == max_retries - 1 and len(s2_items) > 3:
print(f"⚠️ Multiple failures. Attempting with minimal scenes (3)...")
s2_items = s2_items[:3]
continue
else:
raise HTTPException(
status_code=500,
detail=f"Failed to load Sentinel-2 data. This may be due to expired access tokens "
f"or unavailable data. Please try a different time range or area. "
f"Error: {error_msg[:300]}"
)
if s2_data is None:
raise HTTPException(status_code=500, detail="Failed to load Sentinel-2 data")
# Calculate spectral indices
print(f"\n📊 Calculating spectral indices...")
@@ -4476,158 +4849,158 @@ async def ndvi_forecast(config: NDVIForecastConfig):
if len(point_features) == 0:
print("[NDVI FORECAST][FALLBACK] No valid classification points. Falling back to simple seasonal forecast.")
force_simple_forecast = True
print(f"✅ Extracted features for {len(point_features)} valid points")
# Classify points
X = np.array(point_features)
predictions = model.predict(X)
if label_encoder is not None:
try:
predictions = label_encoder.inverse_transform(predictions.astype(int))
except:
pass
# Count land types
unique_types, type_counts = np.unique(predictions, return_counts=True)
print(f"\n📊 Detected land types:")
for land_type, count in zip(unique_types, type_counts):
print(f" Type {land_type}: {count} points ({count/len(predictions)*100:.1f}%)")
# Calculate land-type-specific seasonal patterns
print(f"\n📈 Calculating land-type-specific seasonal patterns...")
land_type_patterns = {}
for time_idx in range(len(s2_data.time)):
time_val = pd.Timestamp(s2_data.time.values[time_idx])
month = time_val.month
# Get valid pixels for this time step
if hasattr(cloud_mask, 'isel'):
mask_t = cloud_mask.isel(time=time_idx)
else:
# Wrap numpy mask to xarray with same coords/dims as ndvi slice
ndvi_slice = ndvi.isel(time=time_idx)
mask_t = xr.DataArray(mask_t := cloud_mask[time_idx], coords=ndvi_slice.coords, dims=ndvi_slice.dims)
ndvi_t = ndvi.isel(time=time_idx).where(mask_t)
ndwi_t = ndwi.isel(time=time_idx).where(mask_t)
ndbi_t = ndbi.isel(time=time_idx).where(mask_t)
evi_t = evi.isel(time=time_idx).where(mask_t)
# Extract values at classified points
for point_idx, (lat, lon) in enumerate(point_coords):
land_type = predictions[point_idx]
try:
ndvi_val = float(ndvi_t.sel(y=lat, x=lon, method='nearest').values)
if not np.isnan(ndvi_val):
ndwi_val = float(ndwi_t.sel(y=lat, x=lon, method='nearest').values)
ndbi_val = float(ndbi_t.sel(y=lat, x=lon, method='nearest').values)
evi_val = float(evi_t.sel(y=lat, x=lon, method='nearest').values)
# Initialize land type if not exists
if land_type not in land_type_patterns:
land_type_patterns[land_type] = {}
if month not in land_type_patterns[land_type]:
land_type_patterns[land_type][month] = {
'ndvi': [], 'ndwi': [], 'ndbi': [], 'evi': []
}
# Append values
land_type_patterns[land_type][month]['ndvi'].append(ndvi_val)
land_type_patterns[land_type][month]['ndwi'].append(ndwi_val)
land_type_patterns[land_type][month]['ndbi'].append(ndbi_val)
land_type_patterns[land_type][month]['evi'].append(evi_val)
except:
continue
# Calculate statistics for each land type and month
land_type_seasonal_stats = {}
for land_type, month_data in land_type_patterns.items():
land_type_seasonal_stats[land_type] = {}
for month, values in month_data.items():
ndvi_vals = values['ndvi']
if len(ndvi_vals) > 0:
land_type_seasonal_stats[land_type][month] = {
'ndvi_mean': float(np.mean(ndvi_vals)),
'ndvi_min': float(np.min(ndvi_vals)),
'ndvi_max': float(np.max(ndvi_vals)),
'ndvi_std': float(np.std(ndvi_vals)),
'ndvi_range': float(np.max(ndvi_vals) - np.min(ndvi_vals)),
'ndwi_mean': float(np.mean(values['ndwi'])),
'ndbi_mean': float(np.mean(values['ndbi'])),
'evi_mean': float(np.mean(values['evi'])),
'n_samples': len(ndvi_vals)
}
print(f"✅ Calculated patterns for {len(land_type_seasonal_stats)} land types")
if len(land_type_seasonal_stats) == 0:
print("[NDVI FORECAST][FALLBACK] No seasonal patterns. Falling back to simple seasonal forecast.")
force_simple_forecast = True
# Generate forecast using land-type-weighted average
if force_simple_forecast:
print("[NDVI FORECAST][FALLBACK] Skipping ML forecast, will use simple seasonal averaging.")
else:
print(f"\n🔮 Generating land-type-specific forecast...")
print(f"✅ Extracted features for {len(point_features)} valid points")
forecast_timeseries = []
current_date = forecast_start
# Classify points
X = np.array(point_features)
predictions = model.predict(X)
# Calculate land type weights
total_points = len(predictions)
if total_points == 0:
print("[NDVI FORECAST][FALLBACK] Zero valid classified points. Switching to simple seasonal forecast.")
force_simple_forecast = True
else:
land_type_weights = {lt: np.sum(predictions == lt) / total_points
for lt in unique_types}
if label_encoder is not None:
try:
predictions = label_encoder.inverse_transform(predictions.astype(int))
except:
pass
if not force_simple_forecast:
while current_date <= forecast_end:
month = current_date.month
# Aggregate forecast across all land types (weighted)
weighted_forecast = {
'ndvi_mean': 0, 'ndvi_min': 0, 'ndvi_max': 0, 'ndvi_std': 0,
'ndvi_range': 0, 'ndwi_mean': 0, 'ndbi_mean': 0, 'evi_mean': 0
}
land_type_contributions = {}
for land_type, weight in land_type_weights.items():
if land_type in land_type_seasonal_stats and month in land_type_seasonal_stats[land_type]:
stats = land_type_seasonal_stats[land_type][month]
land_type_contributions[int(land_type)] = {
**stats,
'weight': float(weight)
}
for key in weighted_forecast:
weighted_forecast[key] += stats[key] * weight
if land_type_contributions:
forecast_data = {
'date': current_date.strftime('%Y-%m-%d'),
'is_forecast': True,
'land_type_specific': land_type_contributions,
**weighted_forecast
}
forecast_timeseries.append(forecast_data)
current_date += relativedelta(months=1)
# Count land types
unique_types, type_counts = np.unique(predictions, return_counts=True)
print(f"\n📊 Detected land types:")
for land_type, count in zip(unique_types, type_counts):
print(f" Type {land_type}: {count} points ({count/len(predictions)*100:.1f}%)")
# Calculate land-type-specific seasonal patterns
print(f"\n📈 Calculating land-type-specific seasonal patterns...")
land_type_patterns = {}
for time_idx in range(len(s2_data.time)):
time_val = pd.Timestamp(s2_data.time.values[time_idx])
month = time_val.month
method_used = "Land-Type-Specific Forecasting (ML-Enhanced)"
# Get valid pixels for this time step
if hasattr(cloud_mask, 'isel'):
mask_t = cloud_mask.isel(time=time_idx)
else:
# Wrap numpy mask to xarray with same coords/dims as ndvi slice
ndvi_slice = ndvi.isel(time=time_idx)
mask_t = xr.DataArray(mask_t := cloud_mask[time_idx], coords=ndvi_slice.coords, dims=ndvi_slice.dims)
ndvi_t = ndvi.isel(time=time_idx).where(mask_t)
ndwi_t = ndwi.isel(time=time_idx).where(mask_t)
ndbi_t = ndbi.isel(time=time_idx).where(mask_t)
evi_t = evi.isel(time=time_idx).where(mask_t)
# Extract values at classified points
for point_idx, (lat, lon) in enumerate(point_coords):
land_type = predictions[point_idx]
try:
ndvi_val = float(ndvi_t.sel(y=lat, x=lon, method='nearest').values)
if not np.isnan(ndvi_val):
ndwi_val = float(ndwi_t.sel(y=lat, x=lon, method='nearest').values)
ndbi_val = float(ndbi_t.sel(y=lat, x=lon, method='nearest').values)
evi_val = float(evi_t.sel(y=lat, x=lon, method='nearest').values)
# Initialize land type if not exists
if land_type not in land_type_patterns:
land_type_patterns[land_type] = {}
if month not in land_type_patterns[land_type]:
land_type_patterns[land_type][month] = {
'ndvi': [], 'ndwi': [], 'ndbi': [], 'evi': []
}
# Append values
land_type_patterns[land_type][month]['ndvi'].append(ndvi_val)
land_type_patterns[land_type][month]['ndwi'].append(ndwi_val)
land_type_patterns[land_type][month]['ndbi'].append(ndbi_val)
land_type_patterns[land_type][month]['evi'].append(evi_val)
except:
continue
# Calculate statistics for each land type and month
land_type_seasonal_stats = {}
for land_type, month_data in land_type_patterns.items():
land_type_seasonal_stats[land_type] = {}
for month, values in month_data.items():
ndvi_vals = values['ndvi']
if len(ndvi_vals) > 0:
land_type_seasonal_stats[land_type][month] = {
'ndvi_mean': float(np.mean(ndvi_vals)),
'ndvi_min': float(np.min(ndvi_vals)),
'ndvi_max': float(np.max(ndvi_vals)),
'ndvi_std': float(np.std(ndvi_vals)),
'ndvi_range': float(np.max(ndvi_vals) - np.min(ndvi_vals)),
'ndwi_mean': float(np.mean(values['ndwi'])),
'ndbi_mean': float(np.mean(values['ndbi'])),
'evi_mean': float(np.mean(values['evi'])),
'n_samples': len(ndvi_vals)
}
print(f"✅ Calculated patterns for {len(land_type_seasonal_stats)} land types")
if len(land_type_seasonal_stats) == 0:
print("[NDVI FORECAST][FALLBACK] No seasonal patterns. Falling back to simple seasonal forecast.")
force_simple_forecast = True
# Generate forecast using land-type-weighted average
if force_simple_forecast:
print("[NDVI FORECAST][FALLBACK] Skipping ML forecast, will use simple seasonal averaging.")
else:
print(f"\n🔮 Generating land-type-specific forecast...")
forecast_timeseries = []
current_date = forecast_start
# Calculate land type weights
total_points = len(predictions)
if total_points == 0:
print("[NDVI FORECAST][FALLBACK] Zero valid classified points. Switching to simple seasonal forecast.")
force_simple_forecast = True
else:
land_type_weights = {lt: np.sum(predictions == lt) / total_points
for lt in unique_types}
if not force_simple_forecast:
while current_date <= forecast_end:
month = current_date.month
# Aggregate forecast across all land types (weighted)
weighted_forecast = {
'ndvi_mean': 0, 'ndvi_min': 0, 'ndvi_max': 0, 'ndvi_std': 0,
'ndvi_range': 0, 'ndwi_mean': 0, 'ndbi_mean': 0, 'evi_mean': 0
}
land_type_contributions = {}
for land_type, weight in land_type_weights.items():
if land_type in land_type_seasonal_stats and month in land_type_seasonal_stats[land_type]:
stats = land_type_seasonal_stats[land_type][month]
land_type_contributions[int(land_type)] = {
**stats,
'weight': float(weight)
}
for key in weighted_forecast:
weighted_forecast[key] += stats[key] * weight
if land_type_contributions:
forecast_data = {
'date': current_date.strftime('%Y-%m-%d'),
'is_forecast': True,
'land_type_specific': land_type_contributions,
**weighted_forecast
}
forecast_timeseries.append(forecast_data)
current_date += relativedelta(months=1)
method_used = "Land-Type-Specific Forecasting (ML-Enhanced)"
if force_simple_forecast or not use_ml_classification:
# Simple seasonal averaging (fallback)