basketballcantho/remote-sensing
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

update 01

Browse Source
This commit is contained in:
Victor Phan
2025-12-21 14:34:18 +07:00
parent 671a6f851b
commit 6104856031
29 changed files with 60073 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files
01.train_ODC.ipynb
+947
View File
@@ -0,0 +1,947 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"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<timed exec>: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": [
"<timed exec>:169: RuntimeWarning: divide by zero encountered in divide\n",
"<timed exec>: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 \"<timed exec>\", line 262, in <module>\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": 7,
"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",
"# 🤖 RANDOM FOREST MODEL TRAINING\n",
"print(\"=\"*70)\n",
"print(\"MODEL TRAINING\")\n",
"print(\"=\"*70)\n",
"\n",
"if train is not None and ndvi is not None:\n",
" print(\"\\n[1] Extracting features from NDVI...\")\n",
" try:\n",
" # Extract NDVI values at training point locations\n",
" X = []\n",
" y = []\n",
" \n",
" for idx, point in train.iterrows():\n",
" try:\n",
" # Get NDVI value at point location (nearest neighbor)\n",
" ndvi_val = float(ndvi.sel(x=point.geometry.x, y=point.geometry.y, method='nearest').values)\n",
" label = label_mapping[point.Hientrang]\n",
" \n",
" X.append([ndvi_val])\n",
" y.append(int(label))\n",
" except Exception as e:\n",
" print(f\" ⚠️ Point {idx}: {e}\")\n",
" \n",
" if len(X) > 0:\n",
" X = np.array(X)\n",
" y = np.array(y)\n",
" print(f\" ✅ Extracted {len(X)} samples\")\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\n",
" )\n",
" print(f\" Train: {len(X_train)}, Test: {len(X_test)}\")\n",
" \n",
" # Train model\n",
" print(f\"\\n[3] Training Random Forest...\")\n",
" from sklearn.ensemble import RandomForestClassifier\n",
" from sklearn.metrics import accuracy_score\n",
" \n",
" model = RandomForestClassifier(n_estimators=100, random_state=42, n_jobs=-1)\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",
" print(f\" ✅ Model trained!\")\n",
" print(f\" Accuracy: {accuracy*100:.2f}%\")\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 NDVI\")\n",
" model = None\n",
"\n",
"print(\"=\"*70)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"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\n",
"print(\"=\"*70)\n",
"print(\"MODEL SAVING\")\n",
"print(\"=\"*70)\n",
"\n",
"if model is not None:\n",
" print(\"\\n🔄 Saving trained model...\")\n",
" try:\n",
" save_model(\"model_rasterio.joblib\", model)\n",
" print(\"✅ Model saved to model_train/model_rasterio.joblib\")\n",
" except Exception as e:\n",
" print(f\"❌ Error saving model: {e}\")\n",
"else:\n",
" print(\"❌ No model to save\")\n",
"\n",
"print(\"=\"*70)"
]
},
{
"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 <a href='https://aka.ms/vscodeJupyterKernelCrash'>here</a> for more info. \n",
"\u001b[1;31mView Jupyter <a href='command:jupyter.viewOutput'>log</a> 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
}