diff --git a/api_server.py b/api_server.py index 6bba65e..b1bd909 100644 --- a/api_server.py +++ b/api_server.py @@ -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("

Change Detection Interface not found.

") +@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)