This commit is contained in:
Victor Phan
2026-02-26 15:54:55 +07:00
parent 0db4148a40
commit 5d6efcf56a
2 changed files with 124 additions and 302 deletions
+121 -302
View File
@@ -497,333 +497,152 @@
],
"source": [
"%%time\n",
"# 💾 LOAD SENTINEL-2 DATA DIRECTLY FROM S3 COGS USING RASTERIO - WITH TEMPORAL FEATURES\n",
"# 💾 LOAD SENTINEL-2 DATA VIA dc.load() - ODC NATIVE APPROACH\n",
"# Dùng dc.load() thay vì rasterio.open(S3) → tránh lỗi 403/VPC\n",
"print(\"=\"*70)\n",
"print(\"LOADING SENTINEL-2 DATA FROM S3 COGs (RASTERIO) - OPTIMAL ACCURACY\")\n",
"print(\"LOADING SENTINEL-2 DATA VIA dc.load() (ODC NATIVE)\")\n",
"print(\"=\"*70)\n",
"\n",
"try:\n",
" import rasterio\n",
" import xarray as xr\n",
" import numpy as np\n",
" from scipy import ndimage\n",
"\n",
" # Cấu hình S3 access cho PROCESS CHÍNH (không phải workers)\n",
" configure_s3_access(aws_unsigned=False, requester_pays=True)\n",
" print(\"✅ S3 access configured for local process\")\n",
"\n",
" # Kiểm tra AWS credentials\n",
" aws_key = os.environ.get('AWS_ACCESS_KEY_ID', '')\n",
" aws_secret = os.environ.get('AWS_SECRET_ACCESS_KEY', '')\n",
" if aws_key and aws_secret:\n",
" print(f\"✅ AWS credentials found (key: {aws_key[:8]}...)\")\n",
" else:\n",
" print(\"⚠️ WARNING: AWS credentials not found in environment!\")\n",
" print(\" Set AWS_ACCESS_KEY_ID and AWS_SECRET_ACCESS_KEY to fix HTTP 403 errors.\")\n",
" print(\" Trying unsigned access as fallback...\")\n",
" configure_s3_access(aws_unsigned=True)\n",
" print(\" Switched to aws_unsigned=True (public buckets only)\")\n",
"\n",
" # ===== CHECK IF SHOULD SKIP DOWNLOAD =====\n",
" if use_cache and data is not None:\n",
" if use_cache and 'data' in dir() and data is not None:\n",
" print(f\"\\n✅ Using cached dataset - skipping download!\")\n",
" print(f\" Variables: {len(data.data_vars)}\")\n",
" print(f\" Shape: {data.dims}\")\n",
" display(data)\n",
" \n",
"\n",
" else:\n",
" # ===== DOWNLOAD FROM S3 =====\n",
" print(f\"\\n📥 Downloading from S3...\")\n",
" \n",
" # Get all scenes from datacube metadata (filtered to Vietnam region)\n",
" datasets = list(dc.find_datasets(\n",
" product='s2_l2a',\n",
" time=date_range,\n",
" lat=latitude_range,\n",
" lon=longtitude_range\n",
" ))\n",
" \n",
" if not datasets:\n",
" raise ValueError(\"No datasets found for date range and region\")\n",
" \n",
" print(f\"\\n📦 Found {len(datasets)} available scenes\")\n",
" print(f\" Date range: {date_range[0]} to {date_range[1]}\")\n",
" print(f\" Region: lon={longtitude_range}, lat={latitude_range}\")\n",
" \n",
" # ===== LOAD ALL SCENES WITH ALL AVAILABLE BANDS (NO MAGNIFICATION) =====\n",
" # num_scenes = len(datasets) # Load ALL scenes\n",
" num_scenes = 1 # Load ALL scenes\n",
" print(f\"\\n[LOADING] Loading {num_scenes} scenes with ALL available bands...\")\n",
" print(f\" (Keeping NATIVE resolution - NO upsampling/magnification)\")\n",
" \n",
" all_data_dict = {}\n",
" failed_scenes = []\n",
" scene_dates = []\n",
" \n",
" # Discover all available bands from first scene\n",
" first_scene = datasets[0]\n",
" all_available_bands = list(first_scene.measurements.keys())\n",
" print(f\" Available bands: {all_available_bands}\")\n",
" \n",
" for scene_idx in range(num_scenes):\n",
" selected = datasets[scene_idx]\n",
" scene_label = selected.metadata.label\n",
" scene_datetime = selected.time.begin if hasattr(selected.time, 'begin') else selected.time\n",
" scene_dates.append(scene_datetime)\n",
" \n",
" # Print progress every 5 scenes\n",
" if scene_idx % 5 == 0 or scene_idx == 0 or scene_idx == num_scenes - 1:\n",
" print(f\"\\n [{scene_idx + 1:2d}/{num_scenes}] {scene_label} ({scene_datetime.date()})\")\n",
" \n",
" # Load ALL available bands from S3 COGs\n",
" scene_data_dict = {}\n",
" \n",
" for band_name in all_available_bands:\n",
" if band_name in selected.measurements:\n",
" band_path = selected.measurements[band_name]['path']\n",
" \n",
" try:\n",
" with rasterio.open(band_path) as src:\n",
" data_band = src.read(1)\n",
" scene_data_dict[band_name] = data_band\n",
" except Exception as e:\n",
" if scene_idx % 5 == 0:\n",
" print(f\" ⚠️ Error loading {band_name}: {str(e)[:50]}\")\n",
" failed_scenes.append((scene_idx, scene_label, band_name, str(e)))\n",
" \n",
" if scene_data_dict:\n",
" all_data_dict[scene_idx] = scene_data_dict\n",
" if scene_idx % 5 == 0 or scene_idx == num_scenes - 1:\n",
" print(f\" ✅ {len(scene_data_dict)} bands loaded\")\n",
" else:\n",
" failed_scenes.append((scene_idx, scene_label, \"all\", \"No bands loaded\"))\n",
" \n",
" if not all_data_dict:\n",
" raise ValueError(\"Could not load any bands from any scene\")\n",
" \n",
" print(f\"\\n✅ Successfully loaded {len(all_data_dict)} scenes!\")\n",
" if failed_scenes:\n",
" print(f\"⚠️ Failed to load {len(failed_scenes)} band instances (will be skipped)\")\n",
" \n",
" # ===== NORMALIZE RESOLUTION (No upsampling - just match to highest) =====\n",
" print(f\"\\n[RESOLUTION NORMALIZATION] Aligning all bands to native resolution (NO magnification)...\")\n",
" \n",
" # Find max resolution\n",
" ref_resolution = None\n",
" max_size = 0\n",
" max_band = None\n",
" \n",
" for scene_idx in all_data_dict.keys():\n",
" for band_name, data_band in all_data_dict[scene_idx].items():\n",
" size = data_band.shape[0]\n",
" if size > max_size:\n",
" max_size = size\n",
" ref_resolution = size\n",
" max_band = band_name\n",
" \n",
" print(f\" Reference resolution: {max_size}×{max_size} pixels (native {max_band})\")\n",
" \n",
" # Resample all bands to match reference resolution (both up and down)\n",
" resampled_count = 0\n",
" for scene_idx in all_data_dict.keys():\n",
" for band_name in list(all_data_dict[scene_idx].keys()):\n",
" band_data_arr = all_data_dict[scene_idx][band_name]\n",
" current_size = band_data_arr.shape[0]\n",
" \n",
" if current_size != ref_resolution:\n",
" scale_factor = ref_resolution / current_size\n",
" \n",
" # Resample to match reference resolution (both up and down)\n",
" if band_name == 'scl':\n",
" resampled_data = ndimage.zoom(band_data_arr, scale_factor, order=0)\n",
" else:\n",
" resampled_data = ndimage.zoom(band_data_arr, scale_factor, order=1)\n",
" \n",
" all_data_dict[scene_idx][band_name] = resampled_data\n",
" new_size = resampled_data.shape[0]\n",
" if scene_idx == 0: # Print for first scene only to reduce clutter\n",
" print(f\" Resampling {band_name}: {current_size}×{current_size} → {new_size}×{new_size}\")\n",
" resampled_count += 1\n",
" \n",
" print(f\"✅ Resolution normalization complete! ({resampled_count} bands resampled)\")\n",
" \n",
" # ===== CALCULATE SPECTRAL INDICES FOR EACH SCENE =====\n",
" print(f\"\\n[SPECTRAL INDICES] Calculating spectral indices for each scene...\")\n",
" \n",
" # Map band names: notebook dùng 'swir16'/'swir22' thay vì 'b11'/'b12'\n",
" indices_count = 0\n",
" for scene_idx in all_data_dict.keys():\n",
" scene_data = all_data_dict[scene_idx]\n",
" \n",
" try:\n",
" # NDVI: (NIR - Red) / (NIR + Red)\n",
" if 'nir' in scene_data and 'red' in scene_data:\n",
" nir = scene_data['nir'].astype(float)\n",
" red = scene_data['red'].astype(float)\n",
" ndvi = (nir - red) / (nir + red + 1e-8)\n",
" scene_data['ndvi'] = ndvi.astype(np.float32)\n",
" indices_count += 1\n",
" \n",
" # NDBI: (SWIR - NIR) / (SWIR + NIR) — dùng swir16 thay b11\n",
" swir_key = 'swir16' if 'swir16' in scene_data else ('b11' if 'b11' in scene_data else None)\n",
" if swir_key and 'nir' in scene_data:\n",
" swir = scene_data[swir_key].astype(float)\n",
" nir = scene_data['nir'].astype(float)\n",
" ndbi = (swir - nir) / (swir + nir + 1e-8)\n",
" scene_data['ndbi'] = ndbi.astype(np.float32)\n",
" indices_count += 1\n",
" \n",
" # NDWI: (NIR - SWIR) / (NIR + SWIR)\n",
" if swir_key and 'nir' in scene_data:\n",
" nir = scene_data['nir'].astype(float)\n",
" swir = scene_data[swir_key].astype(float)\n",
" ndwi = (nir - swir) / (nir + swir + 1e-8)\n",
" scene_data['ndwi'] = ndwi.astype(np.float32)\n",
" indices_count += 1\n",
" \n",
" # EVI: Enhanced Vegetation Index\n",
" if 'nir' in scene_data and 'red' in scene_data and 'blue' in scene_data:\n",
" nir = scene_data['nir'].astype(float)\n",
" red = scene_data['red'].astype(float)\n",
" blue = scene_data['blue'].astype(float)\n",
" evi = 2.5 * (nir - red) / (nir + 6*red - 7.5*blue + 1)\n",
" scene_data['evi'] = evi.astype(np.float32)\n",
" indices_count += 1\n",
" \n",
" except Exception as e:\n",
" pass\n",
" \n",
" print(f\"✅ Calculated {indices_count} spectral indices per scene\")\n",
" \n",
" # ===== STACK SCENES ALONG TIME DIMENSION =====\n",
" print(f\"\\n[STACKING] Stacking all {len(all_data_dict)} scenes to create time-series...\")\n",
" \n",
" data_vars = {}\n",
" band_names = list(all_data_dict[0].keys())\n",
" \n",
" for band_name in band_names:\n",
" band_data_list = []\n",
" for scene_idx in sorted(all_data_dict.keys()):\n",
" if band_name in all_data_dict[scene_idx]:\n",
" band_data_list.append(all_data_dict[scene_idx][band_name])\n",
" \n",
" if band_data_list:\n",
" stacked = np.stack(band_data_list, axis=0)\n",
" data_vars[band_name] = (['time', 'y', 'x'], stacked)\n",
" \n",
" # Create xarray Dataset with time dimension\n",
" first_band_data = list(all_data_dict[0].values())[0]\n",
" y_size, x_size = first_band_data.shape\n",
" \n",
" data = xr.Dataset(\n",
" data_vars,\n",
" coords={\n",
" 'time': np.arange(len(all_data_dict)),\n",
" 'x': np.arange(x_size),\n",
" 'y': np.arange(y_size)\n",
" }\n",
" )\n",
" \n",
" # ===== CALCULATE TEMPORAL FEATURES FOR ACCURACY =====\n",
" print(f\"\\n[TEMPORAL FEATURES] Computing temporal features from time-series...\")\n",
" \n",
" temporal_features_added = 0\n",
" \n",
" # For NDVI: temporal statistics\n",
" if 'ndvi' in data.data_vars:\n",
" ndvi_ts = data['ndvi']\n",
" \n",
" # Min NDVI (vegetation stress indicator)\n",
" data['ndvi_min'] = ndvi_ts.min(dim='time')\n",
" temporal_features_added += 1\n",
" \n",
" # Max NDVI (peak vegetation)\n",
" data['ndvi_max'] = ndvi_ts.max(dim='time')\n",
" temporal_features_added += 1\n",
" \n",
" # Mean NDVI\n",
" data['ndvi_mean'] = ndvi_ts.mean(dim='time')\n",
" temporal_features_added += 1\n",
" \n",
" # NDVI range (variability)\n",
" data['ndvi_range'] = data['ndvi_max'] - data['ndvi_min']\n",
" temporal_features_added += 1\n",
" \n",
" # NDVI std (temporal consistency)\n",
" data['ndvi_std'] = ndvi_ts.std(dim='time')\n",
" temporal_features_added += 1\n",
" \n",
" # For all indices: mean values (aggregate features)\n",
" for band_name in ['ndbi', 'ndwi', 'evi']:\n",
" if band_name in data.data_vars:\n",
" band_ts = data[band_name]\n",
" data[f'{band_name}_mean'] = band_ts.mean(dim='time')\n",
" temporal_features_added += 1\n",
" \n",
" print(f\"✅ Added {temporal_features_added} temporal/aggregate features\")\n",
" \n",
" # ===== SAVE TO CACHE =====\n",
" print(f\"\\n[CACHE] Saving dataset to cache...\")\n",
" # ===== BUILD QUERY =====\n",
" query = {\n",
" 'product': 's2_l2a',\n",
" 'time': date_range,\n",
" 'x': longtitude_range,\n",
" 'y': latitude_range,\n",
" }\n",
"\n",
" # Determine native CRS\n",
" print(f\"\\n[1] Finding native CRS...\")\n",
" try:\n",
" native_crs = notebook_utils.mostcommon_crs(dc, query)\n",
" except Exception:\n",
" native_crs = 'EPSG:32648' # UTM Zone 48N covers Vietnam\n",
" print(f\" CRS: {native_crs}\")\n",
"\n",
" # Bands to load (all main Sentinel-2 bands)\n",
" measurements = ['blue', 'green', 'red', 'nir', 'swir16', 'swir22',\n",
" 'coastal', 'rededge1', 'rededge2', 'rededge3', 'scl']\n",
"\n",
" load_params = {\n",
" 'measurements': measurements,\n",
" 'output_crs': native_crs,\n",
" 'resolution': (-10, 10), # 10m native resolution\n",
" 'group_by': 'solar_day', # Group multiple passes same day\n",
" 'dask_chunks': {'x': 2048, 'y': 2048},\n",
" }\n",
"\n",
" print(f\"\\n[2] Loading data via dc.load() ...\")\n",
" print(f\" Region : lon={longtitude_range}, lat={latitude_range}\")\n",
" print(f\" Time : {date_range[0]} → {date_range[1]}\")\n",
" print(f\" Bands : {measurements}\")\n",
"\n",
" # ---- ODC load — handles S3 credentials internally ----\n",
" data = load_s2l2a_with_offset(dc, query | load_params)\n",
"\n",
" if data is None or len(data.time) == 0:\n",
" raise ValueError(\"dc.load() returned no data for this region/time range\")\n",
"\n",
" print(f\"\\n✅ Loaded {len(data.time)} scenes via dc.load()\")\n",
" print(f\" Dimensions: {dict(data.dims)}\")\n",
"\n",
" # Compute (trigger Dask)\n",
" print(f\"\\n[3] Computing (triggering Dask download)...\")\n",
" data = data.compute()\n",
" print(f\"✅ Data computed!\")\n",
"\n",
" # ===== CALCULATE SPECTRAL INDICES =====\n",
" print(f\"\\n[4] Calculating spectral indices...\")\n",
"\n",
" # NDVI\n",
" if 'nir' in data and 'red' in data:\n",
" nir = data['nir'].astype(float)\n",
" red = data['red'].astype(float)\n",
" data['ndvi'] = ((nir - red) / (nir + red + 1e-8)).astype(np.float32)\n",
" print(f\" ✅ NDVI\")\n",
"\n",
" # NDBI (uses swir16)\n",
" if 'swir16' in data and 'nir' in data:\n",
" swir = data['swir16'].astype(float)\n",
" nir = data['nir'].astype(float)\n",
" data['ndbi'] = ((swir - nir) / (swir + nir + 1e-8)).astype(np.float32)\n",
" print(f\" ✅ NDBI\")\n",
"\n",
" # NDWI (uses swir16)\n",
" if 'swir16' in data and 'nir' in data:\n",
" nir = data['nir'].astype(float)\n",
" swir = data['swir16'].astype(float)\n",
" data['ndwi'] = ((nir - swir) / (nir + swir + 1e-8)).astype(np.float32)\n",
" print(f\" ✅ NDWI\")\n",
"\n",
" # EVI\n",
" if 'nir' in data and 'red' in data and 'blue' in data:\n",
" nir = data['nir'].astype(float)\n",
" red = data['red'].astype(float)\n",
" blue = data['blue'].astype(float)\n",
" data['evi'] = (2.5 * (nir - red) / (nir + 6*red - 7.5*blue + 1)).astype(np.float32)\n",
" print(f\" ✅ EVI\")\n",
"\n",
" # ===== TEMPORAL AGGREGATE FEATURES =====\n",
" print(f\"\\n[5] Computing temporal aggregate features...\")\n",
" added = 0\n",
"\n",
" if 'ndvi' in data:\n",
" ndvi_ts = data['ndvi']\n",
" data['ndvi_min'] = ndvi_ts.min(dim='time'); added += 1\n",
" data['ndvi_max'] = ndvi_ts.max(dim='time'); added += 1\n",
" data['ndvi_mean'] = ndvi_ts.mean(dim='time'); added += 1\n",
" data['ndvi_std'] = ndvi_ts.std(dim='time'); added += 1\n",
" data['ndvi_range'] = data['ndvi_max'] - data['ndvi_min']; added += 1\n",
"\n",
" for idx_name in ['ndbi', 'ndwi', 'evi']:\n",
" if idx_name in data:\n",
" data[f'{idx_name}_mean'] = data[idx_name].mean(dim='time')\n",
" added += 1\n",
"\n",
" print(f\" ✅ {added} aggregate features added\")\n",
"\n",
" # ===== SAVE TO CACHE =====\n",
" print(f\"\\n[6] Saving to cache: {cache_file}\")\n",
" try:\n",
" import os\n",
" os.makedirs(cache_dir, exist_ok=True)\n",
" data.to_netcdf(cache_file, engine='netcdf4')\n",
" cache_size = os.path.getsize(cache_file) / (1024**3)\n",
" print(f\"✅ Dataset saved to cache: {cache_file}\")\n",
" print(f\" Cache size: {cache_size:.2f} GB\")\n",
" size_gb = os.path.getsize(cache_file) / (1024**3)\n",
" print(f\" ✅ Cached ({size_gb:.2f} GB)\")\n",
" except Exception as e:\n",
" print(f\"⚠️ Error saving cache: {e}\")\n",
" \n",
" print(f\"\\n✅ OPTIMAL Dataset with native resolution + temporal features created!\")\n",
" print(f\" {'='*70}\")\n",
" print(f\" 🎬 Total scenes (time steps): {len(all_data_dict)}\")\n",
" print(f\" 📊 Total bands/variables: {len(data.data_vars)}\")\n",
" print(f\" 🖼️ Spatial size: {x_size} × {y_size} pixels (NATIVE resolution)\")\n",
" print(f\" 📏 Native resolution: 10m (Sentinel-2 L2A)\")\n",
" print(f\" ⏰ Temporal range: {scene_dates[0].date()} to {scene_dates[-1].date()}\")\n",
" print(f\" 💾 Total dataset size: {notebook_utils.xarray_object_size(data)}\")\n",
" print(f\" 💿 Cached at: {cache_file}\")\n",
" print(f\" {'='*70}\")\n",
" \n",
" print(f\"\\n Dataset dimensions:\")\n",
" for dim, size in data.dims.items():\n",
" print(f\" {dim}: {size}\")\n",
" \n",
" print(f\"\\n Variables ({len(data.data_vars)}):\")\n",
" spatial_vars = []\n",
" temporal_vars = []\n",
" for var_name in sorted(data.data_vars):\n",
" if len(data[var_name].shape) == 3:\n",
" spatial_vars.append(f\"{var_name} {data[var_name].shape}\")\n",
" else:\n",
" temporal_vars.append(f\"{var_name} {data[var_name].shape}\")\n",
" \n",
" print(f\" Spatial time-series ({len(spatial_vars)}):\")\n",
" for v in spatial_vars:\n",
" print(f\" - {v}\")\n",
" print(f\" Temporal aggregates ({len(temporal_vars)}):\")\n",
" for v in temporal_vars:\n",
" print(f\" - {v}\")\n",
" \n",
" print(f\" {'='*70}\")\n",
" \n",
" print(f\" ⚠️ Cache save failed: {e}\")\n",
"\n",
" print(f\"\\n{'='*70}\")\n",
" print(f\"✅ Dataset ready!\")\n",
" print(f\" Scenes : {len(data.time)}\")\n",
" print(f\" Variables : {len(data.data_vars)}\")\n",
" print(f\" Spatial : {data.dims.get('x',0)} × {data.dims.get('y',0)} px @ 10m\")\n",
" print(f\"{'='*70}\")\n",
" display(data)\n",
" \n",
"\n",
" # ===== EXTRACT NDVI FOR TRAINING =====\n",
" print(f\"\\n[NDVI EXTRACTION] Extracting NDVI for model training...\")\n",
" if 'ndvi_mean' in data.data_vars:\n",
" # Use mean NDVI across time\n",
" ndvi = data['ndvi_mean']\n",
" print(f\"✅ NDVI extracted (mean across time)\")\n",
" print(f\" Shape: {ndvi.shape}\")\n",
" print(f\"\\n✅ NDVI (mean) shape: {ndvi.shape}\")\n",
" elif 'ndvi' in data.data_vars:\n",
" # Use first time step if mean not available\n",
" ndvi = data['ndvi'].isel(time=0)\n",
" print(f\"✅ NDVI extracted (first time step)\")\n",
" print(f\" Shape: {ndvi.shape}\")\n",
" print(f\"\\n✅ NDVI (t=0) shape: {ndvi.shape}\")\n",
" else:\n",
" print(f\"❌ NDVI not found in dataset\")\n",
" ndvi = None\n",
" \n",
" print(f\"\\n❌ NDVI not found in dataset\")\n",
"\n",
"except Exception as e:\n",
" print(f\"❌ Error: {e}\")\n",
" import traceback\n",
Executable
+3
View File
@@ -0,0 +1,3 @@
git add .
git commit -m "update"
git push origin dev_01