273 lines
7.8 KiB
Python
273 lines
7.8 KiB
Python
#!/usr/bin/env python
|
|
# coding: utf-8
|
|
|
|
# In[49]:
|
|
|
|
|
|
get_ipython().run_cell_magic('time', '', '%matplotlib inline\nfrom new_import 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 = datacube.Datacube()\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')
|
|
|
|
|
|
# LOAD VH, VV
|
|
|
|
# In[47]:
|
|
|
|
|
|
## cấu hình thời gian lấy ảnh và tọa độ
|
|
date_range = ('2022-09-01', '2023-10-01')
|
|
longtitude_range = (105.5, 106.4)
|
|
latitude_range = (9.2, 10.0)
|
|
|
|
|
|
# In[3]:
|
|
|
|
|
|
## cấu hình dữ liệu train và vh vv file
|
|
train_path = "train/ST_training data_updated_1130points.shp" # đường dẫn shp file train
|
|
name_vh = "vh-0922_0923-full_ST.tif"
|
|
name_vv = "vv-0922_0923-full_ST.tif"
|
|
|
|
|
|
train = load_train_data(train_path)
|
|
|
|
|
|
# In[4]:
|
|
|
|
|
|
# %%time
|
|
# ## tải về dữ liệu sen1
|
|
# import os
|
|
# if not os.path.exists(name_vh):
|
|
# !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):
|
|
# !aws s3 cp s3://easi-asia-dc-data/staging/ctu/sentinel-1/vv-0922_0923-full_ST.tif vv-0922_0923-full_ST.tif
|
|
|
|
|
|
# In[5]:
|
|
|
|
|
|
|
|
|
|
|
|
# In[38]:
|
|
|
|
|
|
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[43]:
|
|
|
|
|
|
vv_data = ds.vv
|
|
vv_data
|
|
|
|
|
|
# In[44]:
|
|
|
|
|
|
bbox = [105.5, 9.2, 106.4, 10.0]
|
|
time_range = "2022-09-01/2023-10-01"
|
|
dsvh, dsvv = load_sen1(bbox, time_range)
|
|
dsvv
|
|
|
|
|
|
# LOAD SENTINEL 2
|
|
#
|
|
#
|
|
|
|
# In[50]:
|
|
|
|
|
|
data = load_data(dc, date_range, longtitude_range, latitude_range)
|
|
notebook_utils.heading(notebook_utils.xarray_object_size(data))
|
|
display(data)
|
|
|
|
|
|
# In[8]:
|
|
|
|
|
|
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')
|
|
|
|
|
|
# CALCULATING THE MEAN VALUE AND FILL TO NAN POINT
|
|
|
|
# In[9]:
|
|
|
|
|
|
ds1 = calculate_indices(result, index='NDVI', satellite_mission='s2')
|
|
ndvi = ds1["NDVI"]
|
|
average_ndvi = ndvi.resample(time='1M').mean().persist() ## tính mean cho từng tháng -> time = 12
|
|
progress(average_ndvi)
|
|
|
|
|
|
# In[10]:
|
|
|
|
|
|
dsvh.shape
|
|
|
|
|
|
# In[11]:
|
|
|
|
|
|
average_ndvi = average_ndvi.compute()
|
|
average_ndvi = average_ndvi[:, :dsvh.shape[1], :dsvh.shape[2]]
|
|
|
|
|
|
# In[12]:
|
|
|
|
|
|
get_ipython().run_cell_magic('time', '', "filled_ds = average_ndvi.bfill(dim='time')\nfilled_ds = filled_ds.ffill(dim='time')\n")
|
|
|
|
|
|
# FIND NAN POINT AFTER FILLING AND FILLING AGAIN WITH LINEARREGRESSION ALGORITHM
|
|
|
|
# In[13]:
|
|
|
|
|
|
nan_mask = filled_ds.isnull()
|
|
|
|
# Print the NaN mask
|
|
# print(nan_mask)
|
|
|
|
# Count the number of NaNs
|
|
num_nans = nan_mask.sum()
|
|
print(f'Number of NaNs: {num_nans.values}')
|
|
|
|
|
|
# In[14]:
|
|
|
|
|
|
from sklearn.preprocessing import PolynomialFeatures
|
|
from sklearn.linear_model import LinearRegression
|
|
from sklearn.ensemble import RandomForestRegressor
|
|
|
|
mask = ~np.isnan(filled_ds)
|
|
X_train = np.stack([dsvh.values[mask], dsvv.values[mask]], axis=1)
|
|
y_train = filled_ds.values[mask]
|
|
|
|
|
|
# In[15]:
|
|
|
|
|
|
model = LinearRegression()
|
|
model.fit(X_train, y_train)
|
|
|
|
|
|
# In[16]:
|
|
|
|
|
|
X_pred = np.stack([dsvh.values[~mask], dsvv.values[~mask]], axis=1)
|
|
filled_ds.values[~mask] = model.predict(X_pred)
|
|
|
|
|
|
# MATCH LABEL TO DATASET
|
|
|
|
# In[17]:
|
|
|
|
|
|
get_ipython().run_cell_magic('time', '', '\n# Takes 1 minute to complete.\nloaded_datasets = {}\nfor idx, point in train.iterrows():\n key = f"point_{idx + 1}"\n try:\n ndvi_data = filled_ds.sel(x=point.geometry.x, y=point.geometry.y, method=\'nearest\').values\n vh_data = dsvh.sel(x=point.geometry.x, y=point.geometry.y, method=\'nearest\').values\n vv_data = dsvv.sel(x=point.geometry.x, y=point.geometry.y, method=\'nearest\').values\n loaded_datasets[key] = {\n "data": np.concatenate((ndvi_data, vh_data, vv_data)),\n "label": point.HT_code\n }\n except Exception as e:\n # loaded_datasets[key] = None\n print(e)\n')
|
|
|
|
|
|
# In[18]:
|
|
|
|
|
|
label_mapping = {
|
|
"Lua tom": "0",
|
|
"Lua": "1",
|
|
"CHN": "2",
|
|
"CLN": "3",
|
|
"TS": "4",
|
|
"Song": "5",
|
|
"Dat xay dung": "6",
|
|
"Rung": "7"
|
|
}
|
|
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])
|
|
|
|
|
|
# In[19]:
|
|
|
|
|
|
X = []
|
|
x_new = []
|
|
lb_new = []
|
|
for k, v in loaded_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])
|
|
|
|
|
|
# BUILDING DATASETS
|
|
|
|
# In[20]:
|
|
|
|
|
|
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)
|
|
|
|
|
|
# TRAIN MODEL
|
|
|
|
# In[21]:
|
|
|
|
|
|
get_ipython().run_cell_magic('time', '', 'from sklearn.pipeline import Pipeline\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.model_selection import GridSearchCV\nfrom xgboost import XGBClassifier\nfrom sklearn.neighbors import KNeighborsClassifier\nfrom sklearn.naive_bayes import GaussianNB\nfrom sklearn.svm import SVC\nfrom sklearn.metrics import accuracy_score\n\n# Define the models\nrf_model = XGBClassifier(
|
|
n_estimators=200,
|
|
max_depth=30,
|
|
tree_method="hist",
|
|
device="cuda",
|
|
random_state=42,
|
|
n_jobs=-1,
|
|
verbosity=1
|
|
)\nknn_model = KNeighborsClassifier()\nnb_model = GaussianNB()\nsvm_model = SVC()\n\n# Create a pipeline\npipeline = Pipeline([\n (\'scaler\', StandardScaler()), # Apply scaling\n (\'classifier\', rf_model) # Placeholder, will be set by param_grid\n])\n\n# Define the parameter grid for each classifier\nparam_grid = [\n # RandomForest\n {\n \'classifier\': [rf_model],\n \'classifier__n_estimators\': [100, 300, 500, 700],\n \'classifier__max_depth\': [6, 8, 10, 15],\n \'classifier__criterion\': [\'gini\', \'entropy\'],\n },\n # KNeighborsClassifier\n {\n \'classifier\': [knn_model],\n \'classifier__n_neighbors\': [3, 5, 7, 9],\n \'classifier__weights\': [\'uniform\', \'distance\'],\n \'classifier__metric\': [\'euclidean\', \'manhattan\']\n },\n # Naive Bayes (GaussianNB doesn\'t have hyperparameters to tune here)\n {\n \'classifier\': [nb_model],\n },\n # SVM\n {\n \'classifier\': [svm_model],\n \'classifier__C\': [0.1, 1, 10, 100],\n \'classifier__kernel\': [\'linear\', \'rbf\'],\n \'classifier__gamma\': [\'scale\', \'auto\']\n }\n]\n\n# Use GridSearchCV to find the best classifier and hyperparameters\ngrid_search = GridSearchCV(pipeline, param_grid, cv=5, scoring=\'accuracy\', n_jobs=-1)\ngrid_search.fit(X_train, y_train)\n\n# Print out the best parameters and classifier\nbest_params = grid_search.best_params_\nprint("Best Parameters:", best_params)\n\n# Make predictions on the validation set\ny_pred = grid_search.predict(X_val)\n\n# Evaluate the results\naccuracy = accuracy_score(y_val, y_pred)\nprint(f"Accuracy: {round(accuracy, 2)*100} %")\n')
|
|
|
|
|
|
# In[22]:
|
|
|
|
|
|
## check accuracy score
|
|
|
|
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[23]:
|
|
|
|
|
|
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, "model_new2.joblib"))
|
|
|
|
|
|
# In[24]:
|
|
|
|
|
|
client.close()
|
|
cluster.close()
|
|
|