notebooks/A01-fragstats-metrics-comparison.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bbd00868",
   "metadata": {},
   "source": [
    "# FRAGSTATS metrics comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60c74eb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import pylandstats as pls"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc815626",
   "metadata": {},
   "source": [
    "The aim of this notebook is to compare the computed values for the landscape metrics of PyLandStats with those of FRAGSTATS (v4.2). Let us first set the relative tolerance to 0.001, which means that we accept a relative difference of 0.1% between the values computed with PyLandStats and FRAGSTATS. We will also define some utilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2f9ddee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# definitions\n",
    "\n",
    "fragstats_abbrev_dict = pls.settings.fragstats_abbrev_dict\n",
    "tol = 1e-3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f3ee825",
   "metadata": {},
   "outputs": [],
   "source": [
    "# utils\n",
    "\n",
    "\n",
    "def read_fragstats_csv(csv_filepath):\n",
    "    # `na_values` is required because of the leading whitespaces that FRAGSTATS\n",
    "    # leaves when saving CSV files\n",
    "    fragstats_df = pd.read_csv(csv_filepath, na_values=[\" N/A\"])\n",
    "    fragstats_df.columns = fragstats_df.columns.str.strip()\n",
    "    try:\n",
    "        fragstats_df[\"TYPE\"] = (\n",
    "            fragstats_df[\"TYPE\"].str.strip().str.replace(\"cls_\", \"\").astype(int)\n",
    "        )\n",
    "    except KeyError:\n",
    "        pass\n",
    "\n",
    "    return fragstats_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12a6f7e6",
   "metadata": {},
   "source": [
    "Below is the reference basename for the extract of the Canton of Vaud (Switerland) derived from the [CORINE Land Cover dataset](https://land.copernicus.eu/pan-european/corine-land-cover) of the year 2000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b033428",
   "metadata": {},
   "outputs": [],
   "source": [
    "basename = \"vaud_g100_clc00_V18_5\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0277544f",
   "metadata": {},
   "source": [
    "We will instantiate a PyLandStats landscape with the respective raster file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3abf1081",
   "metadata": {},
   "outputs": [],
   "source": [
    "ls = pls.Landscape(\"../data/raw/clc/{}.tif\".format(basename), res=(100, 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d82c0ca4",
   "metadata": {},
   "source": [
    "Now we will use pandas to load the `.patch`, `.class` and `.landscape` files (which are dumps of the metrics computed in FRAGSTATS), and then compare the computed values with those of PyLandStats. If the value for any metric differs more than the relative tolerance defined above, a `RuntimeError` will be raised."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "627f222b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "area: OK\n",
      "perimeter: OK\n",
      "perimeter_area_ratio: OK\n",
      "shape_index: OK\n",
      "fractal_dimension: OK\n",
      "euclidean_nearest_neighbor: OK\n"
     ]
    }
   ],
   "source": [
    "patch_df = read_fragstats_csv(\"../data/raw/fragstats/{}.patch\".format(basename))\n",
    "\n",
    "for patch_metric in pls.Landscape.PATCH_METRICS:\n",
    "    fragstats_abbrev = fragstats_abbrev_dict[patch_metric]\n",
    "    for class_val in ls.classes:\n",
    "        fragstats_ser = patch_df[fragstats_abbrev][patch_df[\"TYPE\"] == class_val]\n",
    "        pls_ser = getattr(ls, patch_metric)(class_val=class_val)\n",
    "        if not np.allclose(fragstats_ser, pls_ser, tol, equal_nan=True):\n",
    "            raise RuntimeError(patch_metric, class_val, fragstats_ser, pls_ser)\n",
    "    print(\"{}: OK\".format(patch_metric))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76ef9ac8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total_edge (class 1): OK\n",
      "total_edge (class 2): OK\n",
      "patch_density (class 1): OK\n",
      "patch_density (class 2): OK\n",
      "proportion_of_landscape (class 1): OK\n",
      "proportion_of_landscape (class 2): OK\n",
      "effective_mesh_size (class 1): OK\n",
      "effective_mesh_size (class 2): OK\n",
      "edge_density (class 1): OK\n",
      "edge_density (class 2): OK\n",
      "landscape_shape_index (class 1): OK\n",
      "landscape_shape_index (class 2): OK\n",
      "number_of_patches (class 1): OK\n",
      "number_of_patches (class 2): OK\n",
      "total_area (class 1): OK\n",
      "total_area (class 2): OK\n",
      "largest_patch_index (class 1): OK\n",
      "largest_patch_index (class 2): OK\n"
     ]
    }
   ],
   "source": [
    "class_df = read_fragstats_csv(\"../data/raw/fragstats/{}.class\".format(basename))\n",
    "\n",
    "for class_metric in set(pls.Landscape.CLASS_METRICS).difference(\n",
    "    set(pls.Landscape.DISTR_METRICS)\n",
    "):\n",
    "    if class_metric == \"total_area\":\n",
    "        fragstats_abbrev = \"CA\"\n",
    "    else:\n",
    "        fragstats_abbrev = fragstats_abbrev_dict[class_metric]\n",
    "    for class_val in ls.classes:\n",
    "        fragstats_val = class_df[fragstats_abbrev][class_df[\"TYPE\"] == class_val].iloc[\n",
    "            0\n",
    "        ]\n",
    "        pls_val = getattr(ls, class_metric)(class_val=class_val)\n",
    "        if not (\n",
    "            np.isclose(fragstats_val, pls_val, tol)\n",
    "            or np.isclose(fragstats_val, pls_val, atol=tol)\n",
    "        ):\n",
    "            raise RuntimeError(\n",
    "                \"{} (class {}): fragstats {}, pylandstats {}\".format(\n",
    "                    class_metric, class_val, fragstats_val, pls_val\n",
    "                )\n",
    "            )\n",
    "\n",
    "        print(\"{} (class {}): OK\".format(class_metric, class_val))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef6b83e1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shannon_diversity_index: OK\n",
      "patch_density: OK\n",
      "number_of_patches: OK\n",
      "effective_mesh_size: OK\n",
      "edge_density: OK\n",
      "landscape_shape_index: OK\n",
      "contagion: OK\n",
      "total_edge: OK\n",
      "total_area: OK\n",
      "largest_patch_index: OK\n"
     ]
    }
   ],
   "source": [
    "landscape_df = read_fragstats_csv(\"../data/raw/fragstats/{}.land\".format(basename))\n",
    "\n",
    "# we exclude the distribution metrics and the entropy metrics, since except SHDI and\n",
    "# CONTAG, they are not implemented in FRAGSTATS\n",
    "for landscape_metric in (\n",
    "    set(pls.Landscape.LANDSCAPE_METRICS)\n",
    "    .difference(set(pls.Landscape.ENTROPY_METRICS).union(pls.Landscape.DISTR_METRICS))\n",
    "    .union([\"shannon_diversity_index\", \"contagion\"])\n",
    "):\n",
    "    fragstats_abbrev = fragstats_abbrev_dict[landscape_metric]\n",
    "    fragstats_val = landscape_df[fragstats_abbrev].iloc[0]\n",
    "    pls_val = getattr(ls, landscape_metric)()\n",
    "    if not np.isclose(fragstats_val, pls_val, rtol=0.01):\n",
    "        raise RuntimeError(\n",
    "            \"{}: fragstats {}, pylandstats {}\".format(\n",
    "                landscape_metric, fragstats_val, pls_val\n",
    "            )\n",
    "        )\n",
    "\n",
    "    print(\"{}: OK\".format(landscape_metric))\n",
    "\n",
    "# # Treat contagion differently: here we use a relative tolerance of .01 (1% relative\n",
    "# # error) because the computed contagion might have greater divergence with FRAGSTATS\n",
    "\n",
    "# landscape_metric = \"contagion\"\n",
    "# fragstats_abbrev = fragstats_abbrev_dict[landscape_metric]\n",
    "# fragstats_val = landscape_df[fragstats_abbrev].iloc[0]\n",
    "# pls_val = getattr(ls, landscape_metric)()\n",
    "\n",
    "# if not np.isclose(fragstats_val, pls_val, rtol=0.01):\n",
    "#     raise RuntimeError(\n",
    "#         \"{}: fragstats {}, pylandstats {}\".format(\n",
    "#             landscape_metric, fragstats_val, pls_val\n",
    "#         )\n",
    "#     )\n",
    "# print(\n",
    "#     \"{}: OK (fragstats: {}, pylandstats: {})\".format(\n",
    "#         landscape_metric, fragstats_val, pls_val\n",
    "#     )\n",
    "# )"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (pylandstats)",
   "language": "python",
   "name": "pylandstats"
  },
  "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}