Malek1ab/ndvi-api
0
1import os2import tempfile3import base644import logging5from typing import Dict, Any, Optional, List, Tuple6from datetime import datetime, timedelta7import io8 9# Configuration pour éviter les problèmes de permissions10import matplotlib11matplotlib.use('Agg')12os.environ['MPLCONFIGDIR'] = '/tmp/matplotlib'13 14from fastapi import FastAPI, File, UploadFile, HTTPException, Query15from fastapi.responses import JSONResponse, StreamingResponse16from fastapi.middleware.cors import CORSMiddleware17import geopandas as gpd18import numpy as np19from sentinelhub import SHConfig, Geometry, CRS, SentinelHubRequest, DataCollection, MimeType20import matplotlib.pyplot as plt21import matplotlib.animation as animation22from matplotlib import ticker23from rasterio.features import rasterize24from rasterio.transform import from_bounds25from shapely.geometry import mapping26import matplotlib.colors as mcolors27from scipy.ndimage import gaussian_filter, median_filter, uniform_filter28from scipy.interpolate import griddata29from matplotlib.patches import Rectangle30 31# Tentative d'importation des bibliothèques optionnelles pour l'amélioration32try:33 import cv234 CV2_AVAILABLE = True35except ImportError:36 CV2_AVAILABLE = False37 logging.warning("OpenCV non disponible, certaines améliorations seront limitées")38 39try:40 from skimage import filters, restoration, morphology41 from skimage.segmentation import slic42 from skimage.transform import resize43 SKIMAGE_AVAILABLE = True44except ImportError:45 SKIMAGE_AVAILABLE = False46 logging.warning("scikit-image non disponible, certaines améliorations seront limitées")47 48# Configuration du logging49logging.basicConfig(level=logging.INFO)50logger = logging.getLogger(__name__)51 52app = FastAPI(53 title="API Indices de Végétation Satellite - Version Améliorée avec Séries Temporelles",54 description="Affichage ultra-clair des indices de végétation avec amélioration des pixels et animations temporelles",55 version="3.0.0"56)57 58app.add_middleware(59 CORSMiddleware,60 allow_origins=["*"],61 allow_credentials=True,62 allow_methods=["*"],63 allow_headers=["*"],64)65 66def get_sentinelhub_config() -> SHConfig:67 config = SHConfig()68 config.instance_id = os.getenv("SENTINELHUB_INSTANCE_ID", "4274ebc8-62a9-4351-a9dc-3f0d1999071b")69 config.sh_client_id = os.getenv("SENTINELHUB_CLIENT_ID", "faf6f181-fea9-4d0e-96f3-d7aae3831e9d")70 config.sh_client_secret = os.getenv("SENTINELHUB_CLIENT_SECRET", "2teo28ufFIHljw2w896bwMgpCgH8E3PF")71 return config72 73def validate_geojson(gdf: gpd.GeoDataFrame) -> None:74 if gdf.empty:75 raise HTTPException(status_code=400, detail="Fichier GeoJSON vide")76 77 if len(gdf) == 0:78 raise HTTPException(status_code=400, detail="Aucune géométrie trouvée dans le GeoJSON")79 80 if gdf.geometry.isna().all():81 raise HTTPException(status_code=400, detail="Aucune géométrie valide trouvée")82 83 valid_geoms = gdf.geometry.dropna()84 if len(valid_geoms) == 0:85 raise HTTPException(status_code=400, detail="Aucune géométrie valide trouvée")86 87def create_mask_from_geometry(geometry, image_shape, transform):88 """Crée un masque à partir de la géométrie pour rendre transparent l'extérieur"""89 try:90 geom_dict = mapping(geometry)91 92 mask = rasterize(93 [geom_dict],94 out_shape=image_shape,95 transform=transform,96 fill=0,97 default_value=1,98 dtype=np.uint899 )100 101 return mask.astype(bool)102 except Exception as e:103 logger.error(f"Erreur lors de la création du masque: {e}")104 return np.ones(image_shape, dtype=bool)105 106# NOUVELLES FONCTIONS DU CODE 2107 108def transparent_cmap(cmap, alpha=0.5):109 """Créer une colormap transparente"""110 colors = cmap(np.arange(cmap.N))111 colors[:, -1] = alpha112 return mcolors.ListedColormap(colors)113 114def s2_to_rgb(bands_data, contrast_factor=3.5):115 """Convertir les bandes Sentinel-2 en image RGB"""116 if bands_data.shape[-1] >= 3:117 rgb = bands_data[:, :, [2, 1, 0]] # R, G, B118 else:119 return bands_data120 121 # Normaliser les valeurs122 rgb = np.clip(rgb * contrast_factor / 10000, 0, 1)123 return rgb124 125def simple_vegetation_segmentation(ndvi, threshold=0.3):126 """Segmentation simple basée sur un seuil NDVI"""127 return np.where(ndvi > threshold, 1, 0)128 129def calculate_masked_average(index_array, mask):130 """Calculer la moyenne de l'indice dans la zone masquée"""131 masked_index = index_array * mask132 return np.nanmean(masked_index[mask == 1])133 134def get_time_series_images(geometry, start_date, end_date, config, interval_days=10, max_cloud_coverage=30):135 """Obtenir une série temporelle d'images Sentinel-2"""136 images = []137 dates = []138 139 current_date = start_date140 while current_date < end_date:141 next_date = min(current_date + timedelta(days=interval_days), end_date)142 143 try:144 evalscript = """145 //VERSION=3146 function setup() {147 return {148 input: ["B02", "B03", "B04", "B08", "CLM"],149 output: [150 {id: "bands", bands: 4, sampleType: "FLOAT32"},151 {id: "cloud_mask", bands: 1, sampleType: "UINT8"}152 ]153 };154 }155 156 function evaluatePixel(sample) {157 return {158 bands: [sample.B02, sample.B03, sample.B04, sample.B08],159 cloud_mask: [sample.CLM]160 };161 }162 """163 164 request = SentinelHubRequest(165 evalscript=evalscript,166 input_data=[SentinelHubRequest.input_data(167 data_collection=DataCollection.SENTINEL2_L2A,168 time_interval=(current_date.strftime('%Y-%m-%d'), next_date.strftime('%Y-%m-%d')),169 maxcc=max_cloud_coverage/100,170 mosaicking_order='leastCC'171 )],172 responses=[173 SentinelHubRequest.output_response('bands', MimeType.TIFF),174 SentinelHubRequest.output_response('cloud_mask', MimeType.TIFF)175 ],176 geometry=geometry,177 size=(512, 512),178 config=config179 )180 181 data = request.get_data()182 183 if data and len(data) > 0:184 bands = data[0]['bands']185 cloud_mask = data[0]['cloud_mask']186 187 # Appliquer le masque de nuages188 bands_masked = np.where(cloud_mask[..., 0] == 0, bands, np.nan)189 images.append(bands_masked)190 dates.append(current_date)191 logger.info(f"Image récupérée pour {current_date.strftime('%Y-%m-%d')}")192 193 except Exception as e:194 logger.error(f"Erreur lors du chargement de l'image {current_date}: {e}")195 196 current_date = next_date197 198 return images, dates199 200# FIN DES NOUVELLES FONCTIONS201 202def calculate_vegetation_indices_improved(bands_img: np.ndarray) -> Dict[str, np.ndarray]:203 """Calcul amélioré des indices de végétation avec gestion des erreurs"""204 logger.info(f"Calcul des indices pour une image de taille: {bands_img.shape}")205 206 blue = bands_img[:, :, 0].astype('float64')207 green = bands_img[:, :, 1].astype('float64') 208 red = bands_img[:, :, 2].astype('float64')209 nir = bands_img[:, :, 3].astype('float64')210 211 # Normaliser les valeurs si nécessaire212 if np.max(bands_img) > 1.5:213 logger.info("Normalisation des bandes detectée nécessaire")214 blue = blue / 10000.0215 green = green / 10000.0216 red = red / 10000.0217 nir = nir / 10000.0218 219 indices = {}220 221 # Calcul avec votre formule exacte222 with np.errstate(divide='ignore', invalid='ignore'):223 # NDVI: (NIR - Red) / (NIR + Red)224 ndvi = (nir - red) / (nir + red)225 ndvi = np.clip(ndvi, -1, 1)226 ndvi = np.where(np.isfinite(ndvi), ndvi, np.nan)227 indices['ndvi'] = ndvi228 229 # NDWI: (NIR - Green) / (NIR + Green)230 ndwi = (nir - green) / (nir + green)231 ndwi = np.clip(ndwi, -1, 1)232 ndwi = np.where(np.isfinite(ndwi), ndwi, np.nan)233 indices['ndwi'] = ndwi234 235 # SAVI: ((NIR - Red) / (NIR + Red + 0.5)) * 1.5236 savi = ((nir - red) / (nir + red + 0.5)) * 1.5237 savi = np.clip(savi, -1, 1)238 savi = np.where(np.isfinite(savi), savi, np.nan)239 indices['savi'] = savi240 241 # EVI: 2.5 * (NIR - Red) / (NIR + 6*Red - 7.5*Blue + 1)242 evi = 2.5 * (nir - red) / (nir + 6 * red - 7.5 * blue + 1)243 evi = np.clip(evi, -1, 1)244 evi = np.where(np.isfinite(evi), evi, np.nan)245 indices['evi'] = evi246 247 logger.info(f"Indices calculés: {list(indices.keys())}")248 return indices249 250# NOUVELLES FONCTIONS D'AMÉLIORATION DES PIXELS251 252def advanced_pixel_enhancement(array: np.ndarray, mask: np.ndarray, enhancement_type: str = 'adaptive') -> np.ndarray:253 """254 Amélioration avancée des pixels avec plusieurs techniques255 """256 enhanced_array = np.where(mask, array, np.nan)257 258 if enhancement_type == 'adaptive':259 enhanced_array = adaptive_smoothing(enhanced_array)260 elif enhancement_type == 'super_resolution':261 enhanced_array = super_resolution_upscale(enhanced_array)262 elif enhancement_type == 'edge_preserving' and CV2_AVAILABLE:263 enhanced_array = edge_preserving_filter(enhanced_array)264 elif enhancement_type == 'bilateral':265 enhanced_array = bilateral_filter_2d(enhanced_array)266 elif enhancement_type == 'segmentation_based' and SKIMAGE_AVAILABLE:267 enhanced_array = segmentation_based_enhancement(enhanced_array)268 else:269 # Fallback vers lissage gaussien amélioré270 enhanced_array = enhanced_gaussian_filter(enhanced_array)271 272 return enhanced_array273 274def adaptive_smoothing(array: np.ndarray, window_size: int = 5) -> np.ndarray:275 """Lissage adaptatif basé sur la variance locale"""276 from scipy.ndimage import generic_filter277 278 def adaptive_mean(values):279 valid_values = values[~np.isnan(values)]280 if len(valid_values) < 3:281 return np.nanmean(values)282 283 variance = np.var(valid_values)284 285 if variance < 0.01: # Zone homogène286 return np.nanmean(values)287 else: # Zone avec détails288 return np.nanmedian(values)289 290 enhanced = generic_filter(array, adaptive_mean, size=window_size, mode='constant', cval=np.nan)291 return enhanced292 293def super_resolution_upscale(array: np.ndarray, scale_factor: int = 2) -> np.ndarray:294 """Super-résolution par interpolation bicubique avancée"""295 h, w = array.shape296 297 y_orig, x_orig = np.mgrid[0:h, 0:w]298 299 h_new, w_new = h * scale_factor, w * scale_factor300 y_new = np.linspace(0, h-1, h_new)301 x_new = np.linspace(0, w-1, w_new)302 x_new_grid, y_new_grid = np.meshgrid(x_new, y_new)303 304 valid_mask = ~np.isnan(array)305 306 if np.sum(valid_mask) < 4:307 return array308 309 try:310 points = np.column_stack((y_orig[valid_mask], x_orig[valid_mask]))311 values = array[valid_mask]312 313 enhanced = griddata(points, values, (y_new_grid, x_new_grid), method='cubic', fill_value=np.nan)314 315 nan_mask = np.isnan(enhanced)316 if np.any(nan_mask):317 enhanced_linear = griddata(points, values, (y_new_grid, x_new_grid), method='linear', fill_value=np.nan)318 enhanced[nan_mask] = enhanced_linear[nan_mask]319 320 # Redimensionner à la taille originale321 if SKIMAGE_AVAILABLE:322 enhanced_downsampled = resize(enhanced, (h, w), order=3, mode='constant', cval=np.nan, anti_aliasing=True)323 else:324 # Fallback simple325 enhanced_downsampled = gaussian_filter(array, sigma=0.8)326 327 return enhanced_downsampled328 329 except Exception as e:330 logger.warning(f"Super-résolution échouée: {e}")331 return gaussian_filter(array, sigma=0.8)332 333def edge_preserving_filter(array: np.ndarray, sigma_spatial: float = 1.0, sigma_intensity: float = 0.1) -> np.ndarray:334 """Filtre préservant les contours avec OpenCV"""335 if not CV2_AVAILABLE:336 return gaussian_filter(array, sigma=1.0)337 338 valid_mask = ~np.isnan(array)339 340 if np.sum(valid_mask) == 0:341 return array342 343 array_norm = array.copy()344 array_norm[~valid_mask] = 0345 346 min_val, max_val = np.nanmin(array), np.nanmax(array)347 if max_val > min_val:348 array_norm = (array_norm - min_val) / (max_val - min_val)349 350 array_uint8 = (array_norm * 255).astype(np.uint8)351 352 try:353 enhanced_uint8 = cv2.edgePreservingFilter(array_uint8, flags=2, sigma_s=sigma_spatial*50, sigma_r=sigma_intensity)354 355 enhanced = enhanced_uint8.astype(np.float64) / 255.0356 357 if max_val > min_val:358 enhanced = enhanced * (max_val - min_val) + min_val359 360 enhanced[~valid_mask] = np.nan361 362 return enhanced363 364 except Exception as e:365 logger.warning(f"Edge-preserving filter échoué: {e}")366 return gaussian_filter(array, sigma=1.0)367 368def bilateral_filter_2d(array: np.ndarray, sigma_spatial: float = 1.5, sigma_intensity: float = 0.2) -> np.ndarray:369 """Filtre bilatéral 2D simplifié"""370 # Version simplifiée pour éviter la complexité excessive371 return gaussian_filter(array, sigma=sigma_spatial)372 373def segmentation_based_enhancement(array: np.ndarray, n_segments: int = 100) -> np.ndarray:374 """Amélioration basée sur la segmentation SLIC"""375 if not SKIMAGE_AVAILABLE:376 return gaussian_filter(array, sigma=1.0)377 378 valid_mask = ~np.isnan(array)379 380 if np.sum(valid_mask) < 100:381 return gaussian_filter(array, sigma=1.0)382 383 array_for_slic = np.stack([array, array, array], axis=2)384 array_for_slic[~valid_mask] = 0385 386 try:387 segments = slic(array_for_slic, n_segments=n_segments, compactness=10, sigma=1, start_label=1)388 389 enhanced = array.copy()390 391 for segment_id in np.unique(segments):392 if segment_id == 0:393 continue394 395 segment_mask = segments == segment_id396 segment_values = array[segment_mask & valid_mask]397 398 if len(segment_values) > 0:399 segment_mean = np.mean(segment_values)400 enhanced[segment_mask] = segment_mean401 402 return enhanced403 404 except Exception as e:405 logger.warning(f"Segmentation échouée: {e}")406 return gaussian_filter(array, sigma=1.0)407 408def enhanced_gaussian_filter(array: np.ndarray, sigma: float = 1.0) -> np.ndarray:409 """Filtre gaussien amélioré par défaut"""410 # Appliquer un filtre gaussien avec préservation des NaN411 valid_mask = ~np.isnan(array)412 413 if np.sum(valid_mask) == 0:414 return array415 416 # Remplacer temporairement les NaN par la moyenne417 temp_array = array.copy()418 temp_array[~valid_mask] = np.nanmean(array)419 420 # Appliquer le filtre421 filtered = gaussian_filter(temp_array, sigma=sigma)422 423 # Remettre les NaN424 filtered[~valid_mask] = np.nan425 426 return filtered427 428def create_smooth_interpretation_map(index_array: np.ndarray, index_name: str) -> Dict[str, Any]:429 """Crée une carte d'interprétation avec classifications claires"""430 431 interpretation_configs = {432 'ndvi': {433 'thresholds': [-1, -0.2, 0.0, 0.2, 0.4, 0.6, 1.0],434 'labels': ['Eau/Roches', 'Sol nu', 'Végétation très faible', 'Végétation clairsemée', 'Végétation modérée', 'Végétation dense'],435 'colors': ['#0066cc', '#8B4513', '#FFFF99', '#ADFF2F', '#32CD32', '#006400'],436 'description': 'Plus la valeur est élevée, plus la végétation est dense et saine'437 },438 'ndwi': {439 'thresholds': [-1, -0.3, -0.1, 0.1, 0.3, 1.0],440 'labels': ['Sol très sec', 'Sol sec', 'Sol modérément humide', 'Sol humide', 'Eau/Zone très humide'],441 'colors': ['#8B4513', '#DEB887', '#98FB98', '#87CEEB', '#0066cc'],442 'description': 'Indique la présence d\'eau et le niveau d\'humidité'443 },444 'savi': {445 'thresholds': [-1, -0.1, 0.1, 0.3, 0.5, 0.8, 1.0],446 'labels': ['Non végétalisé', 'Sol nu/rocailleux', 'Végétation émergente', 'Végétation faible', 'Végétation modérée', 'Végétation dense'],447 'colors': ['#8B4513', '#D2691E', '#FFFF99', '#ADFF2F', '#32CD32', '#006400'],448 'description': 'NDVI ajusté pour réduire l\'influence de la couleur du sol'449 },450 'evi': {451 'thresholds': [-1, -0.1, 0.1, 0.3, 0.5, 0.7, 1.0],452 'labels': ['Aucune végétation', 'Végétation minimale', 'Végétation faible', 'Végétation modérée', 'Végétation dense', 'Végétation très dense'],453 'colors': ['#8B4513', '#DEB887', '#FFFF99', '#ADFF2F', '#32CD32', '#006400'],454 'description': 'Version améliorée du NDVI, plus précise pour la végétation dense'455 }456 }457 458 config = interpretation_configs.get(index_name, interpretation_configs['ndvi'])459 460 classified_map = np.full_like(index_array, np.nan)461 462 for i, threshold in enumerate(config['thresholds'][:-1]):463 next_threshold = config['thresholds'][i + 1]464 mask = (index_array >= threshold) & (index_array < next_threshold)465 classified_map[mask] = i466 467 return {468 'classified_map': classified_map,469 'labels': config['labels'],470 'colors': config['colors'],471 'description': config['description'],472 'thresholds': config['thresholds']473 }474 475def generate_enhanced_visualization(indices: Dict[str, np.ndarray], mask: np.ndarray, geometry_bounds, enhancement_method: str = 'adaptive') -> Dict[str, str]:476 """Génère des visualisations avec amélioration des pixels"""477 478 plots = {}479 480 with tempfile.TemporaryDirectory() as tmpdir:481 for idx, array in indices.items():482 483 # Appliquer l'amélioration des pixels484 enhanced_array = advanced_pixel_enhancement(array, mask, enhancement_method)485 486 # Générer les versions : continue, classifiée, et comparaison487 for version in ['enhanced_continuous', 'enhanced_classified', 'comparison']:488 489 path = os.path.join(tmpdir, f"{idx}_{version}.png")490 491 try:492 if version == 'comparison':493 # Vue comparative494 fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(24, 7), facecolor='white')495 496 # Original497 masked_original = np.where(mask, array, np.nan)498 im1 = ax1.imshow(499 masked_original,500 cmap='RdYlGn' if idx == 'ndvi' else 'viridis',501 extent=[geometry_bounds[0], geometry_bounds[2], geometry_bounds[1], geometry_bounds[3]],502 interpolation='nearest',503 alpha=0.9504 )505 ax1.set_title(f'{idx.upper()} - Original', fontsize=14, fontweight='bold')506 plt.colorbar(im1, ax=ax1, shrink=0.7)507 508 # Amélioré509 im2 = ax2.imshow(510 enhanced_array,511 cmap='RdYlGn' if idx == 'ndvi' else 'viridis',512 extent=[geometry_bounds[0], geometry_bounds[2], geometry_bounds[1], geometry_bounds[3]],513 interpolation='bilinear',514 alpha=0.9515 )516 ax2.set_title(f'{idx.upper()} - Amélioré ({enhancement_method})', fontsize=14, fontweight='bold')517 plt.colorbar(im2, ax=ax2, shrink=0.7)518 519 # Différence520 difference = enhanced_array - masked_original521 im3 = ax3.imshow(522 difference,523 cmap='RdBu',524 extent=[geometry_bounds[0], geometry_bounds[2], geometry_bounds[1], geometry_bounds[3]],525 interpolation='bilinear',526 alpha=0.9527 )528 ax3.set_title(f'{idx.upper()} - Amélioration (Différence)', fontsize=14, fontweight='bold')529 plt.colorbar(im3, ax=ax3, shrink=0.7)530 531 for ax in [ax1, ax2, ax3]:532 ax.set_xlabel('Longitude', fontsize=11)533 ax.set_ylabel('Latitude', fontsize=11)534 ax.grid(True, alpha=0.3)535 536 else:537 # Versions continues et classifiées améliorées538 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 7), facecolor='white')539 540 if version == 'enhanced_continuous':541 # Version continue améliorée542 if idx == 'ndvi':543 cmap = 'RdYlGn'544 vmin, vmax = -0.2, 0.8545 elif idx == 'ndwi':546 cmap = 'Blues'547 vmin, vmax = -0.4, 0.4548 elif idx == 'savi':549 cmap = 'YlGn'550 vmin, vmax = -0.1, 0.7551 else: # evi552 cmap = 'Greens'553 vmin, vmax = -0.1, 0.8554 555 im1 = ax1.imshow(556 enhanced_array,557 cmap=cmap,558 vmin=vmin,559 vmax=vmax,560 extent=[geometry_bounds[0], geometry_bounds[2], geometry_bounds[1], geometry_bounds[3]],561 interpolation='bilinear',562 alpha=0.9563 )564 565 cbar1 = plt.colorbar(im1, ax=ax1, shrink=0.7)566 cbar1.set_label(f'{idx.upper()} Valeurs Améliorées', fontsize=11)567 568 ax1.set_title(f'{idx.upper()} - Continu Amélioré ({enhancement_method})', fontsize=14, fontweight='bold')569 570 # Version classifiée (toujours présente)571 interpretation = create_smooth_interpretation_map(enhanced_array, idx)572 classified_map = interpretation['classified_map']573 574 colors = interpretation['colors']575 n_classes = len(colors)576 cmap_classified = mcolors.ListedColormap(colors)577 bounds = list(range(n_classes + 1))578 norm = mcolors.BoundaryNorm(bounds, cmap_classified.N)579 580 im2 = ax2.imshow(581 classified_map,582 cmap=cmap_classified,583 norm=norm,584 extent=[geometry_bounds[0], geometry_bounds[2], geometry_bounds[1], geometry_bounds[3]],585 interpolation='nearest',586 alpha=0.9587 )588 589 cbar2 = plt.colorbar(im2, ax=ax2, shrink=0.7, ticks=range(n_classes))590 cbar2.ax.set_yticklabels(interpretation['labels'], fontsize=9)591 cbar2.set_label('Classification', fontsize=11)592 593 ax2.set_title(f'{idx.upper()} - Classification Améliorée', fontsize=14, fontweight='bold')594 595 for ax in [ax1, ax2]:596 ax.set_xlabel('Longitude', fontsize=11)597 ax.set_ylabel('Latitude', fontsize=11)598 ax.grid(True, alpha=0.3, linestyle='--', linewidth=0.5)599 ax.tick_params(axis='both', which='major', labelsize=9)600 601 # Ajouter des statistiques améliorées602 valid_data = enhanced_array[np.isfinite(enhanced_array)]603 if len(valid_data) > 0:604 mean_val = np.mean(valid_data)605 std_val = np.std(valid_data)606 min_val = np.min(valid_data)607 max_val = np.max(valid_data)608 609 info_text = f"""Statistiques Améliorées:610Moyenne: {mean_val:.3f}611Min: {min_val:.3f} | Max: {max_val:.3f}612Écart-type: {std_val:.3f}613Pixels valides: {len(valid_data)}614Méthode: {enhancement_method}"""615 616 props = dict(boxstyle='round,pad=0.5', facecolor='lightblue', alpha=0.9, edgecolor='navy')617 fig.text(0.02, 0.98, info_text, transform=fig.transFigure, fontsize=9,618 verticalalignment='top', bbox=props)619 620 plt.tight_layout()621 plt.subplots_adjust(bottom=0.1, top=0.92, left=0.1, right=0.95)622 623 plt.savefig(624 path,625 dpi=300,626 bbox_inches='tight',627 facecolor='white',628 edgecolor='none',629 format='png'630 )631 plt.close()632 633 with open(path, "rb") as f:634 plots[f"{idx}_{version}"] = base64.b64encode(f.read()).decode("utf-8")635 636 logger.info(f"Image {version} générée pour {idx} avec amélioration {enhancement_method}")637 638 except Exception as e:639 logger.error(f"Erreur lors de la génération du graphique {version} pour {idx}: {str(e)}")640 if 'fig' in locals():641 plt.close(fig)642 643 return plots644 645# NOUVELLE FONCTION POUR CRÉER UNE ANIMATION TEMPORELLE646def create_timelapse_animation(images: List[np.ndarray], dates: List[datetime], 647 index_name: str, mask: np.ndarray = None, 648 geometry_bounds: Tuple[float, float, float, float] = None) -> bytes:649 """Créer une animation GIF de la série temporelle"""650 651 # Calculer l'indice pour chaque image652 index_arrays = []653 for img in images:654 if img is not None:655 indices = calculate_vegetation_indices_improved(img)656 index_array = indices.get(index_name, indices['ndvi'])657 if mask is not None:658 index_array = np.where(mask, index_array, np.nan)659 index_arrays.append(index_array)660 661 if not index_arrays:662 logger.error("Aucune donnée pour créer l'animation")663 return None664 665 # Créer l'animation666 fig, ax = plt.subplots(figsize=(10, 10), facecolor='white')667 668 # Déterminer la plage de couleurs669 if index_name == 'ndvi':670 vmin, vmax = -0.2, 0.8671 cmap = 'RdYlGn'672 elif index_name == 'ndwi':673 vmin, vmax = -0.4, 0.4674 cmap = 'Blues'675 elif index_name == 'savi':676 vmin, vmax = -0.1, 0.7677 cmap = 'YlGn'678 else: # evi679 vmin, vmax = -0.1, 0.8680 cmap = 'Greens'681 682 # Première image683 if geometry_bounds:684 extent = [geometry_bounds[0], geometry_bounds[2], geometry_bounds[1], geometry_bounds[3]]685 else:686 extent = None687 688 im = ax.imshow(index_arrays[0], cmap=cmap, vmin=vmin, vmax=vmax, 689 extent=extent, interpolation='bilinear', alpha=0.9)690 691 # Titre et colorbar692 title = ax.set_title(f"{index_name.upper()} - {dates[0].strftime('%Y-%m-%d')}", 693 fontsize=14, fontweight='bold')694 cbar = fig.colorbar(im, ax=ax, shrink=0.7)695 cbar.set_label(f'{index_name.upper()} Valeurs', fontsize=11)696 697 ax.set_xlabel('Longitude', fontsize=11)698 ax.set_ylabel('Latitude', fontsize=11)699 ax.grid(True, alpha=0.3)700 701 def update(frame):702 im.set_array(index_arrays[frame])703 title.set_text(f"{index_name.upper()} - {dates[frame].strftime('%Y-%m-%d')}")704 return [im, title]705 706 ani = animation.FuncAnimation(707 fig, update, frames=len(index_arrays), 708 interval=500, blit=True, repeat=True709 )710 711 # Sauvegarder dans un buffer712 buffer = io.BytesIO()713 ani.save(buffer, writer='pillow', fps=2, format='gif')714 buffer.seek(0)715 716 plt.close()717 718 return buffer.getvalue()719 720# NOUVELLE FONCTION POUR CRÉER UN GRAPHIQUE DE SÉRIE TEMPORELLE721def create_time_series_plot(averages_dict: Dict[str, List[float]], dates: List[datetime]) -> bytes:722 """Créer un graphique de l'évolution temporelle des indices"""723 724 plt.figure(figsize=(12, 8), facecolor='white')725 726 # Couleurs pour chaque indice727 colors = {728 'ndvi': '#2E8B57',729 'ndwi': '#4682B4',730 'savi': '#DAA520',731 'evi': '#228B22'732 }733 734 # Tracer chaque indice735 for index_name, values in averages_dict.items():736 if values:737 plt.plot(dates, values, 'o-', label=index_name.upper(), 738 linewidth=2.5, markersize=8, color=colors.get(index_name, 'gray'))739 740 plt.legend(fontsize=12, loc='best')741 plt.title("Évolution temporelle des indices de végétation", fontsize=16, fontweight='bold')742 plt.xlabel("Date", fontsize=14)743 plt.ylabel("Valeur moyenne de l'indice", fontsize=14)744 plt.grid(True, alpha=0.3, linestyle='--')745 746 # Améliorer l'affichage des dates747 plt.gca().xaxis.set_major_locator(ticker.MaxNLocator(10))748 plt.xticks(rotation=45)749 750 # Ajouter des zones colorées pour les saisons751 if dates:752 for i, date in enumerate(dates):753 month = date.month754 if 3 <= month <= 5: # Printemps755 plt.axvspan(date, dates[min(i+1, len(dates)-1)], alpha=0.1, color='green')756 elif 6 <= month <= 8: # Été757 plt.axvspan(date, dates[min(i+1, len(dates)-1)], alpha=0.1, color='yellow')758 elif 9 <= month <= 11: # Automne759 plt.axvspan(date, dates[min(i+1, len(dates)-1)], alpha=0.1, color='orange')760 else: # Hiver761 plt.axvspan(date, dates[min(i+1, len(dates)-1)], alpha=0.1, color='blue')762 763 plt.tight_layout()764 765 # Sauvegarder dans un buffer766 buffer = io.BytesIO()767 plt.savefig(buffer, format='png', dpi=300, bbox_inches='tight')768 buffer.seek(0)769 770 plt.close()771 772 return buffer.getvalue()773 774@app.get("/")775def health_check():776 return {777 "status": "API en fonctionnement",778 "service": "Calculateur d'Indices avec Amélioration des Pixels et Séries Temporelles",779 "version": "3.0.0",780 "features": [781 "amelioration_pixels_avancee", 782 "affichage_ultra_clair", 783 "comparaison_avant_apres",784 "methodes_multiples",785 "series_temporelles",786 "animations_gif",787 "evolution_temporelle"788 ],789 "methodes_amelioration": [790 "adaptive", 791 "super_resolution", 792 "edge_preserving", 793 "bilateral", 794 "segmentation_based"795 ],796 "bibliotheques_disponibles": {797 "opencv": CV2_AVAILABLE,798 "scikit_image": SKIMAGE_AVAILABLE799 }800 }801 802@app.get("/info")803def api_info():804 return {805 "indices_supportes": {806 "NDVI": "Indice de végétation normalisé",807 "NDWI": "Indice d'eau normalisé", 808 "SAVI": "Indice de végétation ajusté au sol",809 "EVI": "Indice de végétation amélioré"810 },811 "ameliorations_pixels_v2": {812 "adaptive": "Lissage adaptatif basé sur la variance locale",813 "super_resolution": "Amélioration par super-résolution x2",814 "edge_preserving": "Filtre préservant les contours (OpenCV requis)",815 "bilateral": "Filtre bilatéral pour réduction du bruit",816 "segmentation_based": "Amélioration par segmentation SLIC (scikit-image requis)"817 },818 "nouvelles_fonctionnalites_v3": {819 "series_temporelles": "Analyse de l'évolution dans le temps",820 "animations_gif": "Création d'animations temporelles",821 "graphiques_evolution": "Visualisation de l'évolution des indices",822 "export_rgb": "Export des images RGB",823 "segmentation_vegetation": "Masquage automatique de la végétation"824 },825 "endpoints": {826 "/calculate-indices": "Calcul simple des indices",827 "/calculate-time-series": "Calcul de série temporelle",828 "/methods": "Liste des méthodes disponibles",829 "/health": "État de santé de l'API"830 }831 }832 833@app.post("/calculate-indices")834async def calculate_indices(835 geojson: UploadFile = File(...),836 enhancement_method: str = Query(default="adaptive", description="Méthode d'amélioration des pixels")837):838 logger.info(f"Début du calcul avec amélioration '{enhancement_method}' pour: {geojson.filename}")839 840 # Valider la méthode d'amélioration841 valid_methods = ['adaptive', 'super_resolution', 'edge_preserving', 'bilateral', 'segmentation_based']842 if enhancement_method not in valid_methods:843 enhancement_method = 'adaptive'844 logger.warning(f"Méthode non valide, utilisation de 'adaptive' par défaut")845 846 if not geojson.filename.endswith('.geojson') and not geojson.filename.endswith('.json'):847 raise HTTPException(848 status_code=400, 849 detail="Le fichier doit être un GeoJSON"850 )851 852 try:853 # Lire et valider le GeoJSON854 gdf = gpd.read_file(geojson.file)855 validate_geojson(gdf)856 857 aoi_polygon = gdf.geometry.dropna().iloc[0]858 geometry = Geometry(aoi_polygon, CRS.WGS84)859 bounds = aoi_polygon.bounds860 861 # Configurer SentinelHub862 config = get_sentinelhub_config()863 864 end_date = datetime.now()865 start_date = end_date - timedelta(days=180)866 867 evalscript = """868 //VERSION=3869 function setup() {870 return {871 input: ["B02", "B03", "B04", "B08"],872 output: {bands: 4, sampleType: "FLOAT32"}873 };874 }875 876 function evaluatePixel(sample) {877 return [sample.B02, sample.B03, sample.B04, sample.B08];878 }879 """880 881 request_sh = SentinelHubRequest(882 evalscript=evalscript,883 input_data=[SentinelHubRequest.input_data(884 data_collection=DataCollection.SENTINEL2_L2A,885 time_interval=(start_date.strftime('%Y-%m-%d'), end_date.strftime('%Y-%m-%d')),886 mosaicking_order='mostRecent'887 )],888 responses=[SentinelHubRequest.output_response('default', MimeType.TIFF)],889 geometry=geometry,890 size=(512, 512),891 config=config892 )893 894 # Obtenir les données satellite895 logger.info("Récupération des données satellite...")896 bands_data = request_sh.get_data()897 898 if not bands_data or len(bands_data) == 0:899 raise HTTPException(900 status_code=404,901 detail="Aucune imagerie satellite trouvée pour cette zone et période"902 )903 904 bands_img = bands_data[0]905 logger.info(f"Données satellite récupérées: {bands_img.shape}")906 907 # Créer le masque géométrique908 logger.info("Création du masque géométrique...")909 transform = from_bounds(bounds[0], bounds[1], bounds[2], bounds[3], 910 bands_img.shape[1], bands_img.shape[0])911 mask = create_mask_from_geometry(aoi_polygon, bands_img.shape[:2], transform)912 913 # Calculer les indices de végétation914 logger.info("Calcul des indices de végétation...")915 indices = calculate_vegetation_indices_improved(bands_img)916 917 # Générer l'image RGB918 rgb_image = s2_to_rgb(bands_img)919 920 # Générer les visualisations améliorées921 logger.info(f"Génération des visualisations avec méthode: {enhancement_method}")922 plots = generate_enhanced_visualization(indices, mask, bounds, enhancement_method)923 924 # Ajouter l'image RGB aux plots925 with tempfile.TemporaryDirectory() as tmpdir:926 rgb_path = os.path.join(tmpdir, "rgb_image.png")927 928 fig, ax = plt.subplots(figsize=(10, 10), facecolor='white')929 ax.imshow(rgb_image, extent=[bounds[0], bounds[2], bounds[1], bounds[3]])930 ax.set_title("Image RGB Sentinel-2", fontsize=14, fontweight='bold')931 ax.set_xlabel('Longitude', fontsize=11)932 ax.set_ylabel('Latitude', fontsize=11)933 ax.grid(True, alpha=0.3)934 935 plt.tight_layout()936 plt.savefig(rgb_path, dpi=300, bbox_inches='tight', facecolor='white')937 plt.close()938 939 with open(rgb_path, "rb") as f:940 plots["rgb_image"] = base64.b64encode(f.read()).decode("utf-8")941 942 # Calculer les statistiques comparatives943 logger.info("Calcul des statistiques...")944 statistics = {}945 improvement_metrics = {}946 947 for idx, array in indices.items():948 # Données originales949 masked_original = array[mask]950 valid_original = masked_original[np.isfinite(masked_original)]951 952 # Données améliorées953 enhanced_array = advanced_pixel_enhancement(array, mask, enhancement_method)954 masked_enhanced = enhanced_array[mask]955 valid_enhanced = masked_enhanced[np.isfinite(masked_enhanced)]956 957 if len(valid_original) > 0:958 # Statistiques originales959 stats_original = {960 "moyenne": float(np.mean(valid_original)),961 "ecart_type": float(np.std(valid_original)),962 "minimum": float(np.min(valid_original)),963 "maximum": float(np.max(valid_original)),964 "mediane": float(np.median(valid_original)),965 "pixels_valides": int(len(valid_original))966 }967 968 # Statistiques améliorées969 stats_enhanced = {970 "moyenne": float(np.mean(valid_enhanced)),971 "ecart_type": float(np.std(valid_enhanced)),972 "minimum": float(np.min(valid_enhanced)),973 "maximum": float(np.max(valid_enhanced)),974 "mediane": float(np.median(valid_enhanced)),975 "pixels_valides": int(len(valid_enhanced))976 }977 978 # Métriques d'amélioration979 variance_reduction = (stats_original["ecart_type"] - stats_enhanced["ecart_type"]) / stats_original["ecart_type"] * 100980 981 # Calcul de la corrélation entre original et amélioré982 correlation = float(np.corrcoef(valid_original, valid_enhanced)[0, 1]) if len(valid_original) == len(valid_enhanced) else 1.0983 984 # Signal-to-Noise Ratio improvement (approximation)985 snr_original = stats_original["moyenne"] / stats_original["ecart_type"] if stats_original["ecart_type"] > 0 else 0986 snr_enhanced = stats_enhanced["moyenne"] / stats_enhanced["ecart_type"] if stats_enhanced["ecart_type"] > 0 else 0987 snr_improvement = snr_enhanced - snr_original988 989 improvement_metrics[idx] = {990 "reduction_variance_pourcent": float(variance_reduction),991 "correlation_original_ameliore": float(correlation),992 "amelioration_snr": float(snr_improvement),993 "conservation_moyenne": abs(stats_original["moyenne"] - stats_enhanced["moyenne"]) < 0.01,994 "qualite_amelioration": "excellente" if variance_reduction > 20 else "bonne" if variance_reduction > 10 else "moderate"995 }996 997 statistics[idx] = {998 "original": stats_original,999 "ameliore": stats_enhanced,1000 "amelioration": improvement_metrics[idx]1001 }1002 1003 # Calculer les statistiques de classification améliorée1004 classification_stats = {}1005 for idx, array in indices.items():1006 enhanced_array = advanced_pixel_enhancement(array, mask, enhancement_method)1007 interpretation = create_smooth_interpretation_map(enhanced_array, idx)1008 classified_map = interpretation['classified_map']1009 1010 # Distribution des classes1011 class_distribution = {}1012 total_valid_pixels = np.sum(np.isfinite(classified_map))1013 1014 if total_valid_pixels > 0:1015 for i, label in enumerate(interpretation['labels']):1016 count = np.sum(classified_map == i)1017 percentage = (count / total_valid_pixels) * 1001018 class_distribution[label] = {1019 "pixels": int(count),1020 "pourcentage": float(percentage)1021 }1022 1023 classification_stats[idx] = {1024 "distribution_classes": class_distribution,1025 "classe_dominante": max(class_distribution.items(), key=lambda x: x[1]["pourcentage"])[0] if class_distribution else "Inconnue",1026 "diversite_classes": len([k for k, v in class_distribution.items() if v["pourcentage"] > 1.0])1027 }1028 1029 # Informations sur la méthode d'amélioration utilisée1030 method_info = {1031 "adaptive": {1032 "description": "Lissage adaptatif basé sur la variance locale",1033 "avantages": ["Préserve les détails", "Lisse les zones homogènes", "Adaptatif automatiquement"],1034 "cas_usage": "Idéal pour tous types de terrain"1035 },1036 "super_resolution": {1037 "description": "Super-résolution par interpolation bicubique",1038 "avantages": ["Améliore la résolution", "Interpolation avancée", "Détails fins"],1039 "cas_usage": "Terrains avec détails fins importants"1040 },1041 "edge_preserving": {1042 "description": "Filtre préservant les contours (OpenCV)",1043 "avantages": ["Préserve les bordures nettes", "Réduit le bruit", "Contours nets"],1044 "cas_usage": "Zones avec bordures nettes importantes",1045 "disponible": CV2_AVAILABLE1046 },1047 "bilateral": {1048 "description": "Filtre bilatéral pour réduction du bruit",1049 "avantages": ["Réduction du bruit", "Préservation des détails", "Lissage intelligent"],1050 "cas_usage": "Images avec beaucoup de bruit"1051 },1052 "segmentation_based": {1053 "description": "Amélioration par segmentation SLIC",1054 "avantages": ["Segmentation intelligente", "Zones homogènes", "Classification précise"],1055 "cas_usage": "Zones avec régions distinctes",1056 "disponible": SKIMAGE_AVAILABLE1057 }1058 }1059 1060 # Recommandations automatiques1061 recommendations = []1062 1063 # Analyser la qualité d'amélioration pour donner des recommandations1064 avg_variance_reduction = np.mean([metrics["reduction_variance_pourcent"] for metrics in improvement_metrics.values()])1065 1066 if avg_variance_reduction < 5:1067 recommendations.append("L'image semble déjà de bonne qualité. Essayez 'super_resolution' pour plus de détails.")1068 elif avg_variance_reduction > 30:1069 recommendations.append(f"Excellente amélioration avec '{enhancement_method}'. Qualité optimale atteinte.")1070 else:1071 if enhancement_method == 'adaptive':1072 recommendations.append("Bon résultat. Essayez 'edge_preserving' pour préserver davantage les contours.")1073 else:1074 recommendations.append("Amélioration satisfaisante. L'algorithme adaptatif peut aussi donner de bons résultats.")1075 1076 # Recommandations basées sur les données1077 if any(stats["original"]["ecart_type"] > 0.3 for stats in statistics.values()):1078 recommendations.append("Données bruitées détectées. 'bilateral' pourrait donner de meilleurs résultats.")1079 1080 # Préparer la réponse finale1081 response_data = {1082 "status": "succès",1083 "message": f"Indices calculés avec amélioration des pixels (méthode: {enhancement_method})",1084 "data": {1085 "images": plots,1086 "statistiques": statistics,1087 "classification": classification_stats,1088 "metadonnees": {1089 "satellite": "Sentinel-2 L2A",1090 "periode_recherche": f"{start_date.strftime('%Y-%m-%d')} à {end_date.strftime('%Y-%m-%d')}",1091 "resolution_image": "512x512 pixels",1092 "methode_amelioration": enhancement_method,1093 "methode_info": method_info.get(enhancement_method, {}),1094 "types_images_generees": ["enhanced_continuous", "enhanced_classified", "comparison", "rgb_image"],1095 "qualite_export": "300 DPI",1096 "bibliotheques_utilisees": {1097 "opencv": CV2_AVAILABLE,1098 "scikit_image": SKIMAGE_AVAILABLE,1099 "scipy": True,1100 "matplotlib": True1101 }1102 },1103 "recommandations": recommendations,1104 "methodes_disponibles": {1105 method: {1106 "description": info["description"],1107 "disponible": info.get("disponible", True)1108 }1109 for method, info in method_info.items()1110 }1111 }1112 }1113 1114 logger.info(f"Traitement terminé avec succès. Méthode: {enhancement_method}, Images générées: {len(plots)}")1115 return JSONResponse(response_data)1116 1117 except HTTPException:1118 raise1119 except Exception as e:1120 logger.error(f"Erreur lors du traitement: {str(e)}", exc_info=True)1121 raise HTTPException(1122 status_code=500, 1123 detail=f"Erreur lors du traitement: {str(e)}"1124 )1125 1126# NOUVEAU ENDPOINT POUR LES SÉRIES TEMPORELLES1127@app.post("/calculate-time-series")1128async def calculate_time_series(1129 geojson: UploadFile = File(...),1130 start_date: str = Query(..., description="Date de début (YYYY-MM-DD)"),1131 end_date: str = Query(..., description="Date de fin (YYYY-MM-DD)"),1132 interval_days: int = Query(default=10, description="Intervalle en jours entre les images"),1133 enhancement_method: str = Query(default="adaptive", description="Méthode d'amélioration des pixels"),1134 create_animation: bool = Query(default=True, description="Créer des animations GIF"),1135 indices: List[str] = Query(default=["ndvi"], description="Liste des indices à calculer")1136):1137 logger.info(f"Début du calcul de série temporelle pour: {geojson.filename}")1138 1139 try:1140 # Parser les dates1141 start_dt = datetime.strptime(start_date, "%Y-%m-%d")1142 end_dt = datetime.strptime(end_date, "%Y-%m-%d")1143 1144 if start_dt >= end_dt:1145 raise HTTPException(status_code=400, detail="La date de début doit être avant la date de fin")1146 1147 # Lire et valider le GeoJSON1148 gdf = gpd.read_file(geojson.file)1149 validate_geojson(gdf)1150 1151 aoi_polygon = gdf.geometry.dropna().iloc[0]1152 geometry = Geometry(aoi_polygon, CRS.WGS84)1153 bounds = aoi_polygon.bounds1154 1155 # Configurer SentinelHub1156 config = get_sentinelhub_config()1157 1158 # Obtenir la série temporelle d'images1159 logger.info(f"Récupération des images de {start_date} à {end_date} avec intervalle de {interval_days} jours")1160 images, dates = get_time_series_images(geometry, start_dt, end_dt, config, interval_days)1161 1162 if not images:1163 raise HTTPException(1164 status_code=404,1165 detail="Aucune image trouvée pour cette période et cette zone"1166 )1167 1168 logger.info(f"{len(images)} images récupérées")1169 1170 # Créer le masque géométrique1171 transform = from_bounds(bounds[0], bounds[1], bounds[2], bounds[3], 1172 images[0].shape[1], images[0].shape[0])1173 mask = create_mask_from_geometry(aoi_polygon, images[0].shape[:2], transform)1174 1175 # Calculer les indices pour chaque image et collecter les moyennes1176 averages_dict = {idx: [] for idx in indices}1177 1178 for img in images:1179 img_indices = calculate_vegetation_indices_improved(img)1180 for idx in indices:1181 if idx in img_indices:1182 # Appliquer l'amélioration1183 enhanced = advanced_pixel_enhancement(img_indices[idx], mask, enhancement_method)1184 avg = calculate_masked_average(enhanced, mask)1185 averages_dict[idx].append(avg)1186 1187 # Créer les animations si demandé1188 animations = {}1189 if create_animation:1190 logger.info("Création des animations...")1191 for idx in indices:1192 try:1193 animation_data = create_timelapse_animation(images, dates, idx, mask, bounds)1194 if animation_data:1195 animations[f"{idx}_animation"] = base64.b64encode(animation_data).decode("utf-8")1196 except Exception as e:1197 logger.error(f"Erreur lors de la création de l'animation pour {idx}: {e}")1198 1199 # Créer le graphique de série temporelle1200 logger.info("Création du graphique de série temporelle...")