update chức năng ndvi time seriese

api_server.py

@@ -122,6 +122,31 @@ class TrainingStatus(BaseModel):
     end_time: Optional[str]
 
 
+class NDVIConfig(BaseModel):
+    """Cấu hình tính NDVI time series"""
+    bbox: List[float]  # [min_lon, min_lat, max_lon, max_lat]
+    start_date: str
+    end_date: str
+    max_cloud_cover: int = 30
+    resolution: int = 20
+
+
+class PredictionWithNDVIConfig(BaseModel):
+    """Cấu hình predict kết hợp land classification và NDVI"""
+    model_filename: str
+    min_lon: float
+    min_lat: float
+    max_lon: float
+    max_lat: float
+    start_date: str
+    end_date: str
+    max_scenes: int = 12
+    cloud_cover: int = 30
+    resolution: int = 20
+    export_ndvi: bool = True  # Export NDVI raster
+    export_classification: bool = True  # Export classification raster
+
+
 @app.get("/", response_class=HTMLResponse)
 async def root():
     """Serve main index page with tabs"""
@@ -173,6 +198,26 @@ async def dashboard():
         raise HTTPException(status_code=404, detail="Dashboard không tồn tại")
 
 
+@app.get("/batch", response_class=HTMLResponse)
+async def batch_page():
+    """Serve batch processing interface"""
+    html_file = Path(__file__).parent / "batch_interface.html"
+    if html_file.exists():
+        return FileResponse(html_file)
+    else:
+        raise HTTPException(status_code=404, detail="Batch interface không tồn tại")
+
+
+@app.get("/ndvi", response_class=HTMLResponse)
+async def ndvi_page():
+    """Serve NDVI time series interface"""
+    html_file = Path(__file__).parent / "ndvi_interface.html"
+    if html_file.exists():
+        return FileResponse(html_file)
+    else:
+        raise HTTPException(status_code=404, detail="NDVI interface không tồn tại")
+
+
 @app.get("/api/config/presets")
 async def get_presets():
     """Lấy các preset cấu hình sẵn"""
@@ -398,14 +443,39 @@ async def list_reports():
         else:
             report_type = "unknown"
         
-        reports.append({
+        report_info = {
             "filename": report_file.name,
             "type": report_type,
             "created": datetime.fromtimestamp(report_file.stat().st_mtime).isoformat(),
             "size_kb": round(report_file.stat().st_size / 1024, 2),
             "view_url": f"/api/reports/view/{report_file.name}",
-            "download_url": f"/api/reports/download/{report_file.name}"
-        })
+            "download_url": f"/api/reports/download/{report_file.name}",
+            "is_batch_job": False,
+            "batch_metadata": None
+        }
+        
+        # Check if this is a batch job report
+        if report_type == "prediction":
+            predictions_dir = Path("predictions")
+            # Look for batch metadata JSON files that reference this report
+            for json_file in predictions_dir.glob("batch_*.json"):
+                try:
+                    import json
+                    with open(json_file, 'r') as f:
+                        metadata = json.load(f)
+                        if metadata.get("report_filename") == report_file.name or \
+                           (metadata.get("batch_job_id") and report_file.name.endswith('.html')):
+                            report_info["is_batch_job"] = True
+                            report_info["batch_metadata"] = {
+                                "batch_job_id": metadata.get("batch_job_id"),
+                                "batch_name": metadata.get("batch_name"),
+                                "batch_timestamp": metadata.get("batch_timestamp")
+                            }
+                            break
+                except Exception as e:
+                    pass
+        
+        reports.append(report_info)
     
     # Sort by creation time (newest first)
     reports.sort(key=lambda x: x["created"], reverse=True)
@@ -1025,12 +1095,37 @@ async def list_predictions():
     
     predictions = []
     for pred_file in predictions_dir.glob("*.tif"):
-        predictions.append({
+        pred_info = {
             "filename": pred_file.name,
             "created": datetime.fromtimestamp(pred_file.stat().st_mtime).isoformat(),
             "size_mb": round(pred_file.stat().st_size / 1024 / 1024, 2),
-            "download_url": f"/api/predictions/download/{pred_file.name}"
-        })
+            "download_url": f"/api/predictions/download/{pred_file.name}",
+            "is_batch_job": pred_file.name.startswith("batch_"),
+            "batch_metadata": None
+        }
+        
+        # Try to load batch metadata from JSON sidecar if exists
+        json_file = pred_file.with_suffix('.json')
+        if json_file.exists():
+            try:
+                import json
+                with open(json_file, 'r') as f:
+                    metadata = json.load(f)
+                    pred_info["batch_metadata"] = {
+                        "batch_job_id": metadata.get("batch_job_id"),
+                        "batch_name": metadata.get("batch_name"),
+                        "batch_timestamp": metadata.get("batch_timestamp")
+                    }
+            except Exception as e:
+                print(f"[METADATA ERROR] Failed to load {json_file}: {e}")
+        
+        # Check PNG preview
+        png_file = pred_file.with_suffix('.png')
+        pred_info["has_preview"] = png_file.exists()
+        if png_file.exists():
+            pred_info["preview_url"] = f"/api/predictions/preview/{png_file.name}"
+        
+        predictions.append(pred_info)
     
     # Sort by creation time (newest first)
     predictions.sort(key=lambda x: x["created"], reverse=True)
@@ -1366,6 +1461,8 @@ async def process_batch_queue():
     """Process batch prediction queue"""
     global batch_queue, batch_results
     
+    import asyncio
+    
     while batch_queue:
         # Get next job
         job = None
@@ -1379,26 +1476,27 @@ async def process_batch_queue():
         
         # Mark as running
         job["status"] = "running"
+        job["progress"] = 0
         job["started_at"] = datetime.now().isoformat()
         
         try:
             # Create PredictionConfig from job config
             pred_config = PredictionConfig(**job["config"])
             
-            # Run prediction (simplified version)
-            # In real implementation, call the actual prediction function
-            print(f"[BATCH] Processing job: {job['name']}")
+            print(f"[BATCH] Processing job {job['job_id']}: {job['name']}")
+            job["progress"] = 5
             
-            # Simulate prediction (replace with actual prediction call)
-            # await run_prediction(pred_config)
+            # Run prediction synchronously (in the same thread to avoid conflicts)
+            await asyncio.to_thread(run_batch_prediction, job, pred_config)
             
-            # For now, mark as completed
-            job["status"] = "completed"
-            job["completed_at"] = datetime.now().isoformat()
-            job["result"] = {
-                "output_file": f"predictions/batch_{job['job_id']}.tif",
-                "message": "Prediction completed successfully"
-            }
+            # Check if prediction was successful
+            if job.get("result") and not job.get("error"):
+                job["status"] = "completed"
+                job["progress"] = 100
+                job["completed_at"] = datetime.now().isoformat()
+                print(f"[BATCH] Job {job['job_id']} completed successfully")
+            else:
+                raise Exception(job.get("error", "Unknown error during prediction"))
             
         except Exception as e:
             job["error"] = str(e)
@@ -1407,12 +1505,14 @@ async def process_batch_queue():
             if job["retries"] < job["max_retries"]:
                 job["retries"] += 1
                 job["status"] = "queued"  # Retry
-                print(f"[BATCH] Job {job['name']} failed, retrying ({job['retries']}/{job['max_retries']})")
+                job["progress"] = 0
+                print(f"[BATCH] Job {job['job_id']} ({job['name']}) failed, retrying ({job['retries']}/{job['max_retries']}): {e}")
                 continue
             else:
                 job["status"] = "failed"
+                job["progress"] = 0
                 job["completed_at"] = datetime.now().isoformat()
-                print(f"[BATCH] Job {job['name']} failed permanently: {e}")
+                print(f"[BATCH] Job {job['job_id']} ({job['name']}) failed permanently: {e}")
         
         # Move to results
         batch_queue.remove(job)
@@ -1423,6 +1523,540 @@ async def process_batch_queue():
             batch_results = batch_results[-100:]
 
 
+def run_batch_prediction(job: dict, config: PredictionConfig):
+    """Run prediction for a single batch job"""
+    try:
+        job["progress"] = 10
+        
+        # Import required libraries
+        import xarray as xr
+        import numpy as np
+        from datetime import datetime as dt
+        import rioxarray
+        import dask.array as da
+        
+        job["progress"] = 15
+        
+        # Load model
+        model_path = Path("model_train") / config.model_filename
+        if not model_path.exists():
+            raise FileNotFoundError(f"Model không tồn tại: {config.model_filename}")
+        
+        model_data = joblib.load(model_path)
+        
+        if isinstance(model_data, dict):
+            model = model_data.get('model')
+            label_encoder = model_data.get('label_encoder')
+        else:
+            model = model_data
+            label_encoder = None
+        
+        job["progress"] = 20
+        
+        # Check if CNN model
+        is_cnn_model = hasattr(model, '__class__') and 'CNN' in model.__class__.__name__
+        
+        # Load data from Microsoft Planetary Computer
+        import pystac_client
+        import planetary_computer
+        from odc.stac import load
+        
+        catalog = pystac_client.Client.open(
+            "https://planetarycomputer.microsoft.com/api/stac/v1",
+            modifier=planetary_computer.sign_inplace,
+        )
+        
+        bbox = [config.min_lon, config.min_lat, config.max_lon, config.max_lat]
+        time_range = f"{config.start_date}/{config.end_date}"
+        
+        job["progress"] = 25
+        
+        # Search Sentinel-2
+        s2_search = catalog.search(
+            collections=["sentinel-2-l2a"],
+            bbox=bbox,
+            datetime=time_range,
+            query={"eo:cloud_cover": {"lt": config.cloud_cover}}
+        )
+        
+        s2_items = list(s2_search.items())
+        if not s2_items:
+            raise ValueError("Không tìm thấy dữ liệu Sentinel-2")
+        
+        s2_items = s2_items[:config.max_scenes]
+        
+        job["progress"] = 35
+        
+        # Load Sentinel-2 data
+        s2_data = load(
+            s2_items,
+            bbox=bbox,
+            chunks={"time": 1, "x": 2048, "y": 2048},
+            groupby="solar_day",
+            resolution=config.resolution
+        )
+        
+        job["progress"] = 50
+        
+        # Calculate NDVI
+        nir = s2_data["B08"].astype('float32')
+        red = s2_data["B04"].astype('float32')
+        ndvi = (nir - red) / (nir + red + 1e-8)
+        
+        # Mask clouds if SCL available
+        if "SCL" in s2_data:
+            scl = s2_data["SCL"]
+            cloud_mask = (scl == 3) | (scl == 8) | (scl == 9) | (scl == 10)
+            ndvi = ndvi.where(~cloud_mask)
+        
+        # Fill NaN and resample
+        ndvi_filled = ndvi.ffill(dim='time').bfill(dim='time')
+        ndvi_monthly = ndvi_filled.resample(time="1ME").mean().compute()
+        
+        job["progress"] = 70
+        
+        # Prepare features
+        n_times_ndvi = len(ndvi_monthly.time)
+        y_size = len(ndvi_monthly.y)
+        x_size = len(ndvi_monthly.x)
+        n_pixels = y_size * x_size
+        
+        ndvi_features = []
+        for t in range(n_times_ndvi):
+            ndvi_t = ndvi_monthly.isel(time=t).values.flatten()
+            ndvi_features.append(ndvi_t)
+        
+        features = np.column_stack(ndvi_features)
+        features = np.nan_to_num(features, nan=0.0)
+        
+        job["progress"] = 80
+        
+        # Adjust features to match model expectations
+        try:
+            if is_cnn_model:
+                expected_features = model.n_features
+            elif hasattr(model, 'n_features_in_'):
+                expected_features = model.n_features_in_
+            else:
+                try:
+                    expected_features = model.get_booster().num_features()
+                except:
+                    expected_features = features.shape[1]
+            
+            if features.shape[1] > expected_features:
+                features = features[:, :expected_features]
+            elif features.shape[1] < expected_features:
+                n_missing = expected_features - features.shape[1]
+                padding = np.tile(features[:, -1:], (1, n_missing))
+                features = np.column_stack([features, padding])
+        except:
+            pass
+        
+        # Predict
+        if is_cnn_model:
+            predictions = model.predict(features)
+        else:
+            predictions = model.predict(features)
+        
+        # Decode labels
+        if label_encoder is not None:
+            try:
+                predictions = label_encoder.inverse_transform(predictions)
+            except:
+                pass
+        
+        job["progress"] = 90
+        
+        # Reshape and create output
+        pred_shape = (y_size, x_size)
+        predictions_2d = predictions.reshape(pred_shape)
+        
+        prediction_da = xr.DataArray(
+            predictions_2d,
+            coords={"y": ndvi_monthly.y, "x": ndvi_monthly.x},
+            dims=["y", "x"],
+            name="classification"
+        )
+        
+        # Save output
+        output_dir = Path("predictions")
+        output_dir.mkdir(exist_ok=True)
+        
+        output_file = output_dir / f"batch_{job['job_id']}_{job['name'].replace(' ', '_')}.tif"
+        
+        if hasattr(s2_data, 'rio') and s2_data.rio.crs is not None:
+            prediction_da.rio.write_crs(s2_data.rio.crs, inplace=True)
+        else:
+            prediction_da.rio.write_crs("EPSG:4326", inplace=True)
+        
+        prediction_da.rio.to_raster(str(output_file), driver="GTiff")
+        
+        # Generate PNG preview
+        png_file = output_dir / f"batch_{job['job_id']}_{job['name'].replace(' ', '_')}.png"
+        try:
+            import matplotlib
+            matplotlib.use('Agg')
+            import matplotlib.pyplot as plt
+            
+            fig, ax = plt.subplots(figsize=(12, 10), dpi=150)
+            im = ax.imshow(predictions_2d, cmap='tab20', interpolation='nearest')
+            ax.set_title(f"{job['name']} - Batch {job['job_id']}", fontsize=14, fontweight='bold')
+            ax.set_xlabel('X (pixels)', fontsize=10)
+            ax.set_ylabel('Y (pixels)', fontsize=10)
+            
+            cbar = plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
+            cbar.set_label('Class', rotation=270, labelpad=15)
+            ax.grid(True, alpha=0.3, linestyle='--', linewidth=0.5)
+            
+            plt.tight_layout()
+            plt.savefig(str(png_file), dpi=150, bbox_inches='tight')
+            plt.close(fig)
+        except Exception as e:
+            print(f"[BATCH PNG ERROR] {e}")
+            png_file = None
+        
+        # Get unique classes
+        unique_classes = np.unique(predictions_2d)
+        unique_classes = unique_classes[~np.isnan(unique_classes)].tolist()
+        
+        # Store result in job with batch metadata
+        job["result"] = {
+            "output_file": str(output_file),
+            "png_file": str(png_file) if png_file else None,
+            "shape": list(pred_shape),
+            "unique_classes": unique_classes,
+            "bbox": bbox,
+            "time_range": time_range,
+            "n_features": features.shape[1],
+            "n_times_ndvi": n_times_ndvi,
+            "model_used": config.model_filename,
+            "batch_job_id": job["job_id"],
+            "batch_name": job["name"],
+            "batch_timestamp": datetime.now().isoformat()
+        }
+        
+        # Save batch metadata to JSON sidecar file for persistence
+        metadata_file = output_file.with_suffix('.json')
+        try:
+            import json
+            with open(metadata_file, 'w') as f:
+                json.dump(job["result"], f, indent=2, default=str)
+            print(f"[BATCH METADATA] Saved to {metadata_file}")
+        except Exception as e:
+            print(f"[BATCH METADATA ERROR] Failed to save metadata: {e}")
+        
+        # Auto generate prediction report for batch job
+        try:
+            from report_generator import generate_prediction_report
+            report_path, _ = generate_prediction_report(job["result"])
+            job["result"]["report_path"] = report_path
+            job["result"]["report_filename"] = Path(report_path).name
+            print(f"[BATCH REPORT] Generated prediction report: {report_path}")
+        except Exception as e:
+            print(f"[BATCH REPORT ERROR] Failed to generate report: {e}")
+        
+        job["progress"] = 100
+        
+    except Exception as e:
+        job["error"] = str(e)
+        import traceback
+        print(f"[BATCH ERROR] Job {job['job_id']}: {traceback.format_exc()}")
+
+
+# ============ PREDICTION WITH NDVI API ============
+
+@app.post("/api/predict/with-ndvi")
+async def predict_with_ndvi(config: PredictionWithNDVIConfig, background_tasks: BackgroundTasks):
+    """Predict land classification và NDVI cho một khu vực"""
+    try:
+        import numpy as np
+        import xarray as xr
+        from pystac_client import Client
+        import planetary_computer
+        import odc.stac
+        import rasterio
+        from rasterio.transform import from_bounds
+        
+        # Load model
+        model_path = Path(f"model_train/{config.model_filename}")
+        if not model_path.exists():
+            raise HTTPException(status_code=404, detail=f"Model {config.model_filename} không tồn tại")
+        
+        model = joblib.load(model_path)
+        print(f"[PREDICT+NDVI] Loaded model: {config.model_filename}")
+        
+        # Connect to Microsoft Planetary Computer
+        catalog = Client.open(
+            "https://planetarycomputer.microsoft.com/api/stac/v1",
+            modifier=planetary_computer.sign_inplace
+        )
+        
+        bbox = [config.min_lon, config.min_lat, config.max_lon, config.max_lat]
+        time_range = f"{config.start_date}/{config.end_date}"
+        
+        # Search for Sentinel-2 data
+        search = catalog.search(
+            collections=["sentinel-2-l2a"],
+            bbox=bbox,
+            datetime=time_range,
+            query={"eo:cloud_cover": {"lt": config.cloud_cover}}
+        )
+        
+        items = list(search.items())[:config.max_scenes]
+        print(f"[PREDICT+NDVI] Found {len(items)} Sentinel-2 scenes")
+        
+        if len(items) == 0:
+            raise HTTPException(status_code=404, detail="Không tìm thấy dữ liệu vệ tinh")
+        
+        # Load all bands needed for features
+        data = odc.stac.load(
+            items,
+            bbox=bbox,
+            bands=["B02", "B03", "B04", "B08"],  # Blue, Green, Red, NIR
+            resolution=config.resolution,
+            chunks={"x": 2048, "y": 2048}
+        ).compute()
+        
+        print(f"[PREDICT+NDVI] Loaded data shape: {data.dims}")
+        
+        # Calculate NDVI and other indices
+        blue = data["B02"].values
+        green = data["B03"].values
+        red = data["B04"].values
+        nir = data["B08"].values
+        
+        # Calculate indices
+        # NDVI = (NIR - Red) / (NIR + Red)
+        ndvi = (nir - red) / (nir + red + 1e-8)
+        
+        # NDWI = (Green - NIR) / (Green + NIR)
+        ndwi = (green - nir) / (green + nir + 1e-8)
+        
+        # NDBI = (SWIR - NIR) / (SWIR + NIR) - we use Red as proxy
+        ndbi = (red - nir) / (red + nir + 1e-8)
+        
+        # Prepare features for prediction
+        # Assuming model was trained with [NDVI, NDWI, NDBI] features
+        height, width = ndvi.shape[1:3]  # Skip time dimension
+        n_pixels = height * width
+        
+        # Average over time dimension
+        ndvi_mean = np.nanmean(ndvi, axis=0)
+        ndwi_mean = np.nanmean(ndwi, axis=0)
+        ndbi_mean = np.nanmean(ndbi, axis=0)
+        
+        # Reshape for prediction
+        features = np.stack([ndvi_mean.flatten(), ndwi_mean.flatten(), ndbi_mean.flatten()], axis=1)
+        
+        # Handle NaN values
+        valid_mask = ~np.isnan(features).any(axis=1)
+        features_clean = features[valid_mask]
+        
+        print(f"[PREDICT+NDVI] Predicting {features_clean.shape[0]} valid pixels...")
+        
+        # Predict
+        predictions = model.predict(features_clean)
+        
+        # Reshape back to raster
+        prediction_raster = np.full(n_pixels, -1, dtype=np.int16)
+        prediction_raster[valid_mask] = predictions
+        prediction_raster = prediction_raster.reshape(height, width)
+        
+        # Prepare outputs
+        timestamp = datetime.now().strftime('%Y%m%d_%H%M%S')
+        output_dir = Path("predictions")
+        output_dir.mkdir(exist_ok=True)
+        
+        output_files = []
+        
+        # Export NDVI if requested
+        if config.export_ndvi:
+            ndvi_file = output_dir / f"ndvi_{timestamp}.tif"
+            transform = from_bounds(bbox[0], bbox[1], bbox[2], bbox[3], width, height)
+            
+            with rasterio.open(
+                ndvi_file, 'w',
+                driver='GTiff',
+                height=height,
+                width=width,
+                count=1,
+                dtype=ndvi_mean.dtype,
+                crs='EPSG:4326',
+                transform=transform
+            ) as dst:
+                dst.write(ndvi_mean, 1)
+            
+            output_files.append({"type": "ndvi", "path": str(ndvi_file)})
+            print(f"[PREDICT+NDVI] Saved NDVI to {ndvi_file}")
+        
+        # Export classification if requested
+        if config.export_classification:
+            class_file = output_dir / f"classification_{timestamp}.tif"
+            transform = from_bounds(bbox[0], bbox[1], bbox[2], bbox[3], width, height)
+            
+            with rasterio.open(
+                class_file, 'w',
+                driver='GTiff',
+                height=height,
+                width=width,
+                count=1,
+                dtype=prediction_raster.dtype,
+                crs='EPSG:4326',
+                transform=transform
+            ) as dst:
+                dst.write(prediction_raster, 1)
+            
+            output_files.append({"type": "classification", "path": str(class_file)})
+            print(f"[PREDICT+NDVI] Saved classification to {class_file}")
+        
+        # Calculate statistics
+        ndvi_stats = {
+            "mean": float(np.nanmean(ndvi_mean)),
+            "min": float(np.nanmin(ndvi_mean)),
+            "max": float(np.nanmax(ndvi_mean)),
+            "std": float(np.nanstd(ndvi_mean))
+        }
+        
+        # Count classes
+        unique_classes, counts = np.unique(predictions, return_counts=True)
+        class_distribution = {
+            int(cls): int(count) for cls, count in zip(unique_classes, counts)
+        }
+        
+        return {
+            "success": True,
+            "message": "Prediction with NDVI completed",
+            "output_files": output_files,
+            "ndvi_stats": ndvi_stats,
+            "class_distribution": class_distribution,
+            "n_scenes": len(items),
+            "resolution": config.resolution,
+            "bbox": bbox
+        }
+        
+    except Exception as e:
+        print(f"[PREDICT+NDVI ERROR] {str(e)}")
+        import traceback
+        traceback.print_exc()
+        raise HTTPException(status_code=500, detail=str(e))
+
+
+# ============ NDVI TIME SERIES API ============
+
+@app.post("/api/ndvi/timeseries")
+async def calculate_ndvi_timeseries(config: NDVIConfig):
+    """Tính NDVI time series cho một khu vực"""
+    try:
+        import numpy as np
+        import xarray as xr
+        from pystac_client import Client
+        import planetary_computer
+        import odc.stac
+        
+        print(f"[NDVI] Starting calculation for bbox: {config.bbox}, time: {config.start_date} to {config.end_date}")
+        
+        # Connect to Microsoft Planetary Computer STAC API
+        catalog = Client.open(
+            "https://planetarycomputer.microsoft.com/api/stac/v1",
+            modifier=planetary_computer.sign_inplace
+        )
+        
+        bbox = config.bbox
+        time_range = f"{config.start_date}/{config.end_date}"
+        
+        # Search for Sentinel-2 data
+        search = catalog.search(
+            collections=["sentinel-2-l2a"],
+            bbox=bbox,
+            datetime=time_range,
+            query={"eo:cloud_cover": {"lt": config.max_cloud_cover}}
+        )
+        
+        items = list(search.items())
+        print(f"[NDVI] Found {len(items)} Sentinel-2 scenes")
+        
+        if len(items) == 0:
+            raise HTTPException(status_code=404, detail="Không tìm thấy dữ liệu Sentinel-2 cho khu vực và thời gian này")
+        
+        # Load data for each time step
+        ndvi_timeseries = []
+        
+        for item in items:
+            try:
+                # Load NIR (B08) and Red (B04) bands
+                data = odc.stac.load(
+                    [item],
+                    bbox=bbox,
+                    bands=["B04", "B08"],  # Red and NIR
+                    resolution=config.resolution,
+                    chunks={"x": 2048, "y": 2048}
+                ).compute()
+                
+                if data is None or len(data.keys()) == 0:
+                    continue
+                
+                # Calculate NDVI = (NIR - Red) / (NIR + Red)
+                nir = data["B08"].values
+                red = data["B04"].values
+                
+                # Avoid division by zero
+                denominator = nir + red
+                denominator = np.where(denominator == 0, np.nan, denominator)
+                
+                ndvi = (nir - red) / denominator
+                
+                # Calculate mean NDVI (ignore NaN values)
+                mean_ndvi = float(np.nanmean(ndvi))
+                
+                # Get date from item
+                date_str = item.datetime.strftime("%Y-%m-%d")
+                
+                ndvi_timeseries.append({
+                    "date": date_str,
+                    "ndvi": mean_ndvi
+                })
+                
+                print(f"[NDVI] {date_str}: NDVI = {mean_ndvi:.3f}")
+                
+            except Exception as e:
+                print(f"[NDVI WARNING] Failed to process item {item.id}: {e}")
+                continue
+        
+        if len(ndvi_timeseries) == 0:
+            raise HTTPException(status_code=500, detail="Không thể tính NDVI cho bất kỳ ảnh nào")
+        
+        # Sort by date
+        ndvi_timeseries.sort(key=lambda x: x["date"])
+        
+        # Calculate statistics
+        ndvi_values = [item["ndvi"] for item in ndvi_timeseries]
+        mean_ndvi = float(np.mean(ndvi_values))
+        min_ndvi = float(np.min(ndvi_values))
+        max_ndvi = float(np.max(ndvi_values))
+        
+        result = {
+            "timeseries": ndvi_timeseries,
+            "n_images": len(ndvi_timeseries),
+            "mean_ndvi": mean_ndvi,
+            "min_ndvi": min_ndvi,
+            "max_ndvi": max_ndvi,
+            "bbox": bbox,
+            "time_range": time_range
+        }
+        
+        print(f"[NDVI] Calculation complete. Mean NDVI: {mean_ndvi:.3f}, Images: {len(ndvi_timeseries)}")
+        
+        return result
+        
+    except HTTPException:
+        raise
+    except Exception as e:
+        print(f"[NDVI ERROR] {e}")
+        import traceback
+        traceback.print_exc()
+        raise HTTPException(status_code=500, detail=f"Lỗi khi tính NDVI: {str(e)}")
+
+
 if __name__ == "__main__":
     print("=" * 70)
     print("🚀 LAND CLASSIFICATION TRAINING API SERVER")