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

hoàn thành chức năng tính ndvi analysys 2 màn hình

Browse Source
This commit is contained in:
Victor Phan
2025-12-24 13:57:08 +07:00
parent be90d26d24
commit 5342df77b4
25 changed files with 4755 additions and 460 deletions
Download Patch File Download Diff File Expand all files Collapse all files
train_module.py
+139 -39
View File
@@ -118,6 +118,9 @@ import planetary_computer
from pystac_client import Client
from odc.stac import load as stac_load
# Feature extraction
from feature_extractor import get_feature_extractor
def train_model(
bbox=[105.6, 9.3, 106.2, 9.8],
@@ -133,6 +136,7 @@ def train_model(
use_gpu=True,
use_cache=True,
test_size=0.2,
feature_mode='simple',
output_model_path=None,
status_callback=None,
cancel_check=None
@@ -155,6 +159,7 @@ def train_model(
status_callback: Optional callback function to report progress
cancel_check: Optional function that returns True if training should be cancelled
test_size: Fraction of data to use for test set (0-1)
feature_mode: 'simple' (3 features), 'temporal' (39 features), or 'extended' (15 features)
Returns:
Dictionary containing training results
@@ -237,16 +242,26 @@ def train_model(
# Sign and load Sentinel-2 data
update_status("Loading Sentinel-2 data...", 25)
items_s2 = [planetary_computer.sign(item) for item in items_s2]
# Load different bands based on feature mode
if feature_mode == 'simple':
bands_to_load = ["B04", "B08", "SCL"]
else: # temporal or extended
bands_to_load = ["B02", "B03", "B04", "B08", "B11", "SCL"]
ds_s2 = stac_load(
items_s2,
bands=["B04", "B08", "SCL"],
bands=bands_to_load,
crs="EPSG:32648",
resolution=resolution,
bbox=bbox,
patch_url=planetary_computer.sign,
fail_on_error=False,
)
ds_s2 = ds_s2.rename({"B04": "red", "B08": "nir", "SCL": "scl"})
# Rename for compatibility (simple mode)
if "B04" in ds_s2 and "red" not in ds_s2:
ds_s2 = ds_s2.rename({"B04": "red", "B08": "nir", "SCL": "scl"})
check_cancellation()
@@ -285,14 +300,9 @@ def train_model(
check_cancellation()
# Calculate NDVI
update_status("Calculating NDVI...", 50)
ndvi = (ds_s2['nir'] - ds_s2['red']) / (ds_s2['nir'] + ds_s2['red'] + 1e-8)
# Apply cloud mask
cloud_mask = ds_s2['scl'].isin([1, 3, 8, 9, 10])
ndvi_masked = ndvi.where(~cloud_mask)
ndvi_mean = ndvi_masked.mean(dim='time')
# ============ FEATURE EXTRACTION ============
update_status(f"Initializing FeatureExtractor (mode={feature_mode})...", 50)
extractor = get_feature_extractor(mode=feature_mode)
# Load training data
update_status("Loading training data...", 55)
@@ -311,32 +321,106 @@ def train_model(
if label_column is None:
raise ValueError(f"Cannot find label column in shapefile. Available: {list(train_gdf.columns)}")
# Extract features
update_status("Extracting features from training points...", 60)
features = []
labels = []
# Extract features using FeatureExtractor
update_status("Extracting features from satellite data...", 60)
for idx, row in train_gdf.iterrows():
point = row.geometry
x_coord = point.x
y_coord = point.y
label = row[label_column]
if feature_mode == 'simple':
# For simple mode: calculate NDVI first
ndvi = (ds_s2['nir'] - ds_s2['red']) / (ds_s2['nir'] + ds_s2['red'] + 1e-8)
# Apply cloud mask
cloud_mask = ds_s2['scl'].isin([1, 3, 8, 9, 10])
ndvi_masked = ndvi.where(~cloud_mask)
try:
ndvi_val = ndvi_mean.sel(x=x_coord, y=y_coord, method='nearest').values
vh_val = ds_s1['vh_db'].sel(x=x_coord, y=y_coord, method='nearest').mean(dim='time').values
vv_val = ds_s1['vv_db'].sel(x=x_coord, y=y_coord, method='nearest').mean(dim='time').values
# Extract features at training points
features = []
labels = []
for idx, row in train_gdf.iterrows():
point = row.geometry
x_coord = point.x
y_coord = point.y
label = row[label_column]
feature_vec = [ndvi_val, vh_val, vv_val]
if not np.isnan(feature_vec).any():
features.append(feature_vec)
labels.append(label)
except:
continue
try:
ndvi_val = ndvi_masked.sel(x=x_coord, y=y_coord, method='nearest').mean(dim='time').values
vh_val = ds_s1['vh_db'].sel(x=x_coord, y=y_coord, method='nearest').mean(dim='time').values
vv_val = ds_s1['vv_db'].sel(x=x_coord, y=y_coord, method='nearest').mean(dim='time').values
feature_vec = [float(ndvi_val), float(vh_val), float(vv_val)]
if not np.isnan(feature_vec).any():
features.append(feature_vec)
labels.append(label)
except Exception as e:
continue
features = np.array(features)
labels = np.array(labels)
features = np.array(features)
labels = np.array(labels)
else: # temporal or extended mode
# Apply cloud mask for temporal/extended modes
if 'scl' in ds_s2 or 'SCL' in ds_s2:
scl_band = ds_s2['scl'] if 'scl' in ds_s2 else ds_s2['SCL']
cloud_mask = scl_band.isin([1, 3, 8, 9, 10])
for band in ds_s2.data_vars:
if band != 'scl' and band != 'SCL':
ds_s2[band] = ds_s2[band].where(~cloud_mask)
# Extract features at training points
features = []
labels = []
for idx, row in train_gdf.iterrows():
point = row.geometry
x_coord = point.x
y_coord = point.y
label = row[label_column]
try:
# Extract point data from S2
point_s2 = ds_s2.sel(x=x_coord, y=y_coord, method='nearest')
# Extract point data from S1
vh_val = ds_s1['vh_db'].sel(x=x_coord, y=y_coord, method='nearest').mean(dim='time').values
vv_val = ds_s1['vv_db'].sel(x=x_coord, y=y_coord, method='nearest').mean(dim='time').values
# Create minimal dataset for feature extraction
point_data = xr.Dataset({
'B02': point_s2['B02'],
'B03': point_s2['B03'],
'B04': point_s2['B04'],
'B08': point_s2['B08'],
'B11': point_s2['B11']
})
# Create VH/VV DataArrays (without spatial dims, just time if exists)
if 'time' in point_data.dims:
vh_da = xr.DataArray([vh_val] * len(point_data.time), dims=['time'])
vv_da = xr.DataArray([vv_val] * len(point_data.time), dims=['time'])
else:
vh_da = xr.DataArray([vh_val])
vv_da = xr.DataArray([vv_val])
# Extract features using FeatureExtractor
# Note: extractor.extract returns (n_pixels, n_features), we take first row
feature_vec = extractor.extract(
s2_data=point_data,
vh_data=vh_da,
vv_data=vv_da
)
# If feature_vec is 2D, take first row
if len(feature_vec.shape) > 1:
feature_vec = feature_vec[0]
if not np.isnan(feature_vec).any():
features.append(feature_vec)
labels.append(label)
except Exception as e:
continue
features = np.array(features)
labels = np.array(labels)
check_cancellation()
@@ -352,6 +436,7 @@ def train_model(
'bbox': bbox,
'time_range': time_range,
'resolution': resolution,
'feature_mode': feature_mode,
'timestamp': datetime.now().isoformat()
}
joblib.dump(cache_data, cache_file)
@@ -486,17 +571,25 @@ def train_model(
# Confusion matrix
conf_matrix = confusion_matrix(y_test, y_pred).tolist()
# Save model
# Save model using ModelManager
update_status("Saving model...", 95)
os.makedirs(os.path.dirname(output_model_path), exist_ok=True)
joblib.dump({'model': model, 'label_encoder': label_encoder}, output_model_path)
# Save model info
# Get feature names from extractor
if feature_mode == 'temporal':
# Calculate n_timesteps from data
n_timesteps = len(features[0]) // 3 - 1 # (NDVI + NDWI + NDBI) * n_timesteps + 3 radar features
feature_names = extractor.get_feature_names(n_timesteps=n_timesteps)
else:
feature_names = extractor.get_feature_names()
# Prepare metadata
info = {
"timestamp": datetime.now().isoformat(),
"data_source": "Microsoft Planetary Computer STAC",
"collections": ["sentinel-2-l2a", "sentinel-1-rtc"],
"features": ["NDVI_mean", "VH_dB_mean", "VV_dB_mean"],
"features": feature_names,
"feature_mode": feature_mode,
"training_samples": len(X_train),
"testing_samples": len(X_test),
"test_size": test_size,
@@ -518,9 +611,16 @@ def train_model(
"resolution": resolution
}
info_path = output_model_path.replace('.joblib', '_info.json')
with open(info_path, 'w') as f:
json.dump(info, f, indent=2)
# Use ModelManager to save
from model_manager import get_model_manager
model_manager = get_model_manager()
model_filename = os.path.basename(output_model_path)
model_manager.save_model(
model=model,
metadata=info,
model_filename=model_filename,
label_encoder=label_encoder
)
update_status("Training complete!", 100)