03_seagrass_distribution_lacbay.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Mapping seagrass distribution in Lac Bay</b>  \n",
    "Notebook for classifying and analyzing seagrass distribution in Lac Bay, Bonaire with Sentinel-2 images\n",
    "\n",
    "* Decision Tree Classifier (DTC) and Maximum Likelihood Classifier (MLC) are employed\n",
    "* Training sites covering 2 different classes (non-seagrass,seagrass) are used to extract pixel values (training samples) over RGB bands \n",
    "* 80:20 train-test ratio for splitting the training samples\n",
    "* K-Fold cross-validation performed for tuning the DTC model\n",
    "* MLC model developed with 4 different chi-square thresholds: 0% (base), 10%,20%,50%\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import re\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import rasterio as rio\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from glob import glob\n",
    "import geopandas as gpd\n",
    "from joblib import dump,load\n",
    "from tqdm import tqdm,tqdm_notebook\n",
    "\n",
    "#custom functions\n",
    "from Python.prep_raster import stack_bands,clip_raster,pixel_sample,computeIndexStack\n",
    "from Python.spec_analysis import transpose_df,jmd2df\n",
    "from Python.data_viz import ridgePlot,validation_curve_plot\n",
    "from Python.mlc import mlClassifier\n",
    "from Python.calc_acc import calc_acc\n",
    "from Python.pred_raster import stack2pred, dtc_pred_stack\n",
    "\n",
    "#sklearn functions\n",
    "from sklearn.model_selection import train_test_split,validation_curve\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "\n",
    "#setup IO directories\n",
    "parent_dir = os.path.join(os.path.abspath('..'),'objective3')                  #change according to preference\n",
    "sub_dirs = ['fullstack','clippedstack','indexstack','predicted','stack2pred']\n",
    "make_dirs = [os.makedirs(os.path.join(parent_dir,name),exist_ok=True) for name in sub_dirs]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Sentinel-2 data preparation</b>\n",
    "* Resample coarse bands to 10m resolution\n",
    "* Stack multiband images \n",
    "* Calculate spectral indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#dates considered for classification and analysis \n",
    "dates = [20190108,20190128,20190212,20190304,20190821,20191129]\n",
    "\n",
    "#band names\n",
    "bands = ['B02_10m','B03_10m','B04_10m']\n",
    "\n",
    "#get product file paths according to dates and tile ID T19PEP (covers Bonaire)\n",
    "level2_dir = '...' #change according to preference\n",
    "level2_files = glob(level2_dir+'/*.SAFE')\n",
    "\n",
    "scene_paths=[file for date in dates for file in level2_files if str(date) in file and 'T19PEP' in file]\n",
    "\n",
    "#sort multiband image paths according to date\n",
    "image_collection ={}\n",
    "\n",
    "for scene in scene_paths:\n",
    "    date = re.findall(r\"(\\d{8})T\", scene)[0]\n",
    "    \n",
    "    #collect all .jp2 band images in SAFE directory\n",
    "    all_images = [f for f in glob(scene + \"*/**/*.jp2\", recursive=True)]\n",
    "    img_paths = [img_path for band in bands for img_path in all_images if band in img_path]\n",
    "    image_collection[date] = img_paths\n",
    "\n",
    "#check nr. of images per date\n",
    "for key in image_collection.keys():print(f'Date: {key} Images: {len(image_collection[key])}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#polygon for cropping image\n",
    "roi_file = './data/boundaries/objective3/lacbay_roi.geojson'                   \n",
    "cm_20190128 = './data/boundaries/objective3/cloudmask_20190128.geojson'  \n",
    "cm_20190212 = './data/boundaries/objective3/cloudmask_20190212.geojson' \n",
    "\n",
    "#stack multiband images to a geotiff \n",
    "for date in tqdm(image_collection.keys(),position=0, leave=True):\n",
    "    stack_file = os.path.join(parent_dir,'fullstack',f'stack_{date}.tif')\n",
    "    stack_bands(image_collection[date],image_collection[date][1],stack_file)\n",
    "    \n",
    "    clip_outfile = os.path.join(parent_dir,'clippedstack',f'stack_{date}_clipped.tif')\n",
    "    #crop multiband image\n",
    "    if '20190128' in date:\n",
    "        clip_raster(stack_file,cm_20190128,clip_outfile,fill=True,nodat=0)\n",
    "    elif '20190212' in date:\n",
    "        clip_raster(stack_file,cm_20190212,clip_outfile,fill=True,nodat=0)\n",
    "    else:\n",
    "        clip_raster(stack_file,roi_file,clip_outfile,fill=True,nodat=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Sample pixel values from multiband images based on training sites</b>  \n",
    "* Training scenes from 1 and 28 January and 12 February 2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get training sites and corresponding images\n",
    "train_sites = glob(r\".\\data\\training_input\\objective3\\*_lac.geojson\")    \n",
    "dates = [20190108,20190128,20190212]                                 \n",
    "stack_bands = [f for date in dates for f in glob(parent_dir+'/clipped*/*_clipped.tif') if str(date) in f]  \n",
    "\n",
    "#bands\n",
    "band_names = ['B02','B03','B04']\n",
    "\n",
    "dataset = []\n",
    "for i in range(len(train_sites)):\n",
    "\n",
    "    #sample multibands and spectral indices\n",
    "    df_sample= pixel_sample(stack_bands[i],train_sites[i],band_names)\n",
    "    dataset.append(df_sample)\n",
    "    \n",
    "#final dataset\n",
    "dataset=pd.concat(dataset,sort=False).reset_index(drop=True) \n",
    "dataset.to_csv(r'./data/training_input/csv/training_samples_20190108_20190212_seagrass.csv',index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Expore spectral signature</b>  \n",
    "* Jeffries-Matusita distance (JMD) used for feature selection ([reference](https://books.google.nl/books?id=RxHbb3enITYC&pg=PA52&lpg=PA52&dq=for+one+feature+and+two+classes+the+Bhattacharyya+distance+is+given+by&source=bl&ots=sTKLGl1POo&sig=ACfU3U2s7tv0LT9vfSUat98l4L9_dyUgeg&hl=nl&sa=X&ved=2ahUKEwiKgeHYwI7lAhWIIlAKHZfJAC0Q6AEwBnoECAkQAQ#v=onepage&q&f=false))\n",
    "* RGB (bands 4,3,2) are selected as input features for the classifiers (though worst JMD scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load training sample\n",
    "df = pd.read_csv(r'./data/training_input/csv/training_samples_20190108_20190212_seagrass.csv')\n",
    "\n",
    "#plot JMD heatmap for each band\n",
    "jmd_bands = [jmd2df(transpose_df(df,'C',band)) for band in ['B02','B03','B04']]\n",
    "sns.heatmap(pd.concat(jmd_bands,sort=True),annot=True)\n",
    "ridgePlot(df[['C','B02','B03','B04']],'C')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Build classifiers</b>  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load training sample\n",
    "df = pd.read_csv(r'./data/training_input/csv/training_samples_20190108_20190212_seagrass.csv')\n",
    "subset_df = df[['C','B02','B03','B04']]\n",
    "\n",
    "#split into train and test datasets 80:20\n",
    "train,test = train_test_split(subset_df, train_size = 0.8,random_state=1,shuffle=True,stratify=np.array(subset_df['C']))\n",
    "train = train.sort_values(by='C',ascending=True) #sort labels\n",
    "\n",
    "#split pedictors from labels (for DTC)\n",
    "le = LabelEncoder()\n",
    "X_train,y_train = train[['B02','B03','B04']],le.fit_transform(train['C'])\n",
    "X_test,y_test = test[['B02','B03','B04']],le.fit_transform(test['C'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Decision Tree Classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#perform k-fold (=10) cross-validation \n",
    "\n",
    "#parameters considered in this step\n",
    "max_depth = np.arange(1,40,2)                    \n",
    "min_samples_split = list(range(2, 100,10))         \n",
    "max_leaf_nodes= list(range(2, 50,5))              \n",
    "min_samples_leaf= list(range(1, 100,10))           \n",
    "min_impurity_decrease=[0,0.00005,0.0001,0.0002,0.0005,0.001,0.0015,0.002,0.005,0.01,0.02,0.05,0.08]   \n",
    "criterion = ['gini','entropy']\n",
    "\n",
    "#assign parameters to a dictionary\n",
    "params = {'max_depth':max_depth,'min_samples_split':min_samples_split,\n",
    "          'max_leaf_nodes':max_leaf_nodes,'min_samples_leaf':min_samples_leaf,\n",
    "          'min_impurity_decrease':min_impurity_decrease,'criterion':criterion}\n",
    "\n",
    "#plot validation curve\n",
    "fig,axs = plt.subplots(3,2,figsize=(10,8))\n",
    "axs = axs.ravel()\n",
    "dtc = DecisionTreeClassifier(random_state=1,criterion='entropy')                    #default model\n",
    "\n",
    "for (param_name,param_range),i in zip(params.items(),range(len(params.items()))):\n",
    "    train_scores,test_scores = validation_curve(dtc,X_train.values,y_train,cv=10,scoring='accuracy',\n",
    "                                                n_jobs=-1,param_range=param_range,param_name=param_name)\n",
    "    validation_curve_plot(train_scores,test_scores,param_range,param_name,axs[i])\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#train dtc model based on best parameters\n",
    "dtc = DecisionTreeClassifier(max_depth=5,random_state=42,criterion='entropy',\n",
    "                             min_samples_split=50,max_leaf_nodes=10,min_samples_leaf=30,min_impurity_decrease=0.02)\n",
    "dtc = dtc.fit(X_train,y_train)\n",
    "\n",
    "#export model as joblib file\n",
    "dump(dtc,r\".\\data\\models\\dtc_model_seagrass.joblib\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Maximum Likelihood Classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#train mlc model\n",
    "mlc = mlClassifier(train,'C')\n",
    "\n",
    "#export model as joblib file\n",
    "dump(mlc,r\".\\data\\models\\mlc_model_seagrass.joblib\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Compute model accuracies (based on test split)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load models\n",
    "dtc = load(r\".\\data\\models\\dtc_model_seagrass.joblib\")\n",
    "mlc = load(r\".\\data\\models\\mlc_model_seagrass.joblib\")\n",
    "\n",
    "#DTC model accuracy\n",
    "dtc_y_pred = dtc.predict(X_test)\n",
    "con_mat_dtc = calc_acc(le.inverse_transform(y_test),le.inverse_transform(dtc_y_pred))\n",
    "con_mat_dtc['classifier'] = 'DTC'\n",
    "\n",
    "#MLC model accuracies with chi-square threshold\n",
    "chi_table = {'MLC base':None,'MLC 10%':7.78,'MLC 20%':5.99,'MLC 50%':3.36}\n",
    "\n",
    "mlc_conmats = []\n",
    "for key,value in chi_table.items():\n",
    "    con_mat_mlc = mlc.classify_testdata(test,'C',threshold=value)\n",
    "    con_mat_mlc['classifier'] = key\n",
    "    mlc_conmats.append(con_mat_mlc)\n",
    "\n",
    "#export model accuracies\n",
    "mlc_conmats = pd.concat(mlc_conmats)\n",
    "model_acc = pd.concat([con_mat_dtc,mlc_conmats])\n",
    "model_acc.to_csv('./data/output/objective3/dtc_mlc_model_acc_obj3.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Classification</b>  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load models\n",
    "dtc = load(r\".\\data\\models\\dtc_model_seagrass.joblib\")\n",
    "mlc = load(r\".\\data\\models\\mlc_model_seagrass.joblib\")\n",
    "\n",
    "#output dir\n",
    "os.makedirs(os.path.join(parent_dir,'predicted/dtc'),exist_ok=True)\n",
    "os.makedirs(os.path.join(parent_dir,'predicted/mlc'),exist_ok=True)\n",
    "\n",
    "clipped_files = glob(parent_dir+'/clippedstack/*_clipped.tif')\n",
    "dates= [20190108,20190304,20190821,20191129]\n",
    "clipped_files = [path for path in clipped_files for date in dates if str(date) in path]\n",
    "\n",
    "for file in clipped_files:\n",
    "    date = re.findall(r\"(\\d{8})\", file)[0]\n",
    "    chi_probs = [None,7.78,5.99,3.36]\n",
    "    with rio.open(file) as src:\n",
    "        stack2pred_img = src.read()\n",
    "        mlc_imgs = np.array([mlc.classify_raster_gx(stack2pred_img,threshold=prob) for prob in chi_probs])\n",
    "        dtc_img = np.array([dtc_pred_stack(dtc,stack2pred_img)])\n",
    "        \n",
    "        #export results\n",
    "        mlc_profile = src.profile.copy()\n",
    "        mlc_profile.update({'nodata':None,'dtype':rio.uint16,'count':4})\n",
    "        mlc_out = os.path.join(parent_dir,'predicted/mlc',f'mlc_{date}.tif')\n",
    "        \n",
    "        dtc_profile = src.profile.copy()\n",
    "        dtc_profile.update({'nodata':None,'dtype':rio.uint8,'count':1})\n",
    "        dtc_out = os.path.join(parent_dir,'predicted/dtc',f'dtc_{date}.tif')\n",
    "        \n",
    "        with rio.open(mlc_out,'w',**mlc_profile) as mlc_dst, rio.open(dtc_out,'w',**dtc_profile) as dtc_dst:\n",
    "            mlc_dst.write(mlc_imgs.astype(rio.uint16))\n",
    "            dtc_dst.write(dtc_img.astype(rio.uint8))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>External validity</b>  \n",
    "* Classify DTC and MLC results for a scene taken on 2019-03-04\n",
    "* Seagrass pixel value = 2 in the DTC and MLC rasters "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#get file paths\n",
    "val_samples = gpd.read_file(r'./data/training_input/objective3/sg_validation_2019.geojson')\n",
    "dtc_file = glob(parent_dir+'/predicted*/dtc/dtc*20190304*.tif')[0]\n",
    "mlc_file = glob(parent_dir+'/predicted*/mlc/mlc*20190304*.tif')[0]\n",
    "\n",
    "coords =  [(val_samples.geometry[i][0].x,val_samples.geometry[i][0].y) for i in range(len(val_samples))]\n",
    "\n",
    "with rio.open(dtc_file) as dtc_src, rio.open(mlc_file) as mlc_src:\n",
    "    #sample from dtc raster\n",
    "    val_samples['DTC'] = [pt[0] for pt in dtc_src.sample(coords)]\n",
    "    \n",
    "    #sample from multilayer mlc raster\n",
    "    mlc_multi = pd.concat([pd.DataFrame(pt).T for pt in mlc_src.sample(coords)],ignore_index=True)\n",
    "    val_samples[['MLC base','MLC 10%','MLC 20%','MLC 50%']] = mlc_multi\n",
    "    \n",
    "#convert pixel values to 1 if seagrass, else to 0 for others\n",
    "val_samples[val_samples.columns[-5:]] = (val_samples[val_samples.columns[-5:]]==2).astype(int)\n",
    "val_samples.drop(['site','mean_cover'],axis=1,inplace=True)\n",
    "\n",
    "#compute classification (validation) accuracy \n",
    "df_val = pd.DataFrame(val_samples.drop(columns='geometry'))\n",
    "\n",
    "acc_val_dfs = []\n",
    "for pred in df_val.columns[df_val.columns!='label']:\n",
    "    acc = calc_acc(df_val['label'].values, df_val[pred].values)\n",
    "    acc['classifier'] = pred\n",
    "    acc_val_dfs.append(acc)\n",
    "acc_val_dfs = pd.concat(acc_val_dfs)\n",
    "acc_val_dfs.to_csv('./data/output/objective3/dtc_mlc_external_val_obj3.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "model_df = pd.read_csv('./data/output/objective3/dtc_mlc_model_acc_obj3.csv').set_index('Unnamed: 0')\n",
    "val_df  = pd.read_csv('./data/output/objective3/dtc_mlc_external_val_obj3.csv').set_index('Observed')\n",
    "\n",
    "acc2plot = {'Model accuracy (2 classes)':model_df.loc['PA','UA'].str[:4].astype(float),\n",
    "            'Model F1-score (Sg)':model_df.loc['sg','F1-score'].astype(float),\n",
    "            'Validation accuracy (2 classes)':val_df.loc['PA','UA'].str[:4].astype(float),\n",
    "            'Validation F1-score (Sg)':val_df.loc['1','F1-score'].astype(float)}\n",
    "\n",
    "[plt.plot(val_df['classifier'].unique(),value,label=key) for key,value in acc2plot.items()]\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Comparative analysis</b>  \n",
    "* Compare seagrass (Sg classified area across different scenes for each model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get classification result paths\n",
    "dtc_paths = glob(parent_dir+'/predicted/dtc/dtc*.tif')\n",
    "mlc_paths = glob(parent_dir+'/predicted/mlc/mlc*.tif')\n",
    "\n",
    "data = dict.fromkeys(['Date','Sg MLC Base','Sg MLC 10%','Sg MLC 20%','Sg MLC 50%','Sg DTC'], [])\n",
    "\n",
    "for i in range(len(mlc_paths)):\n",
    "    date = re.findall(r\"(\\d{8})\", mlc_paths[i])\n",
    "    data['Date'] = data['Date']+ [str(pd.to_datetime(date)[0].date())]\n",
    "    \n",
    "    with rio.open(dtc_paths[i]) as dtc_src, rio.open(mlc_paths[i]) as mlc_src:\n",
    "        data['Sg DTC'] = data['Sg DTC'] + [np.unique(dtc_src.read(),return_counts=True)[1][1]]\n",
    "        for k,sf_mlc_key in enumerate(list(data.keys())[1:-1]):\n",
    "            data[sf_mlc_key] = data[sf_mlc_key]+ [np.unique(mlc_src.read([k+1]), return_counts=True)[1][1]]\n",
    "            \n",
    "#export data\n",
    "data = pd.DataFrame(data)\n",
    "data.to_csv('./data/output/objective3/classified_area_obj3.csv',index=False)            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Plot seagrass classified area in 2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load data and subset only the 2019 results\n",
    "data = pd.read_csv('./data/output/objective3/classified_area_obj3.csv',index_col='Date')\n",
    "\n",
    "#plot seagrass classified area in Lac Bay\n",
    "plt.ylabel('Classified area (ha)')\n",
    "plt.plot(data/100)\n",
    "plt.legend(data.columns,loc='upper left')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}