3 Commits

195 changed files with 15101 additions and 179622 deletions
+6
View File
@@ -0,0 +1,6 @@
<component name="InspectionProjectProfileManager">
<settings>
<option name="USE_PROJECT_PROFILE" value="false" />
<version value="1.0" />
</settings>
</component>
+8
View File
@@ -0,0 +1,8 @@
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectModuleManager">
<modules>
<module fileurl="file://$PROJECT_DIR$/.idea/remote-sensing.iml" filepath="$PROJECT_DIR$/.idea/remote-sensing.iml" />
</modules>
</component>
</project>
+8
View File
@@ -0,0 +1,8 @@
<?xml version="1.0" encoding="UTF-8"?>
<module type="PYTHON_MODULE" version="4">
<component name="NewModuleRootManager">
<content url="file://$MODULE_DIR$" />
<orderEntry type="inheritedJdk" />
<orderEntry type="sourceFolder" forTests="false" />
</component>
</module>
Generated
+6
View File
@@ -0,0 +1,6 @@
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="VcsDirectoryMappings">
<mapping directory="" vcs="Git" />
</component>
</project>
+48
View File
@@ -0,0 +1,48 @@
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ChangeListManager">
<list default="true" id="8512e5bb-2b73-4d09-a37b-d0b357b3fbe9" name="Changes" comment="" />
<option name="SHOW_DIALOG" value="false" />
<option name="HIGHLIGHT_CONFLICTS" value="true" />
<option name="HIGHLIGHT_NON_ACTIVE_CHANGELIST" value="false" />
<option name="LAST_RESOLUTION" value="IGNORE" />
</component>
<component name="Git.Settings">
<option name="RECENT_GIT_ROOT_PATH" value="$PROJECT_DIR$" />
</component>
<component name="ProjectColorInfo"><![CDATA[{
"associatedIndex": 1
}]]></component>
<component name="ProjectId" id="39nRCQRaBb6bqtrPoIoBjpMe8Bs" />
<component name="ProjectViewState">
<option name="hideEmptyMiddlePackages" value="true" />
<option name="showLibraryContents" value="true" />
</component>
<component name="PropertiesComponent"><![CDATA[{
"keyToString": {
"ModuleVcsDetector.initialDetectionPerformed": "true",
"RunOnceActivity.ShowReadmeOnStart": "true",
"RunOnceActivity.TerminalTabsStorage.copyFrom.TerminalArrangementManager.252": "true",
"RunOnceActivity.git.unshallow": "true",
"git-widget-placeholder": "dev__01",
"last_opened_file_path": "//wsl.localhost/Ubuntu-22.04/home/x79/remote-sensing"
}
}]]></component>
<component name="SharedIndexes">
<attachedChunks>
<set>
<option value="bundled-python-sdk-4762d8aabb82-6d6dccd035ac-com.jetbrains.pycharm.pro.sharedIndexes.bundled-PY-253.30387.173" />
</set>
</attachedChunks>
</component>
<component name="TaskManager">
<task active="true" id="Default" summary="Default task">
<changelist id="8512e5bb-2b73-4d09-a37b-d0b357b3fbe9" name="Changes" comment="" />
<created>1771329720024</created>
<option name="number" value="Default" />
<option name="presentableId" value="Default" />
<updated>1771329720024</updated>
</task>
<servers />
</component>
</project>
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
-171
View File
@@ -1,171 +0,0 @@
#!/usr/bin/env python
# coding: utf-8
# In[1]:
get_ipython().run_cell_magic('time', '', '%matplotlib inline\n\nimport importlib\nimport new_import_ODC \n\nimportlib.reload(new_import_ODC)\n\nfrom new_import_ODC import *\n')
# In[2]:
get_ipython().run_cell_magic('time', '', '# Cấu hình Daskgateway\ncluster, client = notebook_utils.initialize_dask(use_gateway=True, workers=(1, 10))\n# Khai báo 1 Datacube là dc\ndc = None\n\n# Cấu hình truy cập dịch vụ S3\nconfigure_s3_access(aws_unsigned=False, requester_pays=True, client=client)\n\nclient\n')
# In[3]:
## cấu hình thời gian lấy ảnh và tọa độ
date_range = ("2022-09-01", "2022-10-01")
longtitude_range = (105.86, 105.94)
latitude_range = (9.65, 9.69)
coordinates = (longtitude_range, latitude_range)
# In[4]:
## truy vấn ảnh vệ tinh sen2
data = load_data(None, date_range, longtitude_range, latitude_range)
notebook_utils.heading(notebook_utils.xarray_object_size(data))
display(data)
# In[5]:
get_ipython().run_cell_magic('time', '', '# Tiến hành loại bỏ các vị trí bị mây ảnh hưởng\nresult = mask_clean(data)\n# progress(result)\n')
# In[6]:
# Tiến hành tính toán NDVI
ds1 = calculate_indices(result, index="NDVI", satellite_mission="s2")
ndvi = ds1["NDVI"]
display(ndvi)
# In[7]:
## Hiển thị ảnh NDVI chưa điền các giá trị mây (chưa fill nan)
plt.imshow(ndvi.isel(time=0))
# In[8]:
# Thiết lập giá trị trung bình mùa vụ để xử lý các điểm ảnh bị mây dựa vào sự thay đổi theo mùa
time_split = [
slice("2022-09-01", "2023-01-01"),
slice("2023-01-01", "2023-05-01"),
slice("2023-05-01", "2023-07-01"),
slice("2023-07-01", "2022-10-01"),
]
# Điền mây ở các vị trí mang giá trị nan (fill nan)
fill_nan_ndvi = fill_nan(ndvi, time_split)
# In kết quả ảnh NDVI đã điền mây (đã fill nan)
plt.imshow(fill_nan_ndvi.isel(time=0))
# In[9]:
get_ipython().run_cell_magic('time', '', '## tính ndvi theo tháng\naverage_ndvi = fill_nan_ndvi.resample(time="1M").mean().persist()\n# progress(average_ndvi)\n\n# compute average_ndvi\naverage_ndvi = average_ndvi.compute()\n')
# In[10]:
#Load dữ liệu ảnh Sentinel 1
dsvh, dsvv = load_data_sen1(None, date_range, coordinates)
average_vv = calculate_average(dsvv, time_pattern='1M')
average_vh = calculate_average(dsvh, time_pattern='1M')
# In[11]:
## cấu hình bộ dữ liệu điểm huấn luyện mô hình (train file)
train_path = "train/ST_training_data_updated_1130points_new.shp" # đường dẫn shp file train
## load dữ liệu điểm huấn luyện mô hình (train file)
train = load_train_data(train_path)
train.head()
# cấu hình nhãn dữ liệu
label_mapping = {
"Lua tom": "0",
"Lua": "1",
"CHN": "2",
"CLN": "3",
"TS": "4",
"Song": "5",
"Dat xay dung": "6",
"Rung": "7",
}
# xây dựng tập dữ liệu (dataset) chứa dữ liệu VH, VV, NDVI
datasets = get_data_sen1_and_sen2(train, average_ndvi, average_vh, average_vv)
# chia tập dữ liệu thành các phần theo tỉ lệ 80(80-20)-20 tương ứng với tập train, validate, test
X_train, X_val, X_test, y_train, y_val, y_test = split_train_data(
train, label_mapping, datasets
)
# In[ ]:
get_ipython().run_cell_magic('time', '', '# Import XGBoost\nimport xgboost as xgb\nfrom sklearn.metrics import accuracy_score\nimport numpy as np\n\n# Convert to numpy arrays\nX_train_np = np.asarray(X_train, dtype=np.float32)\nX_val_np = np.asarray(X_val, dtype=np.float32)\ny_train_np = np.asarray(y_train, dtype=np.int32)\ny_val_np = np.asarray(y_val, dtype=np.int32)\n\nprint("🚀 Training XGBoost model...")\nprint(f" Train samples: {len(X_train_np)}")\nprint(f" Val samples: {len(X_val_np)}")\nprint(f" Features: {X_train_np.shape[1]}")\nprint(f" Classes: 8\\n")\n\n# XGBoost parameters\nparams = {\n \'objective\': \'multi:softmax\', # Multi-class classification\n \'num_class\': 8, # 8 land use classes\n \'max_depth\': 6, # Maximum tree depth\n \'learning_rate\': 0.1, # Learning rate\n \'n_estimators\': 200, # Number of trees\n \'subsample\': 0.8, # Subsample ratio\n \'colsample_bytree\': 0.8, # Feature sampling ratio\n \'random_state\': 42,\n \'n_jobs\': -1, # Use all CPU cores\n \'eval_metric\': \'mlogloss\' # Multi-class log loss\n}\n\n# Train XGBoost model\nmodel = xgb.XGBClassifier(**params)\n\nmodel.fit(\n X_train_np, y_train_np,\n eval_set=[(X_train_np, y_train_np), (X_val_np, y_val_np)],\n verbose=True\n)\n\n# Validation accuracy\ny_val_pred = model.predict(X_val_np)\nval_accuracy = accuracy_score(y_val_np, y_val_pred)\nprint(f"\\n✅ Training completed!")\nprint(f" Validation Accuracy: {val_accuracy:.4f} ({val_accuracy*100:.2f}%)")\n')
# In[ ]:
get_ipython().run_cell_magic('time', '', '# Evaluate on test set\nX_test_np = np.asarray(X_test, dtype=np.float32)\ny_test_np = np.asarray(y_test, dtype=np.int32)\n\nprint("📊 Evaluating XGBoost model on test set...\\n")\n\n# Predictions\ny_pred_test = model.predict(X_test_np)\n\n# Metrics\nfrom sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, confusion_matrix\n\ntest_accuracy = accuracy_score(y_test_np, y_pred_test)\nprecision = precision_score(y_test_np, y_pred_test, average=\'weighted\', zero_division=0)\nrecall = recall_score(y_test_np, y_pred_test, average=\'weighted\', zero_division=0)\nf1 = f1_score(y_test_np, y_pred_test, average=\'weighted\', zero_division=0)\n\nprint(f"📈 Test Results:")\nprint(f" Accuracy: {test_accuracy:.4f} ({test_accuracy*100:.2f}%)")\nprint(f" Precision: {precision:.4f}")\nprint(f" Recall: {recall:.4f}")\nprint(f" F1-Score: {f1:.4f}\\n")\n\n# Confusion Matrix\nfrom sklearn.metrics import ConfusionMatrixDisplay\nimport matplotlib.pyplot as plt\n\n# Create figure first\nfig, ax = plt.subplots(figsize=(10, 8))\n\nclass_names = list(label_mapping.keys())\ncm = confusion_matrix(y_test_np, y_pred_test)\ndisp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=class_names)\ndisp.plot(cmap=\'Blues\', ax=ax)\nplt.xticks(rotation=45, ha=\'right\')\nplt.title(\'XGBoost Confusion Matrix\')\nplt.tight_layout()\nplt.show()\n')
# In[ ]:
# Lưu mô hình huấn luyện
import json
import joblib
# Save XGBoost model
model_path = "model_xgboost.joblib"
joblib.dump(model, model_path)
print(f"✅ Model saved to {model_path}")
# Save model info
info = {
"model_type": "XGBoost",
"num_classes": 8,
"classes": list(label_mapping.keys()),
"num_features": X_train_np.shape[1],
"params": params,
"accuracy": float(test_accuracy),
"precision": float(precision),
"recall": float(recall),
"f1_score": float(f1),
}
with open("model_xgboost_info.json", "w") as f:
json.dump(info, f, indent=2)
print(f"✅ Model info saved to model_xgboost_info.json")
# In[15]:
# đóng client, cluster
# client.close()
# cluster.close()
@@ -1,56 +0,0 @@
#!/usr/bin/env python
# coding: utf-8
# In[6]:
get_ipython().run_cell_magic('time', '', '%matplotlib inline\n\n# Import Microsoft Planetary Computer libraries\nimport planetary_computer\nfrom pystac_client import Client\nfrom odc.stac import load as stac_load\n\n# Standard imports\nimport xarray as xr\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import accuracy_score, classification_report, confusion_matrix, ConfusionMatrixDisplay\nimport geopandas as gpd\n\n# XGBoost for GPU training\nimport xgboost as xgb\n\nfrom xgboost import XGBClassifier\n\nprint(f" XGBoost version: {xgb.__version__}")\n\nprint("✅ All modules loaded successfully")\n')
# In[7]:
get_ipython().run_cell_magic('time', '', '# Kết nối tới Microsoft Planetary Computer STAC\nfrom pystac_client import Client\n\n# KHÔNG dùng modifier ở catalog level để tránh items bị convert thành dict\ncatalog = Client.open(\n "https://planetarycomputer.microsoft.com/api/stac/v1"\n)\nprint("✅ Connected to Microsoft Planetary Computer")\n\nprint("\\n" + "="*70)\n')
# In[8]:
get_ipython().run_cell_magic('time', '', '# 🌍 Định nghĩa khu vực và thời gian\nprint("="*70)\nprint("CONFIGURATION")\nprint("="*70)\n\n# Khu vực quan tâm (Vietnam - Mekong Delta) - GIẢM DIỆN TÍCH ~40%\nbbox = [105.6, 9.3, 106.2, 9.8] # [min_lon, min_lat, max_lon, max_lat]\n\n# GIẢM THỜI GIAN xuống 3 tháng để giảm kích thước dữ liệu cho PC\ntime_range = "2023-03-01/2023-05-31" # 3 tháng (mùa khô)\n\nprint(f"\\n📍 Area of Interest:")\nprint(f" Longitude: {bbox[0]} to {bbox[2]}")\nprint(f" Latitude: {bbox[1]} to {bbox[3]}")\nprint(f"\\n📅 Time Range: {time_range}")\nprint(f" ⚠️ Optimized for personal computer (3 months, reduced area)")\nprint(f"\\n🗺️ CRS: EPSG:32648")\nprint(f" Resolution: 20m (reduced from 10m for smaller data size)")\n\nprint("="*70)\n')
# In[9]:
get_ipython().run_cell_magic('time', '', '# 📡 LOAD SENTINEL-2 FROM MICROSOFT PLANETARY COMPUTER\nprint("="*70)\nprint("LOADING SENTINEL-2 L2A")\nprint("="*70)\n\nprint("\\n🔍 Searching for Sentinel-2 scenes...")\nquery_s2 = catalog.search(\n collections=["sentinel-2-l2a"],\n bbox=bbox,\n datetime=time_range,\n query={"eo:cloud_cover": {"lt": 30}} # Cloud cover < 30% (giảm từ 50%)\n)\n\nitems_s2 = list(query_s2.item_collection())\nprint(f"✅ Found {len(items_s2)} Sentinel-2 scenes")\n\n# GIỚI HẠN SỐ LƯỢNG SCENES cho PC cá nhân\nmax_scenes = 12 # Giảm xuống 12 scenes để tối ưu cho PC\nif len(items_s2) > max_scenes:\n print(f"⚠️ Limiting to {max_scenes} scenes for personal computer")\n # Chọn scenes đều đặn trong khoảng thời gian\n step = len(items_s2) // max_scenes\n items_s2 = items_s2[::step][:max_scenes]\n print(f" Selected {len(items_s2)} scenes evenly distributed")\n\nif len(items_s2) > 0:\n # Show first few scenes\n print(f"\\n📋 Sample scenes:")\n for i, item in enumerate(items_s2[:5]):\n date = item.datetime.strftime("%Y-%m-%d")\n cloud = item.properties.get("eo:cloud_cover", "N/A")\n print(f" [{i+1}] {date} - Cloud: {cloud}%")\n \n # Re-sign items to ensure fresh URLs (keep as pystac objects)\n print(f"\\n🔑 Signing STAC items...")\n items_s2 = [planetary_computer.sign(item) for item in items_s2]\n \n # Load Sentinel-2 data (without Dask chunks)\n print(f"\\n⏳ Loading Sentinel-2 data...")\n ds_s2 = stac_load(\n items_s2,\n bands=["B04", "B08", "SCL"], # Red (B04), NIR (B08), Scene Classification (SCL)\n crs="EPSG:32648",\n resolution=20, # 20m resolution (4x smaller data than 10m)\n bbox=bbox,\n patch_url=planetary_computer.sign, # Re-sign URLs during loading\n fail_on_error=False, # Skip problematic tiles instead of crashing\n )\n \n # Rename bands to simpler names\n ds_s2 = ds_s2.rename({"B04": "red", "B08": "nir", "SCL": "scl"})\n \n print(f"\\n✅ Sentinel-2 loaded!")\n print(f" Shape: {dict(ds_s2.dims)}")\n print(f" Variables: {list(ds_s2.data_vars)}")\n display(ds_s2)\nelse:\n print(f"❌ No Sentinel-2 scenes found")\n\n ds_s2 = Noneprint("="*70)\n')
# In[10]:
get_ipython().run_cell_magic('time', '', '# 📡 LOAD SENTINEL-1 FROM MICROSOFT PLANETARY COMPUTER\nprint("="*70)\nprint("LOADING SENTINEL-1 RTC")\nprint("="*70)\n\nprint("\\n🔍 Searching for Sentinel-1 scenes...")\nquery_s1 = catalog.search(\n collections=["sentinel-1-rtc"],\n bbox=bbox,\n datetime=time_range,\n)\n\nitems_s1 = list(query_s1.item_collection())\nprint(f"✅ Found {len(items_s1)} Sentinel-1 scenes")\n\n# GIỚI HẠN SỐ LƯỢNG SCENES cho PC cá nhân\nmax_scenes = 12 # Giảm xuống 12 scenes để tối ưu cho PC\nif len(items_s1) > max_scenes:\n print(f"⚠️ Limiting to {max_scenes} scenes for personal computer")\n # Chọn scenes đều đặn trong khoảng thời gian\n step = len(items_s1) // max_scenes\n items_s1 = items_s1[::step][:max_scenes]\n print(f" Selected {len(items_s1)} scenes evenly distributed")\n\nif len(items_s1) > 0:\n # Show first few scenes\n print(f"\\n📋 Sample scenes:")\n for i, item in enumerate(items_s1[:5]):\n date = item.datetime.strftime("%Y-%m-%d")\n orbit = item.properties.get("sat:orbit_state", "N/A")\n print(f" [{i+1}] {date} - Orbit: {orbit}")\n \n # Re-sign items to ensure fresh URLs (keep as pystac objects)\n print(f"\\n🔑 Signing STAC items...")\n items_s1 = [planetary_computer.sign(item) for item in items_s1]\n \n # Load Sentinel-1 data (without Dask chunks)\n print(f"\\n⏳ Loading Sentinel-1 data...")\n ds_s1 = stac_load(\n items_s1,\n bands=["vv", "vh"], # VV and VH polarizations\n crs="EPSG:32648",\n resolution=20, # 20m resolution (4x smaller data than 10m)\n bbox=bbox,\n patch_url=planetary_computer.sign, # Re-sign URLs during loading\n fail_on_error=False, # Skip problematic tiles instead of crashing\n )\n \n # Convert to dB (Microsoft S1 is in linear power)\n print(f"\\n🔄 Converting to dB...")\n ds_s1[\'vv_db\'] = 10 * np.log10(ds_s1[\'vv\'].where(ds_s1[\'vv\'] > 0))\n ds_s1[\'vh_db\'] = 10 * np.log10(ds_s1[\'vh\'].where(ds_s1[\'vh\'] > 0))\n \n print(f"\\n✅ Sentinel-1 loaded!")\n print(f" Shape: {dict(ds_s1.dims)}")\n print(f" Variables: {list(ds_s1.data_vars)}")\n display(ds_s1)\nelse:\n print(f"❌ No Sentinel-1 scenes found")\n\n ds_s1 = Noneprint("="*70)\n')
# In[11]:
get_ipython().run_cell_magic('time', '', '# 🌿 CALCULATE NDVI AND PROCESS DATA\nprint("="*70)\nprint("DATA PROCESSING")\nprint("="*70)\n\nif ds_s2 is not None:\n print("\\n[1] Calculating NDVI...")\n # NDVI = (NIR - Red) / (NIR + Red)\n ndvi = (ds_s2[\'nir\'] - ds_s2[\'red\']) / (ds_s2[\'nir\'] + ds_s2[\'red\'] + 1e-8)\n \n print(f"✅ NDVI calculated")\n print(f" Shape: {ndvi.shape}")\n print(f" Time steps: {len(ndvi.time)}")\n \n # Cloud masking using SCL band\n print(f"\\n[2] Applying cloud mask...")\n # SCL values: 1=defective, 3=cloud shadow, 8=cloud medium, 9=cloud high, 10=cirrus\n cloud_mask = ds_s2[\'scl\'].isin([1, 3, 8, 9, 10])\n ndvi_masked = ndvi.where(~cloud_mask)\n \n print(f"✅ Cloud mask applied")\n \n # Temporal aggregation (mean over time)\n print(f"\\n[3] Computing mean NDVI across time...")\n ndvi_mean = ndvi_masked.mean(dim=\'time\')\n \n # Data already in memory, no need to compute() again\n print(f"✅ Mean NDVI computed")\n print(f" Shape: {ndvi_mean.shape}")\n \nelse:\n print("❌ No Sentinel-2 data to process")\n ndvi_mean = None\n\nprint("="*70)\n')
# In[ ]:
get_ipython().run_cell_magic('time', '', '# 🎯 EXTRACT TRAINING DATA FEATURES\nprint("="*70)\nprint("FEATURE EXTRACTION")\nprint("="*70)\n\n# Check if required data is available\nif \'ndvi_mean\' not in globals() or \'ds_s1\' not in globals():\n print("❌ Error: Please run Cell 6 (DATA PROCESSING) first!")\n print(" Required variables: ndvi_mean, ds_s1")\n raise RuntimeError("Missing required data. Run cells in order: Cell 4 → Cell 5 → Cell 6 → Cell 7")\n\n# Load training shapefile\nimport geopandas as gpd\n\ntrain_path = \'train/ST_training data_updated_1130points_new.shp\'\nprint(f"\\n[1] Loading training data from: {train_path}")\ntrain_gdf = gpd.read_file(train_path)\n\n# Ensure CRS matches\nif train_gdf.crs != \'EPSG:32648\':\n print(f" Reprojecting from {train_gdf.crs} to EPSG:32648...")\n train_gdf = train_gdf.to_crs(\'EPSG:32648\')\n\nprint(f"✅ Loaded {len(train_gdf)} training points")\nprint(f" Available columns: {list(train_gdf.columns)}")\n\n# Auto-detect label column (look for common names)\nlabel_column = None\nfor col in [\'HT_code\', \'Ma_LU\', \'LU2022\', \'class\', \'Class\', \'CLASS\', \'label\', \'Label\', \'LABEL\', \'LU_CODE\', \'LU_code\']:\n if col in train_gdf.columns:\n label_column = col\n break\n\nif label_column is None:\n print(f"❌ Cannot find label column. Available columns: {list(train_gdf.columns)}")\n print(f" Please check your shapefile and update the code.")\nelse:\n print(f" Using label column: \'{label_column}\'")\n print(f" Classes: {sorted(train_gdf[label_column].unique())}")\n \n # Extract features at each training point\n print(f"\\n[2] Extracting features at training points...")\n \n features = []\n labels = []\n skipped = 0\n \n for idx, row in train_gdf.iterrows():\n point = row.geometrychro\n x_coord = point.x\n y_coord = point.y\n label = row[label_column]\n \n # Extract NDVI at this location\n if ndvi_mean is not None and ds_s1 is not None:\n try:\n ndvi_val = ndvi_mean.sel(x=x_coord, y=y_coord, method=\'nearest\').values\n \n # Extract Sentinel-1 VH/VV at this location (mean across time)\n # Data already in memory, no need to compute()\n vh_val = ds_s1[\'vh_db\'].sel(x=x_coord, y=y_coord, method=\'nearest\').mean(dim=\'time\').values\n vv_val = ds_s1[\'vv_db\'].sel(x=x_coord, y=y_coord, method=\'nearest\').mean(dim=\'time\').values\n \n # Create feature vector: [NDVI, VH_dB, VV_dB]\n feature_vec = [ndvi_val, vh_val, vv_val]\n \n # Only add if all features are valid (not NaN)\n if not np.isnan(feature_vec).any():\n features.append(feature_vec)\n labels.append(label)\n else:\n skipped += 1\n except Exception as e:\n # Skip points outside the data extent\n skipped += 1\n continue\n \n features = np.array(features)\n labels = np.array(labels)\n \n print(f"✅ Extracted features for {len(features)} valid points")\n print(f" Skipped {skipped} points (outside extent or NaN values)")\n print(f" Feature shape: {features.shape}")\n print(f" Feature names: [\'NDVI_mean\', \'VH_dB_mean\', \'VV_dB_mean\']")\n print(f"\\n Class distribution:")\n unique, counts = np.unique(labels, return_counts=True)\n for cls, cnt in zip(unique, counts):\n print(f" Class {cls}: {cnt} samples ({cnt/len(labels)*100:.1f}%)")\n\nprint("="*70)\n')
# In[21]:
get_ipython().run_cell_magic('time', '', '# 🤖 TRAIN XGBOOST MODEL ON GPU (RTX 4060)\nprint("="*70)\nprint("MODEL TRAINING - GPU ACCELERATED")\nprint("="*70)\n\nfrom xgboost import XGBClassifier\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.metrics import classification_report, confusion_matrix, ConfusionMatrixDisplay\nimport matplotlib.pyplot as plt\n\n# Encode labels to ensure they are 0, 1, 2, ... n-1\nprint("\\n[1] Encoding labels...")\nlabel_encoder = LabelEncoder()\nlabels_encoded = label_encoder.fit_transform(labels)\nprint(f"✅ Original classes: {label_encoder.classes_}")\nprint(f" Encoded as: {np.unique(labels_encoded)}")\n\n# Split data\nprint("\\n[2] Splitting data (80% train, 20% test)...")\nX_train, X_test, y_train, y_test = train_test_split(\n features, labels_encoded, test_size=0.2, random_state=42, stratify=labels_encoded\n)\nprint(f"✅ Training samples: {len(X_train)}")\nprint(f" Testing samples: {len(X_test)}")\n\n# Train XGBoost on GPU\nprint("\\n[3] Training XGBoost classifier on RTX 4060 GPU...")\nprint(" GPU Settings: device=\'cuda:0\'")\n\nxgb_model = XGBClassifier(\n n_estimators=100,\n max_depth=20,\n learning_rate=0.1,\n device=\'cuda:0\', # Use GPU (updated from deprecated gpu_id)\n tree_method=\'hist\', # Use hist with device for GPU training\n random_state=42,\n eval_metric=\'mlogloss\', # Multi-class log loss\n verbosity=1 # Show GPU training progress\n)\n\nxgb_model.fit(X_train, y_train)\nprint(f"✅ Model trained on GPU")\n\n# Evaluate\nprint("\\n[4] Evaluating model...")\ntrain_score = xgb_model.score(X_train, y_train)\ntest_score = xgb_model.score(X_test, y_test)\nprint(f"✅ Training accuracy: {train_score:.4f}")\nprint(f" Testing accuracy: {test_score:.4f}")\n\n# Classification report\nprint("\\n[5] Classification Report:")\ny_pred = xgb_model.predict(X_test)\nprint(classification_report(y_test, y_pred, target_names=[str(c) for c in label_encoder.classes_]))\n\n# Confusion matrix\nprint("\\n[6] Confusion Matrix:")\nfig, ax = plt.subplots(figsize=(10, 8))\ncm = confusion_matrix(y_test, y_pred)\ndisp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=label_encoder.classes_)\ndisp.plot(ax=ax, cmap=\'Blues\', values_format=\'d\')\nplt.title(\'Confusion Matrix - XGBoost GPU Model (RTX 4060)\')\nplt.tight_layout()\nplt.show()\n\nprint("="*70)\n')
# In[23]:
get_ipython().run_cell_magic('time', '', '# 💾 SAVE MODEL AND CLEANUP\nprint("="*70)\nprint("SAVING MODEL & CLEANUP")\nprint("="*70)\n\nimport joblib\nfrom datetime import datetime\n\n# Save model and label encoder\nmodel_filename = f"model_train/model_xgboost_gpu_{datetime.now().strftime(\'%Y%m%d_%H%M%S\')}.joblib"\nprint(f"\\n[1] Saving model to: {model_filename}")\njoblib.dump({\'model\': xgb_model, \'label_encoder\': label_encoder}, model_filename)\nprint(f"✅ Model and label encoder saved")\n\n# Save model info\ninfo = {\n "timestamp": datetime.now().isoformat(),\n "data_source": "Microsoft Planetary Computer STAC",\n "collections": ["sentinel-2-l2a", "sentinel-1-rtc"],\n "features": ["NDVI_mean", "VH_dB_mean", "VV_dB_mean"],\n "training_samples": len(X_train),\n "testing_samples": len(X_test),\n "train_accuracy": float(train_score),\n "test_accuracy": float(test_score),\n "model_type": "XGBClassifier",\n "device": "cuda:0",\n "gpu_device": "RTX 4060",\n "tree_method": "hist",\n "n_estimators": 100,\n "max_depth": 20,\n "learning_rate": 0.1\n}\n\nimport json\ninfo_filename = model_filename.replace(\'.joblib\', \'_info.json\')\nwith open(info_filename, \'w\') as f:\n json.dump(info, f, indent=2)\nprint(f"✅ Model info saved to: {info_filename}")\n\n# No cleanup needed (Dask removed)\nprint("\\n[2] Cleanup complete")\n\nprint("="*70)\n\nprint("\\n" + "="*70)\n\nprint("🎉 TRAINING COMPLETE!")\n')
File diff suppressed because one or more lines are too long
-191
View File
@@ -1,191 +0,0 @@
#!/usr/bin/env python
# coding: utf-8
# In[1]:
get_ipython().run_cell_magic('time', '', '%matplotlib inline\n\nimport importlib\nimport new_import_ODC \n\nimportlib.reload(new_import_ODC)\n\nfrom new_import_ODC import *\n')
# In[2]:
get_ipython().run_cell_magic('time', '', '# Dask gateway\ncluster, client = notebook_utils.initialize_dask(use_gateway=True, workers=(1,4))\ndc = datacube.Datacube()\n\n# Configure s3 access\nconfigure_s3_access(aws_unsigned=False, requester_pays=True, client=client)\n\nclient\n')
# In[3]:
## cấu hình thời gian lấy ảnh và tọa độ
date_range = ('2022-09-01', '2023-10-01')
longtitude_range = (105.86575, 105.94120)
latitude_range = (9.65070, 9.69850)
# In[4]:
## truy vấn ảnh vệ tinh sen2
data = load_data(dc, date_range, longtitude_range, latitude_range)
notebook_utils.heading(notebook_utils.xarray_object_size(data))
display(data)
# In[5]:
get_ipython().run_cell_magic('time', '', '# Tiến hành loại bỏ các vị trí bị mây ảnh hưởng\nresult = mask_clean(data)\nprogress(result)\n')
# In[6]:
# Tiến hành tính toán NDVI
ds1 = calculate_indices(result, index='NDVI', satellite_mission='s2')
ndvi = ds1["NDVI"]
display(ndvi)
# In[7]:
## ảnh NDVI chưa điền mây (fill nan)
plt.imshow(ndvi.isel(time=50))
# In[8]:
# đặt thời gian các mùa
time_split = [slice('2022-09-01', '2023-01-01'),
slice('2023-01-01', '2023-05-01'),
slice('2023-05-01', '2023-07-01'),
slice('2023-07-01', '2023-10-01')]
# Điền mây ở các vị trí mang giá trị nan (fill nan)
fill_nan_ndvi = fill_nan(ndvi, time_split)
# In kết quả ảnh ndvi đã điền mây (đã fill nan)
plt.imshow(fill_nan_ndvi.isel(time=50))
# In[9]:
get_ipython().run_cell_magic('time', '', "## tính ndvi theo tháng\naverage_ndvi = fill_nan_ndvi.resample(time='1M').mean().persist()\nprogress(average_ndvi)\n")
# In[10]:
# compute average_ndvi
average_ndvi = average_ndvi.compute()
# In[11]:
# load dữ liệu sen1
coordinates = (longtitude_range, latitude_range)
dsvh, dsvv = load_data_sen1(dc, date_range, coordinates)
average_vv = calculate_average(dsvv, time_pattern='1M')
average_vh = calculate_average(dsvh, time_pattern='1M')
# In[12]:
# load model RF
loaded_model = joblib.load(os.path.join("model_train", "model_odc.joblib"))
# dự đoán
data_array = predict(loaded_model, data.rio.crs, average_ndvi, average_vh, average_vv)
# In[13]:
# cấu hình màu cho các loại đất
colors = [
"#abcee9",
"#ffef44",
"#c4ff9e",
"#ffd6a8",
"#93ddda",
"#1aeef7",
"#ffa7f2",
"#33ee33"
]
labels = [
"Lúa tôm",
"Lúa",
"CHN",
"CLN",
"TS",
"Sông",
"Đất xây dựng",
"Rừng"
]
# hiển thị phân loại sử dụng đất
cmap = ListedColormap(colors)
img = data_array.plot(cmap=cmap, add_colorbar=False)
cbar = plt.colorbar(img)
cbar.ax.set_yticklabels(labels)
plt.title("Phân loại sử dụng đất")
plt.axis('off')
plt.show()
# In[14]:
## cấu hình shapefile ranh giới thuận hòa và vh vv file
thuanhoa_path = "ThuanHoa/region/ST_ThuanHoa_Boundaryofficially.shp"
# cắt theo ranh giới xã thuận hòa
region_result = cut_according_shp(thuanhoa_path, average_ndvi, data_array)
# In[15]:
# hiển thị kết quả phân loại sử dụng đất
colorval = list(range(len(colors)))
options = {
'title': 'Phân loại sử dụng đất',
'cmap': colors,
'clim': (0, 8),
'aspect': 'equal',
'colorbar_opts': {
'major_label_overrides': dict(zip(colorval, labels)),
'major_label_text_align': 'left',
'ticker': FixedTicker(ticks=colorval),
},
}
region_result.hvplot(
rasterize = True, # Use Datashader, particularly useful for dask arrays
aggregator = reductions.mode(), # Datashader selects mode value, requires 'hv.Image'
).options(opts.Image(**options))
# In[16]:
# Lưu lại kết quả
region_result.rio.to_raster("KetQuaPhanLoaiDatODC.tif")
# In[17]:
# đóng client, cluster
client.close()
cluster.close()
# In[ ]:
-117
View File
@@ -1,117 +0,0 @@
#!/usr/bin/env python
# coding: utf-8
# In[1]:
# Khai báo các thư viện cần thiết
from new_import_ODC import *
# Khai báo đường dẫn đến kết quả phân loại và dữ liệu của địa phương
KD_path = "ThuanHoa/KhoanhDat/ThuanHoa_TKDD2022.shp"
KetQuaPhanLoaiDat = "KetQuaPhanLoaiDatODC.tif"
# In[2]:
# khai báo các loại đất từ dữ liệu kiểm kê ứng với các hiện trạng được phân loại từ viễn thám
CODE_MAP = {
"BHK": 2,
"CLN": 3,
"DGD": 6,
"DGT": 6,
"DNL": 6,
"DRA": 6,
"DSH": 6,
"DTL": 5,
"DTS": 6,
"DYT": 6,
"LUC": 1,
"NKH": 3,
"NTD": 6,
"NTS": 4,
"ONT": 6,
"SKC": 6,
"SKX": 6,
"SON": 5,
"TMD": 6,
"TON": 6,
"TSC": 6,
}
# Khai báo các nhãn phân loại đất ứng với 3 loại đất chính
HT_MAP = {
"NN": {"name": "Đất Nông Nghiệp", "data": [1, 2, 3, 4]},
"PNN": {"name": "Đất Phi Nông Nghiệp", "data": [6]},
"TQ": {"name": "Đất Thổ Quả", "data": [15]},
}
# In[3]:
# Tiến hành chồng lắp
result = compare(KD_path, KetQuaPhanLoaiDat, CODE_MAP, HT_MAP)
# In[4]:
# cấu hình màu cho các loại sử dụng đất
colors = [
"#abcee9",
"#ffffc0",
"#c4ff9e",
"#ffd6a8",
"#93ddda",
"#1aeef7",
"#ffa7f2",
"#33ee33",
]
labels = ["Lúa tôm", "Lúa", "CHN", "CLN", "TS", "Sông", "Đất xây dựng", "Rừng"]
# In[5]:
# Lưu kết quả
save_result(result, HT_MAP)
# In[6]:
# hiển thị kết quả
xx = []
for k, v in result.items():
rs = merge_arrays(v, nodata=np.nan)
xx.append(rs.squeeze(drop=True))
xx = xr.concat(xx, pd.Index([HT_MAP[x]["name"] for x in HT_MAP], name="name"))
colorval = list(range(len(colors)))
options = {
"cmap": colors,
"clim": (0, 8),
"aspect": "equal",
"height": 400,
"colorbar_opts": {
"major_label_overrides": dict(zip(colorval, labels)),
"major_label_text_align": "left",
"ticker": FixedTicker(ticks=colorval),
},
}
xx.hvplot(
groupby="name",
rasterize=True, # Use Datashader, particularly useful for dask arrays
aggregator=reductions.mode(), # Datashader selects mode value, requires 'hv.Image'
).options(opts.Image(**options))
# In[ ]:
File diff suppressed because one or more lines are too long
-18
View File
@@ -1,18 +0,0 @@
import joblib
import numpy as np
cache_file = "dataset_cache/training_data_2d.joblib"
data = joblib.load(cache_file)
X = np.array(data['X'])
y = np.array(data['y'])
print("X shape:", X.shape)
print("X mean:", np.mean(X))
print("X std:", np.std(X))
print("X min:", np.min(X))
print("X max:", np.max(X))
print("Any NaN:", np.isnan(X).any())
for i in range(6):
print(f"Channel {i} mean: {np.mean(X[:, i, :, :]):.4f}, min: {np.min(X[:, i, :, :]):.4f}, max: {np.max(X[:, i, :, :]):.4f}")
-10
View File
@@ -1,10 +0,0 @@
import joblib
import geopandas as gpd
from shapely.geometry import Point
data = joblib.load('dataset_cache/training_data_2d.joblib')
X, y = data['X'], data['y']
print(f"X shape: {X.shape}, y shape: {y.shape}")
gdf = gpd.read_file("train/ST_training_data_updated_1130points_new.shp")
print("Total points:", len(gdf))
-7
View File
@@ -1,7 +0,0 @@
import joblib
import numpy as np
data = joblib.load('dataset_cache/training_data.joblib')
X, y = data['X'], data['y']
print(f"X shape: {X.shape}")
print(f"y shape: {y.shape}")
-8
View File
@@ -1,8 +0,0 @@
import joblib
import numpy as np
cache_file = "dataset_cache/training_data_2d.joblib"
data = joblib.load(cache_file)
X = np.array(data['X'])
b2 = X[:, 0, :, :]
print("Zeros in B2:", np.sum(b2 == 0) / b2.size)
print("X shape:", X.shape)
-25
View File
@@ -1,25 +0,0 @@
import glob, json
changed_files = []
for file_path in glob.glob('*.ipynb'):
with open(file_path, 'r', encoding='utf-8') as f:
nb = json.load(f)
changed = False
for cell in nb.get('cells', []):
if cell.get('cell_type') == 'code':
source = cell.get('source', [])
for i, line in enumerate(source):
if 'time=50' in line:
source[i] = line.replace('time=50', 'time=0')
changed = True
if 'load_data_sen1(dc,' in line:
source[i] = line.replace('load_data_sen1(dc,', 'load_data_sen1(None,')
changed = True
if changed:
with open(file_path, 'w', encoding='utf-8') as f:
json.dump(nb, f, indent=1)
changed_files.append(file_path)
print('Fixed issues in:', changed_files)
-19
View File
@@ -1,19 +0,0 @@
import joblib
import geopandas as gpd
import numpy as np
cache_file = "dataset_cache/training_data_2d.joblib"
data = joblib.load(cache_file)
X = data['X']
print("X shape:", len(X))
gdf = gpd.read_file("train/ST_training_data_updated_1130points_new.shp")
gdf = gdf.to_crs("EPSG:32648")
print("gdf length:", len(gdf))
if len(X) == len(gdf):
y = [(row['HT_code'] - 1) for idx, row in gdf.iterrows()]
joblib.dump({'X': X, 'y': y}, cache_file)
print("Fixed y in cache! Saved.")
else:
print("Lengths do not match, cannot fix automatically.")
-20
View File
@@ -1,20 +0,0 @@
import json
def fix_import(file_path):
with open(file_path, 'r', encoding='utf-8') as f:
nb = json.load(f)
changed = False
for cell in nb.get('cells', []):
if cell.get('cell_type') == 'code':
source = cell.get('source', [])
if isinstance(source, list):
for i, line in enumerate(source):
if "from new_import import *" in line:
source[i] = line.replace("from new_import import *", "from new_import_ODC import *")
changed = True
if changed:
with open(file_path, 'w', encoding='utf-8') as f:
json.dump(nb, f, indent=1)
print(f"Fixed {file_path}")
fix_import('new_train.ipynb')
-20
View File
@@ -1,20 +0,0 @@
import re
filepath = "01.train_ODC.py"
with open(filepath, 'r', encoding='utf-8') as f:
content = f.read()
# Find the run_cell_magic line
pattern = re.compile(r"get_ipython\(\)\.run_cell_magic\('time', '', '(# 🤖 LAND USE CLASSIFICATION MODEL TRAINING.*?)(?=\n')\n'", re.DOTALL)
def repl(match):
# Get the inner string and escape all actual newlines with \n
inner = match.group(1)
inner = inner.replace('\n', '\\n')
return f"get_ipython().run_cell_magic('time', '', '{inner}')"
content = pattern.sub(repl, content)
with open(filepath, 'w', encoding='utf-8') as f:
f.write(content)
print("Fixed 01.train_ODC.py syntax")
-22
View File
@@ -1,22 +0,0 @@
import json
def fix_filename(file_path):
with open(file_path, 'r', encoding='utf-8') as f:
nb = json.load(f)
changed = False
for cell in nb.get('cells', []):
if cell.get('cell_type') == 'code':
source = cell.get('source', [])
if isinstance(source, list):
for i, line in enumerate(source):
if "ST_training data_updated_1130points.shp" in line:
source[i] = line.replace("ST_training data_updated_1130points.shp", "ST_training_data_updated_1130points.shp")
changed = True
if changed:
with open(file_path, 'w', encoding='utf-8') as f:
json.dump(nb, f, indent=1)
print(f"Fixed typo in {file_path}")
import glob
for nb in glob.glob("*.ipynb"):
fix_filename(nb)
-58
View File
@@ -1,58 +0,0 @@
import nbformat as nbf
nb = nbf.v4.new_notebook()
text_1 = """# Script Tải Dữ liệu Vệ tinh (Cache) qua Google Colab
Mục đích của Notebook này là mượn sức mạnh đường truyền và RAM của Google Colab để tải 270 ảnh Sentinel-2 & Sentinel-1 từ Microsoft Planetary Computer. Sau khi xử lý nội suy, nó sẽ sinh ra một file cache `.joblib` duy nhất chứa toàn bộ mảng dữ liệu.
Bạn chỉ cần tải file `.joblib` đó về máy là xong!"""
code_1 = """!pip install planetary-computer pystac-client odc-stac geopandas rasterio xarray joblib scikit-learn xgboost lightgbm"""
code_2 = """from google.colab import drive
drive.mount('/content/drive')"""
text_2 = """## Hướng dẫn:
1. Nén toàn bộ thư mục `remote-sensing` ở máy tính của bạn thành file `remote-sensing.zip`.
2. Upload file `remote-sensing.zip` đó lên Google Drive (để ngay ngoài cùng).
3. Chạy ô lệnh bên dưới để giải nén và chuyển vào thư mục dự án."""
code_3 = """import os
import shutil
# Giải nén dự án từ Google Drive
!unzip -q /content/drive/MyDrive/remote-sensing.zip -d /content/
os.chdir('/content/remote-sensing')
!ls -la"""
text_3 = """## Bắt đầu tải và Cache Dữ Liệu
Chạy một mô hình CPU đơn giản (Decision Tree) để ép hệ thống gọi hàm `FeatureExtractor`. Hàm này sẽ làm mọi việc nặng nhọc: tìm ảnh, ghép mây, tính trung vị và lưu kết quả vào thư mục `dataset_cache/`."""
code_4 = """# Lệnh này sẽ mất khoảng 5-15 phút để tải toàn bộ ảnh từ Microsoft
!python train_land_decision_tree_gpu.py"""
text_4 = """## Hoàn tất
Bạn hãy kiểm tra xem file `.joblib` lớn (khoảng 40-60MB) đã xuất hiện chưa. Nếu rồi, hãy lưu ngược nó lại Google Drive để tải về máy!"""
code_5 = """# Xem file cache đã được tạo thành công chưa
!ls -lh dataset_cache/
# Copy toàn bộ thư mục cache sang Google Drive để tải về máy dễ dàng
!cp -r dataset_cache/ /content/drive/MyDrive/dataset_cache_finished/
print("Hoàn thành! Bạn hãy mở Google Drive của mình, tìm thư mục 'dataset_cache_finished' và tải file .joblib mới nhất về máy tính.")"""
nb['cells'] = [
nbf.v4.new_markdown_cell(text_1),
nbf.v4.new_code_cell(code_1),
nbf.v4.new_code_cell(code_2),
nbf.v4.new_markdown_cell(text_2),
nbf.v4.new_code_cell(code_3),
nbf.v4.new_markdown_cell(text_3),
nbf.v4.new_code_cell(code_4),
nbf.v4.new_markdown_cell(text_4),
nbf.v4.new_code_cell(code_5)
]
with open('Download_Cache_Colab.ipynb', 'w') as f:
nbf.write(nb, f)
print("Created Download_Cache_Colab.ipynb")
-97
View File
@@ -1,97 +0,0 @@
import os
import glob
import json
from tabulate import tabulate
print("📊 BẢNG SO SÁNH KẾT QUẢ CÁC MÔ HÌNH\n")
# 1. Phân loại đất
print("### 1. Nhóm Phân loại Lớp phủ (Land Classification)")
land_data = []
if os.path.exists("model_xgboost_info.json"):
with open("model_xgboost_info.json", 'r') as f:
data = json.load(f)
params = data.get('params', {})
param_str = f"estimators:{params.get('n_estimators')}, depth:{params.get('max_depth')}" if params else "N/A"
land_data.append([
data.get('model_type', 'XGBoost'),
data.get('accuracy', ''),
data.get('precision', ''),
data.get('recall', ''),
data.get('f1_score', ''),
param_str
])
for info_file in glob.glob("model_train/*_info.json"):
with open(info_file, 'r') as f:
data = json.load(f)
# Support both 'accuracy' and 'test_accuracy'
acc = data.get('accuracy', data.get('test_accuracy', ''))
f1 = data.get('f1_score', '')
precision = data.get('precision', '')
recall = data.get('recall', '')
clf_rep = data.get('classification_report')
if isinstance(clf_rep, dict) and 'macro avg' in clf_rep:
if not f1:
f1 = clf_rep['macro avg'].get('f1-score', '')
if not precision:
precision = clf_rep['macro avg'].get('precision', '')
if not recall:
recall = clf_rep['macro avg'].get('recall', '')
if not acc and not f1:
continue
params = data.get('params', {})
param_str = f"estimators:{params.get('n_estimators')}, depth:{params.get('max_depth')}" if params else "N/A"
if data.get('model_type') == 'RandomForest_RealData':
param_str = "estimators:100, depth:15"
land_data.append([
data.get('model_type', ''),
acc,
precision,
recall,
f1,
param_str
])
if land_data:
print(tabulate(land_data, headers=["Model", "Accuracy", "Precision", "Recall", "F1-Score", "Parameters"], tablefmt="github"))
print("\n")
# 2. Xóa mây
print("### 2. Nhóm Xóa mây (Cloud Removal)")
cloud_data = []
for info_file in glob.glob("cloud_removal_model/*_info.json"):
with open(info_file, 'r') as f:
data = json.load(f)
cloud_data.append([
data.get('model_type', ''),
data.get('epoch', ''),
data.get('train_loss', ''),
data.get('val_loss', '')
])
if cloud_data:
print(tabulate(cloud_data, headers=["Model", "Epochs", "Train Loss", "Val Loss"], tablefmt="github"))
print("\n")
# 3. Dự báo NDVI
print("### 3. Nhóm Dự báo Thực vật (NDVI Forecasting)")
ndvi_data = []
for info_file in glob.glob("ndvi_forecast_model/*_info.json"):
with open(info_file, 'r') as f:
data = json.load(f)
ndvi_data.append([
data.get('model_type', ''),
data.get('rmse', ''),
data.get('mae', ''),
data.get('epoch', 'N/A')
])
if ndvi_data:
print(tabulate(ndvi_data, headers=["Model", "RMSE", "MAE", "Epochs"], tablefmt="github"))
print("\n")
-7
View File
@@ -1,7 +0,0 @@
import json
nb = json.load(open('01.train_ODC.ipynb'))
for idx, cell in enumerate(nb['cells']):
if cell['cell_type'] == 'code':
print(f"Cell {idx}:")
print("".join(cell['source'][:3]))
print("-" * 20)
-18
View File
@@ -1,18 +0,0 @@
🚀 BẮT ĐẦU HUẤN LUYỆN MÔ HÌNH CNN (GPU & CACHE)
Initializing FeatureExtractor (mode=extended)...
📦 Đang load cache: training_data_507bd2ba4ec0d3fe107839cbf73a7a7d.joblib...
✅ Loaded 632 samples từ cache!
⚡ Đã bỏ qua download위성 data (tiết kiệm thời gian)
[CACHE HIT] Using cached dataset with 632 samples
Training CNN model...
Building CNN model on cuda...
Training CNN model with PyTorch...
CNN Epoch 10/15, Loss: 1.2171
Evaluating model...
Generating classification report...
Saving model...
[MODEL MANAGER] Saving model to: model_train/model_cnn_auto.joblib
[MODEL MANAGER] Saving metadata to: model_train/model_cnn_auto_info.json
[MODEL MANAGER] Model saved successfully!
Training complete!
✅ Hoàn thành! Accuracy: 0.5748
-15
View File
@@ -1,15 +0,0 @@
🚀 BẮT ĐẦU HUẤN LUYỆN MÔ HÌNH DECISION TREE (GPU & CACHE)
Initializing FeatureExtractor (mode=extended)...
📦 Đang load cache: training_data_507bd2ba4ec0d3fe107839cbf73a7a7d.joblib...
✅ Loaded 632 samples từ cache!
⚡ Đã bỏ qua download위성 data (tiết kiệm thời gian)
[CACHE HIT] Using cached dataset with 632 samples
Training DECISION_TREE model...
Evaluating model...
Generating classification report...
Saving model...
[MODEL MANAGER] Saving model to: model_train/model_decision_tree_auto.joblib
[MODEL MANAGER] Saving metadata to: model_train/model_decision_tree_auto_info.json
[MODEL MANAGER] Model saved successfully!
Training complete!
✅ Hoàn thành! Accuracy: 0.5906
File diff suppressed because it is too large Load Diff
File diff suppressed because it is too large Load Diff
File diff suppressed because it is too large Load Diff
-30
View File
@@ -1,30 +0,0 @@
🚀 BẮT ĐẦU HUẤN LUYỆN MÔ HÌNH MOBILENET-LRASPP (GPU & CACHE)
Initializing FeatureExtractor (mode=extended)...
📦 Đang load cache: training_data_507bd2ba4ec0d3fe107839cbf73a7a7d.joblib...
✅ Loaded 632 samples từ cache!
⚡ Đã bỏ qua download위성 data (tiết kiệm thời gian)
[CACHE HIT] Using cached dataset with 632 samples
Training MOBILENET-LRASPP model...
Building MobileNetV3 + LR-ASPP model on cuda...
[MOBILENET] Class distribution: [ 48 89 3 86 74 38 117 50]
[MOBILENET] Class weights: [0.37418982 0.20181024 5.98703525 0.20885013 0.24271772 0.47266082
0.15351377 0.35922223]
Training MobileNetV3 + LR-ASPP model with PyTorch...
MobileNet Epoch 5/25, Train Loss: 1.0614, Val Loss: 1.3584, Val Acc: 40.16%, LR: 0.000800
[MOBILENET] Epoch 5/25 - Train Loss: 1.0614, Val Loss: 1.3584, Val Acc: 40.16%
MobileNet Epoch 10/25, Train Loss: 0.9054, Val Loss: 0.8582, Val Acc: 59.06%, LR: 0.000800
[MOBILENET] Epoch 10/25 - Train Loss: 0.9054, Val Loss: 0.8582, Val Acc: 59.06%
MobileNet Epoch 15/25, Train Loss: 0.7728, Val Loss: 0.8231, Val Acc: 61.42%, LR: 0.000800
[MOBILENET] Epoch 15/25 - Train Loss: 0.7728, Val Loss: 0.8231, Val Acc: 61.42%
MobileNet Epoch 20/25, Train Loss: 0.7125, Val Loss: 0.8783, Val Acc: 59.84%, LR: 0.000400
[MOBILENET] Epoch 20/25 - Train Loss: 0.7125, Val Loss: 0.8783, Val Acc: 59.84%
[MOBILENET] Early stopping at epoch 24 (best val loss: 0.8029)
MobileNet early stopped at epoch 24
Evaluating model...
Generating classification report...
Saving model...
[MODEL MANAGER] Saving model to: model_train/model_mobilenet-lraspp_auto.joblib
[MODEL MANAGER] Saving metadata to: model_train/model_mobilenet-lraspp_auto_info.json
[MODEL MANAGER] Model saved successfully!
Training complete!
✅ Hoàn thành! Accuracy: 0.5669
-15
View File
@@ -1,15 +0,0 @@
🚀 BẮT ĐẦU HUẤN LUYỆN MÔ HÌNH RANDOM FOREST (GPU & CACHE)
Initializing FeatureExtractor (mode=extended)...
📦 Đang load cache: training_data_507bd2ba4ec0d3fe107839cbf73a7a7d.joblib...
✅ Loaded 632 samples từ cache!
⚡ Đã bỏ qua download위성 data (tiết kiệm thời gian)
[CACHE HIT] Using cached dataset with 632 samples
Training RANDOM_FOREST model...
Evaluating model...
Generating classification report...
Saving model...
[MODEL MANAGER] Saving model to: model_train/model_random_forest_auto.joblib
[MODEL MANAGER] Saving metadata to: model_train/model_random_forest_auto_info.json
[MODEL MANAGER] Model saved successfully!
Training complete!
✅ Hoàn thành! Accuracy: 0.6142
-15
View File
@@ -1,15 +0,0 @@
🚀 BẮT ĐẦU HUẤN LUYỆN MÔ HÌNH SVM (GPU & CACHE)
Initializing FeatureExtractor (mode=extended)...
📦 Đang load cache: training_data_507bd2ba4ec0d3fe107839cbf73a7a7d.joblib...
✅ Loaded 632 samples từ cache!
⚡ Đã bỏ qua download위성 data (tiết kiệm thời gian)
[CACHE HIT] Using cached dataset with 632 samples
Training SVM model...
Evaluating model...
Generating classification report...
Saving model...
[MODEL MANAGER] Saving model to: model_train/model_svm_auto.joblib
[MODEL MANAGER] Saving metadata to: model_train/model_svm_auto_info.json
[MODEL MANAGER] Model saved successfully!
Training complete!
✅ Hoàn thành! Accuracy: 0.5827
-30
View File
@@ -1,30 +0,0 @@
🚀 BẮT ĐẦU HUẤN LUYỆN MÔ HÌNH SWIN-UNET (GPU & CACHE)
Initializing FeatureExtractor (mode=extended)...
📦 Đang load cache: training_data_507bd2ba4ec0d3fe107839cbf73a7a7d.joblib...
✅ Loaded 632 samples từ cache!
⚡ Đã bỏ qua download위성 data (tiết kiệm thời gian)
[CACHE HIT] Using cached dataset with 632 samples
Training SWIN-UNET model...
Building Swin-UNet model on cuda...
[SWIN-UNET] Class distribution: [ 48 89 3 86 74 38 117 50]
[SWIN-UNET] Class weights: [0.37418982 0.20181024 5.98703525 0.20885013 0.24271772 0.47266082
0.15351377 0.35922223]
Training Swin-UNet model with PyTorch (with class weights)...
Swin-UNet Epoch 5/40, Train Loss: 1.4096, Val Loss: 1.2804, Val Acc: 48.03%, LR: 0.000293
[SWIN-UNET] Epoch 5/40 - Train Loss: 1.4096, Val Loss: 1.2804, Val Acc: 48.03%
Swin-UNet Epoch 10/40, Train Loss: 1.1129, Val Loss: 1.2746, Val Acc: 57.48%, LR: 0.000271
[SWIN-UNET] Epoch 10/40 - Train Loss: 1.1129, Val Loss: 1.2746, Val Acc: 57.48%
Swin-UNet Epoch 15/40, Train Loss: 1.0479, Val Loss: 1.3126, Val Acc: 50.39%, LR: 0.000238
[SWIN-UNET] Epoch 15/40 - Train Loss: 1.0479, Val Loss: 1.3126, Val Acc: 50.39%
Swin-UNet Epoch 20/40, Train Loss: 0.9548, Val Loss: 1.1118, Val Acc: 52.76%, LR: 0.000196
[SWIN-UNET] Epoch 20/40 - Train Loss: 0.9548, Val Loss: 1.1118, Val Acc: 52.76%
[SWIN-UNET] Early stopping at epoch 22 (best val loss: 1.1096)
Swin-UNet early stopped at epoch 22
Evaluating model...
Generating classification report...
Saving model...
[MODEL MANAGER] Saving model to: model_train/model_swin-unet_auto.joblib
[MODEL MANAGER] Saving metadata to: model_train/model_swin-unet_auto_info.json
[MODEL MANAGER] Model saved successfully!
Training complete!
✅ Hoàn thành! Accuracy: 0.5354
-94521
View File
File diff suppressed because it is too large Load Diff
-38
View File
@@ -1,38 +0,0 @@
🚀 BẮT ĐẦU PIPELINE 2D PATCH-BASED & CLOUD REMOVAL
Loading 2D patches from dataset_cache/training_data_2d.joblib...
Training 2D CNN with Data Augmentation...
Epoch 1/150 - Loss: 2.2272 - Test Acc: 0.0752 🌟
Epoch 2/150 - Loss: 2.0843 - Test Acc: 0.0796 🌟
Epoch 4/150 - Loss: 2.0350 - Test Acc: 0.1372 🌟
Epoch 8/150 - Loss: 2.0447 - Test Acc: 0.1637 🌟
Epoch 10/150 - Loss: 2.0397 - Test Acc: 0.1372
Epoch 12/150 - Loss: 2.0353 - Test Acc: 0.2168 🌟
Epoch 20/150 - Loss: 2.0139 - Test Acc: 0.1372
Traceback (most recent call last):
File "/home/x79/remote-sensing/train_land_2d_patch.py", line 298, in <module>
main()
File "/home/x79/remote-sensing/train_land_2d_patch.py", line 294, in main
train_2d_model(X, y)
File "/home/x79/remote-sensing/train_land_2d_patch.py", line 219, in train_2d_model
for batch_X, batch_y in train_loader:
File "/home/x79/miniconda3/envs/env_01/lib/python3.10/site-packages/torch/utils/data/dataloader.py", line 725, in __next__
data = self._next_data()
File "/home/x79/miniconda3/envs/env_01/lib/python3.10/site-packages/torch/utils/data/dataloader.py", line 785, in _next_data
data = self._dataset_fetcher.fetch(index) # may raise StopIteration
File "/home/x79/miniconda3/envs/env_01/lib/python3.10/site-packages/torch/utils/data/_utils/fetch.py", line 54, in fetch
data = [self.dataset[idx] for idx in possibly_batched_index]
File "/home/x79/miniconda3/envs/env_01/lib/python3.10/site-packages/torch/utils/data/_utils/fetch.py", line 54, in <listcomp>
data = [self.dataset[idx] for idx in possibly_batched_index]
File "/home/x79/remote-sensing/train_land_2d_patch.py", line 192, in __getitem__
x = transform(x)
File "/home/x79/miniconda3/envs/env_01/lib/python3.10/site-packages/torchvision/transforms/transforms.py", line 95, in __call__
img = t(img)
File "/home/x79/miniconda3/envs/env_01/lib/python3.10/site-packages/torch/nn/modules/module.py", line 1778, in _wrapped_call_impl
return self._call_impl(*args, **kwargs)
File "/home/x79/miniconda3/envs/env_01/lib/python3.10/site-packages/torch/nn/modules/module.py", line 1789, in _call_impl
return forward_call(*args, **kwargs)
File "/home/x79/miniconda3/envs/env_01/lib/python3.10/site-packages/torchvision/transforms/transforms.py", line 752, in forward
return F.vflip(img)
File "/home/x79/miniconda3/envs/env_01/lib/python3.10/site-packages/torchvision/transforms/functional.py", line 757, in vflip
def vflip(img: Tensor) -> Tensor:
KeyboardInterrupt
File diff suppressed because it is too large Load Diff
File diff suppressed because one or more lines are too long
-39
View File
@@ -1,39 +0,0 @@
🚀 V4: TÍCH HỢP RADAR SENTINEL-1 (32-CHANNELS FUSION)
============================================================
Clean FUSION data: (252, 32, 16, 16), 7 classes, [31, 38, 32, 49, 23, 75, 4]
============================================================
32-CHANNELS FUSION CNN
============================================================
Ep 1 Fusion-Acc=0.1961 🌟
Ep 3 Fusion-Acc=0.2745 🌟
Ep 4 Fusion-Acc=0.4706 🌟
Ep 5 Fusion-Acc=0.5490 🌟
Ep 6 Fusion-Acc=0.7059 🌟
Ep 7 Fusion-Acc=0.7451 🌟
Ep 8 Fusion-Acc=0.8431 🌟
Ep 15 Fusion-Acc=0.8627 🌟
✅ CNN Fusion best: 0.8627
============================================================
HYBRID FUSION: CNN embed + S1/S2 Rich features + XGBoost
============================================================
Extracted 2182 fusion features per sample
Final Feature Vector: (252, 2694)
✅ Hybrid Fusion Acc: 0.8627
Fold 1: 0.9412
Fold 2: 0.9020
Fold 3: 0.9200
Fold 4: 0.8800
Fold 5: 0.9000
✅ CV Mean: 0.9086 ± 0.0206
============================================================
📊 FINAL RESULTS V4 (WITH RADAR)
============================================================
✅ Hybrid Fusion CV: 0.9086
📈 CNN Fusion (32ch): 0.8627
📈 Hybrid Fusion (CNN+XGB): 0.8627
🏆 BEST: 0.9086
-18
View File
@@ -1,18 +0,0 @@
============================================================
HYBRID FUSION ENSEMBLE: CNN embed + S1/S2 Rich features + XGB/LGBM/ETC
============================================================
Final Feature Vector: (443, 2694)
Fold 1: 0.8876
Fold 2: 0.9438
Fold 3: 0.9438
Fold 4: 0.9659
Fold 5: 0.9432
✅ Ensemble CV Mean: 0.9369 ± 0.0261
============================================================
📊 FINAL RESULTS V5 (ENSEMBLE + RADAR)
============================================================
✅ Hybrid Fusion Ensemble CV: 0.9369
🏆 BEST: 0.9369
-21
View File
@@ -1,21 +0,0 @@
🚀 V6: EXHAUSTIVE HYPERPARAMETER TUNING
============================================================
Data: (443, 32, 16, 16), 7 classes, dist=[65, 55, 48, 72, 74, 124, 5]
--- Training Multi-Seed CNN Ensemble ---
Seed 42: CNN Acc = 0.8989
Seed 123: CNN Acc = 0.8652
Seed 777: CNN Acc = 0.8876
Multi-seed CNN embedding: (443, 1536)
Total features: (443, 3940)
============================================================
🔬 EXHAUSTIVE HYPERPARAMETER SEARCH
============================================================
🏆 XGB-deep: 0.9526 ± 0.0110 (folds: ['0.955', '0.933', '0.966', '0.955', '0.955'])
✅ XGB-shallow: 0.9436 ± 0.0173 (folds: ['0.944', '0.910', '0.955', '0.955', '0.955'])
🏆 XGB-balanced: 0.9504 ± 0.0113 (folds: ['0.944', '0.933', '0.955', '0.955', '0.966'])
✅ LGBM-tuned: 0.9458 ± 0.0149 (folds: ['0.955', '0.921', '0.966', '0.943', '0.943'])
✅ LGBM-conservative: 0.9481 ± 0.0152 (folds: ['0.944', '0.921', '0.966', '0.955', '0.955'])
🏆 ETC-deep: 0.9572 ± 0.0082 (folds: ['0.944', '0.955', '0.955', '0.966', '0.966'])
🏆 RF-tuned: 0.9549 ± 0.0099 (folds: ['0.944', '0.955', '0.944', '0.966', '0.966'])
-157
View File
@@ -1,157 +0,0 @@
🚀 BẮT ĐẦU TÌM KIẾM SIÊU THAM SỐ CHO SWIN-UNET
Loading data from dataset_cache/training_data_507bd2ba4ec0d3fe107839cbf73a7a7d.joblib...
Using device: cuda
[1/48] Training with params: {'embed_dim': 64, 'lr': 0.001, 'weight_decay': 0.01, 'epochs': 200}
Test Accuracy: 0.5984
🌟 NEW BEST ACCURACY: 0.5984
[2/48] Training with params: {'embed_dim': 64, 'lr': 0.001, 'weight_decay': 0.01, 'epochs': 500}
Test Accuracy: 0.6063
🌟 NEW BEST ACCURACY: 0.6063
[3/48] Training with params: {'embed_dim': 64, 'lr': 0.001, 'weight_decay': 0.001, 'epochs': 200}
Test Accuracy: 0.6142
🌟 NEW BEST ACCURACY: 0.6142
[4/48] Training with params: {'embed_dim': 64, 'lr': 0.001, 'weight_decay': 0.001, 'epochs': 500}
Test Accuracy: 0.6772
🌟 NEW BEST ACCURACY: 0.6772
[5/48] Training with params: {'embed_dim': 64, 'lr': 0.0005, 'weight_decay': 0.01, 'epochs': 200}
Test Accuracy: 0.5827
[6/48] Training with params: {'embed_dim': 64, 'lr': 0.0005, 'weight_decay': 0.01, 'epochs': 500}
Test Accuracy: 0.5984
[7/48] Training with params: {'embed_dim': 64, 'lr': 0.0005, 'weight_decay': 0.001, 'epochs': 200}
Test Accuracy: 0.5669
[8/48] Training with params: {'embed_dim': 64, 'lr': 0.0005, 'weight_decay': 0.001, 'epochs': 500}
Test Accuracy: 0.6142
[9/48] Training with params: {'embed_dim': 64, 'lr': 0.0001, 'weight_decay': 0.01, 'epochs': 200}
Test Accuracy: 0.6142
[10/48] Training with params: {'embed_dim': 64, 'lr': 0.0001, 'weight_decay': 0.01, 'epochs': 500}
Test Accuracy: 0.5118
[11/48] Training with params: {'embed_dim': 64, 'lr': 0.0001, 'weight_decay': 0.001, 'epochs': 200}
Test Accuracy: 0.5591
[12/48] Training with params: {'embed_dim': 64, 'lr': 0.0001, 'weight_decay': 0.001, 'epochs': 500}
Test Accuracy: 0.6063
[13/48] Training with params: {'embed_dim': 128, 'lr': 0.001, 'weight_decay': 0.01, 'epochs': 200}
Test Accuracy: 0.7008
🌟 NEW BEST ACCURACY: 0.7008
[14/48] Training with params: {'embed_dim': 128, 'lr': 0.001, 'weight_decay': 0.01, 'epochs': 500}
Test Accuracy: 0.6142
[15/48] Training with params: {'embed_dim': 128, 'lr': 0.001, 'weight_decay': 0.001, 'epochs': 200}
Test Accuracy: 0.6772
[16/48] Training with params: {'embed_dim': 128, 'lr': 0.001, 'weight_decay': 0.001, 'epochs': 500}
Test Accuracy: 0.6063
[17/48] Training with params: {'embed_dim': 128, 'lr': 0.0005, 'weight_decay': 0.01, 'epochs': 200}
Test Accuracy: 0.5906
[18/48] Training with params: {'embed_dim': 128, 'lr': 0.0005, 'weight_decay': 0.01, 'epochs': 500}
Test Accuracy: 0.6457
[19/48] Training with params: {'embed_dim': 128, 'lr': 0.0005, 'weight_decay': 0.001, 'epochs': 200}
Test Accuracy: 0.5906
[20/48] Training with params: {'embed_dim': 128, 'lr': 0.0005, 'weight_decay': 0.001, 'epochs': 500}
Test Accuracy: 0.5906
[21/48] Training with params: {'embed_dim': 128, 'lr': 0.0001, 'weight_decay': 0.01, 'epochs': 200}
Test Accuracy: 0.6220
[22/48] Training with params: {'embed_dim': 128, 'lr': 0.0001, 'weight_decay': 0.01, 'epochs': 500}
Test Accuracy: 0.6457
[23/48] Training with params: {'embed_dim': 128, 'lr': 0.0001, 'weight_decay': 0.001, 'epochs': 200}
Test Accuracy: 0.5984
[24/48] Training with params: {'embed_dim': 128, 'lr': 0.0001, 'weight_decay': 0.001, 'epochs': 500}
Test Accuracy: 0.6063
[25/48] Training with params: {'embed_dim': 256, 'lr': 0.001, 'weight_decay': 0.01, 'epochs': 200}
Test Accuracy: 0.7087
🌟 NEW BEST ACCURACY: 0.7087
[26/48] Training with params: {'embed_dim': 256, 'lr': 0.001, 'weight_decay': 0.01, 'epochs': 500}
Test Accuracy: 0.7323
🌟 NEW BEST ACCURACY: 0.7323
[27/48] Training with params: {'embed_dim': 256, 'lr': 0.001, 'weight_decay': 0.001, 'epochs': 200}
Test Accuracy: 0.6457
[28/48] Training with params: {'embed_dim': 256, 'lr': 0.001, 'weight_decay': 0.001, 'epochs': 500}
Test Accuracy: 0.6142
[29/48] Training with params: {'embed_dim': 256, 'lr': 0.0005, 'weight_decay': 0.01, 'epochs': 200}
Test Accuracy: 0.5984
[30/48] Training with params: {'embed_dim': 256, 'lr': 0.0005, 'weight_decay': 0.01, 'epochs': 500}
Test Accuracy: 0.5906
[31/48] Training with params: {'embed_dim': 256, 'lr': 0.0005, 'weight_decay': 0.001, 'epochs': 200}
Test Accuracy: 0.6142
[32/48] Training with params: {'embed_dim': 256, 'lr': 0.0005, 'weight_decay': 0.001, 'epochs': 500}
Test Accuracy: 0.7008
[33/48] Training with params: {'embed_dim': 256, 'lr': 0.0001, 'weight_decay': 0.01, 'epochs': 200}
Test Accuracy: 0.6299
[34/48] Training with params: {'embed_dim': 256, 'lr': 0.0001, 'weight_decay': 0.01, 'epochs': 500}
Test Accuracy: 0.6299
[35/48] Training with params: {'embed_dim': 256, 'lr': 0.0001, 'weight_decay': 0.001, 'epochs': 200}
Test Accuracy: 0.6378
[36/48] Training with params: {'embed_dim': 256, 'lr': 0.0001, 'weight_decay': 0.001, 'epochs': 500}
Test Accuracy: 0.6457
[37/48] Training with params: {'embed_dim': 512, 'lr': 0.001, 'weight_decay': 0.01, 'epochs': 200}
Test Accuracy: 0.6142
[38/48] Training with params: {'embed_dim': 512, 'lr': 0.001, 'weight_decay': 0.01, 'epochs': 500}
Test Accuracy: 0.6457
[39/48] Training with params: {'embed_dim': 512, 'lr': 0.001, 'weight_decay': 0.001, 'epochs': 200}
Test Accuracy: 0.6378
[40/48] Training with params: {'embed_dim': 512, 'lr': 0.001, 'weight_decay': 0.001, 'epochs': 500}
Test Accuracy: 0.6299
[41/48] Training with params: {'embed_dim': 512, 'lr': 0.0005, 'weight_decay': 0.01, 'epochs': 200}
Test Accuracy: 0.6378
[42/48] Training with params: {'embed_dim': 512, 'lr': 0.0005, 'weight_decay': 0.01, 'epochs': 500}
Test Accuracy: 0.6220
[43/48] Training with params: {'embed_dim': 512, 'lr': 0.0005, 'weight_decay': 0.001, 'epochs': 200}
Test Accuracy: 0.6142
[44/48] Training with params: {'embed_dim': 512, 'lr': 0.0005, 'weight_decay': 0.001, 'epochs': 500}
Test Accuracy: 0.7008
[45/48] Training with params: {'embed_dim': 512, 'lr': 0.0001, 'weight_decay': 0.01, 'epochs': 200}
Test Accuracy: 0.6457
[46/48] Training with params: {'embed_dim': 512, 'lr': 0.0001, 'weight_decay': 0.01, 'epochs': 500}
Test Accuracy: 0.7323
[47/48] Training with params: {'embed_dim': 512, 'lr': 0.0001, 'weight_decay': 0.001, 'epochs': 200}
Test Accuracy: 0.6693
[48/48] Training with params: {'embed_dim': 512, 'lr': 0.0001, 'weight_decay': 0.001, 'epochs': 500}
Test Accuracy: 0.6457
✅ Đã lưu mô hình tốt nhất (Acc: 0.7323) vào land_classification_model/model_swin-unet_optimized_95.joblib
Cấu hình tốt nhất: {'embed_dim': 256, 'lr': 0.001, 'weight_decay': 0.01, 'epochs': 500}
-117
View File
@@ -1,117 +0,0 @@
🚀 CHIẾN LƯỢC TOÀN DIỆN ĐẠT >95% ACCURACY
============================================================
Loaded data: X=(706, 24, 16, 16), y=(706,)
Labels unique: [-1 0 1 2 3 4 5 6]
After cleanup: X=(652, 24, 16, 16), y=(652,) (removed 54 bad samples)
Remapped labels: [0 1 2 3 4 5 6]
Class 0: 65 samples
Class 1: 52 samples
Class 2: 48 samples
Class 3: 72 samples
Class 4: 108 samples
Class 5: 219 samples
Class 6: 88 samples
============================================================
STRATEGY 5: Flat pixel features + XGBoost (sanity check)
============================================================
Flat features: (652, 6144)
✅ Flat XGBoost acc: 0.7328
============================================================
STRATEGY 1: Lightweight CNN (no upsampling)
============================================================
Device: cuda
Epoch 1/300 Loss=1.8379 Acc=0.0763 🌟
Epoch 2/300 Loss=1.5825 Acc=0.2824 🌟
Epoch 3/300 Loss=1.4412 Acc=0.5649 🌟
Epoch 4/300 Loss=1.3274 Acc=0.6336 🌟
Epoch 5/300 Loss=1.2893 Acc=0.6870 🌟
Epoch 8/300 Loss=1.2140 Acc=0.7099 🌟
Epoch 11/300 Loss=1.2224 Acc=0.7252 🌟
Epoch 14/300 Loss=1.1087 Acc=0.7328 🌟
Epoch 16/300 Loss=1.1081 Acc=0.7710 🌟
Epoch 18/300 Loss=1.1847 Acc=0.7939 🌟
Epoch 20/300 Loss=1.0316 Acc=0.7786 (patience=2)
Epoch 24/300 Loss=1.0465 Acc=0.8092 🌟
Epoch 26/300 Loss=0.9583 Acc=0.8397 🌟
Epoch 40/300 Loss=0.9361 Acc=0.8626 🌟
Epoch 60/300 Loss=0.9807 Acc=0.8092 (patience=20)
Epoch 80/300 Loss=0.9459 Acc=0.8015 (patience=40)
Epoch 100/300 Loss=0.8443 Acc=0.7939 (patience=60)
Early stop at epoch 100
✅ LightCNN best acc: 0.8626
============================================================
STRATEGY 2: Hybrid CNN embeddings + XGBoost
============================================================
CNN embeddings: (652, 256)
Extracted 316 rich features per sample
Combined features: (652, 572)
✅ Hybrid XGBoost acc: 0.8244
============================================================
STRATEGY 3: Rich Features + Stacking Ensemble
============================================================
Extracted 316 rich features per sample
XGBoost: 0.7939
LightGBM: 0.7786
ExtraTrees: 0.7863
RandomForest: 0.7710
GBM: 0.7710
[06:55:44] WARNING: /__w/xgboost/xgboost/src/learner.cc:782:
Parameters: { "use_label_encoder" } are not used.
[06:55:44] WARNING: /__w/xgboost/xgboost/src/learner.cc:782:
Parameters: { "use_label_encoder" } are not used.
[06:55:44] WARNING: /__w/xgboost/xgboost/src/learner.cc:782:
Parameters: { "use_label_encoder" } are not used.
[06:55:44] WARNING: /__w/xgboost/xgboost/src/learner.cc:782:
Parameters: { "use_label_encoder" } are not used.
[06:55:44] WARNING: /__w/xgboost/xgboost/src/learner.cc:782:
Parameters: { "use_label_encoder" } are not used.
[06:56:17] WARNING: /__w/xgboost/xgboost/src/common/error_msg.cc:62: Falling back to prediction using DMatrix due to mismatched devices. This might lead to higher memory usage and slower performance. XGBoost is running on: cuda:0, while the input data is on: cpu.
Potential solutions:
- Use a data structure that matches the device ordinal in the booster.
- Set the device for booster before call to inplace_predict.
This warning will only be shown once.
Stacking Ensemble: 0.7710
Voting Ensemble: 0.7786
✅ Best ensemble: XGBoost = 0.7939
Extracted 316 rich features per sample
============================================================
STRATEGY 4: 5-Fold Stratified Cross-Validation
============================================================
Fold 1: 0.7939
Fold 2: 0.8244
Fold 3: 0.7769
Fold 4: 0.8154
Fold 5: 0.7615
✅ CV Mean: 0.7944 ± 0.0234
============================================================
📊 TỔNG KẾT KẾT QUẢ
============================================================
📈 LightCNN: 0.8626
📈 Hybrid CNN+XGBoost: 0.8244
📈 CV Mean (XGBoost rich): 0.7944
📈 Ensemble XGBoost: 0.7939
📈 Ensemble ExtraTrees: 0.7863
📈 Ensemble LightGBM: 0.7786
📈 Ensemble Voting: 0.7786
📈 Ensemble RandomForest: 0.7710
📈 Ensemble GBM: 0.7710
📈 Ensemble Stacking: 0.7710
📈 Flat XGBoost (baseline): 0.7328
🏆 BEST: LightCNN = 0.8626
✅ Kết quả đã được lưu vào model_train/ultimate_results.json
⚠️ Chưa đạt 95%. Best = 0.8626. Cần thêm dữ liệu hoặc feature engineering.
-69
View File
@@ -1,69 +0,0 @@
🚀 CHIẾN LƯỢC V2: TOÀN DIỆN ĐẠT >95%
============================================================
Clean data: (652, 24, 16, 16), 7 classes
============================================================
CNN + TTA (Test-Time Augmentation)
============================================================
Ep 1 Loss=1.6500 TTA-Acc=0.1450 🌟
Ep 2 Loss=1.4633 TTA-Acc=0.3511 🌟
Ep 3 Loss=1.3316 TTA-Acc=0.6336 🌟
Ep 4 Loss=1.2101 TTA-Acc=0.7252 🌟
Ep 6 Loss=1.2056 TTA-Acc=0.8015 🌟
Ep 13 Loss=1.1216 TTA-Acc=0.8244 🌟
Ep 23 Loss=0.9355 TTA-Acc=0.8321 🌟
Ep 30 Loss=0.9346 TTA-Acc=0.8092 (pat=7)
Ep 60 Loss=0.9658 TTA-Acc=0.7634 (pat=37)
Ep 69 Loss=0.8988 TTA-Acc=0.8397 🌟
Ep 74 Loss=0.9159 TTA-Acc=0.8550 🌟
Ep 90 Loss=0.8761 TTA-Acc=0.8168 (pat=16)
Ep 120 Loss=0.9473 TTA-Acc=0.8092 (pat=46)
Ep 139 Loss=0.9317 TTA-Acc=0.8626 🌟
Ep 150 Loss=0.7716 TTA-Acc=0.8397 (pat=11)
Ep 175 Loss=0.7432 TTA-Acc=0.8702 🌟
Ep 180 Loss=0.8290 TTA-Acc=0.8702 (pat=5)
Ep 210 Loss=0.8060 TTA-Acc=0.8626 (pat=35)
Ep 240 Loss=0.6980 TTA-Acc=0.8473 (pat=65)
Early stop ep 255
✅ CNN+TTA best: 0.8702
============================================================
RICH FEATURES V2 + ENSEMBLE
============================================================
Extracted 1713 features per sample
XGB: 0.7557
LGBM: 0.7634
ET: 0.7481
RF: 0.7557
Voting: 0.7634
5-Fold CV:
Fold 1: 0.7557
Fold 2: 0.7939
Fold 3: 0.7769
Fold 4: 0.8308
Fold 5: 0.8000
CV: 0.7915 ± 0.0249
============================================================
HYBRID V2: CNN embed + Rich features + XGBoost
============================================================
Extracted 1713 features per sample
Combined: (652, 2097)
✅ Hybrid V2: 0.8244
CV: 0.9142 ± 0.0194
============================================================
📊 KẾT QUẢ TỔNG HỢP V2
============================================================
✅ Hybrid CV: 0.9142
📈 CNN+TTA: 0.8702
📈 Hybrid V2: 0.8244
📈 Ens CV: 0.7915
📈 Ens_LGBM: 0.7634
📈 Ens_Vote: 0.7634
📈 Ens_XGB: 0.7557
📈 Ens_RF: 0.7557
📈 Ens_ET: 0.7481
🏆 BEST: Hybrid CV = 0.9142
-27
View File
@@ -1,27 +0,0 @@
🚀 V3: MULTI-SEED ENSEMBLE + T0-ONLY + SELF-TRAINING
============================================================
Clean: (652, 24, 16, 16), 7 classes, [65, 52, 48, 72, 108, 219, 88]
============================================================
MULTI-SEED CNN ENSEMBLE (10 models)
============================================================
Seed 0: 0.8473
Seed 1: 0.8702
Seed 2: 0.8550
Seed 3: 0.8550
Seed 4: 0.8702
Seed 5: 0.8626
Seed 6: 0.8702
Seed 7: 0.8702
Seed 8: 0.8702
Seed 9: 0.8702
✅ 10-Model Ensemble TTA: 0.8473
============================================================
TIMESTEP-0-ONLY XGBoost (cleanest data)
============================================================
T0 valid: 539/652
Features: (539, 1638)
XGB t0: 0.6852
LGBM t0: 0.6296
ET t0: 0.6852
-15
View File
@@ -1,15 +0,0 @@
🚀 BẮT ĐẦU HUẤN LUYỆN MÔ HÌNH XGBOOST (GPU & CACHE)
Initializing FeatureExtractor (mode=extended)...
📦 Đang load cache: training_data_507bd2ba4ec0d3fe107839cbf73a7a7d.joblib...
✅ Loaded 632 samples từ cache!
⚡ Đã bỏ qua download위성 data (tiết kiệm thời gian)
[CACHE HIT] Using cached dataset with 632 samples
Training XGBOOST model...
Evaluating model...
Generating classification report...
Saving model...
[MODEL MANAGER] Saving model to: model_train/model_xgboost_auto.joblib
[MODEL MANAGER] Saving metadata to: model_train/model_xgboost_auto_info.json
[MODEL MANAGER] Model saved successfully!
Training complete!
✅ Hoàn thành! Accuracy: 0.6378
File diff suppressed because one or more lines are too long
-231
View File
@@ -1,231 +0,0 @@
import re
import os
with open('/home/x79/remote-sensing/new_import_ODC.py', 'r', encoding='utf-8') as f:
content = f.read()
# 1. load_data
load_data_replacement = """def load_data(dc, date_range, longtitude_range, latitude_range):
import os, hashlib
cache_dir = "dataset_cache"
os.makedirs(cache_dir, exist_ok=True)
key_str = f"s2_{date_range}_{longtitude_range}_{latitude_range}"
cache_key = hashlib.md5(key_str.encode()).hexdigest() + ".nc"
cache_path = os.path.join(cache_dir, cache_key)
if os.path.exists(cache_path):
print(f"✅ Loading cached S2 data from {cache_path}")
return xr.open_dataset(cache_path)
product = 's2_l2a'
bbox = [longtitude_range[0], latitude_range[0], longtitude_range[1], latitude_range[1]]
import pystac_client
import planetary_computer
import odc.stac
catalog = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace,
)
search = catalog.search(
collections=["sentinel-2-l2a"],
bbox=bbox,
datetime=f"{date_range[0]}/{date_range[1]}",
)
items = list(search.items())
data = odc.stac.load(
items,
bands=["red", "nir", "SCL"],
bbox=bbox,
crs="EPSG:32648",
resolution=RESOLUTION,
chunks={"x": 2048, "y": 2048, "time": 1},
groupby="solar_day"
)
if "SCL" in data.data_vars:
data = data.rename({"SCL": "scl"})
print(f"💾 Caching S2 data to {cache_path}")
data = data.compute()
data.to_netcdf(cache_path)
return data"""
content = re.sub(r'def load_data\(dc, date_range, longtitude_range, latitude_range\):.*?return data', load_data_replacement, content, flags=re.DOTALL)
# 2. load_sen1
load_sen1_replacement = """def load_sen1(bbox, time_range):
import os, hashlib
cache_dir = "dataset_cache"
os.makedirs(cache_dir, exist_ok=True)
key_str = f"s1_vh_vv_{bbox}_{time_range}"
cache_key_vh = hashlib.md5((key_str + "vh").encode()).hexdigest() + ".nc"
cache_key_vv = hashlib.md5((key_str + "vv").encode()).hexdigest() + ".nc"
cache_path_vh = os.path.join(cache_dir, cache_key_vh)
cache_path_vv = os.path.join(cache_dir, cache_key_vv)
if os.path.exists(cache_path_vh) and os.path.exists(cache_path_vv):
print(f"✅ Loading cached S1 data from {cache_path_vh} and {cache_path_vv}")
return xr.open_dataarray(cache_path_vh), xr.open_dataarray(cache_path_vv)
import pystac_client
import planetary_computer
import odc.stac
# Kết nối STAC Client
catalog = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace,
)
# Tìm kiếm Items
search = catalog.search(
collections=["sentinel-1-rtc"],
bbox=bbox,
datetime=time_range,
)
items = list(search.items())
# Tải dữ liệu thành xarray Dataset
ds_s1 = odc.stac.load(
items,
bands=["vv", "vh"],
bbox=bbox,
crs="EPSG:32648",
resolution=RESOLUTION,
chunks={"x": 2048, "y": 2048, "time": 1}
)
# Tính giá trị trung vị theo thời gian
ds_median = ds_s1.median(dim="time").compute()
vv = ds_median["vv"]
vh = ds_median["vh"]
# Thêm chiều 'band' để giống hệt rioxarray
vv = vv.expand_dims(dim="band")
vh = vh.expand_dims(dim="band")
# Phục hồi metadata về toạ độ
vv = vv.rio.write_crs("EPSG:32648")
vh = vh.rio.write_crs("EPSG:32648")
print(f"💾 Caching S1 data to {cache_dir}")
vh.to_netcdf(cache_path_vh)
vv.to_netcdf(cache_path_vv)
return vh, vv"""
content = re.sub(r'def load_sen1\(bbox, time_range\):.*?return vh, vv', load_sen1_replacement, content, flags=re.DOTALL)
# 3. load_data_sen1
load_data_sen1_replacement = """def load_data_sen1(dc, date_range, coordinates):
import os, hashlib
longtitude_range, latitude_range = coordinates
bbox = [longtitude_range[0], latitude_range[0], longtitude_range[1], latitude_range[1]]
cache_dir = "dataset_cache"
os.makedirs(cache_dir, exist_ok=True)
key_str = f"data_sen1_{date_range}_{bbox}"
cache_key_vh = hashlib.md5((key_str + "vh").encode()).hexdigest() + ".nc"
cache_key_vv = hashlib.md5((key_str + "vv").encode()).hexdigest() + ".nc"
cache_path_vh = os.path.join(cache_dir, cache_key_vh)
cache_path_vv = os.path.join(cache_dir, cache_key_vv)
if os.path.exists(cache_path_vh) and os.path.exists(cache_path_vv):
print(f"✅ Loading cached S1 (coord) data")
return xr.open_dataarray(cache_path_vh), xr.open_dataarray(cache_path_vv)
import pystac_client
import planetary_computer
import odc.stac
catalog = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace,
)
search = catalog.search(
collections=["sentinel-1-rtc"],
bbox=bbox,
datetime=f"{date_range[0]}/{date_range[1]}",
)
items = list(search.items())
data_sen1 = odc.stac.load(
items,
bands=["vv", "vh"],
bbox=bbox,
crs="EPSG:32648",
resolution=RESOLUTION,
chunks={"x": 2048, "y": 2048, "time": 1},
groupby="solar_day"
)
data_sen1 = data_sen1.compute()
dsvh = data_sen1.vh
dsvv = data_sen1.vv
print(f"💾 Caching S1 (coord) data")
dsvh.to_netcdf(cache_path_vh)
dsvv.to_netcdf(cache_path_vv)
return dsvh, dsvv"""
content = re.sub(r'def load_data_sen1\(dc, date_range, coordinates\):.*?return dsvh, dsvv', load_data_sen1_replacement, content, flags=re.DOTALL)
# 4. load_data_sen2
load_data_sen2_replacement = """def load_data_sen2(dc, date_range, coordinates):
import os, hashlib
longtitude_range, latitude_range = coordinates
bbox = [longtitude_range[0], latitude_range[0], longtitude_range[1], latitude_range[1]]
cache_dir = "dataset_cache"
os.makedirs(cache_dir, exist_ok=True)
key_str = f"data_sen2_{date_range}_{bbox}"
cache_key = hashlib.md5(key_str.encode()).hexdigest() + ".nc"
cache_path = os.path.join(cache_dir, cache_key)
if os.path.exists(cache_path):
print(f"✅ Loading cached S2 (coord) data from {cache_path}")
return xr.open_dataset(cache_path)
import pystac_client
import planetary_computer
import odc.stac
catalog = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace,
)
search = catalog.search(
collections=["sentinel-2-l2a"],
bbox=bbox,
datetime=f"{date_range[0]}/{date_range[1]}",
)
items = list(search.items())
data = odc.stac.load(
items,
bands=["red", "nir", "SCL"],
bbox=bbox,
crs="EPSG:32648",
resolution=RESOLUTION,
chunks={"x": 2048, "y": 2048, "time": 1},
groupby="solar_day"
)
if "SCL" in data.data_vars:
data = data.rename({"SCL": "scl"})
data = data.compute()
print(f"💾 Caching S2 (coord) data to {cache_path}")
data.to_netcdf(cache_path)
return data"""
content = re.sub(r'def load_data_sen2\(dc, date_range, coordinates\):.*?return data', load_data_sen2_replacement, content, flags=re.DOTALL)
with open('/home/x79/remote-sensing/new_import_ODC.py', 'w', encoding='utf-8') as f:
f.write(content)
print("Patching successful.")
-29
View File
@@ -1,29 +0,0 @@
import json
import glob
def fix_load_sen1(file_path):
with open(file_path, 'r', encoding='utf-8') as f:
nb = json.load(f)
changed = False
for cell in nb.get('cells', []):
if cell.get('cell_type') == 'code':
source = cell.get('source', [])
for i, line in enumerate(source):
if 'load_sen1(name_vh, name_vv)' in line:
indent = line[:len(line) - len(line.lstrip())]
replacement = (
f"{indent}bbox = [longtitude_range[0], latitude_range[0], longtitude_range[1], latitude_range[1]]\n"
f"{indent}time_range = f'{{date_range[0]}}/{{date_range[1]}}'\n"
f"{indent}{line.lstrip().replace('load_sen1(name_vh, name_vv)', 'load_sen1(bbox, time_range)')}"
)
source[i] = replacement
changed = True
if changed:
with open(file_path, 'w', encoding='utf-8') as f:
json.dump(nb, f, indent=1)
print(f"Patched load_sen1 in {file_path}")
for nb in glob.glob("*.ipynb"):
fix_load_sen1(nb)
-43
View File
@@ -1,43 +0,0 @@
import json
import glob
def patch_notebook(file_path):
with open(file_path, 'r', encoding='utf-8') as f:
nb = json.load(f)
changed = False
for cell in nb.get('cells', []):
if cell.get('cell_type') == 'code':
source = cell.get('source', [])
# Check if this cell should be fully commented out
full_source = ''.join(source)
if 'dc.load(' in full_source or 'ds.vv' in full_source:
for i in range(len(source)):
if not source[i].startswith('#'):
source[i] = '# ' + source[i]
changed = True
continue
# Otherwise, do line-by-line replacements
for i, line in enumerate(source):
if 'ST_training data_updated_1130points.shp' in line:
source[i] = line.replace('ST_training data_updated_1130points.shp', 'ST_training_data_updated_1130points.shp')
changed = True
if 'from new_import import *' in line:
source[i] = line.replace('from new_import import *', 'from new_import_ODC import *')
changed = True
if 'dc = datacube.Datacube()' in line:
source[i] = line.replace('dc = datacube.Datacube()', 'dc = None')
changed = True
if 'load_data(dc,' in line:
source[i] = line.replace('load_data(dc,', 'load_data(None,')
changed = True
if changed:
with open(file_path, 'w', encoding='utf-8') as f:
json.dump(nb, f, indent=1)
print(f"Patched {file_path}")
for nb in glob.glob("*.ipynb"):
patch_notebook(nb)
-75
View File
@@ -1,75 +0,0 @@
import json
import glob
import re
def patch_python_script(filepath):
try:
with open(filepath, 'r', encoding='utf-8') as f:
content = f.read()
original_content = content
# Replace imports
content = re.sub(r'from sklearn\.ensemble import RandomForestClassifier',
'from xgboost import XGBClassifier', content)
# Replace the model instantiations (for 01.train_ODC.py)
rf_pattern = re.compile(r'model\s*=\s*RandomForestClassifier\([^)]+\)', re.DOTALL)
xgb_replacement = """model = XGBClassifier(
n_estimators=200,
max_depth=30,
tree_method="hist",
device="cuda",
random_state=42,
n_jobs=-1,
verbosity=1
)"""
content = rf_pattern.sub(xgb_replacement, content)
if content != original_content:
with open(filepath, 'w', encoding='utf-8') as f:
f.write(content)
print(f"Patched {filepath}")
except Exception as e:
print(f"Error patching {filepath}: {e}")
def patch_notebook(filepath):
try:
with open(filepath, 'r', encoding='utf-8') as f:
nb = json.load(f)
changed = False
import_pattern = re.compile(r'from\s+sklearn\.ensemble\s+import\s+RandomForestClassifier')
inst_pattern = re.compile(r'RandomForestClassifier\([^)]*\)')
for cell in nb.get('cells', []):
if cell.get('cell_type') == 'code':
source = cell.get('source', [])
for i in range(len(source)):
if import_pattern.search(source[i]):
source[i] = import_pattern.sub('from xgboost import XGBClassifier', source[i])
changed = True
if inst_pattern.search(source[i]):
source[i] = inst_pattern.sub("XGBClassifier(tree_method='hist', device='cuda', random_state=42, n_jobs=-1)", source[i])
changed = True
if "'classifier__criterion': ['gini', 'entropy']" in source[i]:
source[i] = source[i].replace("'classifier__criterion': ['gini', 'entropy']",
"'classifier__learning_rate': [0.01, 0.1, 0.2]")
changed = True
if changed:
with open(filepath, 'w', encoding='utf-8') as f:
json.dump(nb, f, indent=1)
print(f"Patched {filepath}")
except Exception as e:
print(f"Error patching {filepath}: {e}")
if __name__ == "__main__":
patch_python_script("01.train_ODC.py")
patch_python_script("new_train.py")
for nb in glob.glob("*.ipynb"):
patch_notebook(nb)
-272
View File
@@ -1,272 +0,0 @@
#!/usr/bin/env python
# coding: utf-8
# In[1]:
get_ipython().run_cell_magic('time', '', '%matplotlib inline\nfrom new_import import *\n')
# In[2]:
get_ipython().run_cell_magic('time', '', '# Dask gateway\ncluster, client = notebook_utils.initialize_dask(use_gateway=True, workers=(1,4))\ndc = datacube.Datacube()\n\n# Configure s3 access\nconfigure_s3_access(aws_unsigned=False, requester_pays=True, client=client)\n\nclient\n')
# In[3]:
## cấu hình thời gian lấy ảnh và tọa độ
# date_range = ('2022-09-01', '2023-10-01')
# longtitude_range = (105.86575, 105.94120)
# latitude_range = (9.65070, 9.69850)
date_range = ('2022-09-01', '2023-10-01')
longtitude_range = (105.5, 106.4)
latitude_range = (9.2, 10.0)
# In[4]:
## truy vấn ảnh vệ tinh sen2
data = load_data(dc, date_range, longtitude_range, latitude_range)
notebook_utils.heading(notebook_utils.xarray_object_size(data))
display(data)
# In[5]:
# Specify the start and end times
min_date = '2022-09-01' # Thời gian bắt đầu lấy data cho quá trình train
max_date = '2023-10-01' # Thời gian kết thúc lấy data cho quá trình train
# Just do 1 month for testing
# max_date = '2022-10-01' # Thời gian kết thúc lấy data cho quá trình train
# Specify a spatail region to search using latitude/longitude cooridinates
min_longitude, max_longitude = (105.5, 106.4)
min_latitude, max_latitude = (9.2, 10.0)
# Specify the product. In this case we want to use Sentinel-2 Level-2A data
product = 's2_l2a'
# Construct the search query dictionary
query = {
'product': product, # Product name
'x': (min_longitude, max_longitude), # "x" axis bounds
'y': (min_latitude, max_latitude), # "y" axis bounds
'time': (min_date, max_date), # Any parsable date strings
}
# In[6]:
# Most common CRS
native_crs = notebook_utils.mostcommon_crs(dc, query)
print(f'Most common native CRS: {native_crs}')
# In[7]:
# Specify the spectral band measurements we want to use for a classification algorithm
measurements = ['red', 'nir', 'scl']
load_params = {
'measurements': measurements, # Selected measurement or alias names
'output_crs': native_crs, # Target EPSG code
'resolution': (-10, 10), # Target resolution
'group_by': 'solar_day', # Scene grouping
'dask_chunks': {'x': 2048, 'y': 2048}, # Dask chunks
}
# In[8]:
get_ipython().run_cell_magic('time', '', '# The replacement "dc.load()" function for this product\ndata = load_s2l2a_with_offset(\n dc,\n query | load_params # Combine the two dicts that contain our search and load parameters\n)\n\n# This line prints the total size of the dataset hat was loaded\nnotebook_utils.heading(notebook_utils.xarray_object_size(data))\n\ndisplay(data)\n')
# In[9]:
# %%time
# # Tiến hành loại bỏ các vị trí bị mây ảnh hưởng
# result = mask_clean(data)
# progress(result)
# In[10]:
# Tiến hành tính toán NDVI
ds1 = calculate_indices(result, index='NDVI', satellite_mission='s2')
ndvi = ds1["NDVI"]
display(ndvi)
# In[11]:
get_ipython().run_cell_magic('time', '', "## tính ndvi theo tháng\naverage_ndvi = ndvi.resample(time='1M').mean().persist()\nprogress(average_ndvi)\n")
# In[12]:
# compute average_ndvi
average_ndvi = average_ndvi.compute()
# In[13]:
# cấu hình vh vv file
# name_vh = "ThuanHoa/ThuanHoa_VH.tif"
# name_vv = "ThuanHoa/ThuanHoa_VV.tif"
# load dữ liệu sen1
bbox = [105.5, 9.2, 106.4, 10.0]
time_range = '2022-09-01/2023-10-01'
# dsvh, dsvv = load_sen1(bbox, time_range)
name_vh = "vh-0922_0923-full_ST.tif"
name_vv = "vv-0922_0923-full_ST.tif"
if not os.path.exists(name_vh):
get_ipython().system('aws s3 cp s3://easi-asia-dc-data/staging/ctu/sentinel-1/vh-0922_0923-full_ST.tif vh-0922_0923-full_ST.tif')
if not os.path.exists(name_vv):
get_ipython().system('aws s3 cp s3://easi-asia-dc-data/staging/ctu/sentinel-1/vv-0922_0923-full_ST.tif vv-0922_0923-full_ST.tif')
bbox = [105.5, 9.2, 106.4, 10.0]
time_range = '2022-09-01/2023-10-01'
dsvh, dsvv = load_sen1(bbox, time_range)
# In[27]:
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
# In[28]:
average_ndvi = average_ndvi[:, :7680, :8687]
mask = ~np.isnan(average_ndvi)
print(average_ndvi.shape)
print(dsvh.shape)
print(dsvv.shape)
print(mask.shape)
X_train = np.stack([dsvh.values[mask], dsvv.values[mask]], axis=1)
y_train = average_ndvi.values[mask]
# In[29]:
model = LinearRegression()
model.fit(X_train, y_train)
# In[30]:
X_pred = np.stack([dsvh.values[~mask], dsvv.values[~mask]], axis=1)
average_ndvi.values[~mask] = model.predict(X_pred)
# In[31]:
average_ndvi_filled = xr.DataArray(average_ndvi, dims=average_ndvi.dims)
# In[32]:
plt.imshow(average_ndvi_filled.isel(time=6))
# In[65]:
plt.imshow(average_ndvi.isel(time=6))
# In[33]:
train_path = "train/ST_training data_updated_1130points.shp"
# In[34]:
train = load_train_data(train_path)
# In[37]:
datasets = get_data_sen1_and_sen2(train, average_ndvi_filled, dsvh, dsvv)
# In[39]:
# cấu hình nhãn dữ liệu
label_mapping = {
"Lua tom": "0",
"Lua": "1",
"CHN": "2",
"CLN": "3",
"TS": "4",
"Song": "5",
"Dat xay dung": "6",
"Rung": "7"
}
# chia tập dữ liệu train, val, test
X_train, X_val, X_test, y_train, y_val, y_test = split_train_data(train, label_mapping, datasets)
# In[40]:
# Huấn luyện mô hình
grid_search = train_with_rf(X_train, X_val, y_train, y_val)
# In[41]:
# kiểm tra độ chính xác với tập test
y_pred_test = grid_search.predict(X_test)
test_accuracy = accuracy_score(y_test, y_pred_test)
print(f"Accuracy for test data {round(test_accuracy, 2)*100} %")
# In[42]:
# Lưu mô hình huấn luyện
save_model("model_new.joblib", grid_search)
# In[43]:
# đóng client, cluster
client.close()
cluster.close()
# In[ ]:
-127
View File
@@ -1,127 +0,0 @@
#!/usr/bin/env python
# coding: utf-8
# In[1]:
get_ipython().run_cell_magic('time', '', '%matplotlib inline\nfrom new_import import *\n')
# In[2]:
get_ipython().run_cell_magic('time', '', '# Dask gateway\ncluster, client = notebook_utils.initialize_dask(use_gateway=True, workers=(1,4))\ndc = datacube.Datacube()\n\n# Configure s3 access\nconfigure_s3_access(aws_unsigned=False, requester_pays=True, client=client)\n\nclient\n')
# In[3]:
## cấu hình thời gian lấy ảnh và tọa độ
date_range = ('2022-09-01', '2023-10-01')
longtitude_range = (105.86575, 105.94120)
latitude_range = (9.65070, 9.69850)
# In[4]:
## truy vấn ảnh vệ tinh sen2
data = load_data(dc, date_range, longtitude_range, latitude_range)
notebook_utils.heading(notebook_utils.xarray_object_size(data))
display(data)
# In[5]:
get_ipython().run_cell_magic('time', '', '# Tiến hành loại bỏ các vị trí bị mây ảnh hưởng\nresult = mask_clean(data)\nprogress(result)\n')
# In[6]:
# Tiến hành tính toán NDVI
ds1 = calculate_indices(result, index='NDVI', satellite_mission='s2')
ndvi = ds1["NDVI"]
display(ndvi)
# In[17]:
get_ipython().run_cell_magic('time', '', "## tính ndvi theo tháng\naverage_ndvi = ndvi.resample(time='1M').mean().persist()\nprogress(average_ndvi)\n")
# In[18]:
# compute average_ndvi
average_ndvi = average_ndvi.compute()
# In[9]:
# cấu hình vh vv file
name_vh = "ThuanHoa/ThuanHoa_VH.tif"
name_vv = "ThuanHoa/ThuanHoa_VV.tif"
# load dữ liệu sen1
bbox = [105.5, 9.2, 106.4, 10.0]
time_range = '2022-09-01/2023-10-01'
dsvh, dsvv = load_sen1(bbox, time_range)
# In[10]:
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
# In[11]:
mask = ~np.isnan(average_ndvi)
X_train = np.stack([dsvh.values[mask], dsvv.values[mask]], axis=1)
y_train = average_ndvi.values[mask]
# In[12]:
model = LinearRegression()
model.fit(X_train, y_train)
# In[13]:
X_pred = np.stack([dsvh.values[~mask], dsvv.values[~mask]], axis=1)
average_ndvi.values[~mask] = model.predict(X_pred)
# In[14]:
average_ndvi_filled = xr.DataArray(average_ndvi, dims=average_ndvi.dims)
# In[16]:
plt.imshow(average_ndvi_filled.isel(time=1))
# In[19]:
plt.imshow(average_ndvi.isel(time=1))
# In[ ]:
-145
View File
@@ -1,145 +0,0 @@
import json
def get_content(filename):
with open(filename, "r", encoding="utf-8") as f:
return f.read()
fe_content = get_content("feature_extractor.py")
tm_content = get_content("train_module.py")
dt_content = get_content("train_land_decision_tree_gpu.py")
notebook = {
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Tải Dữ liệu Vệ tinh qua Colab (Self-contained)\n",
"Notebook này đã được nhúng sẵn toàn bộ mã nguồn xử lý. Bạn không cần upload cả thư mục `remote-sensing` nữa.\n",
"\n",
"## Bước 1: Upload Shapefile (BẮT BUỘC)\n",
"Mô hình cần biết các điểm tọa độ đất để lấy dữ liệu. Hãy nén thư mục `train/` trên máy bạn thành `train.zip` và chạy ô dưới đây để upload nó trực tiếp lên Colab."
]
},
{
"cell_type": "code",
"execution_count": None,
"metadata": {},
"outputs": [],
"source": [
"from google.colab import files\n",
"import os\n",
"\n",
"print(\"Hãy chọn file train.zip từ máy tính của bạn:\")\n",
"uploaded = files.upload()\n",
"\n",
"if \"train.zip\" in uploaded:\n",
" !unzip -q -o train.zip -d /content/train_tmp/\n",
" # Move the extracted files directly to /content/train/\n",
" !mkdir -p /content/train\n",
" !mv /content/train_tmp/*/* /content/train/ 2>/dev/null || mv /content/train_tmp/* /content/train/\n",
" print(\"Đã giải nén shapefile thành công vào thư mục /content/train/\")\n",
"else:\n",
" print(\"LỖI: Bạn chưa upload file có tên là train.zip!\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bước 2: Cài đặt thư viện & Tạo môi trường"
]
},
{
"cell_type": "code",
"execution_count": None,
"metadata": {},
"outputs": [],
"source": [
"!pip install planetary-computer pystac-client odc-stac geopandas rasterio xarray joblib scikit-learn xgboost lightgbm"
]
},
{
"cell_type": "code",
"execution_count": None,
"metadata": {},
"outputs": [],
"source": [
"%%writefile feature_extractor.py\n" + fe_content
]
},
{
"cell_type": "code",
"execution_count": None,
"metadata": {},
"outputs": [],
"source": [
"%%writefile train_module.py\n" + tm_content
]
},
{
"cell_type": "code",
"execution_count": None,
"metadata": {},
"outputs": [],
"source": [
"%%writefile train_land_decision_tree_gpu.py\n" + dt_content
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bước 3: Chạy tiến trình tải ảnh vệ tinh và tạo Cache"
]
},
{
"cell_type": "code",
"execution_count": None,
"metadata": {},
"outputs": [],
"source": [
"!python train_land_decision_tree_gpu.py"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bước 4: Tải file Cache về máy"
]
},
{
"cell_type": "code",
"execution_count": None,
"metadata": {},
"outputs": [],
"source": [
"from google.colab import files\n",
"import glob\n",
"\n",
"cache_files = glob.glob(\"dataset_cache/*.joblib\")\n",
"if cache_files:\n",
" latest_cache = max(cache_files, key=os.path.getctime)\n",
" print(f\"Đang tải file {latest_cache} về máy...\")\n",
" files.download(latest_cache)\n",
"else:\n",
" print(\"Chưa tìm thấy file cache. Hãy chắc chắn bước 3 đã chạy thành công!\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
with open("Download_Cache_Colab.ipynb", "w", encoding="utf-8") as f:
json.dump(notebook, f, indent=1, ensure_ascii=False)
print("Notebook updated successfully!")
Binary file not shown.
-35
View File
@@ -1,35 +0,0 @@
import json
import glob
def fix_notebook(file_path):
with open(file_path, 'r', encoding='utf-8') as f:
nb = json.load(f)
changed = False
for cell in nb.get('cells', []):
if cell.get('cell_type') == 'code':
source = cell.get('source', [])
if isinstance(source, list):
for i, line in enumerate(source):
if "dc = datacube.Datacube()" in line:
source[i] = "dc = None\n"
changed = True
if "ds = dc.load(" in line:
source[i] = "ds = None\n"
changed = True
if "data = dc.load(" in line:
source[i] = "data = None\n"
changed = True
# If ds is None, ds.vv will fail
if "vv_data = ds.vv" in line:
source[i] = "vv_data = None\n"
changed = True
if "notebook_utils.xarray_object_size(ds)" in line:
source[i] = line.replace("notebook_utils.xarray_object_size(ds)", "'ds is None'")
changed = True
if changed:
with open(file_path, 'w', encoding='utf-8') as f:
json.dump(nb, f, indent=1)
print(f"Removed datacube from {file_path}")
for nb in glob.glob("*.ipynb"):
fix_notebook(nb)
-110
View File
@@ -1,110 +0,0 @@
import os
import glob
import json
import subprocess
import time
from tabulate import tabulate
scripts = [
"train_land_randomforest.py",
"train_cloud_cnn.py",
"train_cloud_swin_unet.py",
"train_ndvi_statistical.py",
"train_ndvi_lstm_gru.py",
"train_ndvi_convlstm.py",
"train_ndvi_hybrid_physics.py",
"train_ndvi_ensemble.py"
]
print("🚀 Đang khởi chạy song song tất cả các mô hình...")
processes = []
for script in scripts:
if os.path.exists(script):
cmd = f"source /home/x79/miniconda3/etc/profile.d/conda.sh && conda activate env_01 && python {script}"
p = subprocess.Popen(["bash", "-c", cmd], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
processes.append((script, p))
for script, p in processes:
p.wait()
print("✅ Đã chạy xong tất cả các mô hình!\n")
print("📊 BẢNG SO SÁNH KẾT QUẢ CÁC MÔ HÌNH\n")
# 1. Phân loại đất
print("### 1. Nhóm Phân loại Lớp phủ (Land Classification)")
land_data = []
# Đọc XGBoost từ thư mục gốc
if os.path.exists("model_xgboost_info.json"):
with open("model_xgboost_info.json", 'r') as f:
data = json.load(f)
params = data.get('params', {})
param_str = f"estimators:{params.get('n_estimators')}, depth:{params.get('max_depth')}" if params else "N/A"
land_data.append([
data.get('model_type', 'XGBoost'),
data.get('accuracy', ''),
data.get('precision', ''),
data.get('recall', ''),
data.get('f1_score', ''),
param_str
])
# Đọc các model khác trong model_train
for info_file in glob.glob("model_train/*_info.json"):
with open(info_file, 'r') as f:
data = json.load(f)
# Chỉ lấy các model có độ chính xác (để lọc model rác/cũ)
if 'accuracy' not in data and 'f1_score' not in data:
continue
params = data.get('params', {})
param_str = f"estimators:{params.get('n_estimators')}, depth:{params.get('max_depth')}" if params else "N/A"
# Fallback for Random Forest
if data.get('model_type') == 'RandomForest_RealData':
param_str = "estimators:100, depth:15"
land_data.append([
data.get('model_type', ''),
data.get('accuracy', ''),
data.get('precision', ''),
data.get('recall', ''),
data.get('f1_score', ''),
param_str
])
if land_data:
print(tabulate(land_data, headers=["Model", "Accuracy", "Precision", "Recall", "F1-Score", "Parameters"], tablefmt="github"))
print("\n")
# 2. Xóa mây
print("### 2. Nhóm Xóa mây (Cloud Removal)")
cloud_data = []
for info_file in glob.glob("cloud_removal_model/*_info.json"):
with open(info_file, 'r') as f:
data = json.load(f)
cloud_data.append([
data.get('model_type', ''),
data.get('epoch', ''),
data.get('train_loss', ''),
data.get('val_loss', '')
])
if cloud_data:
print(tabulate(cloud_data, headers=["Model", "Epochs", "Train Loss", "Val Loss"], tablefmt="github"))
print("\n")
# 3. Dự báo NDVI
print("### 3. Nhóm Dự báo Thực vật (NDVI Forecasting)")
ndvi_data = []
for info_file in glob.glob("ndvi_forecast_model/*_info.json"):
with open(info_file, 'r') as f:
data = json.load(f)
ndvi_data.append([
data.get('model_type', ''),
data.get('rmse', ''),
data.get('mae', ''),
data.get('epoch', 'N/A')
])
if ndvi_data:
print(tabulate(ndvi_data, headers=["Model", "RMSE", "MAE", "Epochs"], tablefmt="github"))
print("\n")
-72
View File
@@ -1,72 +0,0 @@
#!/usr/bin/env python
# coding: utf-8
# In[2]:
get_ipython().run_line_magic('matplotlib', 'inline')
from new_import import *
# Dask gateway
cluster, client = notebook_utils.initialize_dask(use_gateway=True, workers=(1,4))
dc = datacube.Datacube()
# Configure s3 access
configure_s3_access(aws_unsigned=False, requester_pays=True, client=client)
# In[3]:
ds = dc.load(
product="sentinel1_grd_gamma0_20m",
x=(105.5, 106.4),
y=(9.2, 10.0),
time=("2022-09-01", "2023-10-01"),
measurements=["vv", "vh"],
output_crs="EPSG:32648",
resolution=(-10,10),
dask_chunks={"x":2048, "y":2048},
skip_broken_datasets=True,
group_by="solar_day"
)
notebook_utils.heading(notebook_utils.xarray_object_size(ds))
ds
# In[18]:
vh = ds.vh.resample(time='1M').mean().persist()
vh = vh.compute()
vv = ds.vv.resample(time='1M').mean().persist()
vv = vv.compute()
# In[28]:
vv.min()
# In[33]:
import matplotlib.pyplot as plt
# Plot the data
plt.imshow(vh.isel(time=0), cmap='viridis', vmin=0, vmax=1)
plt.colorbar() # Add colorbar for reference
plt.show()
# In[ ]:
-33
View File
@@ -1,33 +0,0 @@
import torch
import numpy as np
device = "cpu"
input_array = np.zeros((4, 16, 16), dtype=np.float32)
input_tensor = torch.from_numpy(input_array).unsqueeze(0).to(device)
print("Before pad:", input_tensor.shape)
from train_cloud_removal import UNet
model = UNet(in_channels=6, out_channels=4).to(device)
if hasattr(model, 'inc') and hasattr(model.inc.double_conv[0], 'in_channels'):
expected_channels = model.inc.double_conv[0].in_channels
elif hasattr(model, 'conv1') and hasattr(model.conv1, 'in_channels'):
expected_channels = model.conv1.in_channels
else:
expected_channels = list(model.parameters())[0].shape[1]
print("Expected channels:", expected_channels)
if expected_channels > input_tensor.shape[1]:
pad_channels = expected_channels - input_tensor.shape[1]
padding = torch.zeros(1, pad_channels, *input_tensor.shape[2:]).to(device)
input_tensor = torch.cat([input_tensor, padding], dim=1)
print("After pad:", input_tensor.shape)
try:
model(input_tensor)
print("Success!")
except Exception as e:
print("Error:", e)
-25
View File
@@ -1,25 +0,0 @@
import os
import sys
sys.path.insert(0, os.getcwd())
import new_import_ODC
import time
import xarray as xr
# Mock minimal params to test load_data
date_range = ('2023-01-01', '2023-01-31')
longtitude_range = (105.0, 105.1)
latitude_range = (9.5, 9.6)
print("--- First Call (Downloading & Caching) ---")
start = time.time()
data1 = new_import_ODC.load_data(None, date_range, longtitude_range, latitude_range)
end = time.time()
print(f"Time taken: {end - start:.2f}s")
print("--- Second Call (Loading from Cache) ---")
start = time.time()
data2 = new_import_ODC.load_data(None, date_range, longtitude_range, latitude_range)
end = time.time()
print(f"Time taken: {end - start:.2f}s")
print("✅ Test completed")
-7
View File
@@ -1,7 +0,0 @@
import torch
checkpoint = torch.load('cloud_removal_model/cloud_removal_unet_best.pth', map_location='cpu')
print(checkpoint.keys())
print("in_channels in checkpoint:", 'in_channels' in checkpoint)
if 'in_channels' in checkpoint:
print(checkpoint['in_channels'])
print("Shape of inc.double_conv.0.weight:", checkpoint['model_state_dict']['inc.double_conv.0.weight'].shape)
-4
View File
@@ -1,4 +0,0 @@
import torch
from cloud_removal import DeepInpaintingStrategy
cloud_remover = DeepInpaintingStrategy()
print("Model channels:", list(cloud_remover.model.parameters())[0].shape[1])
-10
View File
@@ -1,10 +0,0 @@
from cloud_removal import DeepInpaintingStrategy
import torch
import numpy as np
cr = DeepInpaintingStrategy(model_path="cloud_removal_model/cloud_removal_unet_best.pth")
if cr.model is not None:
expected = list(cr.model.parameters())[0].shape[1]
print("Expected channels:", expected)
else:
print("Failed to load model")
-53
View File
@@ -1,53 +0,0 @@
import geopandas as gpd
import planetary_computer
import pystac_client
import odc.stac
import numpy as np
from shapely.geometry import Point, shape
from pyproj import Transformer
bbox = [105.5, 9.2, 106.3, 10.0]
time_range = "2023-01-01/2023-04-30"
catalog = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace,
)
items = list(catalog.search(collections=["sentinel-2-l2a"], bbox=bbox, datetime=time_range, query={"eo:cloud_cover": {"lt": 30}}).items())
items = sorted(items, key=lambda x: x.properties["eo:cloud_cover"])
gdf = gpd.read_file("train/ST_training_data_updated_1130points_new.shp")
gdf = gdf.to_crs("EPSG:32648")
# Find a point that fails. Let's just test a few points.
for idx, row in gdf.head(20).iterrows():
x_coord = row['geometry'].x
y_coord = row['geometry'].y
transformer = Transformer.from_crs("epsg:32648", "epsg:4326", always_xy=True)
lon, lat = transformer.transform(x_coord, y_coord)
point = Point(lon, lat)
filtered = []
for item in items:
if shape(item.geometry).contains(point):
filtered.append(item)
filtered = [planetary_computer.sign(item) for item in filtered]
if not filtered:
print(f"Point {idx}: NO ITEMS CONTAINS POINT!")
continue
ds = odc.stac.load(
filtered,
bands=["B02"],
x=(x_coord - 80, x_coord + 80),
y=(y_coord - 80, y_coord + 80),
crs="EPSG:32648",
resolution=10,
patch_url=planetary_computer.sign,
fail_on_error=False
).compute()
sums = ds["B02"].sum(dim=["x", "y"]).values
non_zero = (sums > 0).sum()
print(f"Point {idx}: {len(filtered)} items, {non_zero} non-zero time steps")
-37
View File
@@ -1,37 +0,0 @@
import geopandas as gpd
import planetary_computer
import pystac_client
import odc.stac
import sys
bbox = [105.5, 9.2, 106.3, 10.0]
time_range = "2023-01-01/2023-04-30"
catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1", modifier=planetary_computer.sign_inplace)
search = catalog.search(collections=["sentinel-2-l2a"], bbox=bbox, datetime=time_range, query={"eo:cloud_cover": {"lt": 30}})
items = list(search.items())
items = sorted(items, key=lambda x: x.properties.get("eo:cloud_cover", 100))[:4]
items = [planetary_computer.sign(item) for item in items]
gdf = gpd.read_file("train/ST_training_data_updated_1130points_new.shp")
gdf = gdf.to_crs("EPSG:32648")
row = gdf.iloc[0]
x, y_coord = row.geometry.x, row.geometry.y
point_bbox = [x - 80, y_coord - 80, x + 80, y_coord + 80]
patch_s2 = odc.stac.load(
items,
bands=["B02", "B03", "B04", "B08", "SCL"],
x=(x - 80, x + 80),
y=(y_coord - 80, y_coord + 80),
crs="EPSG:32648",
resolution=10,
patch_url=planetary_computer.sign,
fail_on_error=False
).compute()
print("patch_s2 vars:", patch_s2.data_vars)
if patch_s2.dims['x'] < 16 or patch_s2.dims['y'] < 16:
print("Too small:", patch_s2.dims)
else:
print("Success dimension:", patch_s2.dims)
-3
View File
@@ -1,3 +0,0 @@
import geopandas as gpd
gdf = gpd.read_file("train/ST_training_data_updated_1130points_new.shp")
print(gdf.head(1)['HT_code'])
-3
View File
@@ -1,3 +0,0 @@
import geopandas as gpd
gdf = gpd.read_file("train/ST_training_data_updated_1130points_new.shp")
print(gdf.columns)
-23
View File
@@ -1,23 +0,0 @@
import numpy as np
from xgboost import XGBClassifier
from sklearn.datasets import make_classification
from sklearn.metrics import accuracy_score
print("🚀 Testing XGBoost with CUDA GPU...")
try:
X, y = make_classification(n_samples=10000, n_features=20, n_classes=2, random_state=42)
model = XGBClassifier(
n_estimators=100,
max_depth=10,
tree_method="hist",
device="cuda",
random_state=42,
verbosity=1
)
print("Training model...")
model.fit(X, y)
y_pred = model.predict(X)
acc = accuracy_score(y, y_pred)
print(f"✅ Training successful! Accuracy: {acc*100:.2f}%")
except Exception as e:
print(f"❌ Error during training: {e}")
-3
View File
@@ -1,3 +0,0 @@
import torch
checkpoint = torch.load('cloud_removal_model/cloud_removal_unet_best.pth', map_location='cpu')
print(list(checkpoint['model_state_dict'].keys())[:5])
-13
View File
@@ -1,13 +0,0 @@
import torch
from pathlib import Path
model = torch.load('cloud_removal_model/cloud_removal_unet_best.pth', map_location='cpu')
print(type(model))
print("hasattr inc:", hasattr(model, 'inc'))
if hasattr(model, 'inc'):
print("hasattr double_conv:", hasattr(model.inc, 'double_conv'))
if hasattr(model.inc, 'double_conv'):
print("in_channels:", model.inc.double_conv[0].in_channels)
else:
for name, param in model.named_parameters():
print(name, param.shape)
break
-19
View File
@@ -1,19 +0,0 @@
import geopandas as gpd
import planetary_computer
import pystac_client
import odc.stac
bbox = [105.5, 9.2, 106.3, 10.0]
time_range = "2023-01-01/2023-01-30"
catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1", modifier=planetary_computer.sign_inplace)
items = list(catalog.search(collections=["sentinel-2-l2a"], bbox=bbox, datetime=time_range).items())[:1]
items = [planetary_computer.sign(item) for item in items]
x = 561609
y = 1024183
try:
ds = odc.stac.load(items, bands=["B02"], crs="EPSG:32648", resolution=10, x=(x-80, x+80), y=(y-80, y+80))
print("Success with x/y:", ds.dims)
except Exception as e:
print("Error with x/y:", e)
-18
View File
@@ -1,18 +0,0 @@
import geopandas as gpd
import planetary_computer
import pystac_client
import odc.stac
import numpy as np
bbox = [105.5, 9.2, 106.3, 10.0]
time_range = "2023-01-01/2023-04-30"
catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1", modifier=planetary_computer.sign_inplace)
# Get ALL items
items = list(catalog.search(collections=["sentinel-2-l2a"], bbox=bbox, datetime=time_range, query={"eo:cloud_cover": {"lt": 30}}).items())
items = [planetary_computer.sign(item) for item in items]
print(f"Total items: {len(items)}")
x = 561609
y = 1024183
ds = odc.stac.load(items, bands=["B02"], x=(x-80, x+80), y=(y-80, y+80), crs="EPSG:32648", resolution=10, patch_url=planetary_computer.sign).compute()
print("Time dimension size:", ds.dims['time'])
-20
View File
@@ -1,20 +0,0 @@
import geopandas as gpd
import planetary_computer
import pystac_client
import odc.stac
import numpy as np
bbox = [105.5, 9.2, 106.3, 10.0]
time_range = "2023-01-01/2023-04-30"
catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1", modifier=planetary_computer.sign_inplace)
items = list(catalog.search(collections=["sentinel-2-l2a"], bbox=bbox, datetime=time_range, query={"eo:cloud_cover": {"lt": 30}}).items())
items = [planetary_computer.sign(item) for item in items]
x = 561609
y = 1024183
ds = odc.stac.load(items, bands=["B02"], x=(x-80, x+80), y=(y-80, y+80), crs="EPSG:32648", resolution=10, patch_url=planetary_computer.sign, fail_on_error=False).compute()
print("Original time size:", len(ds.time))
ds2 = ds.dropna(dim="time", how="all")
print("After dropna time size:", len(ds2.time))
print("B02 mean:", np.nanmean(ds2["B02"].values))
print("B02 non-nan count:", np.sum(~np.isnan(ds2["B02"].values)))
-25
View File
@@ -1,25 +0,0 @@
import geopandas as gpd
import planetary_computer
import pystac_client
import odc.stac
import numpy as np
bbox = [105.5, 9.2, 106.3, 10.0]
time_range = "2023-01-01/2023-04-30"
catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1", modifier=planetary_computer.sign_inplace)
items = list(catalog.search(collections=["sentinel-2-l2a"], bbox=bbox, datetime=time_range, query={"eo:cloud_cover": {"lt": 30}}).items())
items = [planetary_computer.sign(item) for item in items]
x = 561609
y = 1024183
ds = odc.stac.load(items, bands=["B02", "B03", "B04", "B08", "SCL"], x=(x-80, x+80), y=(y-80, y+80), crs="EPSG:32648", resolution=10, patch_url=planetary_computer.sign, fail_on_error=False).compute()
print("Original shape:", ds["B02"].shape)
ds2 = ds.dropna(dim="time", how="all")
print("After dropna time size:", len(ds2.time))
if len(ds2.time) > 0:
ds2 = ds2.isel(time=slice(0, 4))
median = ds2["B02"].median(dim="time", skipna=True).values
print("Median shape:", median.shape)
print("Zeros in median:", np.sum(median == 0) / median.size)
print("NaNs in median:", np.sum(np.isnan(median)) / median.size)
-15
View File
@@ -1,15 +0,0 @@
import sys
# Thêm đường dẫn hiện tại vào PYTHONPATH để import được new_import_ODC nếu cần
sys.path.append('.')
import warnings
warnings.filterwarnings('ignore')
from new_import_ODC import load_sen1
print("Testing load_sen1 with a short time range to speed up Dask compute...")
bbox = [105.5, 9.2, 106.4, 10.0]
time_range = "2023-01-01/2023-01-31" # Short time range for fast testing
vh, vv = load_sen1(bbox, time_range)
print("VH shape:", vh.shape)
print("VV shape:", vv.shape)
print("VH CRS:", vh.rio.crs)
print("Success!")
-49
View File
@@ -1,49 +0,0 @@
import warnings
warnings.filterwarnings('ignore')
def load_sen1(bbox, time_range):
import pystac_client
import planetary_computer
import odc.stac
catalog = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace,
)
search = catalog.search(
collections=["sentinel-1-rtc"],
bbox=bbox,
datetime=time_range,
)
items = list(search.items())
print("Found items:", len(items))
ds_s1 = odc.stac.load(
items,
bands=["vv", "vh"],
bbox=bbox,
crs="EPSG:32648",
resolution=10,
chunks={"x": 2048, "y": 2048, "time": 1}
)
ds_median = ds_s1.median(dim="time").compute()
vv = ds_median["vv"]
vh = ds_median["vh"]
vv = vv.expand_dims(dim="band")
vh = vh.expand_dims(dim="band")
vv = vv.rio.write_crs("EPSG:32648")
vh = vh.rio.write_crs("EPSG:32648")
return vh, vv
print("Testing load_sen1...")
bbox = [105.5, 9.2, 106.4, 10.0]
time_range = "2022-09-01/2023-10-01"
vh, vv = load_sen1(bbox, time_range)
print("VH shape:", vh.shape)
print("VV shape:", vv.shape)
print("Success!")
-17
View File
@@ -1,17 +0,0 @@
import geopandas as gpd
import planetary_computer
import pystac_client
import odc.stac
import numpy as np
bbox = [105.5, 9.2, 106.3, 10.0]
time_range = "2023-01-01/2023-04-30"
catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1", modifier=planetary_computer.sign_inplace)
items = list(catalog.search(collections=["sentinel-2-l2a"], bbox=bbox, datetime=time_range, query={"eo:cloud_cover": {"lt": 30}}).items())[:4]
items = [planetary_computer.sign(item) for item in items]
x = 561609
y = 1024183
ds = odc.stac.load(items, bands=["B02", "B03", "B04", "B08"], x=(x-80, x+80), y=(y-80, y+80), crs="EPSG:32648", resolution=10, patch_url=planetary_computer.sign).compute()
print("B04 nanmean:", np.nanmean(ds["B04"].values))
print("B04 nanmax:", np.nanmax(ds["B04"].values))
-41
View File
@@ -1,41 +0,0 @@
import new_import_ODC
importlib = __import__('importlib')
importlib.reload(new_import_ODC)
from new_import_ODC import *
import numpy as np
date_range = ("2022-09-01", "2022-10-01")
longtitude_range = (105.86, 105.94)
latitude_range = (9.65, 9.69)
coordinates = (longtitude_range, latitude_range)
print("Loading S2...")
data = load_data(None, date_range, longtitude_range, latitude_range)
result = mask_clean(data)
ds1 = calculate_indices(result, index="NDVI", satellite_mission="s2")
ndvi = ds1["NDVI"]
time_split = [
slice("2022-09-01", "2023-01-01"),
slice("2023-01-01", "2023-05-01"),
slice("2023-05-01", "2023-07-01"),
slice("2023-07-01", "2022-10-01"),
]
fill_nan_ndvi = fill_nan(ndvi, time_split)
average_ndvi = fill_nan_ndvi.resample(time="1M").mean().compute()
print("Loading S1...")
dsvh, dsvv = load_data_sen1(None, date_range, coordinates)
average_vv = calculate_average(dsvv, time_pattern='1M')
average_vh = calculate_average(dsvh, time_pattern='1M')
train = load_train_data("train/ST_training_data_updated_1130points_new.shp")
point = train.iloc[0]
ndvi_val = average_ndvi.sel(x=point.geometry.x, y=point.geometry.y, method='nearest').values
vh_val = average_vh.sel(x=point.geometry.x, y=point.geometry.y, method='nearest').values
vv_val = average_vv.sel(x=point.geometry.x, y=point.geometry.y, method='nearest').values
print("NDVI shape:", ndvi_val.shape, "ndim:", ndvi_val.ndim)
print("VH shape:", vh_val.shape, "ndim:", vh_val.ndim)
print("VV shape:", vv_val.shape, "ndim:", vv_val.ndim)
-42
View File
@@ -1,42 +0,0 @@
import geopandas as gpd
import planetary_computer
import pystac_client
import odc.stac
import numpy as np
import time
from shapely.geometry import Point, box, shape
bbox = [105.5, 9.2, 106.3, 10.0]
time_range = "2023-01-01/2023-04-30"
catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1", modifier=planetary_computer.sign_inplace)
items = list(catalog.search(collections=["sentinel-2-l2a"], bbox=bbox, datetime=time_range, query={"eo:cloud_cover": {"lt": 30}}).items())
x = 561609
y = 1024183
start = time.time()
# Filter items by spatial intersection
from pyproj import Transformer
# The items geometry are in EPSG:4326 (lon, lat)
# Our x, y are in EPSG:32648
transformer = Transformer.from_crs("epsg:32648", "epsg:4326", always_xy=True)
lon, lat = transformer.transform(x, y)
point = Point(lon, lat)
filtered_items = []
for item in items:
geom = shape(item.geometry)
if geom.contains(point):
filtered_items.append(item)
filtered_items = sorted(filtered_items, key=lambda x: x.properties["eo:cloud_cover"])
print("Original items:", len(items))
print("Filtered items:", len(filtered_items))
print("Time to filter:", time.time() - start)
start = time.time()
filtered_items = [planetary_computer.sign(item) for item in filtered_items]
ds = odc.stac.load(filtered_items[:4], bands=["B02"], x=(x-80, x+80), y=(y-80, y+80), crs="EPSG:32648", resolution=10, patch_url=planetary_computer.sign, fail_on_error=False).compute()
print("Time to load 4 items:", time.time() - start)
print(ds["B02"].shape)
-5
View File
@@ -1,5 +0,0 @@
from train_cloud_removal import UNet
model = UNet(in_channels=6, out_channels=4)
print(hasattr(model, 'inc'))
print(hasattr(model, 'conv1'))
print(list(model.parameters())[0].shape)
-13
View File
@@ -1,13 +0,0 @@
import numpy as np
y = []
# simulate appending 1130 labels
for i in range(1130):
y.append(i % 5)
y = np.array(y)
unique_labels = sorted(list(np.unique(y)))
label_map = {lbl: i for i, lbl in enumerate(unique_labels)}
y_mapped = np.array([label_map[l] for l in y])
print(len(y), len(y_mapped))
+1135
View File
File diff suppressed because it is too large Load Diff
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
+4465
View File
File diff suppressed because one or more lines are too long
+638
View File
@@ -0,0 +1,638 @@
# GEMINI PROJECT CONTEXT - Land Classification & Remote Sensing System
**Last Updated**: March 26, 2026
**Project Location**: `/home/x79/remote-sensing`
**Purpose**: Complete land classification and environmental monitoring system using satellite remote sensing for Vietnam
---
## 📋 PROJECT OVERVIEW
### High-Level Purpose & Problem Domain
- **Core Task**: Classify land use/land cover (8 land classes) in Vietnam using multispectral Sentinel-2 and radar Sentinel-1 data from Microsoft Planetary Computer
- **Geographic Focus**: Vietnam provinces/regions with bounding-box (bbox) based Area-of-Interest (AOI) selection
- **Key Capabilities**:
- Dynamic training with user-selected regions and time periods
- Pixel-wise inference (prediction) on new regions
- Cloud removal using 7 different strategies
- NDVI time-series forecasting and change detection workflows
- Auto-generated HTML reports with visualizations
- Batch processing of multiple regions
- Model lifecycle management (save, load, validate, delete)
### Data Pipeline
```
Sentinel-2 (optical) + Sentinel-1 (SAR)
[Feature Extraction: 4 modes - simple (3) / temporal (39) / extended (15) / odc (8)]
[Model Training: XGBoost, RF, SVM, CNN, Swin-UNet, MobileNet-LRASPP]
[Prediction: Pixel-wise classification]
[Output: GeoTIFF + PNG preview + HTML report + JSON metadata]
```
---
## 🏗️ SYSTEM ARCHITECTURE
### Core Technology Stack
- **Backend**: FastAPI (~4200 lines in `api_server.py`)
- **ML Training**: scikit-learn (XGBoost, RF, SVM, DT) + PyTorch (CNN, Swin-UNet, MobileNet)
- **Geospatial**: rasterio, rioxarray, geopandas, xarray, odc.stac
- **Data Access**: Microsoft Planetary Computer STAC API (Sentinel-2 L2A, Sentinel-1 RTC)
- **Frontend**: HTML + Leaflet.js (map drawing) + Fetch API + Chart.js
- **GPU Support**: PyTorch with CUDA 12.x (optional fallback to CPU)
### Folder Structure
```
remote-sensing/
├── Core Backend
│ ├── api_server.py # FastAPI app (~4200 LOC, 70+ endpoints)
│ ├── train_module.py # Training pipeline engine
│ ├── feature_extractor.py # Unified feature extraction (4 modes)
│ ├── model_manager.py # Model lifecycle management
│ ├── cloud_removal.py # 7 cloud removal strategies
│ ├── report_generator.py # Auto HTML/PNG report generation
│ ├── generate_previews.py # GeoTIFF → PNG conversion
│ │
├── Utilities & Lookup
│ ├── vietnam_provinces.py # Province bboxes & metadata
│ ├── vietnam_provinces_merged.py # 32-province variant
│ ├── utils.py # Geospatial helper functions
│ ├── create_odc_metadata.py # Metadata generator utility
│ │
├── Frontend Pages (HTML)
│ ├── index.html # Main dashboard hub
│ ├── training_interface.html # Training UI
│ ├── prediction_interface.html # Prediction UI
│ ├── batch_interface.html # Batch processing UI
│ ├── ndvi_interface.html # NDVI time-series UI
│ ├── dashboard.html # Analytics dashboard
│ ├── reports_interface.html # Reports management
│ ├── change_detection_interface.html # Change detection UI
│ ├── cloud_training_interface.html # Cloud removal training UI
│ │
├── Tests & Notebooks
│ ├── test_*.py # Unit & integration tests
│ ├── 01.train_ODC*.ipynb # Training notebooks
│ ├── 02.predict_ODC.ipynb # Prediction notebooks
│ ├── cloud_removal_train.ipynb # Cloud removal training
│ │
├── Model Storage & Caches
│ ├── model_train/ # Trained models (*.joblib, *.pth)
│ │ ├── model_odc.joblib # Legacy GridSearchCV model
│ │ ├── model_*_info.json # Metadata sidecar files
│ ├── cloud_removal_model/ # Cloud removal U-Net models (.pth)
│ ├── predictions/ # Prediction output (GeoTIFF + PNG)
│ ├── reports/ # Generated HTML reports
│ ├── dataset_cache/ # Cached Sentinel data (optional)
│ │
├── Config & Documentation
│ ├── requirement.txt # Python dependencies
│ ├── requirements_api.txt # API-specific deps
│ ├── IMPLEMENTATION_SUMMARY.md # Model manager summary
│ ├── MODEL_MANAGER_GUIDE.md # Full model management guide
│ ├── NDVI_FORECAST_METHODOLOGY.md # NDVI algorithm docs
│ ├── CLOUD_TRAINING_GUIDE.md # Cloud removal training guide
│ └── [Other guides & docs]
```
---
## 🔧 MAIN MODULES & RESPONSIBILITIES
| **Module** | **File(s)** | **Key Responsibility** |
|---|---|---|
| **API Server** | `api_server.py` | FastAPI app with 70+ endpoints; routes all training, prediction, batch, cloud removal, dashboard, reports, model management tasks |
| **Training Engine** | `train_module.py` | Complete training pipeline: fetch data → feature extraction → train/test split → model training → evaluation → save with metadata |
| **Feature Extraction** | `feature_extractor.py` | Standardized feature extraction with 4 modes: simple, temporal, extended, odc; used by both training and prediction |
| **Model Manager** | `model_manager.py` | Lifecycle management: list, load, save, validate, delete models; handles metadata JSON; auto-detects CNN/PyTorch models |
| **Cloud Removal** | `cloud_removal.py` | 7 cloud removal strategies: classic (3-step), temporal_only, median_composite, none, speckle filter, ML inpainting, deep learning U-Net |
| **Report Generator** | `report_generator.py` | Auto-generates HTML/PNG reports with confusion matrices, class distributions, accuracy trends |
| **Preview Generator** | `generate_previews.py` | Converts GeoTIFF outputs to PNG previews (NDVI or classification rasters) |
| **Province Lookup** | `vietnam_provinces*.py` | Lookup tables for 32+ Vietnamese provinces with bboxes and region grouping |
| **Utilities** | `utils.py` | Geospatial helper functions (load GeoDataFrames, etc.) |
---
## 📊 END-TO-END WORKFLOWS
### 1. TRAINING WORKFLOW
```
User Input → Training Configuration
API Endpoint: POST /api/training/start
train_module.py: train_model()
1. Fetch Sentinel-2 & Sentinel-1 from Planetary Computer STAC
2. Apply cloud mask (SCL band: clouds, shadows, cirrus masked)
3. Extract features via FeatureExtractor (mode: simple/temporal/extended/odc)
4. Train/test split (default 0.2)
5. Train selected model type (XGBoost, RF, CNN, Swin-UNet, MobileNet)
6. Evaluate: accuracy, precision, recall, F1, confusion matrix
model_manager.py: Save model + JSON metadata
report_generator.py: Auto-generate HTML training report
Return: {model_filename, accuracy_metrics, training_time}
```
**Key Metadata Saved**:
```json
{
"timestamp": "2026-03-26T14:30:00",
"model_type": "xgboost",
"feature_mode": "temporal",
"n_features": 39,
"n_classes": 8,
"features": ["NDVI_t1", "NDVI_t2", ..., "NDWI_t1", ...],
"test_accuracy": 0.85,
"train_accuracy": 0.92,
"bbox": [105.6, 9.3, 106.2, 9.8],
"time_range": "2023-03-01/2023-05-31",
"resolution": 20,
"data_source": "Microsoft Planetary Computer STAC"
}
```
### 2. PREDICTION WORKFLOW
```
User Input → Prediction Configuration (model_filename, bbox, time_range, cloud_strategy)
API Endpoint: POST /api/predict or POST /api/predict/with-ndvi
run_prediction() function:
1. Load model via model_manager.py (retrieves metadata, feature requirements)
2. Fetch Sentinel-2 & Sentinel-1 for new region
3. Apply chosen cloud_removal_method (classic/temporal_only/median_composite/none/deep_learning)
4. Extract features matching model's metadata requirements
5. Auto-adjust if feature count mismatch (pad/trim)
6. Predict land class for each pixel
7. (Optional) Calculate NDVI: (NIR - Red) / (NIR + Red)
8. Save outputs: GeoTIFF + PNG preview
generate_previews.py: Create PNG from GeoTIFF
report_generator.py: Generate prediction report
Return: {prediction_file, ndvi_file, class_distribution, statistics}
```
### 3. BATCH PROCESSING WORKFLOW
```
User uploads CSV with multiple regions:
(name, min_lon, min_lat, max_lon, max_lat, start_date, end_date, max_scenes, cloud_cover, resolution)
API Endpoint: POST /api/batch/start
Enqueue all regions; process sequentially
For each region: Run same prediction workflow
Track status per region: Queued → Running → Completed/Failed
UI shows progress bar, auto-retry on failure (max 3 retries)
Return: Bulk results with per-region status & output files
```
### 4. CLOUD REMOVAL WORKFLOW
```
User selects cloud_removal_method in prediction config:
cloud_removal.py: process_cloud_removal()
Strategy Selection:
• 'classic': temporal interpolation → median composite → spatial interpolation (3-step)
• 'temporal_only': ffill + bfill across time dimension (fast, good for many scenes)
• 'median_composite': Prioritize median across scenes (best for noise reduction)
• 'none': Keep original, just fill NaN with 0
• 'deep': Use trained U-Net model (S2 cloudy + S1 → clean S2)
• 'ml_inpainting': KNN or Random Forest based inpainting
• 'speckle_filter': Reduce radar noise
Return cleaned Sentinel-2 data for subsequent feature extraction
```
### 5. NDVI TIME-SERIES WORKFLOW
```
User requests NDVI calculation (bbox + time_range + aggregation)
API Endpoint: POST /api/ndvi/timeseries or /api/ndvi/predict-timeseries
Load Sentinel-2 (B04 Red, B08 NIR)
Calculate NDVI = (NIR - Red) / (NIR + Red + 0.00001)
Resample to monthly or user-defined aggregation
Export as GeoTIFF + PNG visualization
Show time-series graph & statistics (mean, min, max, std, trend)
```
### 6. CHANGE DETECTION WORKFLOW
```
User selects: model + current_period + prediction_period
API Endpoint: POST /api/change-detection/compare-periods
Run prediction for both time periods
Compute difference map (current - prediction)
Classify changes: increased vegetation, decreased vegetation, stable
Generate change map GeoTIFF + report with statistics
```
---
## 🌐 API ENDPOINTS SUMMARY (70+ endpoints)
### Model Management
- `GET /api/models/list` - List all trained models with metadata
- `GET /api/models/{filename}/info` - Get model details
- `GET /api/models/{filename}/validate` - Validate model integrity
- `DELETE /api/models/{filename}` - Delete model file
### Training APIs
- `POST /api/training/start` - Start land classification training
- `GET /api/training/status` - Get training progress
- `POST /api/training/stop` - Cancel ongoing training
- `POST /api/cloud-removal/train` - Train cloud removal U-Net
### Prediction APIs
- `POST /api/predict` - Standard prediction (classification only)
- `POST /api/predict/with-ndvi` - Prediction with NDVI export
- `POST /api/change-detection/compare-periods` - Change detection
- `GET /api/prediction/status` - Check prediction progress
- `GET /api/predictions/list` - List prediction outputs
- `GET /api/predictions/download/{filename}` - Download prediction file
- `GET /api/predictions/preview/{filename}` - View PNG preview
### Batch Processing
- `POST /api/batch/start` - Enqueue multiple predictions from CSV
- `GET /api/batch/status` - Check batch queue
- `GET /api/batch/results/{batch_id}` - Retrieve batch results
- `POST /api/batch/cancel/{batch_id}` - Cancel batch job
### Cloud Removal
- `GET /api/cloud-removal/methods` - List available strategies
- `GET /api/cloud-removal/models` - List trained .pth models
- `POST /api/cloud-removal/upload` - Upload .pth cloud removal model
- `DELETE /api/cloud-removal/models/{filename}` - Delete cloud removal model
### Dashboard & Reports
- `GET /api/dashboard/statistics` - Overall system stats
- `GET /api/dashboard/accuracy-trends` - Accuracy over time
- `GET /api/dashboard/class-distribution/{model_filename}` - Class distribution
- `GET /api/reports/list` - List generated reports
- `GET /api/reports/view/{filename}` - View HTML report
- `GET /api/reports/download/{filename}` - Download report
- `DELETE /api/reports/delete/{filename}` - Delete report
### Provinces & Utilities
- `GET /api/provinces/list` - List all Vietnamese provinces
- `GET /api/provinces/by-region` - Group provinces by region
- `GET /api/provinces/{province_name}/bbox` - Get province bbox
- `GET /api/provinces/search/{query}` - Search province by name
- `GET /api/provinces-32/*` - Alternative 32-province variant
- `GET /api/network/check` - Check connectivity to Planetary Computer
- `GET /api/cache/info` - Show cache statistics
- `POST /api/cache/clear` - Clear local cache
### NDVI & Time-Series
- `POST /api/ndvi/timeseries` - Calculate NDVI time-series
- `POST /api/ndvi/predict-timeseries` - NDVI prediction/forecast
- `POST /api/ndvi/forecast` - NDVI forecasting
### File Management
- `GET /api/training/files` - List training files
- `GET /api/overlay/shapefiles` - List available shapefiles
- `GET /api/training/shapefile/{filename}/labels` - Get shapefile labels
- `POST /api/land-classification/upload` - Upload custom model
- `POST /api/cloud-removal/upload` - Upload cloud removal model
### Frontend Routes (Serve HTML)
- `GET /` - Main dashboard
- `GET /training` - Training interface
- `GET /prediction` - Prediction interface
- `GET /dashboard` - Analytics dashboard
- `GET /batch` - Batch processing UI
- `GET /ndvi` - NDVI time-series UI
- `GET /reports` - Reports management
- `GET /cloud-training` - Cloud removal training
- `GET /change-detection` - Change detection UI
---
## 💾 DATA INPUTS / OUTPUTS & FOLDER CONVENTIONS
### Input Data Sources
- **Sentinel-2 L2A** from Microsoft Planetary Computer STAC API
- Bands: B02 (blue), B03 (green), B04 (red), B08 (NIR), B11 (SWIR), SCL (cloud mask)
- Resolution: 10m or 20m (user selectable)
- Collection: `sentinel-2-l2a`
- **Sentinel-1 RTC** from Planetary Computer
- Bands: VH, VV (radar polarizations)
- Converted to dB scale: `10 * log10(intensity)`
- Collection: `sentinel-1-rtc`
- **Training Labels**: User-provided shapefiles with pixel-level class labels
### Output File Structure
```
predictions/
├── prediction_YYYYMMDD_HHMMSS.tif # Classification GeoTIFF
├── prediction_YYYYMMDD_HHMMSS.png # PNG preview
├── ndvi_YYYYMMDD_HHMMSS.tif # NDVI raster
├── ndvi_YYYYMMDD_HHMMSS.png # NDVI preview
reports/
├── training_report_*.html # Auto training reports
├── prediction_report_*.html # Auto prediction reports
model_train/
├── model_odc.joblib # Legacy model
├── model_odc_info.json # Metadata
├── model_xgboost_*.joblib # XGBoost models
├── model_xgboost_*_info.json # Metadata
├── model_cnn_*.joblib # CNN models
├── model_cnn_*_info.json # Metadata
cloud_removal_model/
├── cloud_removal_unet_best.pth # Trained U-Net
├── *.pth # Custom models
├── *.json # Model metadata
```
---
## 🔌 EXTERNAL DEPENDENCIES & PLATFORMS
### Critical External Services
- **Microsoft Planetary Computer** (STAC API)
- Hosts Sentinel-2 L2A and Sentinel-1 RTC archives
- URL: `https://planetarycomputer.microsoft.com/api/stac/v1`
- Auto-signed access tokens via `planetary_computer.sign_inplace`
- Network connectivity check: `GET /api/network/check`
### Key Python Libraries
- **Geospatial**: rasterio, rioxarray, geopandas, shapely, Cartopy, folium, ipyleaflet
- **Data Processing**: numpy, pandas, xarray, dask
- **ML**: scikit-learn, xgboost
- **Deep Learning**: torch, torchvision
- **Web**: fastapi, uvicorn, pydantic
- **Visualization**: matplotlib, Pillow (PIL)
- **Document Gen**: markdown, Pillow
### GPU Support
- PyTorch with CUDA 12.x (optional; falls back to CPU)
- Benefits Swin-UNet and CNN models (10-100x speedup)
- CPU training for XGBoost/RF typically <1 hour; deep models need GPU for reasonable speed
---
## ⚙️ FEATURE EXTRACTION MODES (CRITICAL)
Train and prediction **MUST** use same feature mode and dimension; metadata auto-detects this.
| Mode | # Features | Description | Best For | Training Time |
|---|---|---|---|---|
| **simple** | 3 | NDVI_mean, VH_db_mean, VV_db_mean | Fast iteration, baseline | ~5-10 min |
| **temporal** | 39 | NDVI/NDWI/NDBI across 13 months + radar stats | High accuracy (~85%+) | ~30-60 min |
| **extended** | 15 | NDVI/NDWI/NDBI stats (mean/std/min/max) + radar | Balanced speed/accuracy | ~15-30 min |
| **odc** | 8 | NDVI stats + NDWI/NDBI/EVI mean (legacy ODC mode) | Legacy compatibility | ~10-20 min |
**Critical**: If feature mode = "temporal" (39 features) at training, prediction MUST extract 39 features. System auto-detects from metadata but will fail if mismatched.
---
## 🎯 OPERATIONAL NOTES & CONSTRAINTS
### Performance Limits
1. **Planetary Computer Timeout Issues**
- Large bbox (>10km × 10km) + long time range (>1 month) + high max_scenes → timeouts
- **Solution**: Progressive loading (subdivide bbox), reduce time window, reduce max_scenes
- **Safe Settings**: bbox ≤ 10km × 10km, time ≤ 1 month, max_scenes ≤ 12
2. **Memory Usage**
- Temporal mode (39 features) requires ~2-3x RAM vs simple mode
- Large regions: reduce resolution (10m → 20m) or split into sub-tiles
- Batch processing: sequential (one region at a time due to API limits)
3. **GPU Training**
- Swin-UNet: ~15-60 min on GPU vs ~2-4 hours on CPU
- CNN: ~10-30 min on GPU vs ~1-2 hours on CPU
- XGBoost/RF: CPU-bound; GPU not beneficial
### Data Quality Issues
1. **Cloud Cover**
- SCL band values: 3=cloud shadow, 8=cloud medium, 9=cloud high, 10=cirrus
- Recommend multiple scenes (≥5) for temporal aggregation
- Cloud removal strategy critical—test different approaches
2. **Radar Data (Sentinel-1)**
- Not always available for all regions/dates
- System gracefully falls back to zeros if unavailable
- Safe for "extended" & "odc" modes that have radar fallback
3. **Feature Mode Mismatch**
- Model trained with "temporal" (39 features) needs 39-dim input
- System auto-adjusts (pads/trims) from metadata but may degrade accuracy
- **Best Practice**: Align feature mode explicitly; don't mix
### Known Caveats
1. **Legacy Model (model_odc.joblib)**: Hardcoded 39 temporal features; auto-detected via `model_odc_info.json`
2. **Metadata Consistency**: Old models may lack `.json` sidecar; system generates default (may be incorrect)
3. **Batch Processing**: Sequential only; large batches (100+ regions) take hours
4. **Change Detection**: Simple differencing approach; requires same model & feature mode for both periods
5. **Rate Limiting**: Planetary Computer may rate-limit if too many concurrent requests
### Recommended Best Practices
- Test model on small bbox first (2km × 2km, 1 week, 3 scenes)
- Use "simple" mode for fast iteration, "temporal" for best accuracy (85%+)
- Store metadata JSON alongside model file (sidecar pattern)
- Version control: record feature_mode & n_features in every training
- Monitor training accuracy; retrain if <70% accuracy
- Cache Sentinel data locally to avoid repeated downloads
- Use "median_composite" cloud strategy if >5 scenes; "temporal_only" if 3-4 scenes
---
## 🔍 TEST COVERAGE MAP
| Test File | Coverage | Status |
|---|---|---|
| `test_model_manager.py` | ModelManager lifecycle (list, load, validate) | ✅ Well-tested |
| `test_feature_extractor.py` | All 4 feature extraction modes | ✅ Well-tested |
| `test_training_api.py` | Training API endpoints | ✅ Partial |
| `test_cloud_removal.py` | 7 cloud removal strategies | ✅ Well-tested |
| `test_cloud_training.py` | U-Net cloud removal training | ✅ Partial |
| `test_shapefile_api.py` | Shapefile overlay feature | ✅ Partial |
| `test_planetary_computer.py` | Planetary Computer STAC access | ✅ Well-tested |
| `test_new_features.py` | Recent feature releases | ✅ Partial |
| Jupyter Notebooks | Training & prediction workflows | ✅ Mix of unit/integration/notebooks |
**Coverage Notes**: Model management, feature extraction, and cloud removal well-tested; Dashboard UI, change detection, NDVI time-series mostly tested via notebooks.
---
## 📚 FILE REFERENCE MAP
### Core Execution
- `api_server.py` — Main FastAPI application (~4200 LOC)
- `train_module.py` — Training logic (data fetch → feature extraction → training)
- `run_prediction_new.py` — Prediction execution function
- `feature_extractor.py` — Unified feature extraction (4 modes)
- `model_manager.py` — Model lifecycle (load/save/validate/list)
- `cloud_removal.py` — Cloud removal strategies (7 methods)
- `report_generator.py` — HTML/PNG report auto-generation
- `generate_previews.py` — GeoTIFF → PNG conversion
### Data & Config
- `vietnam_provinces.py` — 32+ province lookup tables & bboxes
- `vietnam_provinces_merged.py` — Alternative 32-province variant
- `utils.py` — Geospatial utility functions
- `create_odc_metadata.py` — Legacy metadata generator
### Frontend
- `index.html` — Main dashboard hub (tab navigation)
- `training_interface.html` — Training configuration UI
- `prediction_interface.html` — Prediction configuration UI
- `batch_interface.html` — Batch processing (CSV upload)
- `ndvi_interface.html` — NDVI time-series visualization
- `dashboard.html` — Analytics & model performance dashboard
- `reports_interface.html` — Report management & viewing
- `change_detection_interface.html` — Change detection visualization
- `cloud_training_interface.html` — Cloud removal U-Net training
### Documentation
- `IMPLEMENTATION_SUMMARY.md` — Model manager & system overview
- `MODEL_MANAGER_GUIDE.md` — Complete model management guide
- `NDVI_FORECAST_METHODOLOGY.md` — NDVI algorithm documentation
- `CLOUD_TRAINING_GUIDE.md` — Cloud removal training guide
- `NDVI_PREDICTION_GUIDE.md` — NDVI prediction workflow
- `CLOUD_PROCESSING.md` — Cloud processing notes
- `UPDATE_SUMMARY.md` — Recent updates & features
---
## 🚀 BOOTSTRAP PROMPT FOR GEMINI
### System Context (Copy & Paste for Gemini)
```
You are assisting a remote-sensing land-classification project for Vietnam.
## ARCHITECTURE SNAPSHOT
- **Backend**: FastAPI (~4200 LOC, 70+ endpoints) for orchestrating training, prediction, batch, cloud removal, reporting
- **Data Source**: Microsoft Planetary Computer STAC API (Sentinel-2 L2A + Sentinel-1 RTC)
- **Training**: scikit-learn (XGBoost/RF/SVM/DT) + PyTorch (CNN/Swin-UNet/MobileNet)
- **Feature Extraction**: 4 modes (simple 3-feat / temporal 39-feat / extended 15-feat / odc 8-feat)
- **Cloud Removal**: 7 strategies (classic, temporal_only, median_composite, none, ML inpainting, deep U-Net)
- **Output**: GeoTIFF + PNG + HTML report + JSON metadata
## CORE FILES TO UNDERSTAND (Priority Order)
1. api_server.py — Main API server (training, prediction, batch, models, reports)
2. train_module.py — Training pipeline (data fetch → feature extraction → train → save)
3. feature_extractor.py — Unified feature extraction with auto mode detection
4. model_manager.py — Model lifecycle (load/save/validate/list)
5. cloud_removal.py — Cloud removal strategies (7 methods)
6. report_generator.py — Auto-generate HTML reports
7. run_prediction_new.py — Prediction execution
8. vietnam_provinces.py — Province lookup & bbox tables
## CRITICAL CONSTRAINTS & GOTCHAS
1. **Feature Mode Consistency**: Training & prediction MUST use same mode (simple/temporal/extended/odc)
→ Auto-detected from metadata JSON
→ Mismatch causes dimension error or accuracy degradation
2. **Planetary Computer Limits**:
→ Timeout if bbox >10km×10km OR time range >1 month OR max_scenes >12
→ Solution: subdivide bbox, reduce time window, limit scenes
3. **Cloud Strategy Selection**:
→ ≥5 scenes → use "median_composite" (best noise reduction)
→ 3-4 scenes → use "temporal_only" (fast temporal interp)
→ <3 scenes → use "none" (skip cloud removal)
4. **Radar Data Fallback**:
→ Sentinel-1 may be unavailable for some regions
→ System gracefully falls back to zeros (safe for all modes)
5. **Model Metadata**:
→ Always stored as `model_name_info.json` sidecar file
→ Contains: n_features, feature_mode, features list, accuracy, bbox, time_range
→ Missing metadata → system uses defaults (may be incorrect)
6. **Legacy Model (model_odc.joblib)**:
→ Hardcoded 39 temporal features
→ Metadata in model_odc_info.json
## REASONING CHECKLIST (before answering)
□ Is feature_mode consistent between train and prediction?
□ Is metadata.json present and correct?
□ Does bbox exceed 10km×10km? (Planetary Computer timeout risk)
□ Is cloud_removal_strategy appropriate for # of scenes?
□ Is Sentinel-1 data available for this region/date?
□ Is model a joblib (scikit-learn) or .pth (PyTorch) file?
□ Is GPU available for deep models (CNN, Swin-UNet)?
□ Does memory allow temporal feature extraction (39-feat)?
## RESPONSE FORMAT
- Always cite api_server.py endpoint, function name, or module being discussed
- Verify feature_mode & n_features from metadata JSON
- Suggest cloud_removal_strategy based on # of scenes available
- For unknown issues: offer alternative approaches (reduce bbox, cache results, use simpler model)
- Explain reasoning using checklist above
## DATA FLOW SUMMARY
Sentinel-2/S1 → [Cloud Remove] → [Feature Extract] → [Train/Predict] → [GeoTIFF + PNG + Report]
```
---
## 📞 QUICK REFERENCE CHECKLIST
### Before Troubleshooting Any Issue
- [ ] Check feature_mode consistency (metadata JSON)
- [ ] Verify metadata.json exists for the model
- [ ] Check Planetary Computer connectivity (`GET /api/network/check`)
- [ ] Review cloud_removal_method choice (≥5 scenes = median_composite)
- [ ] Confirm Sentinel-1 availability (or fallback to zeros if missing)
- [ ] Validate bbox size (≤10km×10km for safety)
- [ ] Check memory usage for temporal feature mode
- [ ] Verify GPU if using CNN/Swin-UNet models
### Common Issues & Solutions
| Issue | Likely Cause | Solution |
|---|---|---|
| Training timeout | Large bbox / long time / many scenes | Subdivide bbox, reduce time window, max_scenes ≤ 12 |
| Feature dimension mismatch | Different feature_mode between train & predict | Check metadata.json, ensure same mode |
| Low prediction accuracy | Cloud cover, poor training data, feature mode too simple | Use "temporal" mode, increase training data, try cloud removal |
| Out of memory | Temporal features + large region | Reduce resolution (20m), split into sub-tiles, increase RAM |
| Model not found | Wrong filename or model_train/ path issue | `GET /api/models/list` to verify, check file path |
| Planetary Computer error | Network issue or API rate limit | Check DNS, retry later, reduce concurrent requests |
| Cloud removal failing | Strategy not suitable for scene count | Try "none" or "median_composite" depending on scenes |
---
## 🎓 LEARNING RESOURCES IN REPO
- **Notebooks**: `01.train_ODC.ipynb`, `02.predict_ODC.ipynb`, `cloud_removal_train.ipynb`
- **Tests**: `test_*.py` files for unit test patterns
- **Docs**: All `*.md` files for detailed guides and methodology
- **Code Comments**: API server and modules heavily commented
---
**Generated**: March 26, 2026
**For Use By**: Gemini, Claude, GPT, or any AI system needing project context
**Maintainer**: Remote-Sensing Project Team
+296 -57
View File
@@ -26,20 +26,20 @@ import socket
import urllib.request
# Import report generator
from scripts.inference.report_generator import generate_training_report, generate_prediction_report
from report_generator import generate_training_report, generate_prediction_report
# Import Model Manager
from core.model_manager import ModelManager, get_model_manager
from model_manager import ModelManager, get_model_manager
# Import Vietnam provinces data
from core.vietnam_provinces import get_all_provinces, get_provinces_by_region, get_province_bbox, search_province
from core.vietnam_provinces_merged import (
from vietnam_provinces import get_all_provinces, get_provinces_by_region, get_province_bbox, search_province
from vietnam_provinces_merged import (
get_all_provinces_32, get_provinces_by_region_32, get_province_bbox_32,
search_province_32, get_merged_info, get_provinces_statistics
)
# Import cloud removal module
from core.cloud_removal import process_cloud_removal, get_available_methods
from cloud_removal import process_cloud_removal, get_available_methods
# Import planetary computer libraries (conditional)
try:
@@ -225,6 +225,9 @@ class PredictionConfig(BaseModel):
cloud_removal_method: str = "classic"
cloud_removal_model: Optional[str] = None # Optional: .pth model filename for deep learning cloud removal
# Shapefile overlay for visualization
shapefile_overlay: Optional[str] = None # Path to shapefile for overlaying boundaries
class TrainingStatus(BaseModel):
"""Trạng thái training"""
@@ -286,6 +289,7 @@ class PredictionWithNDVIConfig(BaseModel):
export_classification: bool = True # Export classification raster
cloud_removal_method: str = "classic"
cloud_removal_model: Optional[str] = None # Optional: .pth model filename for deep learning cloud removal
shapefile_overlay: Optional[str] = None # Path to shapefile for overlaying boundaries
class CloudRemovalTrainingConfig(BaseModel):
@@ -507,57 +511,6 @@ async def get_cloud_removal_methods():
}
@app.get("/api/ndvi-forecast/models")
async def list_ndvi_forecast_models():
"""Liệt kê các NDVI forecast models đã train"""
model_dir = Path("ndvi_forecast_model")
if not model_dir.exists():
return {"models": [], "count": 0}
models = []
# Search for all models
for model_file in list(model_dir.rglob("*.pth")) + list(model_dir.rglob("*.joblib")):
try:
import json
# Try to load metadata from .json sidecar file first
metadata_file = model_file.with_name(model_file.stem + "_info.json")
if metadata_file.exists():
try:
with open(metadata_file, 'r') as f:
metadata = json.load(f)
models.append({
"filename": model_file.name,
"path": str(model_file),
"model_type": metadata.get('model_type', 'Unknown'),
"target": metadata.get('target', 'NDVI'),
"rmse": metadata.get('rmse', 0),
"mae": metadata.get('mae', 0),
"epoch": metadata.get('epoch', 0),
"created": model_file.stat().st_mtime,
"size_mb": model_file.stat().st_size / (1024 * 1024),
})
continue
except Exception as e:
print(f"Error reading JSON {metadata_file}: {e}")
# Fallback for models without metadata
models.append({
"filename": model_file.name,
"path": str(model_file),
"model_type": "Unknown",
"created": model_file.stat().st_mtime,
"size_mb": model_file.stat().st_size / (1024 * 1024)
})
except Exception as e:
print(f"Error loading model info for {model_file}: {e}")
# Sort by creation time (newest first)
models.sort(key=lambda x: x['created'], reverse=True)
return {"models": models, "count": len(models)}
@app.get("/api/cloud-removal/models")
async def list_cloud_removal_models():
"""Liệt kê các cloud removal models đã train"""
@@ -1279,6 +1232,87 @@ async def list_training_files():
}
@app.get("/api/overlay/shapefiles")
async def list_overlay_shapefiles():
"""Liệt kê các shapefile có sẵn cho overlay trên prediction"""
overlay_dirs = ["region", "ChauThanh", "ThuanHoa"]
shapefiles = []
for overlay_dir in overlay_dirs:
dir_path = Path(overlay_dir)
if not dir_path.exists():
continue
# Find all .shp files in this directory and subdirectories
for shp_file in dir_path.rglob("*.shp"):
try:
file_size = shp_file.stat().st_size
file_modified = datetime.fromtimestamp(shp_file.stat().st_mtime).isoformat()
# Try to read shapefile to get feature count and bbox
import geopandas as gpd
gdf = gpd.read_file(str(shp_file))
feature_count = len(gdf)
# Calculate bbox (always in EPSG:4326 for consistency)
bbox = None
if not gdf.empty and gdf.crs:
try:
# Reproject to EPSG:4326 if needed
if gdf.crs != "EPSG:4326":
gdf_4326 = gdf.to_crs("EPSG:4326")
else:
gdf_4326 = gdf
# Get total bounds [minx, miny, maxx, maxy]
bounds = gdf_4326.total_bounds
if len(bounds) == 4:
bbox = [
float(bounds[0]), # min_lon
float(bounds[1]), # min_lat
float(bounds[2]), # max_lon
float(bounds[3]) # max_lat
]
except Exception as bbox_error:
print(f"[WARNING] Cannot calculate bbox for {shp_file}: {bbox_error}")
# Get relative path from workspace root
relative_path = str(shp_file)
shapefiles.append({
"filename": shp_file.name,
"path": relative_path,
"directory": overlay_dir,
"size_bytes": file_size,
"size_mb": round(file_size / 1024 / 1024, 2),
"modified": file_modified,
"feature_count": feature_count,
"crs": str(gdf.crs) if gdf.crs else "Unknown",
"bbox": bbox # [min_lon, min_lat, max_lon, max_lat] in EPSG:4326
})
except Exception as e:
# If cannot read shapefile, just add basic info
file_size = shp_file.stat().st_size
file_modified = datetime.fromtimestamp(shp_file.stat().st_mtime).isoformat()
relative_path = str(shp_file)
shapefiles.append({
"filename": shp_file.name,
"path": relative_path,
"directory": overlay_dir,
"size_bytes": file_size,
"size_mb": round(file_size / 1024 / 1024, 2),
"modified": file_modified,
"error": f"Cannot read shapefile: {str(e)}"
})
return {
"shapefiles": shapefiles,
"count": len(shapefiles),
"directories": overlay_dirs
}
@app.get("/api/training/shapefile/{filename}/labels")
async def get_shapefile_labels(filename: str):
"""Lấy các label từ một shapefile cụ thể"""
@@ -1780,6 +1814,85 @@ def update_progress(message: str):
print(f"[PROGRESS] {message}")
def rasterize_shapefile_overlay(shapefile_path, reference_raster, boundary_value=255):
"""
Rasterize shapefile boundaries to overlay on prediction result.
Args:
shapefile_path: Path to shapefile
reference_raster: xarray DataArray to match dimensions and CRS
boundary_value: Value to use for boundaries (default 255 for white)
Returns:
numpy array with boundaries, same shape as reference_raster
"""
try:
import geopandas as gpd
from rasterio.features import rasterize
import numpy as np
# Read shapefile
gdf = gpd.read_file(shapefile_path)
print(f"[OVERLAY] Loaded shapefile with {len(gdf)} features, CRS: {gdf.crs}")
# Ensure CRS matches
target_crs = reference_raster.rio.crs
if gdf.crs != target_crs:
print(f"[OVERLAY] Reprojecting from {gdf.crs} to {target_crs}")
gdf = gdf.to_crs(target_crs)
# Get raster dimensions and transform
height, width = reference_raster.shape
transform = reference_raster.rio.transform()
print(f"[OVERLAY] Raster dimensions: {height}x{width}")
print(f"[OVERLAY] Transform: {transform}")
# Calculate appropriate buffer size based on pixel resolution
# Get pixel size from transform (transform[0] is x resolution)
pixel_size = abs(transform[0]) # in CRS units
# Very thin boundary - only 0.2 pixels wide for 1px line
buffer_distance = pixel_size * 0.2
print(f"[OVERLAY] Pixel size: {pixel_size}, Buffer distance: {buffer_distance} (thin 1px line)")
# Create boundary geometries with minimal buffering
boundary_geoms = []
for idx, geom in enumerate(gdf.geometry):
if geom is not None and geom.is_valid:
# Get boundary of each polygon
boundary = geom.boundary
if boundary is not None:
# Minimal buffer for 1-pixel thin line
buffered = boundary.buffer(buffer_distance)
boundary_geoms.append((buffered, boundary_value))
if not boundary_geoms:
print(f"[WARNING] No valid boundary geometries found in {shapefile_path}")
return np.zeros((height, width), dtype=np.uint8)
print(f"[OVERLAY] Rasterizing {len(boundary_geoms)} boundaries...")
# Rasterize boundaries
boundary_mask = rasterize(
shapes=boundary_geoms,
out_shape=(height, width),
transform=transform,
fill=0, # Background
dtype=np.uint8
)
boundary_count = np.count_nonzero(boundary_mask)
print(f"[OVERLAY] Boundary pixels: {boundary_count} / {height*width} ({boundary_count/(height*width)*100:.2f}%)")
if boundary_count == 0:
print(f"[OVERLAY WARNING] No boundary pixels were rasterized! Check CRS and geometry overlap.")
return boundary_mask
except Exception as e:
print(f"[ERROR] Failed to rasterize shapefile {shapefile_path}: {e}")
return None
def update_prediction_progress(message: str):
"""Cập nhật prediction progress message"""
global prediction_status
@@ -2134,6 +2247,33 @@ async def run_prediction(config: PredictionConfig):
prediction_da.rio.to_raster(str(output_file), driver="GTiff")
# ============ SHAPEFILE OVERLAY ============
overlay_mask = None
if config.shapefile_overlay:
prediction_status["progress"] = f"Đang overlay shapefile: {config.shapefile_overlay}..."
print(f"[OVERLAY] Shapefile overlay requested: {config.shapefile_overlay}")
# Validate shapefile path exists
shapefile_path = Path(config.shapefile_overlay)
if not shapefile_path.exists():
print(f"[OVERLAY ERROR] Shapefile not found: {shapefile_path}")
print(f"[OVERLAY ERROR] Absolute path: {shapefile_path.absolute()}")
print(f"[OVERLAY ERROR] Current working directory: {Path.cwd()}")
else:
print(f"[OVERLAY] Shapefile exists: {shapefile_path.absolute()}")
try:
overlay_mask = rasterize_shapefile_overlay(str(shapefile_path), prediction_da)
if overlay_mask is not None and np.count_nonzero(overlay_mask) > 0:
print(f"[OVERLAY] Successfully rasterized shapefile boundaries ({np.count_nonzero(overlay_mask)} pixels)")
else:
print(f"[OVERLAY WARNING] Shapefile rasterized but no boundary pixels found")
except Exception as overlay_error:
print(f"[OVERLAY ERROR] Exception: {overlay_error}")
import traceback
traceback.print_exc()
else:
print(f"[OVERLAY] No shapefile overlay requested")
# Generate PNG preview for web display
prediction_status["progress"] = "Đang tạo PNG preview..."
png_file = output_dir / f"prediction_{timestamp}.png"
@@ -2148,6 +2288,42 @@ async def run_prediction(config: PredictionConfig):
# Plot prediction with colormap
im = ax.imshow(predictions_2d, cmap='tab20', interpolation='nearest')
# Overlay shapefile boundaries if available
if overlay_mask is not None and np.count_nonzero(overlay_mask) > 0:
print(f"[PNG OVERLAY] Overlaying {np.count_nonzero(overlay_mask)} boundary pixels")
# Create a mask for boundaries (where overlay_mask > 0)
boundary_mask = overlay_mask > 0
# Method: Direct overlay with high-contrast colors
# Create RGBA overlay image
overlay_rgba = np.zeros((*predictions_2d.shape, 4))
overlay_rgba[boundary_mask, 0] = 1.0 # Red = 1.0 (white)
overlay_rgba[boundary_mask, 1] = 1.0 # Green = 1.0 (white)
overlay_rgba[boundary_mask, 2] = 1.0 # Blue = 1.0 (white)
overlay_rgba[boundary_mask, 3] = 1.0 # Alpha = 1.0 (fully opaque)
# Overlay on top of prediction
ax.imshow(overlay_rgba, interpolation='nearest')
# Also add a black outline for better contrast
from scipy import ndimage
boundary_dilated = ndimage.binary_dilation(boundary_mask, iterations=1)
boundary_outline = boundary_dilated & ~boundary_mask
outline_rgba = np.zeros((*predictions_2d.shape, 4))
outline_rgba[boundary_outline, 0] = 0.0 # Black outline
outline_rgba[boundary_outline, 1] = 0.0
outline_rgba[boundary_outline, 2] = 0.0
outline_rgba[boundary_outline, 3] = 0.8
ax.imshow(outline_rgba, interpolation='nearest')
print(f"[PNG OVERLAY] Added shapefile boundaries to visualization (direct overlay method)")
else:
print(f"[PNG OVERLAY] No overlay mask or empty mask (pixels: {np.count_nonzero(overlay_mask) if overlay_mask is not None else 0})")
ax.set_title(f'Land Classification - {timestamp}', fontsize=16, fontweight='bold', pad=20)
ax.set_xlabel('X (pixels)', fontsize=11)
ax.set_ylabel('Y (pixels)', fontsize=11)
@@ -2227,7 +2403,9 @@ async def run_prediction(config: PredictionConfig):
"n_features": features.shape[1],
"feature_mode": feature_mode,
"used_radar": use_radar,
"model_used": config.model_filename
"model_used": config.model_filename,
"shapefile_overlay": config.shapefile_overlay,
"overlay_applied": overlay_mask is not None
}
# Auto generate prediction report
@@ -4406,8 +4584,69 @@ async def predict_with_ndvi(config: PredictionWithNDVIConfig, background_tasks:
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
# ============ SHAPEFILE OVERLAY ============
overlay_mask = None
if config.shapefile_overlay:
print(f"[OVERLAY] Shapefile overlay requested: {config.shapefile_overlay}")
# Validate shapefile path exists
shapefile_path = Path(config.shapefile_overlay)
if not shapefile_path.exists():
print(f"[OVERLAY ERROR] Shapefile not found: {shapefile_path}")
print(f"[OVERLAY ERROR] Absolute path: {shapefile_path.absolute()}")
else:
print(f"[OVERLAY] Shapefile exists: {shapefile_path.absolute()}")
try:
# Import rioxarray for rio accessor
import rioxarray
# Create temporary DataArray for rasterization
temp_da = xr.DataArray(
prediction_raster,
coords={
"y": np.linspace(bbox[3], bbox[1], height),
"x": np.linspace(bbox[0], bbox[2], width)
},
dims=["y", "x"]
)
temp_da.rio.write_crs("EPSG:4326", inplace=True)
temp_da.rio.write_transform(transform, inplace=True)
overlay_mask = rasterize_shapefile_overlay(str(shapefile_path), temp_da)
if overlay_mask is not None and np.count_nonzero(overlay_mask) > 0:
print(f"[OVERLAY] Successfully rasterized shapefile boundaries ({np.count_nonzero(overlay_mask)} pixels)")
else:
print(f"[OVERLAY WARNING] Shapefile rasterized but no boundary pixels found")
except Exception as overlay_error:
print(f"[OVERLAY ERROR] Exception: {overlay_error}")
import traceback
traceback.print_exc()
else:
print(f"[OVERLAY] No shapefile overlay requested")
fig, ax = plt.subplots(figsize=(14, 10), dpi=150)
im = ax.imshow(prediction_raster, cmap='tab20', interpolation='nearest')
# Overlay shapefile boundaries if available
if overlay_mask is not None and np.count_nonzero(overlay_mask) > 0:
print(f"[PNG OVERLAY] Overlaying {np.count_nonzero(overlay_mask)} boundary pixels (1px thin line)")
# Create a mask for boundaries
boundary_mask = overlay_mask > 0
# Create RGBA overlay image - thin 1px white line only
overlay_rgba = np.zeros((*prediction_raster.shape, 4))
overlay_rgba[boundary_mask, 0] = 1.0 # White (R=1)
overlay_rgba[boundary_mask, 1] = 1.0 # White (G=1)
overlay_rgba[boundary_mask, 2] = 1.0 # White (B=1)
overlay_rgba[boundary_mask, 3] = 1.0 # Fully opaque
ax.imshow(overlay_rgba, interpolation='nearest')
print(f"[PNG OVERLAY] Added thin 1px shapefile boundaries to visualization")
else:
print(f"[PNG OVERLAY] No overlay mask or empty mask")
ax.set_title(f'Land Classification - {timestamp}', fontsize=16, fontweight='bold', pad=20)
ax.set_xlabel('X (pixels)', fontsize=11)
ax.set_ylabel('Y (pixels)', fontsize=11)
+354
View File
@@ -0,0 +1,354 @@
import matplotlib.pyplot as plt
# Common imports and settings
import os, sys
os.environ['USE_PYGEOS'] = '0'
from IPython.display import Markdown
import pandas as pd
pd.set_option("display.max_rows", None)
import xarray as xr
# Datacube
import datacube
from datacube.utils.rio import configure_s3_access
from datacube.utils import masking
from datacube.utils.cog import write_cog
# https://github.com/GeoscienceAustralia/dea-notebooks/tree/develop/Tools
from dea_tools.plotting import display_map, rgb
from dea_tools.datahandling import mostcommon_crs
# EASI defaults
easinotebooksrepo = '/home/jovyan/easi-notebooks'
if easinotebooksrepo not in sys.path: sys.path.append(easinotebooksrepo)
from easi_tools import EasiDefaults, xarray_object_size, notebook_utils, unset_cachingproxy
from easi_tools.load_s2l2a import load_s2l2a_with_offset
from dask.distributed import progress
# Data tools
import numpy as np
from datetime import datetime
# Datacube
from datacube.utils import masking # https://github.com/opendatacube/datacube-core/blob/develop/datacube/utils/masking.py
from odc.algo import enum_to_bool # https://github.com/opendatacube/odc-algo/blob/main/odc/algo/_masking.py
from odc.algo import xr_reproject # https://github.com/opendatacube/odc-algo/blob/main/odc/algo/_warp.py
from datacube.utils.geometry import GeoBox, box # https://github.com/opendatacube/datacube-core/blob/develop/datacube/utils/geometry/_base.py
# Holoviews, Datashader and Bokeh
import hvplot.pandas
import hvplot.xarray
import holoviews as hv
import panel as pn
import colorcet as cc
import cartopy.crs as ccrs
from datashader import reductions
from holoviews import opts
from utils import load_data_geo
import rasterio
import rioxarray
# import geoviews as gv
# from holoviews.operation.datashader import rasterize
hv.extension('bokeh', logo=False)
from deafrica_tools.bandindices import calculate_indices
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, classification_report
from sklearn.preprocessing import LabelEncoder
from sklearn.pipeline import Pipeline
from sklearn.ensemble import RandomForestClassifier
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from shapely.geometry import Point, Polygon
import geopandas as gpd
from pyproj import CRS
from matplotlib.colors import ListedColormap
from holoviews import opts
from datashader import reductions
from bokeh.models.tickers import FixedTicker
from rioxarray.merge import merge_arrays
import joblib
def load_data(dc, date_range, longtitude_range, latitude_range):
product = 's2_l2a'
query = {
'product': product, # Product name
'x': longtitude_range, # "x" axis bounds
'y': latitude_range, # "y" axis bounds
'time': date_range, # Any parsable date strings
}
native_crs = notebook_utils.mostcommon_crs(dc, query)
print(f'Most common native CRS: {native_crs}')
measurements = [ 'red', 'nir', 'scl']
load_params = {
'measurements': measurements, # Selected measurement or alias names
'output_crs': native_crs, # Target EPSG code
'resolution': (-10, 10), # Target resolution
'group_by': 'solar_day', # Scene grouping
'dask_chunks': {'x': 2048, 'y': 2048}, # Dask chunks
}
data = load_s2l2a_with_offset(
dc,
query | load_params # Combine the two dicts that contain our search and load parameters
)
return data
def mask_clean(data):
flag_name = 'scl'
flag_desc = masking.describe_variable_flags(data[flag_name]) # Pandas dataframe
display(flag_desc)
display(flag_desc.loc['qa'].values[1])
# Create a "data quality" Mask layer
flags_def = flag_desc.loc['qa'].values[1]
good_pixel_flags = [flags_def[str(i)] for i in [2, 4, 5, 6]] # To pass strings to enum_to_bool()
# enum_to_bool calculates the pixel-wise "or" of each set of pixels given by good_pixel_flags
# 1 = good data
# 0 = "bad" data
good_pixel_mask = enum_to_bool(data[flag_name], good_pixel_flags)
data_layer_names = [x for x in data.data_vars if x != 'scl']
# Apply good pixel mask to blue, green, red and nir.
result = data[data_layer_names].where(good_pixel_mask).persist()
return result
def fill_nan(ndvi, time_split):
rs = []
for times in time_split:
tmp = ndvi.sel(time=times)
fill_ds = tmp.sel(time=times).bfill(dim='time')
fill_ds = fill_ds.sel(time=times).ffill(dim='time')
rs.append(fill_ds)
merged_ndvi = xr.concat([i for i in rs], dim="time")
fill_m = merged_ndvi.bfill(dim="time")
fill_m = fill_m.ffill(dim="time")
return fill_m
def load_train_data(train_path):
train = load_data_geo(train_path)
return train
def load_sen1(name_vh, name_vv):
dsvv = rioxarray.open_rasterio(name_vv)
dsvh = rioxarray.open_rasterio(name_vh)
return dsvh, dsvv
def get_data_sen1_and_sen2(train, average_ndvi, dsvh, dsvv):
loaded_datasets = {}
for idx, point in train.iterrows():
key = f"point_{idx + 1}"
try:
ndvi_data = average_ndvi.sel(x=point.geometry.x, y=point.geometry.y, method='nearest').values
vh_data = dsvh.sel(x=point.geometry.x, y=point.geometry.y, method='nearest').values
vv_data = dsvv.sel(x=point.geometry.x, y=point.geometry.y, method='nearest').values
loaded_datasets[key] = {
"data": np.concatenate((ndvi_data, vh_data, vv_data)),
"label": point.HT_code
}
except Exception as e:
# loaded_datasets[key] = None
print(e)
return loaded_datasets
def split_train_data(train, label_mapping, datasets):
label_encoder = LabelEncoder()
# Fit and transform the labels
labels = train.Hientrang.values
numeric_labels = label_encoder.fit_transform([label_mapping[label] for label in labels])
X = []
x_new = []
lb_new = []
for k, v in datasets.items():
X.append(v)
for i in range(len(X)):
if X[i] is not None:
x_new.append(X[i]["data"])
lb_new.append(numeric_labels[i])
X_train, X_temp, y_train, y_temp= train_test_split(x_new, lb_new, test_size=0.4, random_state=42)
X_val, X_test, y_val, y_test = train_test_split(X_temp, y_temp, test_size=0.5, random_state=42)
return X_train, X_val, X_test, y_train, y_val, y_test
def train_with_rf(X_train, X_val, y_train, y_val):
# Takes 1-2 minutes to complete
# Tạo RandomForestClassifier mặc định để sử dụng làm mô hình ban đầu trong pipeline
base_model = RandomForestClassifier(random_state=42, n_jobs=-1)
# Tạo pipeline
pipeline = Pipeline([
# ('imputer', SimpleImputer(strategy='mean')),
('scaler', StandardScaler()),
('classifier', base_model),
])
# Thiết lập các tham số bạn muốn tối ưu hóa
param_grid = {
'classifier__n_estimators': [100, 300, 500, 700, 1000],
'classifier__max_depth': [6, 8, 10, 15, 20],
'classifier__criterion': ['gini', 'entropy'],
}
# Sử dụng GridSearchCV để tìm bộ tham số tốt nhất
grid_search = GridSearchCV(pipeline, param_grid, cv=5, scoring='accuracy', n_jobs=-1)
grid_search.fit(X_train, y_train)
# In ra bộ tham số tốt nhất
best_params = grid_search.best_params_
print("Best Parameters:", best_params)
# Dự đoán trên tập kiểm tra
y_pred = grid_search.predict(X_val)
# Đánh giá kết quả
accuracy = accuracy_score(y_val, y_pred)
print(f"Accuracy: {round(accuracy, 2)*100} %")
return grid_search
def save_model(name_file, grid_search):
dir_save_model = "model_train"
if not os.path.exists(dir_save_model):
os.mkdir(dir_save_model)
joblib.dump(grid_search, os.path.join(dir_save_model, name_file))
print("Done!")
def predict(model, data_crs, ndvi, vh, vv):
data_predict = []
for i in range(ndvi.shape[1]):
ndvi_tmp = ndvi.isel(y=i).values
vh_data = vh.sel(y=ndvi.y.values[i], method='nearest').values
vv_data = vv.sel(y=ndvi.y.values[i], method='nearest').values
all_tmp = np.concatenate((ndvi_tmp, vh_data, vv_data), axis=0)
data_predict.extend(all_tmp.T)
y_pred = model.predict(data_predict)
final_label = y_pred.reshape(ndvi.y.shape[0], ndvi.x.shape[0])
final_xarray_save = xr.DataArray(final_label, dims=("y", "x"))
final_xarray_save = final_xarray_save.rio.write_crs(data_crs)
x_values = ndvi.x.values
y_values = ndvi.y.values
data_array = xr.DataArray(final_xarray_save,
coords={'x': x_values, 'y': y_values},
dims=['y', 'x'])
data_array = data_array.rio.write_crs(ndvi.rio.crs)
return data_array
def cut_according_shp(thuanhoa_path, average_ndvi, data_array):
gdf = gpd.read_file(thuanhoa_path)
gdf = gdf.to_crs(average_ndvi.rio.crs)
polygon_coords = list(gdf.geometry.values[0].exterior.coords)
polygon_coordinates = [(x, y) for x, y in polygon_coords]
geometries = [
{
'type': 'Polygon',
'coordinates': [polygon_coordinates]
}
]
region_result = data_array.rio.clip(geometries, data_array.rio.crs, drop=False)
region_result = region_result.where(region_result >= 0, float('nan'))
return region_result
def compare(KD_path, KetQuaPhanLoaiDat, CODE_MAP, HT_MAP):
gdf = gpd.read_file(KD_path, crs="EPSG:9209")
polygon = gdf.geometry.values
label = gdf.tenchu.values
ouput_image = rioxarray.open_rasterio(KetQuaPhanLoaiDat)
code_tq = HT_MAP["TQ"]["data"][0]
code_pnn = HT_MAP["PNN"]["data"][0]
result = {}
for key, values in HT_MAP.items():
print(f"process {key}")
array_list = []
for i in range(len(polygon)):
po = polygon[i]
lb = label[i]
code_lb = CODE_MAP.get(lb, code_tq)
try:
qr = ouput_image.rio.clip([po], "EPSG:9209")
if code_lb in values["data"]:
if code_lb == code_pnn:
qr = qr.where((qr != float(code_pnn)), np.nan)
# qr = qr.where((qr != 3.0), np.nan)
elif code_lb == code_tq:
qr = qr.where((qr != float(code_pnn)), np.nan)
qr = qr.where((qr != 3.0), np.nan)
else:
qr = qr.where(qr != float(code_lb), np.nan)
else:
qr.values[:, :, :] = np.nan
array_list.append(qr)
except Exception as e:
pass
result.update({key: array_list})
return result
def save_result(result, HT_MAP):
# cmap = ListedColormap(colors)
save_path = "ThuanHoa/KetQua"
if not os.path.exists(save_path):
os.mkdir(save_path)
for k, v in result.items():
rs = merge_arrays(v, nodata = np.nan)
rs.rio.to_raster(f"{save_path}/{k}.tif")
print(f"save {save_path}/{k}.tif")
# img = rs.plot(cmap=cmap, add_colorbar=False)
# cbar = plt.colorbar(img)
# cbar.ax.set_yticklabels(labels)
# plt.title(f'{HT_MAP[k]["name"]}')
# plt.axis('off')
# plt.show()
def accuracy_test(test, data_array):
# cấu hình nhãn dữ liệu
label_mapping = {
"Lua tom": "0",
"Lua": "1",
"CHN": "2",
"CLN": "3",
"TS": "4",
"Song": "5",
"Dat xay dung": "6",
"Rung": "7"
}
chk = []
pred = []
dd = []
for idx, point in test.iterrows():
label = point.LULC
predict = data_array.sel(x=point.geometry.x, y=point.geometry.y, method='nearest').values
pred.append(label_mapping[label])
dd.append(str(predict))
chk.append(predict == int(label_mapping[label]))
test["code"] = pred
test["dd"] = dd
test["check"] = chk
path = "ThuanHoa/TestAccuracy"
if not os.path.exists(path):
os.mkdir(path)
test.to_file(f"{path}/result.shp")
percentage_true = np.mean(chk) * 100
print(f"độ chính xác: {percentage_true:.2f}%")
Submodule backup_code_training/CSIROBoeingPhase5-Vietnam added at 8f0cb55cba
@@ -9,7 +9,6 @@ from typing import Tuple, Optional, Dict
from sklearn.neighbors import KNeighborsRegressor
from sklearn.ensemble import RandomForestRegressor
import warnings
from pathlib import Path
warnings.filterwarnings('ignore')
@@ -376,22 +375,6 @@ class DeepInpaintingStrategy(CloudRemovalStrategy):
# Convert to tensor and add batch dimension
input_tensor = torch.from_numpy(input_array).unsqueeze(0).to(self.device)
if hasattr(self.model, 'encoder'):
# Custom UNet from train_cloud_removal.py
expected_channels = self.model.encoder[0].double_conv[0].in_channels
elif hasattr(self.model, 'inc') and hasattr(self.model.inc.double_conv[0], 'in_channels'):
expected_channels = self.model.inc.double_conv[0].in_channels
elif hasattr(self.model, 'conv1') and hasattr(self.model.conv1, 'in_channels'):
expected_channels = self.model.conv1.in_channels
else:
expected_channels = 6
if expected_channels > input_tensor.shape[1]:
pad_channels = expected_channels - input_tensor.shape[1]
padding = torch.zeros(1, pad_channels, *input_tensor.shape[2:]).to(self.device)
input_tensor = torch.cat([input_tensor, padding], dim=1)
# Run through U-Net
with torch.no_grad():
output_tensor = self.model(input_tensor)
View File
-93
View File
@@ -1,93 +0,0 @@
import os
import torch
import numpy as np
import xarray as xr
from torch.utils.data import Dataset
import glob
class NDVITimeSeriesDataset(Dataset):
def __init__(self, sequence_length=3, spatial=False):
"""
Đọc dữ liệu S2 từ cache, tính NDVI và tạo Time-Series.
spatial=False -> Output 1D cho LSTM/ARIMA
spatial=True -> Output 2D cho ConvLSTM
"""
self.sequence_length = sequence_length
self.spatial = spatial
self.data_seqs = []
self.targets = []
# Load from cache
cache_files = glob.glob("dataset_cache/*.nc")
s2_files = [f for f in cache_files if len(os.path.basename(f)) == 35] # S2 cache filenames usually have length 32 + 3 (.nc)
if not s2_files:
print("[WARNING] Không tìm thấy dữ liệu S2 trong cache! Dùng dummy data.")
self._create_dummy()
return
try:
print(f"[DATA] Loading real data from {s2_files[0]}")
ds = xr.open_dataset(s2_files[0], engine='netcdf4')
if 'time' not in ds.dims or len(ds.time) < sequence_length + 1:
self._create_dummy()
return
# Tính NDVI: (B08 - B04) / (B08 + B04)
b8 = ds['B08'].astype(np.float32)
b4 = ds['B04'].astype(np.float32)
ndvi = (b8 - b4) / (b8 + b4 + 1e-8)
ndvi = ndvi.fillna(0).values # shape: (time, y, x)
# Lấy 1 pixel trung tâm hoặc toàn bộ ảnh
if not self.spatial:
# Average pooling over space for 1D time series
ndvi = ndvi.mean(axis=(1, 2)) # shape: (time,)
for i in range(len(ndvi) - sequence_length):
self.data_seqs.append(ndvi[i:i+sequence_length])
self.targets.append(ndvi[i+sequence_length])
else:
# Spatial data for ConvLSTM
# Downsample to 64x64 to avoid OOM
from skimage.transform import resize
T = len(ndvi)
ndvi_resized = np.zeros((T, 64, 64))
for t in range(T):
ndvi_resized[t] = resize(ndvi[t], (64, 64))
for i in range(T - sequence_length):
self.data_seqs.append(ndvi_resized[i:i+sequence_length]) # (seq, 64, 64)
self.targets.append(ndvi_resized[i+sequence_length]) # (64, 64)
except Exception as e:
print(f"[ERROR] {e}. Dùng dummy data.")
self._create_dummy()
def _create_dummy(self):
T = 20
if not self.spatial:
ndvi = np.random.rand(T).astype(np.float32)
for i in range(T - self.sequence_length):
self.data_seqs.append(ndvi[i:i+self.sequence_length])
self.targets.append(ndvi[i+self.sequence_length])
else:
ndvi = np.random.rand(T, 64, 64).astype(np.float32)
for i in range(T - self.sequence_length):
self.data_seqs.append(ndvi[i:i+self.sequence_length])
self.targets.append(ndvi[i+self.sequence_length])
def __len__(self):
return len(self.data_seqs)
def __getitem__(self, idx):
x = torch.tensor(self.data_seqs[idx], dtype=torch.float32)
y = torch.tensor(self.targets[idx], dtype=torch.float32)
if not self.spatial:
x = x.unsqueeze(1) # (seq_len, features=1)
y = y.unsqueeze(0) # (1,)
else:
x = x.unsqueeze(1) # (seq_len, channels=1, H, W)
y = y.unsqueeze(0) # (1, H, W)
return x, y
Binary file not shown.

Some files were not shown because too many files have changed in this diff Show More