Files
remote-sensing/new_import_ODC.py
T

589 lines
19 KiB
Python

TEST_MODE = True
RESOLUTION = 1000 if TEST_MODE else 10
import matplotlib.pyplot as plt
# Common imports and settings
import os, sys
os.environ['USE_PYGEOS'] = '0'
os.environ["GDAL_HTTP_MAX_RETRY"] = "5"
os.environ["GDAL_HTTP_RETRY_DELAY"] = "2"
os.environ["GDAL_HTTP_CONNECTION_TIMEOUT"] = "10"
os.environ["GDAL_HTTP_TIMEOUT"] = "30"
os.environ["CPL_VSIL_CURL_ALLOWED_EXTENSIONS"] = ".tif,.tiff"
os.environ["GDAL_DISABLE_READDIR_ON_OPEN"] = "YES"
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
# removed deafrica_tools imports to avoid ipyleaflet error
# EASI defaults
easinotebooksrepo = '/home/x79/CSIROBoeingPhase4-Vietnam'
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
# removed xr_reproject
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
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score
import joblib
def load_data(dc, date_range, longtitude_range, latitude_range):
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"})
return data
def mask_clean(data):
# For Sentinel-2 L2A SCL:
# 2: Dark Area Pixels, 4: Vegetation, 5: Not Vegetated, 6: Water
good_pixel_mask = data['scl'].isin([2, 4, 5, 6])
data_layer_names = [x for x in data.data_vars if x != 'scl']
# Apply good pixel mask
result = data[data_layer_names].where(good_pixel_mask).persist()
return result
def fill_nan(ndvi, time_split):
if len(ndvi.time) == 0:
return ndvi
# If the total time duration is less than 90 days, skip seasonal splitting
try:
total_days = (ndvi.time[-1] - ndvi.time[0]).dt.days.item()
if total_days < 90:
return ndvi.bfill(dim="time").ffill(dim="time")
except Exception:
pass
rs = []
for times in time_split:
try:
tmp = ndvi.sel(time=times)
if len(tmp.time) == 0:
continue
fill_ds = tmp.bfill(dim='time').ffill(dim='time')
rs.append(fill_ds)
except Exception:
continue
if len(rs) == 0:
return ndvi.bfill(dim="time").ffill(dim="time")
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(bbox, time_range):
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")
return vh, vv
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, model, metadata=None, label_encoder=None):
"""
Save model với metadata để tương thích với ModelManager
Args:
name_file: Tên file model
model: Model object
metadata: Dict chứa thông tin về model (optional)
label_encoder: Label encoder (optional)
"""
from model_manager import get_model_manager
dir_save_model = "model_train"
if not os.path.exists(dir_save_model):
os.mkdir(dir_save_model)
# Nếu có metadata, sử dụng ModelManager
if metadata is not None:
model_manager = get_model_manager()
model_manager.save_model(
model=model,
metadata=metadata,
model_filename=name_file,
label_encoder=label_encoder
)
else:
# Legacy mode: save trực tiếp (backward compatibility)
model_data = {
'model': model,
'label_encoder': label_encoder
} if label_encoder is not None else model
joblib.dump(model_data, os.path.join(dir_save_model, name_file))
print(f"✅ Model saved: {name_file}")
if metadata:
print(f" - Type: {metadata.get('model_type', 'N/A')}")
print(f" - Features: {metadata.get('n_features', 'N/A')}")
print(f" - Accuracy: {metadata.get('test_accuracy', 'N/A')}")
def predict(model, data_crs, ndvi, vh, vv):
# Unpack model if it is wrapped in a dictionary (from ModelManager)
if isinstance(model, dict) and 'model' in model:
model = model['model']
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 load_data_sen1(dc, date_range, coordinates):
longtitude_range, latitude_range = coordinates
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-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"
)
# notebook_utils.heading(notebook_utils.xarray_object_size(data_sen1))
# display(data_sen1)
dsvh = data_sen1.vh
dsvv = data_sen1.vv
return dsvh, dsvv
def calculate_average(data, time_pattern='1M'):
return data.resample(time=time_pattern).mean().persist()
def load_data_sen2(dc, date_range, coordinates):
longtitude_range, latitude_range = coordinates
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"})
return data
def mask_cloud(data):
# For Sentinel-2 L2A SCL:
# 2: Dark Area Pixels, 4: Vegetation, 5: Not Vegetated, 6: Water
good_pixel_mask = data['scl'].isin([2, 4, 5, 6])
data_layer_names = [x for x in data.data_vars if x != 'scl']
# Apply good pixel mask
result = data[data_layer_names].where(good_pixel_mask).persist()
return result
def find_best_model(dataset):
X_train, X_val, y_train, y_val = dataset
# 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_result_new(result, save_path, HT_MAP):
# cmap = ListedColormap(colors)
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}%")