{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Tree crown delineation using detectreeRGB\n", "\n", "{opticon}`tag`\n", "{bdg-primary}`Forest`\n", "{bdg-secondary}`Modelling`\n", "{bdg-warning}`Standard`\n", "{bdg-info}`Python`\n", "\n", "

\n", " \n", " \"license\"\n", " \n", " \n", " \"render\"\n", " \n", " \n", " \"review\"\n", " \n", "
\n", "

\n", "\n", "

\n", " \n", " \"binder\"\n", " \n", " \n", " \"binder\"\n", " \n", "
\n", "

\n", "\n", "

\n", " \n", " \"RoHub\"\n", " \n", " \n", " \"doi\"\n", " \n", "

\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Context\n", "### Purpose\n", "Accurately delineating trees using `detectron2`, a library that provides state-of-the-art deep learning detection and segmentation algorithms.\n", "\n", "### Modelling approach\n", "An established deep learning model, **Mask R-CNN** was deployed from `detectron2` library to delineate tree crowns accurately. A pre-trained model, named `detectreeRGB`, is provided to predict the location and extent of tree crowns from a top-down RGB image, captured by drone, aircraft or satellite. `detectreeRGB` was implemented in `python` 3.8 using `pytorch` v1.7.1 and `detectron2` v0.5. Further details can be found in the [repository documentation](https://github.com/shmh40/detectreeRGB).\n", "\n", "### Highlights\n", "* [detectreeRGB](https://github.com/shmh40/detectreeRGB) advances the state-of-the-art in tree identification from RGB images by delineating exactly the extent of the tree crown.\n", "* We demonstrate how to apply the pretrained model to a sample image fetched from a Zenodo repository.\n", "* Our pre-trained model was developed using aircraft images of tropical forests in Malaysia.\n", "* The model can be further trained using the user's own images.\n", "\n", "### Contributions\n", "\n", "#### Notebook\n", "* Sebastian H. M. Hickman (author), University of Cambridge, [@shmh40](https://github.com/shmh40)\n", "* Alejandro Coca-Castro (reviewer), The Alan Turing Institute, [@acocac](https://github.com/acocac)\n", "\n", "#### Modelling codebase\n", "* Sebastian H. M. Hickman (author), University of Cambridge [@shmh40](https://github.com/shmh40)\n", "* James G. C. Ball (contributor), University of Cambridge [@PatBall1](https://github.com/PatBall1)\n", "* David A. Coomes (contributor), University of Cambridge\n", "* Toby Jackson (contributor), University of Cambridge\n", "\n", "#### Modelling funding\n", "The project was supported by the UKRI Centre for Doctoral Training in Application of Artificial Intelligence to the study of Environmental Risks ([AI4ER](https://ai4er-cdt.esc.cam.ac.uk/)) (EP/S022961/1).\n", "\n", ":::{note}\n", "The authors acknowledge the authors of the [Detectron2](https://github.com/facebookresearch/detectron2) package which provides the Mask R-CNN architecture.\n", ":::" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Load libraries" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "import cv2\n", "from PIL import Image\n", "import os\n", "import numpy as np\n", "import urllib.request\n", "import glob\n", "\n", "# intake library and plugin\n", "import intake\n", "from intake_zenodo_fetcher import download_zenodo_files_for_entry\n", "\n", "# geospatial libraries\n", "import geopandas as gpd\n", "\n", "from rasterio.transform import from_origin\n", "import rasterio.features\n", "\n", "import fiona\n", "\n", "from shapely.geometry import shape, mapping, box\n", "from shapely.geometry.multipolygon import MultiPolygon\n", "\n", "# machine learning libraries\n", "from detectron2 import model_zoo\n", "from detectron2.engine import DefaultPredictor\n", "from detectron2.utils.visualizer import Visualizer, ColorMode\n", "from detectron2.config import get_cfg\n", "from detectron2.engine import DefaultTrainer\n", "\n", "# visualisation\n", "import holoviews as hv\n", "import geoviews.tile_sources as gts\n", "import matplotlib.pyplot as plt\n", "\n", "import hvplot.pandas\n", "import hvplot.xarray\n", "\n", "import warnings\n", "warnings.filterwarnings(action='ignore')\n", "\n", "hv.extension('bokeh', width=100)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Set folder structure" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Define the project main folder\n", "notebook_folder = './notebook'\n", "\n", "# Set the folder structure\n", "config = {\n", " 'in_geotiff': os.path.join(notebook_folder, 'input','tiff'),\n", " 'in_png': os.path.join(notebook_folder, 'input','png'),\n", " 'model': os.path.join(notebook_folder, 'model'),\n", " 'out_geotiff': os.path.join(notebook_folder, 'output','raster'),\n", " 'out_shapefile': os.path.join(notebook_folder, 'output','vector'),\n", "}\n", "\n", "# List comprehension for the folder structure code\n", "[os.makedirs(val) for key, val in config.items() if not os.path.exists(val)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Load and prepare input image" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fetch a GeoTIFF file of aerial forest imagery using `intake`\n", "\n", "Let's fetch a sample aerial image from a [Zenodo repository](https://zenodo.org/record/5494629#.YWQQetnMKjA)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# write a catalog YAML file\n", "catalog_file = os.path.join(notebook_folder, 'catalog.yaml')\n", "\n", "with open(catalog_file, 'w') as f:\n", " f.write('''\n", "sources:\n", " sepilok_rgb:\n", " driver: rasterio\n", " description: 'NERC RGB images of Sepilok, Sabah, Malaysia (collection)'\n", " metadata:\n", " zenodo_doi: \"10.5281/zenodo.5494629\"\n", " args:\n", " urlpath: \"{{ CATALOG_DIR }}/input/tiff/Sep_2014_RGB_602500_646600.tif\"\n", " ''')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cat_tc = intake.open_catalog(catalog_file)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for catalog_entry in list(cat_tc):\n", " download_zenodo_files_for_entry(\n", " cat_tc[catalog_entry],\n", " force_download=False\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tc_rgb = cat_tc[\"sepilok_rgb\"].to_dask()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Inspect the aerial image\n", "\n", "Let's investigate the `data-array`, what is the shape? Bounds? Bands? CRS?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "jupyter": { "outputs_hidden": true }, "tags": [] }, "outputs": [], "source": [ "print('shape =', tc_rgb.shape,',', 'and number of bands =', tc_rgb.count, ', crs =', tc_rgb.crs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Save the RGB bands of the GeoTIFF file as a PNG\n", "\n", "**Mask R-CNN** requires images in `png` format. Let's export the RGB bands to a `png` file." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "minx = 602500\n", "miny = 646600\n", "\n", "R = tc_rgb[0]\n", "G = tc_rgb[1]\n", "B = tc_rgb[2]\n", " \n", "# stack up the bands in an order appropriate for saving with cv2, then rescale to the correct 0-255 range for cv2\n", "\n", "# you will have to change the rescaling depending on the values of your tiff!\n", "rgb = np.dstack((R,G,B)) # BGR for cv2\n", "rgb_rescaled = 255*rgb/65535 # scale to image\n", " \n", "# save this as png, named with the origin of the specific tile - change the filepath!\n", "filepath = config['in_png'] + '/' + 'tile_' + str(minx) + '_' + str(miny) + '.png'\n", "cv2.imwrite(filepath, rgb_rescaled)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Read in and display the PNG file" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "im = cv2.imread(filepath)\n", "plot_input = plt.figure(figsize=(15,15))\n", "plt.imshow(Image.fromarray(im))\n", "plt.title('Input image',fontsize='xx-large')\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Modelling\n", "\n", "### Download the pretrained model" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# define the URL to retrieve the model\n", "fn = 'model_final.pth'\n", "url = f'https://zenodo.org/record/5515408/files/{fn}?download=1'\n", "\n", "urllib.request.urlretrieve(url, config['model'] + '/' + fn)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Settings of `detectron2` config\n", "\n", "The following lines allow configuring the main settings for predictions and load them into a `DefaultPredictor` object." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cfg = get_cfg()\n", "\n", "# if you want to make predictions using a CPU, run the following line. If using GPU, hash it out.\n", "cfg.MODEL.DEVICE='cpu'\n", "\n", "# model and hyperparameter selection\n", "cfg.merge_from_file(model_zoo.get_config_file(\"COCO-InstanceSegmentation/mask_rcnn_R_101_FPN_3x.yaml\"))\n", "cfg.DATALOADER.NUM_WORKERS = 2\n", "cfg.SOLVER.IMS_PER_BATCH = 2\n", "cfg.MODEL.ROI_HEADS.NUM_CLASSES = 1\n", "\n", "### path to the saved pre-trained model weights\n", "cfg.MODEL.WEIGHTS = config['model'] + '/model_final.pth'\n", "\n", "# set confidence threshold at which we predict\n", "cfg.MODEL.ROI_HEADS.SCORE_THRESH_TEST = 0.15\n", "\n", "#### Settings for predictions using detectron config\n", "\n", "predictor = DefaultPredictor(cfg)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Outputs\n", "\n", "### Showing the predictions from `detectreeRGB`" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "outputs = predictor(im)\n", "v = Visualizer(im[:, :, ::-1], scale=1.5, instance_mode=ColorMode.IMAGE_BW) # remove the colors of unsegmented pixels\n", "v = v.draw_instance_predictions(outputs[\"instances\"].to(\"cpu\"))\n", "image = cv2.cvtColor(v.get_image()[:, :, :], cv2.COLOR_BGR2RGB)\n", "\n", "plot_predictions = plt.figure(figsize=(15,15))\n", "plt.imshow(Image.fromarray(image))\n", "plt.title('Predictions',fontsize='xx-large')\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Convert predictions into geospatial files\n", "\n", "#### To GeoTIFF" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mask_array = outputs['instances'].pred_masks.cpu().numpy()\n", "\n", "# get confidence scores too \n", "mask_array_scores = outputs['instances'].scores.cpu().numpy()\n", "\n", "num_instances = mask_array.shape[0]\n", "mask_array_instance = []\n", "output = np.zeros_like(mask_array) \n", "\n", "mask_array_instance.append(mask_array)\n", "output = np.where(mask_array_instance[0] == True, 255, output)\n", "fresh_output = output.astype(float)\n", "x_scaling = 140/fresh_output.shape[1]\n", "y_scaling = 140/fresh_output.shape[2]\n", "# this is an affine transform. This needs to be altered significantly.\n", "transform = from_origin(int(filepath[-17:-11])-20, int(filepath[-10:-4])+120, y_scaling, x_scaling)\n", "\n", "output_raster = config['out_geotiff'] + '/' + 'predicted_rasters_' + filepath[-17:-4]+ '.tif'\n", "\n", "new_dataset = rasterio.open(output_raster, 'w', driver='GTiff',\n", " height = fresh_output.shape[1], width = fresh_output.shape[2], count = fresh_output.shape[0],\n", " dtype=str(fresh_output.dtype),\n", " crs='+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs', \n", " transform=transform)\n", "\n", "new_dataset.write(fresh_output)\n", "new_dataset.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### To shapefile" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Read input band with Rasterio\n", " \n", "with rasterio.open(output_raster) as src:\n", " shp_schema = {'geometry': 'MultiPolygon','properties': {'pixelvalue': 'int', 'score': 'float'}} \n", "\n", " crs = src.crs\n", " for i in range(src.count):\n", " src_band = src.read(i+1)\n", " src_band = np.float32(src_band)\n", " conf = mask_array_scores[i]\n", " # Keep track of unique pixel values in the input band\n", " unique_values = np.unique(src_band)\n", " # Polygonize with Rasterio. `shapes()` returns an iterable\n", " # of (geom, value) as tuples\n", " shapes = list(rasterio.features.shapes(src_band, transform=src.transform))\n", "\n", " if i == 0:\n", " with fiona.open(config['out_shapefile'] + '/predicted_polygons_' + filepath[-17:-4] + '_' + str(0) + '.shp', 'w', 'ESRI Shapefile',\n", " shp_schema) as shp:\n", " polygons = [shape(geom) for geom, value in shapes if value == 255.0] \n", " multipolygon = MultiPolygon(polygons)\n", " # simplify not needed here\n", " #multipolygon = multipolygon_a.simplify(0.1, preserve_topology=False) \n", " shp.write({\n", " 'geometry': mapping(multipolygon),\n", " 'properties': {'pixelvalue': int(unique_values[1]), 'score': float(conf)} \n", " })\n", " else:\n", " with fiona.open(config['out_shapefile'] + '/predicted_polygons_' + filepath[-17:-4] + '_' + str(0)+'.shp', 'a', 'ESRI Shapefile',\n", " shp_schema) as shp:\n", " polygons = [shape(geom) for geom, value in shapes if value == 255.0] \n", " multipolygon = MultiPolygon(polygons)\n", " # simplify not needed here\n", " #multipolygon = multipolygon_a.simplify(0.1, preserve_topology=False) \n", " shp.write({\n", " 'geometry': mapping(multipolygon),\n", " 'properties': {'pixelvalue': int(unique_values[1]), 'score': float(conf)} \n", " })" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Interactive map to visualise the exported files\n", "\n", "### Plot tree crown delineation predictions and scores" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# load and plot polygons\n", "in_shp = glob.glob(config['out_shapefile'] + '/*.shp')\n", "\n", "poly_df = gpd.read_file(in_shp[0])\n", "\n", "plot_vector = poly_df.hvplot(hover_cols=['score'], legend=False).opts(fill_color=None,line_color=None,alpha=0.5, width=800, height=600, xaxis=None, yaxis=None)\n", "plot_vector" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot the exported GeoTIFF file" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# load and plot RGB image\n", "r = tc_rgb.sel(band=[1,2,3])\n", "\n", "normalized = r/(r.quantile(.99,skipna=True)/255)\n", "\n", "mask = normalized.where(normalized < 255)\n", "\n", "int_arr = mask.astype(int)\n", "\n", "plot_rgb = int_arr.astype('uint8').hvplot.rgb(\n", " x='x', y='y', bands='band', data_aspect=0.8, hover=False, legend=False, rasterize=True, xaxis=None, yaxis=None, title='Tree crown delineation by detectreeRGB'\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Overlay the prediction labels and image\n", "\n", "Note we have some artifacts in the RGB image due to the transformations using the normalization procedure." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "plot_predictions_interactive = plot_rgb * plot_vector\n", "plot_predictions_interactive" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Save plot" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hvplot.save(plot_predictions_interactive, notebook_folder + '/interactive_predictions.html')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary\n", "\n", "We have read in a raster, chosen a tile and made predictions on it. These predictions can then be transformed to shapefiles and examined in GIS software!\n", "\n", "* We made the predictions on the `png` using a pretrained **Mask R-CNN** model, `detectreeRGB`.\n", "* The outputs showed the model capability to detect and delineate tree crowns from aerial imagery.\n", "* We then extracted our predictions, added the geospatial location back in, and exported them as shapefiles including the confidence score assigned to each prediction by the model.\n", "* Visualised the predictions on an interactive map!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Citing this Notebook" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Please see [CITATION.cff](https://github.com/eds-book-gallery/94486a7f-e046-461f-bbb9-334ec7b57040/blob/main/CITATION.cff) for the full citation information. The citation file can be exported to APA or BibTex formats (learn more [here](https://docs.github.com/en/repositories/managing-your-repositorys-settings-and-features/customizing-your-repository/about-citation-files))." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Additional information\n", "**Codebase**: version 1.0.0 with commit [16a5a1c](https://github.com/shmh40/detectreeRGB/compare/main...dev)\n", "\n", "**License**: The code in this notebook is licensed under the MIT License. The Environmental Data Science book is licensed under the Creative Commons by Attribution 4.0 license. See further details [here](https://github.com/alan-turing-institute/environmental-ds-book/blob/master/LICENSE.md).\n", "\n", "**Contact**: If you have any suggestion or report an issue with this notebook, feel free to [create an issue](https://github.com/alan-turing-institute/environmental-ds-book/issues/new/choose) or send a direct message to [environmental.ds.book@gmail.com](mailto:environmental.ds.book@gmail.com)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "remove-input" ] }, "outputs": [], "source": [ "from datetime import date\n", "\n", "print('Notebook repository version: v1.0.2')\n", "print(f'Last tested: {date.today()}')" ] }, { "cell_type": "markdown", "metadata": { "tags": [ "remove-cell" ] }, "source": [ "## Outputs registration\n", "The cell below is dedicated to save the notebook outputs for registering them into a Zenodo repository curated by the Environmental DS book." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "outputs = {\n", " 'static_figures': {\n", " 'filenames': ['static_input','static_predictions'],\n", " 'data':[plot_input, plot_predictions]},\n", " 'interactive_figures': {\n", " 'filenames': ['interactive_vector','interactive_predictions'],\n", " 'data':[plot_vector, plot_predictions_interactive]},\n", "}\n", "\n", "#save static figures\n", "if len(outputs['static_figures']['filenames']) > 0:\n", " [data.savefig(os.path.join(notebook_folder,outputs['static_figures']['filenames'][x] + '.png')) for x, data in enumerate(outputs['static_figures']['data'])]\n", "\n", "#save interactive figures\n", "if len(outputs['interactive_figures']['filenames']) > 0:\n", " [hvplot.save(data, os.path.join(notebook_folder,outputs['interactive_figures']['filenames'][x] + '.html')) for x, data in enumerate(outputs['interactive_figures']['data'])]" ] } ], "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" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": {}, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 4 }