Victor Phan
948 lines
39 KiB
Plaintext
948 lines
39 KiB
Plaintext
{
|
||
"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
|
||
}
|