02_sargassum_detection_sea.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Detection of Sargassum in open sea</b>  \n",
    "Notebook for classifying and analyzing Sargassum in Bonaire with Sentinel-2 images\n",
    "\n",
    "* Density slicing (GNDVI),PCA, Decision Tree Classifier (DTC) and Maximum Likelihood Classifier (MLC) are employed\n",
    "* 7 different subsets in open sea are investigated\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",
    "from rasterio.merge import merge\n",
    "from rasterio.features import sieve\n",
    "import matplotlib.pyplot as plt\n",
    "from glob import glob\n",
    "import geopandas as gpd\n",
    "from joblib import load\n",
    "from tqdm import tqdm,tqdm_notebook\n",
    "from tqdm.contrib import tzip\n",
    "import rasterio.features \n",
    "from shapely.geometry import shape\n",
    "from shapely.geometry.multipolygon import MultiPolygon\n",
    "\n",
    "#custom functions\n",
    "from Python.prep_raster import stack_bands,clip_raster,pixel_sample,computeIndexStack,compute_index\n",
    "from Python.pca_slice import compute_pca_score,compute_pca_image,pca_triangle,density_slice\n",
    "from Python.mlc import mlClassifier\n",
    "from Python.pred_raster import stack2pred, dtc_pred_stack\n",
    "from Python.misc import get_feat_layer_order\n",
    "\n",
    "#setup IO directories\n",
    "parent_dir = os.path.join(os.path.abspath('..'),'objective2')                  #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": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#dates considered for classification and analysis \n",
    "dates_tiles = {\"T19PEP\":[20180304, 20190309],\n",
    "               \"T19PFP\":[20180304, 20190304, 20190428]}\n",
    "\n",
    "#band names\n",
    "bands = ['B01_60m','B02_10m','B03_10m','B04_10m','B05_20m','B06_20m',\n",
    "         'B07_20m','B08_10m','B8A_20m','B09_60m','B11_20m','B12_20m']\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 = []\n",
    "for file in level2_files:\n",
    "    for key,value in dates_tiles.items():\n",
    "        if key in file and any(str(v) in file for v in value):\n",
    "            scene_paths.append(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",
    "    tileid = re.findall(r\"(T\\d{2}P..)\", 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[f'{tileid}_{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": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#stack multiband images to a geotiff (!computationaly intensive)\n",
    "for date in tqdm(image_collection.keys(),position=0, leave=True):\n",
    "    ref10m= image_collection[date][1]                                   #use band B02 (10m) as reference metadata\n",
    "    outfile = os.path.join(parent_dir,'fullstack',f'stack_{date}.tif')\n",
    "    stack_bands(image_collection[date],ref10m,outfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#crop multiband image stack \n",
    "roi_file = gpd.read_file('./data/boundaries/objective2/sea_rois.geojson') #polygon for cropping image\n",
    "\n",
    "sub_geoms = dict(zip(roi_file['subset'],\n",
    "                     [[roi_file.__geo_interface__['features'][i]['geometry']] for i in range(len(roi_file))]))\n",
    "\n",
    "stack_files = glob(parent_dir + \"/fullstack/*.tif\")\n",
    "clip_dir = os.path.join(parent_dir,'clippedstack')\n",
    "\n",
    "for file in tqdm(stack_files,position=0, leave=True):\n",
    "    date = re.findall(r\"(\\d{8})\",file)[0]\n",
    "    subset = os.path.basename(file)\n",
    "    if 'T19PEP_20180304' in subset:\n",
    "        clip_raster(file,sub_geoms['20180304_A'],\n",
    "                    os.path.join(clip_dir,subset.replace(date,'20180304_A_clipped')),fill=True,nodat=0)\n",
    "        clip_raster(file,sub_geoms['20180304_B'],\n",
    "                    os.path.join(clip_dir,subset.replace(date,'20180304_B_clipped')),fill=True,nodat=0)\n",
    "    elif 'T19PEP_20190309' in subset:\n",
    "        clip_raster(file,sub_geoms['20190309'],\n",
    "                    os.path.join(clip_dir,subset.replace(date,'20190309_clipped')),fill=True,nodat=0)\n",
    "    elif 'T19PFP_20180304' in subset:\n",
    "        clip_raster(file,sub_geoms['20180304_C'],\n",
    "                    os.path.join(clip_dir,subset.replace(date,'20180304_C_clipped')),fill=True,nodat=0)\n",
    "    elif 'T19PFP_20190304' in subset:\n",
    "        clip_raster(file,sub_geoms['20190304'],\n",
    "                    os.path.join(clip_dir,subset.replace(date,'20190304_clipped')),fill=True,nodat=0)\n",
    "    elif 'T19PFP_20190428' in subset:\n",
    "        clip_raster(file,sub_geoms['20190428_A'],\n",
    "                    os.path.join(clip_dir,subset.replace(date,'20190428_A_clipped')),fill=True,nodat=0)\n",
    "        clip_raster(file,sub_geoms['20190428_B'],\n",
    "                    os.path.join(clip_dir,subset.replace(date,'20190428_B_clipped')),fill=True,nodat=0)\n",
    "        \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>PCA</b>\n",
    "* Sparsely distributed Sargassum slicks (aka windrows) in rough waters are more visible on PCA 3-5\n",
    "* Large seaweed at dark calm waters are visible on PCA 1-2 \n",
    "* Additional stripe/ detector artifacts may cause false/ noisy pixels \n",
    "* PCA performance depends on the scene and image scale "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#crop pca border (fix)\n",
    "roi_file = gpd.read_file('./data/boundaries/objective2/sea_rois.geojson') #polygon for cropping image\n",
    "sub_geoms = dict(zip(roi_file['subset'],\n",
    "                     [[roi_file.__geo_interface__['features'][i]['geometry']] for i in range(len(roi_file))]))\n",
    "\n",
    "clip_files = glob(parent_dir+'/clippedstack/*_clipped.tif' )\n",
    "\n",
    "#setup pca parameters and dir\n",
    "os.makedirs(os.path.join(parent_dir,'predicted\\\\pca'),exist_ok=True)\n",
    "bands = ['B02','B03','B04','B05','B06','B07','B08','B8A','B11','B12']\n",
    "band_order = [2,3,4,5,6,7,8,9,11,12]\n",
    "\n",
    "#collect pca score\n",
    "pca_score_dfs = []\n",
    "\n",
    "for file in tqdm(clip_files,position=0, leave=True):\n",
    "    subset = re.search(r'stack_(.*?)_clipped', os.path.basename(file)).group(1)\n",
    "    \n",
    "    #compute pca raster with components = 5, also crop to fix nodata boundary\n",
    "    outfile = os.path.join(parent_dir,'predicted/pca',f'pca_{subset}_x.tif')\n",
    "    compute_pca_image(file,band_order,nr_components=5,outfile=outfile)\n",
    "    clip_raster(outfile,sub_geoms[subset[7:]],outfile.replace('_x',\"\"),fill=True,nodat=0)\n",
    "    os.remove(outfile)\n",
    "    \n",
    "    #compute pca score and variance\n",
    "    pca_score_dfs.append(compute_pca_score(file,band_order,bands,nr_components=5,subset))\n",
    "\n",
    "#export pca score and variance\n",
    "pd.concat(pca_score_dfs,axis=0).to_csv(r'./data/output/objective2/pca_score_obj2.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Triangle thresholding for unimodal image segmentation ([reference](https://imagej.net/Auto_Threshold))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#get PC images and create new output dir\n",
    "pca_files = glob(parent_dir+'/predicted/pca/pca_*.tif' )\n",
    "pc_sel = dict(zip(pca_files,[1,1,5,2,4,4,5]))\n",
    "\n",
    "#perform triangle thresholding \n",
    "for file,pc in pc_sel.items():\n",
    "    subset = re.search(r'pca_(.*?).tif', os.path.basename(file)).group(1)\n",
    "    outfile = os.path.join(parent_dir,'predicted/pca',f'pca2_{subset}.tif')\n",
    "    pca_triangle(file,pc,outfile)\n",
    "    os.remove(file) #remove this line if you like to keep the original component images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Density Slicing of GNDVI image </b>\n",
    "* GNDVI value between -0.25 and -0.30 suitable for detecting Sargassum patches in calm dark waters (no clouds)\n",
    "* GNDVI value > 0.05 more suitable for noisy scenes (rough waters/ cloudy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get clipped rasters and create output dir\n",
    "clip_files = glob(parent_dir+'/clippedstack/*_clipped.tif')\n",
    "os.makedirs(os.path.join(parent_dir,'predicted/slice/'),exist_ok=True)\n",
    "\n",
    "index_sel = dict(zip(clip_files,[-0.25,-0.3]+[0.05]*5))\n",
    "\n",
    "for clip_file,thr in index_sel.items():\n",
    "    \n",
    "    #create GNDVI raster\n",
    "    subset = re.search(r'stack_(.*?)_clipped.tif', os.path.basename(clip_file)).group(1)\n",
    "    gndvi_outfile = os.path.join(parent_dir,f'indexstack/gndvi_{subset}.tif')\n",
    "    computeIndexStack(clip_file,['GNDVI'],gndvi_outfile)\n",
    "    \n",
    "    #perform density slicing\n",
    "    slice_outfile = os.path.join(parent_dir,f'predicted/slice/slice_{subset}.tif')\n",
    "    density_slice(gndvi_outfile,1,thr,slice_outfile)\n",
    "    \n",
    "    os.remove(gndvi_outfile) #remove this line if you like to keep the original GNDVI image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>DTC and MLC classification</b>\n",
    "* Only the MLC base model (no threshold) is used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load models\n",
    "dtc = load(r\".\\data\\models\\dtc_model_sargassum.joblib\")\n",
    "mlc = load(r\".\\data\\models\\mlc_model_sargassum.joblib\")\n",
    "\n",
    "make_dirs = [os.makedirs(os.path.join(parent_dir,f'predicted/{name}'),exist_ok=True) for name in ['dtc','mlc']]\n",
    "\n",
    "#clipped files\n",
    "clip_files = glob(parent_dir+'/clippedstack/*_clipped.tif')\n",
    "\n",
    "for file in tqdm(clip_files,position=0, leave=True):\n",
    "    subset = re.search(r'stack_(.*?)_clipped.tif',file).group(1)\n",
    "\n",
    "    with rio.open(file) as src:\n",
    "        #DTC and MLC classifications\n",
    "        stack2pred_img = np.concatenate((computeIndexStack(src.read(),['NDVI','REP']), src.read([5,11])))\n",
    "        mlc_img = np.array([mlc.classify_raster_gx(stack2pred_img)])\n",
    "        dtc_img = np.array([dtc_pred_stack(dtc,stack2pred_img)])\n",
    "        \n",
    "        #export results\n",
    "        profile = src.profile.copy()\n",
    "        profile.update({'nodata':0,'dtype':rio.uint8,'count':1})\n",
    "        dtc_out,mlc_out = f'{parent_dir}/predicted/dtc/dtc_{subset}.tif',f'{parent_dir}/predicted/mlc/mlc_{subset}.tif'\n",
    "        with rio.open(dtc_out,'w',**profile) as dtc_dst, rio.open(mlc_out,'w',**profile) as mlc_dst:\n",
    "            dtc_dst.write(dtc_img.astype(rio.uint8))\n",
    "            mlc_dst.write(mlc_img.astype(rio.uint8))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Comparative analysis</b>  \n",
    "* Compare Sargassum classified area in open sea across classifications (PCA, density slicing, DTC, MLC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "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",
    "pca_paths = glob(parent_dir+'/predicted*/pca/pca2*.tif')\n",
    "slice_paths = glob(parent_dir+'/predicted*/slice/slice*.tif')\n",
    "\n",
    "#collect Sargassum pixel count \n",
    "data = {}\n",
    "for dtc_file,mlc_file,pca_file,slice_file in zip(dtc_paths,mlc_paths,pca_paths,slice_paths):\n",
    "    \n",
    "    data.setdefault('Subset',[]).append(re.search(r'dtc_(.*?).tif', dtc_file).group(1))\n",
    "    \n",
    "    with rio.open(dtc_file) as dtc_src, rio.open(mlc_file) as mlc_src,rio.open(pca_file) as pca_src, rio.open(slice_file) as slice_src:\n",
    "\n",
    "        data.setdefault('DTC',[]).append(np.count_nonzero(np.where(dtc_src.read(1)==3,1,0)))\n",
    "        data.setdefault('MLC Base',[]).append(np.count_nonzero(np.where(mlc_src.read(1)==3,1,0)))\n",
    "        data.setdefault('PCA Triangle',[]).append(np.count_nonzero(pca_src.read(1)))\n",
    "        data.setdefault('GNDVI Slice',[]).append(np.count_nonzero(slice_src.read(1))) \n",
    "        \n",
    "#export data\n",
    "pd.DataFrame(data).to_csv('./data/output/objective2/sargassum_count_obj2.csv',index=False)   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Plot Sargassum classified area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load data and subset only the 2019 results\n",
    "data = pd.read_csv('./data/output/objective2/sargassum_count_obj2.csv',index_col='Subset')\n",
    "data = data.T[[data.T.columns[i] for i in [0,1,3,4,2,5,6]]].T \n",
    "\n",
    "#plot Sargassum \n",
    "plots = [plt.plot(data[col]/100,label=col) for col in data.columns]\n",
    "plt.ylabel('Classified area (ha)')\n",
    "plt.xticks(rotation=45)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Large scene classification</b>  \n",
    "* Classify a larger scale image with the DTC, MLC, PCA and density slicing methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#crop roi file and create output dir\n",
    "roi_file = './data/boundaries/objective2/sea_extent_small.geojson'\n",
    "os.makedirs(os.path.join(parent_dir,'predicted/large_scene'),exist_ok=True)\n",
    "\n",
    "#load models\n",
    "dtc = load(r\".\\data\\models\\dtc_model_sargassum.joblib\")\n",
    "mlc = load(r\".\\data\\models\\mlc_model_sargassum.joblib\")\n",
    "\n",
    "stack_files = glob(parent_dir+'/fullstack/*20180304.tif')\n",
    "\n",
    "for file,pc in tzip(stack_files,[3,2]): \n",
    "    scene_name = re.search(r'stack_(.*?).tif',file).group(1)\n",
    "    outfile = os.path.join(parent_dir,'predicted/large_scene',f'{scene_name}_clipped.tif')\n",
    "    clip_raster(file,roi_file,outfile,fill=True)\n",
    "        \n",
    "    if os.path.exists(outfile):\n",
    "        \n",
    "        with rio.open(outfile) as src:\n",
    "            \n",
    "            #DTC and MLC classifications\n",
    "            stack2pred_img = np.concatenate((computeIndexStack(src.read(),['NDVI','REP']), src.read([5,11])))\n",
    "            mlc_img = np.where(np.array([mlc.classify_raster_gx(stack2pred_img)])==3,1,0)\n",
    "            dtc_img = np.where(np.array([dtc_pred_stack(dtc,stack2pred_img)])==3,1,0)\n",
    "            \n",
    "            #PCA and GNDVI density slicing classifications\n",
    "            pca_img = np.array([pca_triangle(compute_pca_image(outfile,[2,3,4,5,6,7,8,9,11,12],nr_components=pc),pc-1)])\n",
    "            slice_img = np.array([np.where(compute_index(src.read(),'GNDVI')>=0,1,0)])\n",
    "            \n",
    "            #stack all results and apply cloud mask\n",
    "            results = np.concatenate((mlc_img,dtc_img,pca_img,slice_img))\n",
    "            cloud_mask = np.where(src.read([2])>=0.09,1,0)\n",
    "            results_masked = np.where(cloud_mask!=1,results,0)\n",
    "\n",
    "            #export results\n",
    "            profile = src.profile.copy()\n",
    "            profile.update({'nodata':None,'dtype':rio.uint16,'count':4})\n",
    "            \n",
    "            results_out = outfile.replace('_clipped','_multi')\n",
    "            with rio.open(results_out,'w',**profile) as dst:\n",
    "                dst.write(results_masked.astype(rio.uint16))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Merge classification results over two tiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img_files = glob(parent_dir+'/predicted/large_scene/T*_multi.tif')\n",
    "\n",
    "#open image files and merge them\n",
    "data = [rio.open(file) for file in img_files]\n",
    "merged_data,data_affine = merge(data)\n",
    "\n",
    "#sieve filter (spatial filter)\n",
    "sieved_data= []\n",
    "for img in merged_data:\n",
    "    sieved_data.append(sieve(img,size=8))\n",
    "sieved_data = np.array(sieved_data)\n",
    "\n",
    "#export merged and sieve image\n",
    "profile = data[0].profile.copy()\n",
    "profile.update({'transform':data_affine,'width':sieved_data.shape[2],'height':sieved_data.shape[1],\n",
    "                'nodata':0,'dtype':rio.uint16,'count':4})\n",
    "outfile = os.path.join(parent_dir,'predicted/large_scene/merged_20180304_multi.tif')\n",
    "with rio.open(outfile,'w',**profile) as dst:\n",
    "    dst.write(sieved_data.astype(rio.uint16))  \n",
    "    \n",
    "#close image files\n",
    "list(map(lambda x:x.close(),data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Convert merged results into vector data to reduce data size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#merged image\n",
    "img_file = glob(parent_dir+'/predicted/large_scene/merged*_multi.tif')[0]\n",
    "\n",
    "with rasterio.open(img_file) as src:\n",
    "    imgs = src.read()\n",
    "    data = []\n",
    "    \n",
    "    #iterate over the results in the multiband image\n",
    "    for img,method in tzip(imgs,['MLC Base','DTC','PCA Triangle','GNDVI Slice']):\n",
    "    \n",
    "        #get coordinates of every pixel and filter only those with value =1 (Sargassum)\n",
    "        shapes = list(rasterio.features.shapes(img, transform=src.transform))\n",
    "        shapes = list(filter(lambda x: x[1]==1, shapes))\n",
    "        \n",
    "        #convert into geodataframe\n",
    "        multipolygon = [MultiPolygon([shape(geom[0]) for geom in shapes])]\n",
    "        gdf = gpd.GeoDataFrame({'Result': [method]},geometry=multipolygon,crs=src.crs) \n",
    "        data.append(gdf)\n",
    "        \n",
    "    #export polygonize features as geojson\n",
    "    data = pd.concat(data,ignore_index=True)\n",
    "    data.to_file('./data/output/objective2/4classifications_20180304.geojson',driver='GeoJSON')"
   ]
  },
  {
   "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
}