diff --git a/01.train_ODC.ipynb b/01.train_ODC.ipynb index 533460c..96cf2b0 100755 --- a/01.train_ODC.ipynb +++ b/01.train_ODC.ipynb @@ -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", diff --git a/git.sh b/git.sh new file mode 100755 index 0000000..59e3236 --- /dev/null +++ b/git.sh @@ -0,0 +1,3 @@ +git add . +git commit -m "update" +git push origin dev_01