{ "cells": [ { "cell_type": "code", "execution_count": null, "id": "912ed572-1658-406b-976c-cd6de2d4e89e", "metadata": { "tags": [] }, "outputs": [ { "ename": "ModuleNotFoundError", "evalue": "No module named 'easi_tools'", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", "File \u001b[0;32m:4\u001b[0m\n", "File \u001b[0;32m~/CSIROBoeingPhase5-Vietnam/new_import_ODC.py:23\u001b[0m\n\u001b[1;32m 21\u001b[0m easinotebooksrepo \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m/home/jovyan/easi-notebooks\u001b[39m\u001b[38;5;124m'\u001b[39m\n\u001b[1;32m 22\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m easinotebooksrepo \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m sys\u001b[38;5;241m.\u001b[39mpath: sys\u001b[38;5;241m.\u001b[39mpath\u001b[38;5;241m.\u001b[39mappend(easinotebooksrepo)\n\u001b[0;32m---> 23\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21;01measi_tools\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m EasiDefaults, xarray_object_size, notebook_utils, unset_cachingproxy\n\u001b[1;32m 24\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21;01measi_tools\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mload_s2l2a\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m load_s2l2a_with_offset\n\u001b[1;32m 25\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21;01mdask\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mdistributed\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m progress\n", "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'easi_tools'" ] } ], "source": [ "%%time\n", "%matplotlib inline\n", "\n", "import importlib\n", "import new_import_ODC \n", "\n", "importlib.reload(new_import_ODC)\n", "\n", "from new_import_ODC import *\n", "\n", "print(\"✅ All modules loaded successfully\")" ] }, { "cell_type": "code", "execution_count": 2, "id": "d824dc4f-994b-4d1c-8d24-ce6674da141c", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "✅ AWS credentials loaded from environment variables\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/home/x79/miniconda/envs/env_01/lib/python3.10/site-packages/distributed/node.py:187: UserWarning: Port 8787 is already in use.\n", "Perhaps you already have a cluster running?\n", "Hosting the HTTP server on port 41709 instead\n", " warnings.warn(\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "✅ Dask cluster initialized\n", " Cluster: LocalCluster(9f2167a3, 'tcp://127.0.0.1:41233', workers=4, threads=24, memory=31.26 GiB)\n", "✅ Datacube connected (metadata only)\n", "\n", "======================================================================\n", "CPU times: user 4.98 s, sys: 831 ms, total: 5.81 s\n", "Wall time: 7.69 s\n" ] } ], "source": [ "%%time\n", "import os\n", "import sys\n", "\n", "print(\"✅ AWS credentials loaded from environment variables\")\n", "\n", "# Cấu hình Dask local\n", "from dask.distributed import Client, LocalCluster\n", "\n", "cluster = LocalCluster(n_workers=4)\n", "client = Client(cluster)\n", "print(\"✅ Dask cluster initialized\")\n", "print(f\" Cluster: {cluster}\")\n", "\n", "# Khai báo Datacube (chỉ để lấy metadata, không dùng load())\n", "import datacube\n", "try:\n", " dc = datacube.Datacube()\n", " print(\"✅ Datacube connected (metadata only)\")\n", "except Exception as e:\n", " print(f\"⚠️ Datacube connection not critical: {e}\")\n", " dc = None\n", "\n", "print(\"\\n\" + \"=\"*70)" ] }, { "cell_type": "code", "execution_count": 3, "id": "1e113730", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "======================================================================\n", "GETTING SENTINEL-2 SCENE METADATA\n", "======================================================================\n", "\n", "[1] Loading metadata from datacube...\n", " ✅ Found 40 scenes\n", "\n", "[2] Selected scene: S2A_48PWR_20231226_0_L2A\n", " Date: 2023-12-26 03:35:26.919000+00:00\n", "\n", "[3] Available bands:\n", " - nir: https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/48/P/WR/20\n", " - red: https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/48/P/WR/20\n", " - scl: https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/48/P/WR/20\n", " - blue: https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/48/P/WR/20\n", " - green: https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/48/P/WR/20\n", " - nir08: https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/48/P/WR/20\n", " - nir09: https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/48/P/WR/20\n", " - swir16: https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/48/P/WR/20\n", " - swir22: https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/48/P/WR/20\n", " - coastal: https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/48/P/WR/20\n", " - rededge1: https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/48/P/WR/20\n", " - rededge2: https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/48/P/WR/20\n", " - rededge3: https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/48/P/WR/20\n", "======================================================================\n", "CPU times: user 3.34 s, sys: 76.3 ms, total: 3.41 s\n", "Wall time: 3.23 s\n" ] } ], "source": [ "%%time\n", "# 🔧 Get Sentinel-2 scene metadata from datacube\n", "print(\"=\"*70)\n", "print(\"GETTING SENTINEL-2 SCENE METADATA\")\n", "print(\"=\"*70)\n", "\n", "date_range = (\"2023-03-01\", \"2023-12-31\")\n", "longtitude_range = (105.5, 106.4)\n", "latitude_range = (9.2, 10.0)\n", "\n", "try:\n", " print(f\"\\n[1] Loading metadata from datacube...\")\n", " datasets = list(dc.find_datasets(product='s2_l2a', time=date_range))\n", " print(f\" ✅ Found {len(datasets)} scenes\")\n", " \n", " if datasets:\n", " selected = datasets[0]\n", " print(f\"\\n[2] Selected scene: {selected.metadata.label}\")\n", " scene_datetime = selected.time.begin if hasattr(selected.time, 'begin') else selected.time\n", " print(f\" Date: {scene_datetime}\")\n", " \n", " # Display measurement paths\n", " print(f\"\\n[3] Available bands:\")\n", " for name, measurement in selected.measurements.items():\n", " print(f\" - {name}: {measurement['path'][:80]}\")\n", " \n", "except Exception as e:\n", " print(f\"❌ Error: {e}\")\n", " import traceback\n", " traceback.print_exc()\n", "\n", "print(\"=\"*70)" ] }, { "cell_type": "code", "execution_count": 4, "id": "3cd69645", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "======================================================================\n", "CHECKING FOR CACHED DATASET\n", "======================================================================\n", "\n", "⏳ Cache file not found: dataset_cache/sentinel2_timeseries_40scenes.nc\n", " Will download from S3 and save cache\n", " (Next run will use cache automatically)\n", "======================================================================\n", "CPU times: user 4.36 ms, sys: 3.7 ms, total: 8.06 ms\n", "Wall time: 7.22 ms\n" ] } ], "source": [ "%%time\n", "# 🔍 CHECK IF DATASET CACHE EXISTS (Skip download if available)\n", "print(\"=\"*70)\n", "print(\"CHECKING FOR CACHED DATASET\")\n", "print(\"=\"*70)\n", "\n", "import os\n", "import xarray as xr\n", "\n", "cache_dir = \"dataset_cache\"\n", "cache_file = f\"{cache_dir}/sentinel2_timeseries_40scenes.nc\"\n", "\n", "use_cache = False\n", "\n", "if os.path.exists(cache_file):\n", " print(f\"\\n✅ Cache file found: {cache_file}\")\n", " \n", " # Get file info\n", " file_size_gb = os.path.getsize(cache_file) / (1024**3)\n", " print(f\" File size: {file_size_gb:.2f} GB\")\n", " \n", " # Try to load\n", " try:\n", " print(f\"\\n🔄 Loading dataset from cache...\")\n", " data = xr.open_dataset(cache_file)\n", " \n", " print(f\"✅ Dataset loaded from cache!\")\n", " print(f\" Total scenes: {len(data['time'])}\")\n", " print(f\" Variables: {len(data.data_vars)}\")\n", " print(f\" Dimensions: {dict(data.dims)}\")\n", " print(f\"\\n ⏭️ Skipping S3 download (using cached data)\")\n", " \n", " use_cache = True\n", " \n", " except Exception as e:\n", " print(f\"❌ Error loading cache: {e}\")\n", " print(f\" Will download fresh data from S3\")\n", " use_cache = False\n", "else:\n", " print(f\"\\n⏳ Cache file not found: {cache_file}\")\n", " print(f\" Will download from S3 and save cache\")\n", " print(f\" (Next run will use cache automatically)\")\n", "\n", "print(\"=\"*70)" ] }, { "cell_type": "code", "execution_count": 5, "id": "435f9f78-a9a4-4226-86ca-d4bec42d454e", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "======================================================================\n", "LOADING SENTINEL-2 DATA FROM S3 COGs (RASTERIO) - OPTIMAL ACCURACY\n", "======================================================================\n", "\n", "📥 Downloading from S3...\n", "\n", "📦 Found 40 available scenes\n", " Date range: 2023-03-01 to 2023-12-31\n", "\n", "[LOADING] Loading ALL 40 scenes with ALL available bands...\n", " (Keeping NATIVE resolution - NO upsampling/magnification)\n", " Available bands: ['nir', 'red', 'scl', 'blue', 'green', 'nir08', 'nir09', 'swir16', 'swir22', 'coastal', 'rededge1', 'rededge2', 'rededge3']\n", "\n", " [ 1/1] S2A_48PWR_20231226_0_L2A (2023-12-26)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " ✅ 13 bands loaded\n", "\n", "✅ Successfully loaded 1 scenes!\n", "\n", "[RESOLUTION NORMALIZATION] Aligning all bands to native resolution (NO magnification)...\n", " Reference resolution: 10980×10980 pixels (native nir)\n", " Resampling scl: 5490×5490 → 10980×10980\n", " Resampling nir08: 5490×5490 → 10980×10980\n", " Resampling nir09: 1830×1830 → 10980×10980\n", " Resampling swir16: 5490×5490 → 10980×10980\n", " Resampling swir22: 5490×5490 → 10980×10980\n", " Resampling coastal: 1830×1830 → 10980×10980\n", " Resampling rededge1: 5490×5490 → 10980×10980\n", " Resampling rededge2: 5490×5490 → 10980×10980\n", " Resampling rededge3: 5490×5490 → 10980×10980\n", "✅ Resolution normalization complete! (9 bands resampled)\n", "\n", "[SPECTRAL INDICES] Calculating spectral indices for each scene...\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ ":169: RuntimeWarning: divide by zero encountered in divide\n", ":169: RuntimeWarning: invalid value encountered in divide\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "✅ Calculated 2 spectral indices per scene\n", "\n", "[STACKING] Stacking all 1 scenes to create time-series...\n", "\n", "[TEMPORAL FEATURES] Computing temporal features from time-series...\n", "✅ Added 6 temporal/aggregate features\n", "\n", "[CACHE] Saving dataset to cache...\n", "✅ Dataset saved to cache: dataset_cache/sentinel2_timeseries_40scenes.nc\n", " Cache size: 6.40 GB\n", "\n", "✅ OPTIMAL Dataset with native resolution + temporal features created!\n", " ======================================================================\n", " 🎬 Total scenes (time steps): 1\n", " 📊 Total bands/variables: 21\n", " 🖼️ Spatial size: 10980 × 10980 pixels (NATIVE resolution)\n", " 📏 Native resolution: 10m (Sentinel-2 L2A)\n", " ⏰ Temporal range: 2023-12-26 to 2023-12-26\n", "❌ Error: name 'notebook_utils' is not defined\n", "======================================================================\n", "CPU times: user 9min, sys: 4min 16s, total: 13min 16s\n", "Wall time: 16min 42s\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Traceback (most recent call last):\n", " File \"\", line 262, in \n", "NameError: name 'notebook_utils' is not defined\n" ] } ], "source": [ "%%time\n", "# 💾 LOAD SENTINEL-2 DATA DIRECTLY FROM S3 COGS USING RASTERIO - WITH TEMPORAL FEATURES\n", "print(\"=\"*70)\n", "print(\"LOADING SENTINEL-2 DATA FROM S3 COGs (RASTERIO) - OPTIMAL ACCURACY\")\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", " # ===== CHECK IF SHOULD SKIP DOWNLOAD =====\n", " if use_cache 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", " else:\n", " # ===== DOWNLOAD FROM S3 =====\n", " print(f\"\\n📥 Downloading from S3...\")\n", " \n", " # Get all scenes from datacube metadata\n", " datasets = list(dc.find_datasets(\n", " product='s2_l2a',\n", " time=date_range\n", " ))\n", " \n", " if not datasets:\n", " raise ValueError(\"No datasets found for date range\")\n", " \n", " print(f\"\\n📦 Found {len(datasets)} available scenes\")\n", " print(f\" Date range: {date_range[0]} to {date_range[1]}\")\n", " \n", " # ===== LOAD ALL SCENES WITH ALL AVAILABLE BANDS (NO MAGNIFICATION) =====\n", " print(f\"\\n[LOADING] Loading ALL {len(datasets)} scenes with ALL available bands...\")\n", " print(f\" (Keeping NATIVE resolution - NO upsampling/magnification)\")\n", " \n", " # num_scenes = len(datasets) # Load ALL scenes\n", " num_scenes = 1 # Load ALL scenes\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)[:30]}\")\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", " 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)\n", " if 'b11' in scene_data and 'nir' in scene_data:\n", " swir = scene_data['b11'].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 'nir' in scene_data and 'b11' in scene_data:\n", " nir = scene_data['nir'].astype(float)\n", " swir = scene_data['b11'].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", " try:\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", " 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", " display(data)\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", " 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", " else:\n", " print(f\"❌ NDVI not found in dataset\")\n", " ndvi = None\n", " \n", "except Exception as e:\n", " print(f\"❌ Error: {e}\")\n", " import traceback\n", " traceback.print_exc()\n", " data = None\n", " ndvi = None\n", "\n", "print(\"=\"*70)" ] }, { "cell_type": "code", "execution_count": 6, "id": "d2585562-88aa-4c7d-bf70-1f6affcf65d4", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "======================================================================\n", "TRAINING DATA SETUP\n", "======================================================================\n", "\n", "[1] Loading training data: train/ST_training data_updated_1130points_new.shp\n", " ❌ Error: name 'load_train_data' is not defined\n", "\n", "[2] Label mapping:\n", " 0: Lua tom\n", " 1: Lua\n", " 2: CHN\n", " 3: CLN\n", " 4: TS\n", " 5: Song\n", " 6: Dat xay dung\n", " 7: Rung\n", "\n", "======================================================================\n" ] } ], "source": [ "# 🎯 LOAD TRAINING DATA & EXTRACT FEATURES\n", "print(\"=\"*70)\n", "print(\"TRAINING DATA SETUP\")\n", "print(\"=\"*70)\n", "\n", "# Load training points\n", "train_path = \"train/ST_training data_updated_1130points_new.shp\"\n", "print(f\"\\n[1] Loading training data: {train_path}\")\n", "\n", "try:\n", " train = load_train_data(train_path)\n", " print(f\" ✅ Loaded {len(train)} training points\")\n", " print(f\" Columns: {list(train.columns)}\")\n", " train.head()\n", "except Exception as e:\n", " print(f\" ❌ Error: {e}\")\n", " train = None\n", "\n", "# Label mapping\n", "label_mapping = {\n", " \"Lua tom\": \"0\",\n", " \"Lua\": \"1\",\n", " \"CHN\": \"2\",\n", " \"CLN\": \"3\",\n", " \"TS\": \"4\",\n", " \"Song\": \"5\",\n", " \"Dat xay dung\": \"6\",\n", " \"Rung\": \"7\",\n", "}\n", "\n", "print(f\"\\n[2] Label mapping:\")\n", "for label, code in label_mapping.items():\n", " print(f\" {code}: {label}\")\n", "\n", "print(\"\\n\" + \"=\"*70)" ] }, { "cell_type": "code", "execution_count": null, "id": "2e955884-d4af-422d-a8e6-d436199540e0", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "======================================================================\n", "MODEL TRAINING\n", "======================================================================\n", "❌ Missing training data or NDVI\n", "======================================================================\n", "CPU times: user 700 μs, sys: 0 ns, total: 700 μs\n", "Wall time: 638 μs\n" ] } ], "source": [ "%%time\n", "# 🤖 LAND USE CLASSIFICATION MODEL TRAINING (MỤC TIÊU CHÍNH)\n", "print(\"=\"*70)\n", "print(\"LAND USE CLASSIFICATION TRAINING\")\n", "print(\"=\"*70)\n", "print(\"\\n🎯 Mục tiêu: Dự đoán phân loại sử dụng đất (8 lớp)\")\n", "print(\" - NDVI/NDWI/NDBI/EVI là INPUT FEATURES\")\n", "print(\" - Sau khi predict xong → có thể hiển thị NDVI map như chỉ số phụ\")\n", "print(\"=\"*70)\n", "\n", "if train is not None and data is not None:\n", " print(\"\\n[1] Extracting MULTIPLE features from satellite data...\")\n", " print(\" (Sử dụng nhiều spectral indices để cải thiện accuracy)\")\n", " \n", " try:\n", " # Extract features at training point locations\n", " X = []\n", " y = []\n", " \n", " # Available features from data\n", " available_features = ['ndvi_mean', 'ndvi_min', 'ndvi_max', 'ndvi_std', 'ndvi_range',\n", " 'ndwi_mean', 'ndbi_mean', 'evi_mean']\n", " \n", " # Check which features are actually available\n", " features_to_use = [f for f in available_features if f in data.data_vars]\n", " \n", " if not features_to_use:\n", " print(\" ❌ No spectral features found in dataset!\")\n", " print(\" Available variables:\", list(data.data_vars))\n", " model = None\n", " else:\n", " print(f\" Using {len(features_to_use)} features: {features_to_use}\")\n", " \n", " for idx, point in train.iterrows():\n", " try:\n", " # Extract all available features at this point\n", " feature_vec = []\n", " for feat_name in features_to_use:\n", " feat_val = float(data[feat_name].sel(\n", " x=point.geometry.x, \n", " y=point.geometry.y, \n", " method='nearest'\n", " ).values)\n", " feature_vec.append(feat_val)\n", " \n", " # Get label\n", " label = label_mapping[point.Hientrang]\n", " \n", " # Only add if no NaN values\n", " if not np.isnan(feature_vec).any():\n", " X.append(feature_vec)\n", " y.append(int(label))\n", " except Exception as e:\n", " # Skip points with errors\n", " continue\n", " \n", " if len(X) > 0:\n", " X = np.array(X)\n", " y = np.array(y)\n", " print(f\" ✅ Extracted {len(X)} samples with {X.shape[1]} features each\")\n", " \n", " # Show feature statistics\n", " print(f\"\\n Feature statistics:\")\n", " for i, feat_name in enumerate(features_to_use):\n", " print(f\" {feat_name:15s}: mean={X[:,i].mean():.3f}, std={X[:,i].std():.3f}\")\n", " \n", " # Split data\n", " print(f\"\\n[2] Splitting data (80-20)...\")\n", " from sklearn.model_selection import train_test_split\n", " X_train, X_test, y_train, y_test = train_test_split(\n", " X, y, test_size=0.2, random_state=42, stratify=y\n", " )\n", " print(f\" Train: {len(X_train)}, Test: {len(X_test)}\")\n", " \n", " # Show class distribution\n", " unique, counts = np.unique(y_train, return_counts=True)\n", " print(f\"\\n Class distribution in training set:\")\n", " for cls, count in zip(unique, counts):\n", " cls_name = [k for k, v in label_mapping.items() if v == str(cls)][0]\n", " print(f\" {cls}: {cls_name:15s} - {count:4d} samples ({count/len(y_train)*100:.1f}%)\")\n", " \n", " # Train model\n", " print(f\"\\n[3] Training Random Forest for LAND USE CLASSIFICATION...\")\n", " from sklearn.ensemble import RandomForestClassifier\n", " from sklearn.metrics import accuracy_score, classification_report\n", " \n", " model = RandomForestClassifier(\n", " n_estimators=200, # More trees for better accuracy\n", " max_depth=30,\n", " min_samples_split=5,\n", " random_state=42,\n", " n_jobs=-1,\n", " verbose=1\n", " )\n", " model.fit(X_train, y_train)\n", " \n", " # Evaluate\n", " y_pred = model.predict(X_test)\n", " accuracy = accuracy_score(y_test, y_pred)\n", " \n", " print(f\"\\n ✅ Model trained!\")\n", " print(f\" Training accuracy: {model.score(X_train, y_train)*100:.2f}%\")\n", " print(f\" Testing accuracy: {accuracy*100:.2f}%\")\n", " \n", " # Show feature importance\n", " print(f\"\\n Feature importance:\")\n", " importances = model.feature_importances_\n", " indices = np.argsort(importances)[::-1]\n", " for i, idx in enumerate(indices):\n", " print(f\" {i+1}. {features_to_use[idx]:15s}: {importances[idx]:.4f}\")\n", " \n", " # Classification report\n", " print(f\"\\n[4] Classification Report:\")\n", " class_names = [k for k, v in sorted(label_mapping.items(), key=lambda x: x[1])]\n", " print(classification_report(y_test, y_pred, target_names=class_names, zero_division=0))\n", " \n", " else:\n", " print(f\" ❌ No samples extracted\")\n", " model = None\n", " \n", " except Exception as e:\n", " print(f\" ❌ Error: {e}\")\n", " import traceback\n", " traceback.print_exc()\n", " model = None\n", "else:\n", " print(\"❌ Missing training data or satellite data\")\n", " model = None\n", "\n", "print(\"\\n\" + \"=\"*70)\n", "print(\"📝 NOTE: Model này dự đoán PHÂN LOẠI SỬ DỤNG ĐẤT (8 lớp)\")\n", "print(\" NDVI là một trong các features đầu vào, không phải mục tiêu dự đoán\")\n", "print(\" Sau khi predict → có thể hiển thị NDVI map như chỉ số phụ\")\n", "print(\"=\"*70)" ] }, { "cell_type": "code", "execution_count": null, "id": "f1a14379-ed6e-4897-9ca4-2669743fab40", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "======================================================================\n", "MODEL SAVING\n", "======================================================================\n", "❌ No model to save\n", "======================================================================\n" ] } ], "source": [ "# 💾 SAVE MODEL WITH METADATA\n", "print(\"=\"*70)\n", "print(\"MODEL SAVING\")\n", "print(\"=\"*70)\n", "\n", "if model is not None:\n", " print(\"\\n🔄 Saving trained LAND USE CLASSIFICATION model with metadata...\")\n", " try:\n", " from datetime import datetime\n", " \n", " # Prepare metadata for ModelManager\n", " metadata = {\n", " \"timestamp\": datetime.now().isoformat(),\n", " \"data_source\": \"Local S3 ODC (Open Data Cube)\",\n", " \"collections\": [\"sentinel-2-l2a\"],\n", " \"features\": features_to_use, # All features used\n", " \"feature_mode\": \"extended\", # Using extended aggregate features\n", " \"training_samples\": len(X_train),\n", " \"testing_samples\": len(X_test),\n", " \"test_size\": 0.2,\n", " \"train_accuracy\": float(model.score(X_train, y_train)),\n", " \"test_accuracy\": float(accuracy),\n", " \"model_type\": \"random_forest\",\n", " \"device\": \"cpu\",\n", " \"n_estimators\": 200,\n", " \"max_depth\": 30,\n", " \"learning_rate\": None,\n", " \"cnn_epochs\": None,\n", " \"n_features\": X_train.shape[1],\n", " \"n_classes\": len(np.unique(y)),\n", " \"class_names\": list(label_mapping.keys()),\n", " \"classification_report\": classification_report(y_test, y_pred, \n", " target_names=class_names, \n", " output_dict=True,\n", " zero_division=0),\n", " \"bbox\": None,\n", " \"time_range\": f\"{date_range[0]}/{date_range[1]}\",\n", " \"resolution\": 10,\n", " \"notes\": \"LAND USE CLASSIFICATION model trained from 01.train_ODC.ipynb. Predicts 8 land use classes using multiple spectral indices. NDVI is one of the input features, not the prediction target.\"\n", " }\n", " \n", " # Save model with metadata using updated save_model function\n", " save_model(\"model_land_use_odc.joblib\", model, metadata=metadata, label_encoder=None)\n", " \n", " print(\"✅ Model saved to model_train/model_land_use_odc.joblib\")\n", " print(f\" - Purpose: Land Use Classification (8 classes)\")\n", " print(f\" - Features: {len(features_to_use)} ({', '.join(features_to_use[:3])}...)\")\n", " print(f\" - Train Accuracy: {metadata['train_accuracy']*100:.2f}%\")\n", " print(f\" - Test Accuracy: {metadata['test_accuracy']*100:.2f}%\")\n", " print(f\" - Classes: {metadata['n_classes']}\")\n", " print(f\"\\n📝 NDVI là một trong các features, không phải prediction target\")\n", " print(f\" Sau khi predict → có thể tính NDVI map riêng để hiển thị\")\n", " except Exception as e:\n", " print(f\"❌ Error saving model: {e}\")\n", " import traceback\n", " traceback.print_exc()\n", "else:\n", " print(\"❌ No model to save\")\n", "\n", "print(\"=\"*70)" ] }, { "cell_type": "markdown", "id": "4a8579f4", "metadata": {}, "source": [ "# 📖 Hướng dẫn sử dụng Model\n", "\n", "## Mục đích của Model\n", "\n", "Model này được train để **DỰ ĐOÁN PHÂN LOẠI SỬ DỤNG ĐẤT** với 8 lớp:\n", "\n", "1. **Lua tom** (0) - Lúa tôm\n", "2. **Lua** (1) - Lúa\n", "3. **CHN** (2) - Cây hàng năm\n", "4. **CLN** (3) - Cây lâu năm \n", "5. **TS** (4) - Thủy sản\n", "6. **Song** (5) - Sông\n", "7. **Dat xay dung** (6) - Đất xây dựng\n", "8. **Rung** (7) - Rừng\n", "\n", "## Features đầu vào\n", "\n", "Model sử dụng **nhiều spectral indices** làm features:\n", "- NDVI (mean, min, max, std, range)\n", "- NDWI (mean)\n", "- NDBI (mean)\n", "- EVI (mean)\n", "\n", "## NDVI là gì trong hệ thống này?\n", "\n", "⚠️ **QUAN TRỌNG**: NDVI **KHÔNG PHẢI** là mục tiêu dự đoán!\n", "\n", "- **NDVI là INPUT FEATURE**: Một trong các chỉ số dùng để train model\n", "- **Mục tiêu dự đoán**: Phân loại sử dụng đất (8 lớp)\n", "- **NDVI map**: Có thể hiển thị NDVI map như chỉ số phụ sau khi predict xong\n", "\n", "## Workflow Prediction\n", "\n", "```python\n", "# 1. Load model\n", "model, label_encoder, metadata = model_manager.load_model(\"model_land_use_odc.joblib\")\n", "\n", "# 2. Extract features từ satellite data\n", "features = extract_features(satellite_data) # NDVI, NDWI, NDBI, EVI\n", "\n", "# 3. Predict land use classification\n", "land_use_prediction = model.predict(features)\n", "# → Kết quả: Mảng với giá trị 0-7 (8 lớp sử dụng đất)\n", "\n", "# 4. (Optional) Tính NDVI map riêng để hiển thị\n", "ndvi_map = (NIR - Red) / (NIR + Red)\n", "# → NDVI map chỉ để visualize, không phải prediction target\n", "```\n", "\n", "## So sánh với approach cũ\n", "\n", "| Approach | Features | Target | NDVI Role |\n", "|----------|----------|--------|-----------|\n", "| ❌ Cũ (sai) | Chỉ NDVI | 8 lớp đất | Input duy nhất |\n", "| ✅ Mới (đúng) | NDVI + NDWI + NDBI + EVI | 8 lớp đất | Một trong nhiều features |\n", "\n", "## Test Model\n", "\n", "```python\n", "# Test với website\n", "# 1. Upload model_land_use_odc.joblib lên server\n", "# 2. Chọn model trong prediction interface\n", "# 3. Chọn vùng và thời gian\n", "# 4. System sẽ tự động:\n", "# - Extract features (NDVI, NDWI, NDBI, EVI)\n", "# - Predict land use classification\n", "# - (Optional) Generate NDVI visualization map\n", "```" ] }, { "cell_type": "code", "execution_count": null, "id": "33dd516d-9824-499e-96b9-5cd9224c194c", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "======================================================================\n", "CLEANUP\n", "======================================================================\n", "\n", "🔄 Closing Dask client and cluster...\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "✅ Cleanup complete\n", "\n", "======================================================================\n", "✅ PIPELINE COMPLETE\n", "======================================================================\n" ] }, { "ename": "", "evalue": "", "output_type": "error", "traceback": [ "\u001b[1;31mThe Kernel crashed while executing code in the current cell or a previous cell. \n", "\u001b[1;31mPlease review the code in the cell(s) to identify a possible cause of the failure. \n", "\u001b[1;31mClick here for more info. \n", "\u001b[1;31mView Jupyter log for further details." ] } ], "source": [ "# 🛑 CLEANUP\n", "print(\"=\"*70)\n", "print(\"CLEANUP\")\n", "print(\"=\"*70)\n", "\n", "print(\"\\n🔄 Closing Dask client and cluster...\")\n", "try:\n", " client.close()\n", " cluster.close()\n", " print(\"✅ Cleanup complete\")\n", "except Exception as e:\n", " print(f\"⚠️ Error during cleanup: {e}\")\n", "\n", "print(\"\\n\" + \"=\"*70)\n", "print(\"✅ PIPELINE COMPLETE\")\n", "print(\"=\"*70)" ] } ], "metadata": { "kernelspec": { "display_name": "env_01", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.18" } }, "nbformat": 4, "nbformat_minor": 5 }