remote-sensing/scripts/training/new_train.py

#!/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()