3_Land_Cover_Classification.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "03ef128d-6251-4df2-88ad-8b06e57d8ca2",
   "metadata": {},
   "source": [
    "# Land Cover Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12b36467-a13f-44e4-83eb-212d6e5e97a0",
   "metadata": {},
   "source": [
    "## Description\n",
    "\n",
    "Having succesfully run the `2_Evaluate_Fit_classifier.ipynb` notebook, we can now use our classification model to predict values on new satellite data. This notebook will guide you through loading satellite data from the ODC, computing the same feature layers as we did in the first notebook when we extracted training data from the ODC, and using our model to classify the satellite data.  We will classify a region of interest, visualize how well our model is performing, and save the results to disk as a Cloud-Optimized GeoTIFF (COG).\n",
    "\n",
    "The steps are as follows:\n",
    "1. Import the model we output in the previous notebook, `2_Evaluate_Fit_classifier.ipynb`\n",
    "2. Redefine the feature layer function that we used to extract training data from the ODC in the first notebook, `0_Extract_Training_Features.ipynb`\n",
    "3. Loop through a set of locations to extract satellite data from the ODC for a selected region, then compute the feature layers and classify the data using our model\n",
    "4. Predict land cover for the a user selected district\n",
    "5. Visualise and save the classification results\n",
    "\n",
    "To run this analysis, run all the cells in the notebook, starting with the \"Load packages\" cell."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54da22e6-867b-4ea9-b47e-b23f088b3273",
   "metadata": {},
   "source": [
    "### Load Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62ec3cc7-47d6-4194-bb5c-40d6fad0f0ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import datacube\n",
    "import warnings\n",
    "import time\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from odc.algo import xr_geomedian\n",
    "import xarray as xr\n",
    "from joblib import load\n",
    "from deafrica_tools.classification import predict_xr\n",
    "from deafrica_tools.dask import create_local_dask_cluster\n",
    "from deafrica_tools.datahandling import load_ard\n",
    "from deafrica_tools.bandindices import calculate_indices\n",
    "from deafrica_tools.plotting import rgb, display_map\n",
    "from datacube.utils.cog import write_cog\n",
    "from matplotlib.colors import ListedColormap,BoundaryNorm\n",
    "from matplotlib.patches import Patch\n",
    "import geopandas as gpd\n",
    "from datacube.utils import geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6aac7805-f36b-4221-889b-651e439a3ef6",
   "metadata": {},
   "source": [
    "### Set up a dask cluster\n",
    "This will help keep our memory use down and conduct the analysis in parallel. If you'd like to view the dask dashboard, click on the hyperlink that prints below the cell. You can use the dashboard to monitor the progress of calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e56e268-e7bc-48df-b38a-90a9a0a50110",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up a dask cluster\n",
    "create_local_dask_cluster()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cefa4d8-a035-4b50-b176-3b78c072aa0a",
   "metadata": {},
   "source": [
    "## Analysis parameters\n",
    "\n",
    "* `rf_model_path`: The path to the location where the model exported from the previous notebook is stored.\n",
    "* `testing_locations`: A dictionary with values containing latitude and longitude points, and keys representing a unique ID to identify the locations. The `lat` and `lon` points define the center of the satellite images we will load for running small test classifications.\n",
    "* `buffer`: The size, in decimal degrees, to load around the central latitude and longitude points in `locations`. This number here will depend on the compute/RAM available on the Sandbox instance, and the type and number of feature layers calculated.  \n",
    "* `input_polygons_path`: The path to the shapefile containing polygons of Rwanda administrative boundaries.  \n",
    "* `district_name_attribute`: This is the name of column in your shapefile attribute table that identifies level 2 names, i.e. district names.  \n",
    "* `district_name`: This is the district name that you selected for prediction.  \n",
    "* `dask_chunks`: Dask works by breaking up large datasets into chunks, which can be read individually. This parameter specifies the size of the chunks in numbers of pixels, e.g. `{'x':200,'y':200}`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b00bc773-17d9-49d5-b866-4162a63d7c4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "rf_model_path='Results/Mozambique_RF_model.joblib'\n",
    "testing_locations = {\n",
    "    '1': (-17.869, 36.855),\n",
    "    '2': (-17.827, 36.939),\n",
    "    '3': (-17.888, 36.919),\n",
    "}\n",
    "buffer = 0.01\n",
    "input_polygons_path = 'Data/moz_admbnda_adm2_ine_20190607.shp'\n",
    "district_name_attribute = 'ADM2_PT'\n",
    "district_name='Cidade De Quelimane'\n",
    "dask_chunks = {'x':1500,'y':1500}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8622f85-e723-41bc-92f1-060d0f9acec7",
   "metadata": {},
   "source": [
    "## View the selected location\n",
    "The next cell will display the first location as example on an interactive map. Feel free to zoom in and out to get a better understanding of the area you’ll be predicting. Clicking on any point of the map will reveal the latitude and longitude coordinates of that point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c058d650-7852-4af3-b575-6795b7cf8cff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# choose the fist testing location\n",
    "testing_location=testing_locations['1']\n",
    "\n",
    "# define spatial extent\n",
    "lon_range=(testing_location[1]-buffer, testing_location[1]+buffer)\n",
    "lat_range=(testing_location[0]+buffer, testing_location[0]-buffer)\n",
    "\n",
    "# display basemap\n",
    "display_map(x=lon_range, y=lat_range)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3778873-9e8f-4def-aebe-a658ba29d3b6",
   "metadata": {},
   "source": [
    "## Import the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a433c9ed-f4a6-43a6-bf0a-ac7edb203821",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load trained classifier\n",
    "model = load(rf_model_path).set_params(n_jobs=1)\n",
    "print('loaded random forest model:\\n',model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebf2d6f1-1b62-4425-a261-c8cbe30f8b46",
   "metadata": {},
   "source": [
    "## Making a prediction\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62bd1d76-0e03-462a-b429-00b2db246d67",
   "metadata": {},
   "source": [
    "### Redefining the feature layer function\n",
    "\n",
    "Because we elected to use all the features extracted in `0_Extract_Training_Features.ipynb`, we can simply copy-and-paste the `feature_layers` function from the [notebook](0_Extract_Training_Features.ipynb) into the cell below (this has already been done for you). \n",
    "\n",
    "If you're using this notebook to run your own classifications (i.e. not running the default example), then you'll need to redefine your feature layer function below, taking care to match the features in the trained model. For example, if you conducted feature selection and removed features from the model, then you'll need to mimic that process here by removing features in the prediction data. In short, the features in the model must precisely match those in the data you're classifying."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43b6f9ed-9657-4c9e-99f6-97d7ec7a2d89",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define a function to feature layers\n",
    "def feature_layers(query): \n",
    "    # connect to the datacube so we can access DE Africa data\n",
    "    dc = datacube.Datacube(app='feature_layers')\n",
    "    \n",
    "    # load Sentinel-2 analysis ready data\n",
    "    ds = load_ard(dc=dc,\n",
    "                  products=['s2_l2a'],\n",
    "                  group_by='solar_day',\n",
    "                  verbose=False,\n",
    "                  **query)\n",
    "    \n",
    "    # calculate NDVI\n",
    "    ds = calculate_indices(ds,\n",
    "                           index=['NDVI'],\n",
    "                           drop=False,\n",
    "                           satellite_mission='s2')\n",
    "    \n",
    "    # interpolate nodata using mean of previous and next observation\n",
    "#     ds=ds.interpolate_na(dim='time',method='linear',use_coordinate=False,fill_value='extrapolate')\n",
    "\n",
    "    # calculate bi-monthly geomedian\n",
    "    ds=ds.resample(time='2MS').map(xr_geomedian)\n",
    "    \n",
    "    # stack multi-temporal measurements and rename them\n",
    "    n_time=ds.dims['time']\n",
    "    list_measurements=list(ds.keys())\n",
    "    list_stack_measures=[]\n",
    "    for j in range(len(list_measurements)):\n",
    "        for k in range(n_time):\n",
    "            variable_name=list_measurements[j]+'_'+str(k)\n",
    "            measure_single=ds[list_measurements[j]].isel(time=k).rename(variable_name)\n",
    "            list_stack_measures.append(measure_single)\n",
    "    ds_stacked=xr.merge(list_stack_measures,compat='override')\n",
    "    return ds_stacked"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97e18412-2c9a-4593-b602-a935912317ce",
   "metadata": {},
   "source": [
    "### Set up datacube query\n",
    "\n",
    "These query options should match the query params in `0_Extract_Training_Features.ipynb`, unless there are measurements that no longer need to be loaded because they were dropped during a feature selection process (which has not been done in the default example)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39287909-be32-4579-9498-9ad21ef129a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "time = ('2021')\n",
    "measurements = ['blue','green','red','red_edge_1','red_edge_2', 'red_edge_3','nir_1','nir_2','swir_1','swir_2']\n",
    "resolution = (-10,10)\n",
    "output_crs='epsg:32736' # WGS84/UTM Zone 36S"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a4e2dce-7961-4913-8730-1abb25dc0152",
   "metadata": {},
   "source": [
    "### Loop through test locations and predict\n",
    "\n",
    "For every location we listed in the `test_locations` dictionary, we calculate the feature layers, and then use the DE Africa function [predict_xr](https://github.com/digitalearthafrica/deafrica-sandbox-notebooks/blob/83116e80ebb4f8744e3de74e7a713aadd0a7577a/Tools/deafrica_tools/classification.py#L237) to classify the data.\n",
    "\n",
    "The `predict_xr` function is an xarray wrapper around the sklearn estimator `.predict()` and `.predict_proba()` methods, and relies on [dask-ml](https://ml.dask.org/) [ParallelPostfit](https://ml.dask.org/modules/generated/dask_ml.wrappers.ParallelPostFit.html) to run the predictions with dask. `Predict_xr` can compute predictions, prediction probabilites, and return the input feature layers. Read the [documentation](https://github.com/digitalearthafrica/deafrica-sandbox-notebooks/blob/83116e80ebb4f8744e3de74e7a713aadd0a7577a/Tools/deafrica_tools/classification.py#L247) for more insights into this function's capabilities.\n",
    "\n",
    "As the feature number is large and the geomedians can take some time to calculate, you may expect a few minutes for the prediction to complete in this example. While waiting, you can check the dask dashboard to view the progress of the processing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8fbe7d3-0128-4593-a179-1c20b5ac9168",
   "metadata": {},
   "outputs": [],
   "source": [
    "predictions = []\n",
    "\n",
    "for key, value in testing_locations.items():\n",
    "    \n",
    "    print('working on: ' + key)\n",
    "\n",
    "     # generate a datacube query object\n",
    "    query = {\n",
    "        'x': (value[1]-buffer, value[1]+buffer),\n",
    "        'y': (value[0]+buffer, value[0]-buffer),\n",
    "        'time': time,\n",
    "        'measurements': measurements,\n",
    "        'resolution': resolution,\n",
    "        'output_crs': output_crs,\n",
    "        'dask_chunks' : dask_chunks,\n",
    "    }\n",
    "\n",
    "    #calculate features\n",
    "    data = feature_layers(query)\n",
    "\n",
    "    #predict using the imported model\n",
    "    predicted = predict_xr(model,\n",
    "                           data,\n",
    "                           proba=True,\n",
    "                           persist=False,\n",
    "                           clean=True,\n",
    "                           return_input=True\n",
    "                          ).compute()\n",
    "        \n",
    "    predictions.append(predicted)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67088c1b-edea-49db-9dc9-3016ee2a5a96",
   "metadata": {},
   "source": [
    "### Plotting results\n",
    "\n",
    "In the plots below you'll see on the left a true-colour image of the region and on the middle the classified image. We've pre-defined a set of colours for each class value through the `colours` dictionary but you can change them as you like. We also display here the prediction probability on the right, where the higher the values indicate higher confidence on the prediction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a68e62ca-8cab-4c27-b929-ed8bebbe7eec",
   "metadata": {},
   "outputs": [],
   "source": [
    "dict_map={11:'Tree crops',12:'Field crops',21:'Forest plantations',31:'Grassland',\n",
    "                 41:'Aquatic or regularly flooded herbaceous vegetation',44:'Water body',\n",
    "                 51:'Settlements',61:'Bare soils',70:'Mangrove',71:'Mecrusse',\n",
    "                72:'Broadleaved (Semi-) evergreen forest',74:'Broadleaved (Semi-) deciduous forest',75:'Mopane'} # dictionary of class names corresponding to pixel values\n",
    "colours = {11:'gold', 12:'yellow', 21:'darkred',31:'bisque',41:'lightgreen',44:'blue',51:'black',61:'gray',\n",
    "           70:'red',71:'steelblue',72:'blueviolet',74:'green',75:'chocolate'}\n",
    "for i in range(0, len(predictions)):\n",
    "    fig, axes = plt.subplots(1, 3, figsize=(24, 8))\n",
    "\n",
    "    # Plot true colour image\n",
    "    rgb(predictions[i], bands=['red_3','green_3','blue_3'],\n",
    "        ax=axes[0], percentile_stretch=(0.01, 0.99))\n",
    "\n",
    "    # Plot classified image\n",
    "    prediction_values=np.unique(predictions[i].Predictions)\n",
    "    cmap=ListedColormap([colours[k] for k in prediction_values])\n",
    "    norm = BoundaryNorm(list(prediction_values)+[np.max(prediction_values)+1], cmap.N)\n",
    "    predictions[i].Predictions.plot.imshow(ax=axes[1],\n",
    "                                    cmap=cmap,\n",
    "                                    norm=norm,\n",
    "                                    add_labels=True, \n",
    "                                    add_colorbar=False,\n",
    "                                    interpolation='none')\n",
    "    # Plot probability\n",
    "    predictions[i].Probabilities.plot.imshow(ax=axes[2],\n",
    "                                      cmap='RdYlGn',\n",
    "                                    add_labels=True,\n",
    "                                      add_colorbar=True)\n",
    "    # Remove axis on middle and right plot\n",
    "    axes[1].get_yaxis().set_visible(False)\n",
    "    axes[2].get_yaxis().set_visible(False)\n",
    "    # add colour legend\n",
    "    patches_list=[Patch(facecolor=colour) for colour in colours.values()]\n",
    "    axes[1].legend(patches_list, list(dict_map.values()),\n",
    "        loc='upper center', ncol =3, bbox_to_anchor=(0.5, -0.1))\n",
    "    # Add plot titles\n",
    "    axes[0].set_title('True Colour Geomedian')\n",
    "    axes[1].set_title('Classified Image')\n",
    "    axes[2].set_title('Prediction Probability')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea7a27ad-3c34-446d-9760-01e7d202b0e0",
   "metadata": {},
   "source": [
    "The areas with low probabilities indicate locations where the model has low confidence on the predictions. One potential way of using the probability maps in your case studies is to filter out predictions lower than a pre-set threshold."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91c65e7c-bbee-4eb9-8d83-70eb36cc15e2",
   "metadata": {},
   "source": [
    "## Classification on the entire district\n",
    "\n",
    "To produce the classification map for the entire district, we only need to update the query to the extent of the district boundary and then run the prediction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e41de225-5392-4385-9459-71d86aa1b49a",
   "metadata": {},
   "outputs": [],
   "source": [
    "polygons=gpd.read_file(input_polygons_path)\n",
    "polygon=polygons.loc[polygons[district_name_attribute]==district_name].iloc[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d729f80-1573-4529-b6f2-9727cb137893",
   "metadata": {},
   "outputs": [],
   "source": [
    "geom=geometry.Geometry(geom=polygon.geometry, crs=polygons.crs)\n",
    "query = {\n",
    "    'geopolygon': geom,\n",
    "    'time': time,\n",
    "    'measurements': measurements,\n",
    "    'resolution': resolution,\n",
    "    'output_crs': output_crs,\n",
    "    'dask_chunks' : dask_chunks,\n",
    "}\n",
    "data = feature_layers(query)\n",
    "#predict using the imported model\n",
    "predicted = predict_xr(model,\n",
    "                       data,\n",
    "                       proba=True,\n",
    "                       persist=False,\n",
    "                       clean=True,\n",
    "                       return_input=True\n",
    "                      ).compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42241553-cb2b-44a9-9e97-65385604d7ab",
   "metadata": {},
   "source": [
    "**Note**: Depending on which sandbox instance you are using, the prediciton may take up to a few hours to finish, e.g. more than three hours if you are using the 32GB sandbox instance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39db5452-ea2d-4c02-b37d-81049d179576",
   "metadata": {},
   "source": [
    "In real projects, you may need to produce the land cover map at a national or even global scale. It is suggested that you divide your large study area into smaller tiles/regions so that the computation for a single tile can fit in the RAM limit. Using the spatial extent of each tile, you can adjust the query and loop through all the tiles to produce the classification maps. You may need to adjust the `dask_chunks` size to optimize for the larger region. While you adjust your region size and dask chunk sizes, try to control the size to avoid the computation running out of the RAM limit of your sandbox instance, which will interrupt the calculation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3601e796-fd22-44bc-b37d-706ba1690378",
   "metadata": {},
   "source": [
    "## Display the predicted map of the entire district"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30ec52a0-6516-468c-91b3-20efe7892d1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3,figsize=(30,8))\n",
    "\n",
    "# Plot true colour image\n",
    "rgb(predicted, bands=['red_3','green_3','blue_3'],\n",
    "    ax=axes[0], percentile_stretch=(0.01, 0.99))\n",
    "\n",
    "# Plot classified image\n",
    "prediction_values=np.unique(predicted.Predictions)\n",
    "cmap=ListedColormap([colours[k] for k in prediction_values])\n",
    "norm = BoundaryNorm(list(prediction_values)+[np.max(prediction_values)+1], cmap.N)\n",
    "predicted.Predictions.plot.imshow(ax=axes[1],\n",
    "                                cmap=cmap,\n",
    "                                norm=norm,\n",
    "                                add_labels=True, \n",
    "                                add_colorbar=False,\n",
    "                                interpolation='none')\n",
    "\n",
    "# Plot probability\n",
    "predicted.Probabilities.plot.imshow(ax=axes[2],\n",
    "                                  cmap='RdYlGn',\n",
    "                                add_labels=True,\n",
    "                                  add_colorbar=True)\n",
    "# Remove axis on middle and right plot\n",
    "axes[1].get_yaxis().set_visible(False)\n",
    "axes[2].get_yaxis().set_visible(False)\n",
    "\n",
    "# add colour legend\n",
    "patches_list=[Patch(facecolor=colour) for colour in colours.values()]\n",
    "axes[1].legend(patches_list, list(dict_map.values()),\n",
    "               loc='upper center', ncol =4, bbox_to_anchor=(0.5, -0.1))\n",
    "# Add plot titles\n",
    "axes[0].set_title('True Colour Geomedian')\n",
    "axes[1].set_title('Classified Image')\n",
    "axes[2].set_title('Prediction Probability')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d606614-62f7-458b-a06e-62d722c67421",
   "metadata": {},
   "source": [
    "## Write the results to GeoTIFFs\n",
    "\n",
    "We can export our predictions at the testing locations and the entire district to sandbox disk as Cloud-Optimised GeoTIFFs. Here we also export the true colour images for the testing locations which we can load and compare the post-processed results with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f322b240-c5fb-457f-9461-65c6b4ea85f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "outname_prediction='Results/Land_cover_prediction_entire_district.tif'\n",
    "outname_probability='Results/Land_cover_probability_entire_district.tif'\n",
    "write_cog(predicted.Predictions, outname_prediction, overwrite=True)\n",
    "write_cog(predicted.Probabilities, outname_probability, overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0cc29d3-078e-48ae-b5c5-b74a0af1dfe4",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(0, len(predictions)):\n",
    "    outname_prediction='Results/Land_cover_prediction_location_'+str(i)+'.tif'\n",
    "    outname_rgb='Results/S2_RGB_image_location_'+str(i)+'.tif'\n",
    "    write_cog(predictions[i].Predictions, outname_prediction, overwrite=True)\n",
    "    write_cog(predictions[i][['red_3','green_3','blue_3']].to_array(), outname_rgb, overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "912db701-fe8e-47b5-ad0b-71090891b22e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  },
  "vscode": {
   "interpreter": {
    "hash": "01a9736267bf300689e8e016092cd01f0c67384dd94651ae6e139a291bc8cc97"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}