Species_Distribution_Modeling_Maxent.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "8kdsGkYJXXKc"
   },
   "outputs": [],
   "source": [
    "#@title Copyright 2023 The Earth Engine Community Authors { display-mode: \"form\" }\n",
    "#\n",
    "# Licensed under the Apache License, Version 2.0 (the \"License\");\n",
    "# you may not use this file except in compliance with the License.\n",
    "# You may obtain a copy of the License at\n",
    "#\n",
    "# https://www.apache.org/licenses/LICENSE-2.0\n",
    "#\n",
    "# Unless required by applicable law or agreed to in writing, software\n",
    "# distributed under the License is distributed on an \"AS IS\" BASIS,\n",
    "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
    "# See the License for the specific language governing permissions and\n",
    "# limitations under the License."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "l18M9_r5XmAQ"
   },
   "source": [
    "# Species Distribution Modeling\n",
    "Author: osgeokr\n",
    "\n",
    "In this tutorial, the methodology of Species Distribution Modeling using Google Earth Engine will be introduced. A brief overview of Species Distribution Modeling  will be provided, followed by the process of predicting and analyzing the habitat of an endangered bird species known as the Fairy pitta (scientific name: *Pitta nympha*)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "U7i55vr_aKCB"
   },
   "source": [
    "### Run me first\n",
    "\n",
    "Run the following cell to initialize the API. The output will contain instructions on how to grant this notebook access to Earth Engine using your account."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "XeFsiSp2aDL6"
   },
   "outputs": [],
   "source": [
    "import ee\n",
    "\n",
    "# Trigger the authentication flow.\n",
    "ee.Authenticate()\n",
    "\n",
    "# Initialize the library.\n",
    "ee.Initialize(project='my-project')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VOf_UnIcZKBJ"
   },
   "source": [
    "## A brief overview of Species Distribution Modeling\n",
    "\n",
    "Let's explore what species distribution models are, the advantages of using Google Earth Engine for their processing, the required data for the models, and how the workflow is structured.\n",
    "\n",
    "### What is Species Distribution Modeling?\n",
    "\n",
    "Species Distribution Modeling (SDM below) is the most common methodology used to estimate the actual or potential geographic distribution of a species. It involves characterizing the environmental conditions suitable for a particular species and then identifying where these suitable conditions are distributed geographically.\n",
    "\n",
    "SDM has emerged as a crucial component of conservation planning in recent years, and various modeling techniques have been developed for this purpose. Implementing SDM in Google Earth Engine (GEE below) provides easy access to large-scale environmental data, along with powerful computing capabilities and support for machine learning algorithms, allowing for rapid modeling.\n",
    "\n",
    "  > Note: Conservation biologist Dr. Ramiro D. Crego implemented SDM using the GEE JavaScript Code Editor and published his research findings [(Crego et al, 2022)](https://onlinelibrary.wiley.com/doi/10.1111/ddi.13491). The methodology of SDM introduced here has been translated and modified from the [JavaScript source code](https://smithsonian.github.io/SDMinGEE/) he shared into the Python language.\n",
    "\n",
    "### Data Required for SDM\n",
    "\n",
    "SDM typically utilizes the relationship between known species occurrence records and environmental variables to identify the conditions under which a population can sustain. In other words, two types of model input data are required:\n",
    "\n",
    "1. Occurrence records of known species\n",
    "1. Various environmental variables\n",
    "\n",
    "These data are input into algorithms to identify environmental conditions associated with the presence of species.\n",
    "\n",
    "### Workflow of SDM using GEE\n",
    "\n",
    "The workflow for SDM using GEE is as follows:\n",
    "\n",
    "1. Collection and preprocessing of species occurrence data\n",
    "1. Definition of the Area of Interest\n",
    "1. Addition of GEE environmental variables\n",
    "1. Generation of pseudo-absence data\n",
    "1. Model fitting and prediction\n",
    "1. Variable importance and accuracy assessment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "tjomxWfVcTmN"
   },
   "source": [
    "## Habitat Prediction and Analysis Using GEE\n",
    "\n",
    "The [Fairy pitta (*Pitta nympha*)](https://datazone.birdlife.org/species/factsheet/22698684) will be used as a case study to demonstrate the application of GEE-based SDM. While this specific species has been selected for one example, researchers can apply the methodology to any target species of interest with slight modifications to the provided source code.\n",
    "\n",
    "The Fairy pitta is a rare summer migrant and passage migrant in South Korea, whose distribution area is expanding due to recent climate warming on the Korean Peninsula. It is classified as a rare species, endangered wildlife of class II, Natural Monument No. 204, evaluated as Regionally Extinct (RE) in the National Red List, and Vulnerable (VU) according to the IUCN categories.\n",
    "\n",
    "Conducting SDM for the conservation planning of the Fairy pitta appears to be quite valuable. Now, let's proceed with habitat prediction and analysis through GEE."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "pViK9PM-gLjh"
   },
   "source": [
    "First, the Python libraries are imported.The `import` statement brings in the entire contents of a module, while the `from import` statement allows for the importation of specific objects from a module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "4jbM03uIrjST"
   },
   "outputs": [],
   "source": [
    "# Import libraries\n",
    "import geemap\n",
    "\n",
    "import geemap.colormaps as cm\n",
    "import pandas as pd, geopandas as gpd\n",
    "import numpy as np, matplotlib.pyplot as plt\n",
    "import os, requests, math, random\n",
    "\n",
    "from ipyleaflet import TileLayer\n",
    "from statsmodels.stats.outliers_influence import variance_inflation_factor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "SRrwa4ROghr9"
   },
   "source": [
    "### Collection and Preprocessing of Species Occurrence Data\n",
    "\n",
    "Now, let's collect occurrence data for the Fairy pitta. Even if you don't currently have access to occurrence data for the species of interest, you can obtain observational data about specific species through the GBIF API. The [GBIF API](https://techdocs.gbif.org/en/openapi/) is an interface that allows access to the species distribution data provided by GBIF, enabling users to search, filter, and download data, as well as acquire various information related to species.\n",
    "\n",
    "In the code below, the `species_name` variable is assigned the scientific name of the species (e.g., *Pitta nympha* for Fairy pitta), and the `country_code` variable is assigned the country code (e.g., KR for South Korea). The `base_url` variable stores the address of the GBIF API. `params` is a dictionary containing parameters to be used in the API request:\n",
    "\n",
    "* `scientificName`: Sets the scientific name of the species to be searched.\n",
    "* `country`: Limits the search to a specific country.\n",
    "* `hasCoordinate`: Ensures only data with coordinates (true) are searched.\n",
    "* `basisOfRecord`: Chooses only records of human observation (`HUMAN_OBSERVATION`).\n",
    "* `limit`: Sets the maximum number of results returned to 10000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "oHtKaH0FgXTz"
   },
   "outputs": [],
   "source": [
    "def get_gbif_species_data(species_name, country_code):\n",
    "    \"\"\"\n",
    "    Retrieves observational data for a specific species using the GBIF API and returns it as a pandas DataFrame.\n",
    "\n",
    "    Parameters:\n",
    "    species_name (str): The scientific name of the species to query.\n",
    "    country_code (str): The country code of the where the observation data will be queried.\n",
    "\n",
    "    Returns:\n",
    "    pd.DataFrame: A pandas DataFrame containing the observational data.\n",
    "    \"\"\"\n",
    "    base_url = \"https://api.gbif.org/v1/occurrence/search\"\n",
    "    params = {\n",
    "        \"scientificName\": species_name,\n",
    "        \"country\": country_code,\n",
    "        \"hasCoordinate\": \"true\",\n",
    "        \"basisOfRecord\": \"HUMAN_OBSERVATION\",\n",
    "        \"limit\": 10000,\n",
    "    }\n",
    "\n",
    "    try:\n",
    "        response = requests.get(base_url, params=params)\n",
    "        response.raise_for_status()  # Raises an exception for a response error.\n",
    "        data = response.json()\n",
    "        occurrences = data.get(\"results\", [])\n",
    "\n",
    "        if occurrences:  # If data is present\n",
    "            df = pd.json_normalize(occurrences)\n",
    "            return df\n",
    "        else:\n",
    "            print(\"No data found for the given species and country code.\")\n",
    "            return pd.DataFrame()  # Returns an empty DataFrame\n",
    "    except requests.RequestException as e:\n",
    "        print(f\"Request failed: {e}\")\n",
    "        return pd.DataFrame()  # Returns an empty DataFrame in case of an exception"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Zs5ZUfZUnjZ2"
   },
   "source": [
    "Using the parameters set previously, we query the GBIF API for observational records of the Fairy pitta (*Pitta nympha*), and load the results into a DataFrame to check the first row. A DataFrame is a data structure for handling table-formatted data, consisting of rows and columns. If necessary, the DataFrame can be saved as a CSV file and read back in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Mx-DjtGNnUXk"
   },
   "outputs": [],
   "source": [
    "# Retrieve Fairy Pitta data\n",
    "df = get_gbif_species_data(\"Pitta nympha\", \"KR\")\n",
    "\"\"\"\n",
    "# Save DataFrame to CSV and read back in.\n",
    "df.to_csv(\"pitta_nympha_data.csv\", index=False)\n",
    "df = pd.read_csv(\"pitta_nympha_data.csv\")\n",
    "\"\"\"\n",
    "df.head(1)  # Display the first row of the DataFrame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "TEjSEmK3pfe0"
   },
   "source": [
    "Next, we convert the DataFrame into a GeoDataFrame that includes a column for geographic information (`geometry`) and check the first row. A GeoDataFrame can be saved as a GeoPackage file (*.gpkg) and read back in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "qjt0jgJCpALg"
   },
   "outputs": [],
   "source": [
    "# Convert DataFrame to GeoDataFrame\n",
    "gdf = gpd.GeoDataFrame(\n",
    "    df,\n",
    "    geometry=gpd.points_from_xy(df.decimalLongitude,\n",
    "                                df.decimalLatitude),\n",
    "    crs=\"EPSG:4326\"\n",
    ")[[\"species\", \"year\", \"month\", \"geometry\"]]\n",
    "\"\"\"\n",
    "# Convert GeoDataFrame to GeoPackage (requires pycrs module)\n",
    "%pip install -U -q pycrs\n",
    "gdf.to_file(\"pitta_nympha_data.gpkg\", driver=\"GPKG\")\n",
    "gdf = gpd.read_file(\"pitta_nympha_data.gpkg\")\n",
    "\"\"\"\n",
    "gdf.head(1)  # Display the first row of the GeoDataFrame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5Lj919AaqmUq"
   },
   "source": [
    "This time, we have created a function to visualize the distribution of data by year and month from the GeoDataFrame and display it as a graph, which can then be saved as an image file. The use of a heatmap allows us to quickly grasp the frequency of species occurrence by year and month, providing an intuitive visualization of the temporal changes and patterns within the data. This allows for the identification of temporal patterns and seasonal variations in species occurrence data, as well as the rapid detection of outliers or quality issues within the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Iju-dNZErzkJ"
   },
   "outputs": [],
   "source": [
    "# Yearly and monthly data distribution heatmap\n",
    "def plot_heatmap(gdf, h_size=8):\n",
    "\n",
    "    statistics = gdf.groupby([\"month\", \"year\"]).size().unstack(fill_value=0)\n",
    "\n",
    "    # Heatmap\n",
    "    plt.figure(figsize=(h_size, h_size - 6))\n",
    "    heatmap = plt.imshow(\n",
    "        statistics.values, cmap=\"YlOrBr\", origin=\"upper\", aspect=\"auto\"\n",
    "    )\n",
    "\n",
    "    # Display values above each pixel\n",
    "    for i in range(len(statistics.index)):\n",
    "        for j in range(len(statistics.columns)):\n",
    "            plt.text(\n",
    "                j, i, statistics.values[i, j], ha=\"center\", va=\"center\", color=\"black\"\n",
    "            )\n",
    "\n",
    "    plt.colorbar(heatmap, label=\"Count\")\n",
    "    plt.title(\"Monthly Species Count by Year\")\n",
    "    plt.xlabel(\"Year\")\n",
    "    plt.ylabel(\"Month\")\n",
    "    plt.xticks(range(len(statistics.columns)), statistics.columns)\n",
    "    plt.yticks(range(len(statistics.index)), statistics.index)\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(\"heatmap_plot.png\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "dlW1nIQZrjx3"
   },
   "outputs": [],
   "source": [
    "plot_heatmap(gdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vIj9CO1RrCSs"
   },
   "source": [
    "The data from 1995 is very sparse, with significant gaps compared to other years, and the months of August and September also have limited samples and exhibit different seasonal characteristics compared to other periods. Excluding this data could contribute to improving the stability and predictive power of the model.\n",
    "\n",
    "However, it's important to note that excluding data may enhance the model's generalization ability, but it could also lead to the loss of valuable information relevant to the research objectives. Therefore, such decisions should be made with careful consideration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "k1PpbNsGq9qh"
   },
   "outputs": [],
   "source": [
    "# Filtering data by year and month\n",
    "filtered_gdf = gdf[\n",
    "    (~gdf['year'].eq(1995)) &\n",
    "    (~gdf['month'].between(8, 9))\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jy0BBd-6sdLz"
   },
   "source": [
    "Now, the filtered GeoDataFrame is converted into a Google Earth Engine object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Fxu0OsBksMKM"
   },
   "outputs": [],
   "source": [
    "# Convert GeoDataFrame to Earth Engine object\n",
    "data_raw = geemap.geopandas_to_ee(filtered_gdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "55p_GfB6sv3U"
   },
   "source": [
    "Next, we will define the raster pixel size of the SDM results as 1km resolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "FsTbNQ17s1l-"
   },
   "outputs": [],
   "source": [
    "# Spatial resolution setting (meters)\n",
    "grain_size = 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "A20gAGUZtN6S"
   },
   "source": [
    "When multiple occurrence points are present within the same 1km resolution raster pixel, there is a high likelihood that they share the same environmental conditions at the same geographic location. Using such data directly in the analysis can introduce bias into the results.\n",
    "\n",
    "In other words, we need to limit the potential impact of geographic sampling bias. To achieve this, we will retain only one location within each 1km pixel and remove all others, allowing the model to more objectively reflect the environmental conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "dHtyQyMQs82v"
   },
   "outputs": [],
   "source": [
    "def remove_duplicates(data, grain_size):\n",
    "    # Select one occurrence record per pixel at the chosen spatial resolution\n",
    "    random_raster = ee.Image.random().reproject(\"EPSG:4326\", None, grain_size)\n",
    "    rand_point_vals = random_raster.sampleRegions(\n",
    "        collection=ee.FeatureCollection(data), geometries=True\n",
    "    )\n",
    "    return rand_point_vals.distinct(\"random\")\n",
    "\n",
    "\n",
    "data = remove_duplicates(data_raw, grain_size)\n",
    "\n",
    "# Before selection and after selection\n",
    "print(\"Original data size:\", data_raw.size().getInfo())\n",
    "print(\"Final data size:\", data.size().getInfo())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Rhu4b4BxuMHE"
   },
   "source": [
    "The visualization comparing geographic sampling bias before preprocessing (in blue) and after preprocessing (in red) is shown below. To facilitate comparison, the map has been centered on the area with a high concentration of Fairy pitta occurrence coordinates in Hallasan National Park."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "d9pFgpUztsB-"
   },
   "outputs": [],
   "source": [
    "# Visualization of geographic sampling bias before (blue) and after (red) preprocessing\n",
    "Map = geemap.Map(layout={\"height\": \"400px\", \"width\": \"800px\"})\n",
    "\n",
    "# Add the random raster layer\n",
    "random_raster = ee.Image.random().reproject(\"EPSG:4326\", None, grain_size)\n",
    "Map.addLayer(\n",
    "    random_raster,\n",
    "    {\"min\": 0, \"max\": 1, \"palette\": [\"black\", \"white\"], \"opacity\": 0.5},\n",
    "    \"Random Raster\",\n",
    ")\n",
    "\n",
    "# Add the original data layer in blue\n",
    "Map.addLayer(data_raw, {\"color\": \"blue\"}, \"Original data\")\n",
    "\n",
    "# Add the final data layer in red\n",
    "Map.addLayer(data, {\"color\": \"red\"}, \"Final data\")\n",
    "\n",
    "# Set the center of the map to the coordinates\n",
    "Map.setCenter(126.712, 33.516, 14)\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3oPtTj_O6d3H"
   },
   "source": [
    "### Definition of the Area of Interest\n",
    "\n",
    "Defining the Area of Interest (AOI below) refers to the term used by researchers to denote the geographical area they want to analyze. It has a similar meaning to the term Study Area.\n",
    "\n",
    "In this context, we obtained the bounding box of the occurrence point layer geometry and created a 50-kilometer buffer around it (with a maximum tolerance of 1,000 meters) to define the AOI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "XIyhhdzUvyZw"
   },
   "outputs": [],
   "source": [
    "# Define the AOI\n",
    "aoi = data.geometry().bounds().buffer(distance=50000, maxError=1000)\n",
    "\n",
    "# Add the AOI to the map\n",
    "outline = ee.Image().byte().paint(\n",
    "    featureCollection=aoi, color=1, width=3)\n",
    "\n",
    "Map.remove_layer(\"Random Raster\")\n",
    "Map.addLayer(outline, {'palette': 'FF0000'}, \"AOI\")\n",
    "Map.centerObject(aoi, 6)\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fWsplwHD8ZFd"
   },
   "source": [
    "### Addition of GEE environmental variables\n",
    "\n",
    "Now, let's add environmental variables to the analysis. GEE provides a wide range of datasets for environmental variables such as temperature, precipitation, elevation, land cover, and terrain. These datasets enable us to comprehensively analyze various factors that may influence the habitat preferences of the Fairy pitta.\n",
    "\n",
    "The selection of GEE environmental variables in SDM should reflect the habitat preference characteristics of the species. To do this, prior research and literature review on the Fairy pitta's habitat preferences should be conducted. This tutorial primarily focuses on the workflow of SDM using GEE, so some in-depth details are omitted.\n",
    "\n",
    "[**WorldClim V1 Bioclim**](https://developers.google.com/earth-engine/datasets/catalog/WORLDCLIM_V1_BIO): This dataset provides 19 bioclimatic variables derived from monthly temperature and precipitation data. It covers the period from 1960 to 1991 and has a resolution of 927.67 meters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "TNW23qpX7shy"
   },
   "outputs": [],
   "source": [
    "# WorldClim V1 Bioclim\n",
    "bio = ee.Image(\"WORLDCLIM/V1/BIO\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AzrCwWu39iLw"
   },
   "source": [
    "[**NASA SRTM Digital Elevation 30m**](https://developers.google.com/earth-engine/datasets/catalog/USGS_SRTMGL1_003): This dataset contains digital elevation data from the Shuttle Radar Topography Mission (SRTM). The data was primarily collected around the year 2000 and is provided at a resolution of approximately 30 meters (1 arc-second). The following code calculates elevation, slope, aspect, and hillshade layers from the SRTM data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "O8lPyhWv9hHn"
   },
   "outputs": [],
   "source": [
    "# NASA SRTM Digital Elevation 30m\n",
    "terrain = ee.Algorithms.Terrain(ee.Image(\"USGS/SRTMGL1_003\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "H-kFHdWP9_3N"
   },
   "source": [
    "[**Global Forest Cover Change (GFCC) Tree Cover Multi-Year Global 30m**](https://developers.google.com/earth-engine/datasets/catalog/NASA_MEASURES_GFCC_TC_v3): The Vegetation Continuous Fields (VCF) dataset from Landsat estimates the proportion of vertically projected vegetation cover when the vegetation height is greater than 5 meters. This dataset is provided for four time periods centered around the years 2000, 2005, 2010, and 2015, with a resolution of 30 meters. Here, the median values from these four time periods are used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "-ymsifA098H6"
   },
   "outputs": [],
   "source": [
    "# Global Forest Cover Change (GFCC) Tree Cover Multi-Year Global 30m\n",
    "tcc = ee.ImageCollection(\"NASA/MEASURES/GFCC/TC/v3\")\n",
    "median_tcc = (\n",
    "    tcc.filterDate(\"2000-01-01\", \"2015-12-31\")\n",
    "    .select([\"tree_canopy_cover\"], [\"TCC\"])\n",
    "    .median()\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "CMrmqQ5X-uOB"
   },
   "source": [
    "`bio` (Bioclimatic variables), `terrain` (topography), and `median_tcc` (tree canopy cover) are combined into a single multiband image. The `elevation` band is selected from `terrain`, and a `watermask` is created for locations where `elevation` is greater than `0`. This masks regions below sea level (e.g. the ocean) and prepares the researcher to analyze various environmental factors for the AOI comprehensively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "gSbuqe6k-k6B"
   },
   "outputs": [],
   "source": [
    "# Combine bands into a multi-band image\n",
    "predictors = bio.addBands(terrain).addBands(median_tcc)\n",
    "\n",
    "# Create a water mask\n",
    "watermask = terrain.select('elevation').gt(0)\n",
    "\n",
    "# Mask out ocean pixels and clip to the area of interest\n",
    "predictors = predictors.updateMask(watermask).clip(aoi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WmuXiScSAlxX"
   },
   "source": [
    "When highly correlated predictor variables are included together in a model, multicollinearity issues can arise. Multicollinearity is a phenomenon that occurs when there are strong linear relationships among independent variables in a model, leading to instability in the estimation of the model's coefficients (weights). This instability can reduce the model's reliability and make predictions or interpretations for new data challenging. Therefore, we will consider multicollinearity and proceed with the process of selecting predictor variables.\n",
    "\n",
    "First, we will generate 5,000 random points and then extract the predictor variable values of the single multiband image at those points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "vYLvzPDwAeyk"
   },
   "outputs": [],
   "source": [
    "# Generate 5,000 random points\n",
    "data_cor = predictors.sample(scale=grain_size, numPixels=5000, geometries=True)\n",
    "\n",
    "# Extract predictor variable values\n",
    "pvals = predictors.sampleRegions(collection=data_cor, scale=grain_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oHWmXtTxCmfV"
   },
   "source": [
    "We will convert the extracted predictor values for each point into a DataFrame and then check the first row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "x1Bv0knrBbWn"
   },
   "outputs": [],
   "source": [
    "# Converting predictor values from Earth Engine to a DataFrame\n",
    "pvals_df = geemap.ee_to_df(pvals)\n",
    "pvals_df.head(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "H3K-knJnC65R"
   },
   "outputs": [],
   "source": [
    "# Displaying the columns\n",
    "columns = pvals_df.columns\n",
    "print(columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "f3YiwTb8DOzb"
   },
   "source": [
    "Calculating Spearman correlation coefficients between the given predictor variables and visualizing them in a heatmap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "J72HX3kRPMXH"
   },
   "outputs": [],
   "source": [
    "def plot_correlation_heatmap(dataframe, h_size=10, show_labels=False):\n",
    "    # Calculate Spearman correlation coefficients\n",
    "    correlation_matrix = dataframe.corr(method=\"spearman\")\n",
    "\n",
    "    # Create a heatmap\n",
    "    plt.figure(figsize=(h_size, h_size-2))\n",
    "    plt.imshow(correlation_matrix, cmap='coolwarm', interpolation='nearest')\n",
    "\n",
    "    # Optionally display values on the heatmap\n",
    "    if show_labels:\n",
    "        for i in range(correlation_matrix.shape[0]):\n",
    "            for j in range(correlation_matrix.shape[1]):\n",
    "                plt.text(j, i, f\"{correlation_matrix.iloc[i, j]:.2f}\",\n",
    "                         ha='center', va='center', color='white', fontsize=8)\n",
    "\n",
    "    columns = dataframe.columns.tolist()\n",
    "    plt.xticks(range(len(columns)), columns, rotation=90)\n",
    "    plt.yticks(range(len(columns)), columns)\n",
    "    plt.title(\"Variables Correlation Matrix\")\n",
    "    plt.colorbar(label=\"Spearman Correlation\")\n",
    "    plt.savefig('correlation_heatmap_plot.png')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "DP5ect8cDb7C"
   },
   "outputs": [],
   "source": [
    "# Plot the correlation heatmap of variables\n",
    "plot_correlation_heatmap(pvals_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "07B62CNfDyGz"
   },
   "source": [
    "Spearman correlation coefficient is useful for understanding the general associations among predictor variables but does not directly assess how multiple variables interact, specifically detecting multicollinearity.\n",
    "\n",
    "The **Variance Inflation Factor (VIF below)** is a statistical metric used to evaluate multicollinearity and guide variable selection. It indicates the degree of linear relationship of each independent variable with the other independent variables, and high VIF values can be evidence of multicollinearity.\n",
    "\n",
    "Typically, when VIF values exceed 5 or 10, it suggests that the variable has a strong correlation with other variables, potentially compromising the stability and interpretability of the model. In this tutorial, a criterion of VIF values less than 10 was used for variable selection. The following 6 variables were selected based on VIF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "TJdGzd3SDisO"
   },
   "outputs": [],
   "source": [
    "# Filter variables based on Variance Inflation Factor (VIF)\n",
    "def filter_variables_by_vif(dataframe, threshold=10):\n",
    "\n",
    "    original_columns = dataframe.columns.tolist()\n",
    "    remaining_columns = original_columns[:]\n",
    "\n",
    "    while True:\n",
    "        vif_data = dataframe[remaining_columns]\n",
    "        vif_values = [\n",
    "            variance_inflation_factor(vif_data.values, i)\n",
    "            for i in range(vif_data.shape[1])\n",
    "        ]\n",
    "\n",
    "        max_vif_index = vif_values.index(max(vif_values))\n",
    "        max_vif = max(vif_values)\n",
    "\n",
    "        if max_vif < threshold:\n",
    "            break\n",
    "\n",
    "        print(f\"Removing '{remaining_columns[max_vif_index]}' with VIF {max_vif:.2f}\")\n",
    "\n",
    "        del remaining_columns[max_vif_index]\n",
    "\n",
    "    filtered_data = dataframe[remaining_columns]\n",
    "    bands = filtered_data.columns.tolist()\n",
    "    print(\"Bands:\", bands)\n",
    "\n",
    "    return filtered_data, bands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "e3iiKaK5Eg0q"
   },
   "outputs": [],
   "source": [
    "filtered_pvals_df, bands = filter_variables_by_vif(pvals_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "TEFmil6WErkl"
   },
   "outputs": [],
   "source": [
    "# Variable Selection Based on VIF\n",
    "predictors = predictors.select(bands)\n",
    "\n",
    "# Plot the correlation heatmap of variables\n",
    "plot_correlation_heatmap(filtered_pvals_df, h_size=6, show_labels=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "j0dc6evyFFhw"
   },
   "source": [
    "Next, let's visualize the 6 selected predictor variables on the map.\n",
    "![Predictor Variables for Analysis](https://github.com/google/earthengine-community/blob/master/tutorials/species-distribution-modeling/predictor_variables.png?raw=1)\n",
    "\n",
    "You can explore the available palettes for map visualization using the following code. For example, the `terrain` palette looks like this.\n",
    "```\n",
    "cm.plot_colormaps(width=8.0, height=0.2)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "H8uHX-83E8x7"
   },
   "outputs": [],
   "source": [
    "cm.plot_colormap('terrain', width=8.0, height=0.2, orientation='horizontal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "RswVTmLFFUb3"
   },
   "outputs": [],
   "source": [
    "# Elevation layer\n",
    "Map = geemap.Map(layout={'height':'400px', 'width':'800px'})\n",
    "\n",
    "vis_params = {'bands':['elevation'], 'min': 0, 'max': 1800, 'palette': cm.palettes.terrain}\n",
    "Map.addLayer(predictors, vis_params, 'elevation')\n",
    "Map.add_colorbar(vis_params, label=\"Elevation (m)\", orientation=\"vertical\", layer_name=\"elevation\")\n",
    "Map.centerObject(aoi, 6)\n",
    "Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "DOSMZpbsF8sH"
   },
   "outputs": [],
   "source": [
    "# Calculate the minimum and maximum values for bio09\n",
    "min_max_val = (\n",
    "    predictors.select(\"bio09\")\n",
    "    .multiply(0.1)\n",
    "    .reduceRegion(reducer=ee.Reducer.minMax(), scale=1000)\n",
    "    .getInfo()\n",
    ")\n",
    "\n",
    "# bio09 (Mean temperature of driest quarter) layer\n",
    "Map = geemap.Map(layout={\"height\": \"400px\", \"width\": \"800px\"})\n",
    "\n",
    "vis_params = {\n",
    "    \"min\": math.floor(min_max_val[\"bio09_min\"]),\n",
    "    \"max\": math.ceil(min_max_val[\"bio09_max\"]),\n",
    "    \"palette\": cm.palettes.hot,\n",
    "}\n",
    "Map.addLayer(predictors.select(\"bio09\").multiply(0.1), vis_params, \"bio09\")\n",
    "Map.add_colorbar(\n",
    "    vis_params,\n",
    "    label=\"Mean temperature of driest quarter (℃)\",\n",
    "    orientation=\"vertical\",\n",
    "    layer_name=\"bio09\",\n",
    ")\n",
    "Map.centerObject(aoi, 6)\n",
    "Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "cQpG8qRzIQcw"
   },
   "outputs": [],
   "source": [
    "# Slope layer\n",
    "Map = geemap.Map(layout={'height':'400px', 'width':'800px'})\n",
    "\n",
    "vis_params = {'bands':['slope'], 'min': 0, 'max': 25, 'palette': cm.palettes.RdYlGn_r}\n",
    "Map.addLayer(predictors, vis_params, 'slope')\n",
    "Map.add_colorbar(vis_params, label=\"Slope\", orientation=\"vertical\", layer_name=\"slope\")\n",
    "Map.centerObject(aoi, 6)\n",
    "Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "u1hW_wAXIvvG"
   },
   "outputs": [],
   "source": [
    "# Aspect layer\n",
    "Map = geemap.Map(layout={'height':'400px', 'width':'800px'})\n",
    "\n",
    "vis_params = {'bands':['aspect'], 'min': 0, 'max': 360, 'palette': cm.palettes.rainbow}\n",
    "Map.addLayer(predictors, vis_params, 'aspect')\n",
    "Map.add_colorbar(vis_params, label=\"Aspect\", orientation=\"vertical\", layer_name=\"aspect\")\n",
    "Map.centerObject(aoi, 6)\n",
    "Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "C49-BnnxJGUD"
   },
   "outputs": [],
   "source": [
    "# Calculate the minimum and maximum values for bio14\n",
    "min_max_val = (\n",
    "    predictors.select(\"bio14\")\n",
    "    .reduceRegion(reducer=ee.Reducer.minMax(), scale=1000)\n",
    "    .getInfo()\n",
    ")\n",
    "\n",
    "# bio14 (Precipitation of driest month) layer\n",
    "Map = geemap.Map(layout={\"height\": \"400px\", \"width\": \"800px\"})\n",
    "\n",
    "vis_params = {\n",
    "    \"bands\": [\"bio14\"],\n",
    "    \"min\": math.floor(min_max_val[\"bio14_min\"]),\n",
    "    \"max\": math.ceil(min_max_val[\"bio14_max\"]),\n",
    "    \"palette\": cm.palettes.Blues,\n",
    "}\n",
    "Map.addLayer(predictors, vis_params, \"bio14\")\n",
    "Map.add_colorbar(\n",
    "    vis_params,\n",
    "    label=\"Precipitation of driest month (mm)\",\n",
    "    orientation=\"vertical\",\n",
    "    layer_name=\"bio14\",\n",
    ")\n",
    "Map.centerObject(aoi, 6)\n",
    "Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "zKVnKhLYJicg"
   },
   "outputs": [],
   "source": [
    "# TCC layer\n",
    "Map = geemap.Map(layout={\"height\": \"400px\", \"width\": \"800px\"})\n",
    "\n",
    "vis_params = {\n",
    "    \"bands\": [\"TCC\"],\n",
    "    \"min\": 0,\n",
    "    \"max\": 100,\n",
    "    \"palette\": [\"ffffff\", \"afce56\", \"5f9c00\", \"0e6a00\", \"003800\"],\n",
    "}\n",
    "Map.addLayer(predictors, vis_params, \"TCC\")\n",
    "Map.add_colorbar(\n",
    "    vis_params, label=\"Tree Canopy Cover (%)\", orientation=\"vertical\", layer_name=\"TCC\"\n",
    ")\n",
    "Map.centerObject(aoi, 6)\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xIFAutp6M16q"
   },
   "source": [
    "### Generation of pseudo-absence data\n",
    "\n",
    "In the process of SDM, the selection of input data for a species is mainly approached using two methods:\n",
    "\n",
    "1. **Presence-Background Method**: This method compares the locations where a particular species has been observed (presence) with other locations where the species has not been observed (background). Here, the background data does not necessarily mean areas where the species does not exist but rather is set up to reflect the overall environmental conditions of the study area. It is used to distinguish suitable environments where the species could exist from less suitable ones.\n",
    "\n",
    "2. **Presence-Absence Method**: This method compares locations where the species has been observed (presence) with locations where it has definitively not been observed (absence). Here, absence data represents specific locations where the species is known not to exist. It does not reflect the overall environmental conditions of the study area but rather points to locations where the species is estimated not to exist.\n",
    "\n",
    "In practice, it is often difficult to collect true absence data, so pseudo-absence data generated artificially is frequently used. However, it's important to acknowledge the limitations and potential errors of this method, as artificially generated pseudo-absence points may not accurately reflect true absence areas.\n",
    "\n",
    "The choice between these two methods depends on data availability, research objectives, model accuracy and reliability, as well as time and resources. Here, we will use occurrence data collected from GBIF and artificially generated pseudo-absence data to model using the \"Presence-Absence\" method.\n",
    "\n",
    "The generation of pseudo-absence data will be done through the \"environmental profiling approach\", and the specific steps are as follows:\n",
    "\n",
    "1. Environmental Classification Using k-means Clustering: The k-means clustering algorithm, based on Euclidean distance, will be used to divide the pixels within the study area into two clusters. One cluster will represent areas with similar environmental characteristics to randomly selected 100 presence locations, while the other cluster will represent areas with different characteristics.\n",
    "\n",
    "2. Generation of Pseudo-Absence Data within Dissimilar Clusters: Within the second cluster identified in the first step (which has different environmental characteristics from the presence data), randomly generated pseudo-absence points will be created. These pseudo-absence points will represent locations where the species is not expected to exist."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "hMnMjWpOYG_B"
   },
   "outputs": [],
   "source": [
    "# Randomly select 100 locations for occurrence\n",
    "pvals = predictors.sampleRegions(\n",
    "    collection=data.randomColumn().sort('random').limit(100),\n",
    "    properties=[],\n",
    "    scale=grain_size\n",
    ")\n",
    "\n",
    "# Perform k-means clustering\n",
    "clusterer = ee.Clusterer.wekaKMeans(\n",
    "    nClusters=2,\n",
    "    distanceFunction=\"Euclidean\"\n",
    ").train(pvals)\n",
    "\n",
    "cl_result = predictors.cluster(clusterer)\n",
    "\n",
    "# Get cluster ID for locations similar to occurrence\n",
    "cl_id = cl_result.sampleRegions(\n",
    "    collection=data.randomColumn().sort('random').limit(200),\n",
    "    properties=[],\n",
    "    scale=grain_size\n",
    ")\n",
    "\n",
    "# Define non-occurrence areas in dissimilar clusters\n",
    "cl_id = ee.FeatureCollection(cl_id).reduceColumns(ee.Reducer.mode(),['cluster'])\n",
    "cl_id = ee.Number(cl_id.get('mode')).subtract(1).abs()\n",
    "cl_mask = cl_result.select(['cluster']).eq(cl_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Z6sGrkHGYWat"
   },
   "outputs": [],
   "source": [
    "# Presence location mask\n",
    "presence_mask = data.reduceToImage(properties=['random'],\n",
    "reducer=ee.Reducer.first()\n",
    ").reproject('EPSG:4326', None,\n",
    "            grain_size).mask().neq(1).selfMask()\n",
    "\n",
    "# Masking presence locations in non-occurrence areas and clipping to AOI\n",
    "area_for_pa = presence_mask.updateMask(cl_mask).clip(aoi)\n",
    "\n",
    "# Area for Pseudo-absence\n",
    "Map = geemap.Map(layout={'height':'400px', 'width':'800px'})\n",
    "Map.addLayer(area_for_pa, {'palette': 'black'}, 'AreaForPA')\n",
    "Map.centerObject(aoi, 6)\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "iTWaV5fQZPgA"
   },
   "source": [
    "### Model fitting and prediction\n",
    "\n",
    "We will now divide the data into training data and test data. The training data will be used to find the optimal parameters by training the model, while the test data will be used to evaluate the model trained beforehand. An important concept to consider in this context is spatial autocorrelation.\n",
    "\n",
    "**Spatial autocorrelation** is an essential element in SDM, associated with Tobler's law. It embodies the concept that \"everything is related to everything else, but near things are more related than distant things\". Spatial autocorrelation represents the significant relationship between the location of species and environmental variables. However, if spatial autocorrelation exists between the training and test data, the independence between the two data sets can be compromised. This significantly impacts the evaluation of the model's generalization ability.\n",
    "\n",
    "One method to address this issue is the spatial block cross-validation technique, which involves dividing the data into training and testing datasets. This technique involves dividing the data into multiple blocks, using each block independently as training and test datasets to reduce the impact of spatial autocorrelation. This enhances the independence between datasets, allowing for a more accurate evaluation of the model's generalization ability.\n",
    "\n",
    "The specific procedure is as follows:\n",
    "1. Creation of spatial blocks: Divide the entire dataset into spatial blocks of equal size (e.g., 50x50 km).\n",
    "2. Assignment of training and testing sets: Each spatial block is randomly assigned to either the training set (70%) or the test set (30%). This prevents the model from overfitting to data from specific areas and aims to achieve more generalized results.\n",
    "3. Iterative cross-validation: The entire process is repeated n times (e.g., 10 times). In each iteration, the blocks are randomly divided into training and test sets again, which is intended to improve the model's stability and reliability.\n",
    "4. Generation of pseudo-absence data: In each iteration, pseudo-absence data are randomly generated to evaluate the model's performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "HASln4gfb1pZ"
   },
   "outputs": [],
   "source": [
    "Scale = 50000\n",
    "grid = watermask.reduceRegions(\n",
    "    collection=aoi.coveringGrid(scale=Scale, proj='EPSG:4326'),\n",
    "    reducer=ee.Reducer.mean()).filter(ee.Filter.neq('mean', None))\n",
    "\n",
    "Map = geemap.Map(layout={'height':'400px', 'width':'800px'})\n",
    "Map.addLayer(grid, {}, \"Grid for spatial block cross validation\")\n",
    "Map.addLayer(outline, {'palette': 'FF0000'}, \"Study Area\")\n",
    "Map.centerObject(aoi, 6)\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EsPbyMZ-cJg2"
   },
   "source": [
    "Now we can fit the model. Fitting a model involves understanding the patterns in the data and adjusting the model's parameters (weights and biases) accordingly. This process enables the model to make more accurate predictions when presented with new data. For this purpose, we have defined a function called SDM() to fit the model.\n",
    "\n",
    "We will use the **Random Forest** algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "89R79JvGb6c8"
   },
   "outputs": [],
   "source": [
    "def sdm(x):\n",
    "    seed = ee.Number(x)\n",
    "\n",
    "    # Random block division for training and validation\n",
    "    rand_blk = ee.FeatureCollection(grid).randomColumn(seed=seed).sort(\"random\")\n",
    "    training_grid = rand_blk.filter(ee.Filter.lt(\"random\", split))  # Grid for training\n",
    "    testing_grid = rand_blk.filter(ee.Filter.gte(\"random\", split))  # Grid for testing\n",
    "\n",
    "    # Presence points\n",
    "    presence_points = ee.FeatureCollection(data)\n",
    "    presence_points = presence_points.map(lambda feature: feature.set(\"PresAbs\", 1))\n",
    "    tr_presence_points = presence_points.filter(\n",
    "        ee.Filter.bounds(training_grid)\n",
    "    )  # Presence points for training\n",
    "    te_presence_points = presence_points.filter(\n",
    "        ee.Filter.bounds(testing_grid)\n",
    "    )  # Presence points for testing\n",
    "\n",
    "    # Pseudo-absence points for training\n",
    "    tr_pseudo_abs_points = area_for_pa.sample(\n",
    "        region=training_grid,\n",
    "        scale=grain_size,\n",
    "        numPixels=tr_presence_points.size().add(300),\n",
    "        seed=seed,\n",
    "        geometries=True,\n",
    "    )\n",
    "    # Same number of pseudo-absence points as presence points for training\n",
    "    tr_pseudo_abs_points = (\n",
    "        tr_pseudo_abs_points.randomColumn()\n",
    "        .sort(\"random\")\n",
    "        .limit(ee.Number(tr_presence_points.size()))\n",
    "    )\n",
    "    tr_pseudo_abs_points = tr_pseudo_abs_points.map(lambda feature: feature.set(\"PresAbs\", 0))\n",
    "\n",
    "    te_pseudo_abs_points = area_for_pa.sample(\n",
    "        region=testing_grid,\n",
    "        scale=grain_size,\n",
    "        numPixels=te_presence_points.size().add(100),\n",
    "        seed=seed,\n",
    "        geometries=True,\n",
    "    )\n",
    "    # Same number of pseudo-absence points as presence points for testing\n",
    "    te_pseudo_abs_points = (\n",
    "        te_pseudo_abs_points.randomColumn()\n",
    "        .sort(\"random\")\n",
    "        .limit(ee.Number(te_presence_points.size()))\n",
    "    )\n",
    "    te_pseudo_abs_points = te_pseudo_abs_points.map(lambda feature: feature.set(\"PresAbs\", 0))\n",
    "\n",
    "    # Merge training and pseudo-absence points\n",
    "    training_partition = tr_presence_points.merge(tr_pseudo_abs_points)\n",
    "    testing_partition = te_presence_points.merge(te_pseudo_abs_points)\n",
    "\n",
    "    # Extract predictor variable values at training points\n",
    "    train_pvals = predictors.sampleRegions(\n",
    "        collection=training_partition,\n",
    "        properties=[\"PresAbs\"],\n",
    "        scale=grain_size,\n",
    "        geometries=True,\n",
    "    )\n",
    "\n",
    "    # Random Forest classifier\n",
    "    classifier = ee.Classifier.gmoMaxEnt(\n",
    "        numberOfTrees=500,\n",
    "        variablesPerSplit=None,\n",
    "        minLeafPopulation=10,\n",
    "        bagFraction=0.5,\n",
    "        maxNodes=None,\n",
    "        seed=seed,\n",
    "    )\n",
    "    # Presence probability: Habitat suitability map\n",
    "    classifier_pr = classifier.setOutputMode(\"PROBABILITY\").train(\n",
    "        train_pvals, \"PresAbs\", bands\n",
    "    )\n",
    "    classified_img_pr = predictors.select(bands).classify(classifier_pr)\n",
    "\n",
    "    # Binary presence/absence map: Potential distribution map\n",
    "    classifier_bin = classifier.setOutputMode(\"CLASSIFICATION\").train(\n",
    "        train_pvals, \"PresAbs\", bands\n",
    "    )\n",
    "    classified_img_bin = predictors.select(bands).classify(classifier_bin)\n",
    "\n",
    "    return [\n",
    "        classified_img_pr,\n",
    "        classified_img_bin,\n",
    "        training_partition,\n",
    "        testing_partition,\n",
    "    ], classifier_pr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "yCGB0Y6meW2B"
   },
   "source": [
    "Spatial blocks are divided into 70% for model training and 30% for model testing, respectively. Pseudo-absence data are randomly generated within each training and testing set in every iteration. As a result, each execution yields different sets of presence and pseudo-absence data for model training and testing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "KV6Jg50AeA7B"
   },
   "outputs": [],
   "source": [
    "split = 0.7\n",
    "numiter = 10\n",
    "\n",
    "# Random Seed\n",
    "runif = lambda length: [random.randint(1, 1000) for _ in range(length)]\n",
    "items = runif(numiter)\n",
    "\n",
    "# Fixed seed\n",
    "# items = [287, 288, 553, 226, 151, 255, 902, 267, 419, 538]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "BHCnkZvveq_Q"
   },
   "outputs": [],
   "source": [
    "results_list = [] # Initialize SDM results list\n",
    "importances_list = [] # Initialize variable importance list\n",
    "\n",
    "for item in items:\n",
    "    result, trained = sdm(item)\n",
    "    # Accumulate SDM results into the list\n",
    "    results_list.extend(result)\n",
    "\n",
    "    # Accumulate variable importance into the list\n",
    "    importance = ee.Dictionary(trained.explain()).get('importance')\n",
    "    importances_list.extend(importance.getInfo().items())\n",
    "\n",
    "# Flatten the SDM results list\n",
    "results = ee.List(results_list).flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UHQ6Hx07gG2m"
   },
   "source": [
    "Now we can visualize the **habitat suitability map** and **potential distribution map** for the Fairy pitta. In this case, the habitat suitability map is created by using the `mean()` function to calculate the average for each pixel location across all images, and the potential distribution map is generated by using the `mode()` function to determine the most frequently occurring value at each pixel location across all images.\n",
    "\n",
    "![SDM results](https://github.com/google/earthengine-community/blob/master/tutorials/species-distribution-modeling/sdm_results.png?raw=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "J_QUmonhfuSM"
   },
   "outputs": [],
   "source": [
    "# Habitat suitability map\n",
    "images = ee.List.sequence(\n",
    "    0, ee.Number(numiter).multiply(4).subtract(1), 4).map(\n",
    "    lambda x: results.get(x))\n",
    "model_average = ee.ImageCollection.fromImages(images).mean()\n",
    "\n",
    "Map = geemap.Map(layout={'height':'400px', 'width':'800px'}, basemap='Esri.WorldImagery')\n",
    "\n",
    "vis_params = {\n",
    "    'min': 0,\n",
    "    'max': 1,\n",
    "    'palette': cm.palettes.viridis_r}\n",
    "Map.addLayer(model_average, vis_params, 'Habitat suitability')\n",
    "Map.add_colorbar(vis_params, label=\"Habitat suitability\",\n",
    "                 orientation=\"horizontal\",\n",
    "                 layer_name=\"Habitat suitability\")\n",
    "Map.addLayer(data, {'color':'red'}, 'Presence')\n",
    "Map.centerObject(aoi, 6)\n",
    "Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "cBZa-J43jAL2"
   },
   "outputs": [],
   "source": [
    "# Potential distribution map\n",
    "images2 = ee.List.sequence(1, ee.Number(numiter).multiply(4).subtract(1), 4).map(\n",
    "    lambda x: results.get(x)\n",
    ")\n",
    "distribution_map = ee.ImageCollection.fromImages(images2).mode()\n",
    "\n",
    "Map = geemap.Map(\n",
    "    layout={\"height\": \"400px\", \"width\": \"800px\"}, basemap=\"Esri.WorldImagery\"\n",
    ")\n",
    "\n",
    "vis_params = {\"min\": 0, \"max\": 1, \"palette\": [\"white\", \"green\"]}\n",
    "Map.addLayer(distribution_map, vis_params, \"Potential distribution\")\n",
    "Map.addLayer(data, {\"color\": \"red\"}, \"Presence\")\n",
    "Map.add_colorbar(\n",
    "    vis_params,\n",
    "    label=\"Potential distribution\",\n",
    "    discrete=True,\n",
    "    orientation=\"horizontal\",\n",
    "    layer_name=\"Potential distribution\",\n",
    ")\n",
    "Map.centerObject(data.geometry(), 6)\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "B5cjWyzAmd4Z"
   },
   "source": [
    "### Variable importance and accuracy assessment\n",
    "\n",
    "Random Forest (`ee.Classifier.smileRandomForest`) is one of the ensemble learning methods, which operates by constructing multiple decision trees to make predictions. Each decision tree independently learns from different subsets of the data, and their results are aggregated to enable more accurate and stable predictions.\n",
    "\n",
    "Variable importance is a measure that evaluates the impact of each variable on the predictions within the Random Forest model. We will use the previously defined `importances_list` to calculate and print the average variable importance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "uoJ0mBeRme-E"
   },
   "outputs": [],
   "source": [
    "def plot_variable_importance(importances_list):\n",
    "    # Extract each variable importance value into a list\n",
    "    variables = [item[0] for item in importances_list]\n",
    "    importances = [item[1] for item in importances_list]\n",
    "\n",
    "    # Calculate the average importance for each variable\n",
    "    average_importances = {}\n",
    "    for variable in set(variables):\n",
    "        indices = [i for i, var in enumerate(variables) if var == variable]\n",
    "        average_importance = np.mean([importances[i] for i in indices])\n",
    "        average_importances[variable] = average_importance\n",
    "\n",
    "    # Sort the importances in descending order of importance\n",
    "    sorted_importances = sorted(average_importances.items(),\n",
    "                                key=lambda x: x[1], reverse=False)\n",
    "    variables = [item[0] for item in sorted_importances]\n",
    "    avg_importances = [item[1] for item in sorted_importances]\n",
    "\n",
    "    # Adjust the graph size\n",
    "    plt.figure(figsize=(8, 4))\n",
    "\n",
    "    # Plot the average importance as a horizontal bar chart\n",
    "    plt.barh(variables, avg_importances)\n",
    "    plt.xlabel('Importance')\n",
    "    plt.ylabel('Variables')\n",
    "    plt.title('Average Variable Importance')\n",
    "\n",
    "    # Display values above the bars\n",
    "    for i, v in enumerate(avg_importances):\n",
    "        plt.text(v + 0.02, i, f\"{v:.2f}\", va='center')\n",
    "\n",
    "    # Adjust the x-axis range\n",
    "    plt.xlim(0, max(avg_importances) + 5)  # Adjust to the desired range\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.savefig('variable_importance.png')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "2nhhAUydurwm"
   },
   "outputs": [],
   "source": [
    "plot_variable_importance(importances_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "m889_YCzwNhR"
   },
   "source": [
    "Using the Testing Datasets, we calculate AUC-ROC and AUC-PR for each run. Then, we compute the average AUC-ROC and AUC-PR over n iterations.\n",
    "\n",
    "**AUC-ROC** represents the area under the curve of the 'Sensitivity (Recall) vs. 1-Specificity' graph, illustrating the relationship between sensitivity and specificity as the threshold changes. Specificity is based on all observed non-occurrences. Therefore, AUC-ROC encompasses all quadrants of the confusion matrix.\n",
    "\n",
    "**AUC-PR** represents the area under the curve of the 'Precision vs. Recall (Sensitivity)' graph, showing the relationship between precision and recall as the threshold varies. Precision is based on all predicted occurrences. Hence, AUC-PR does not include the true negatives (TN).\n",
    "\n",
    "> Note: It's important to ensure that each run has a sufficient number of points for model validation. The final number of points may vary due to the random partitioning of spatial blocks, so it's crucial to verify if there are enough presence and pseudo-absence points for model validation. In the case of endangered or rare species, there might be a shortage of occurrence data, leading to an insufficient test dataset. In such cases, alternatives may include additional data collection based on expert knowledge and experience or utilizing relevant alternative data sources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "sdyvrG7swBUp"
   },
   "outputs": [],
   "source": [
    "def print_pres_abs_sizes(TestingDatasets, numiter):\n",
    "    # Check and print the sizes of presence and pseudo-absence coordinates\n",
    "    def get_pres_abs_size(x):\n",
    "        fc = ee.FeatureCollection(TestingDatasets.get(x))\n",
    "        presence_size = fc.filter(ee.Filter.eq(\"PresAbs\", 1)).size()\n",
    "        pseudo_absence_size = fc.filter(ee.Filter.eq(\"PresAbs\", 0)).size()\n",
    "        return ee.List([presence_size, pseudo_absence_size])\n",
    "\n",
    "    sizes_info = (\n",
    "        ee.List.sequence(0, ee.Number(numiter).subtract(1), 1)\n",
    "        .map(get_pres_abs_size)\n",
    "        .getInfo()\n",
    "    )\n",
    "\n",
    "    for i, sizes in enumerate(sizes_info):\n",
    "        presence_size = sizes[0]\n",
    "        pseudo_absence_size = sizes[1]\n",
    "        print(\n",
    "            f\"Iteration {i + 1}: Presence Size = {presence_size}, Pseudo-absence Size = {pseudo_absence_size}\"\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "MccWAsdTxjUS"
   },
   "outputs": [],
   "source": [
    "# Extracting the Testing Datasets\n",
    "testing_datasets = ee.List.sequence(\n",
    "    3, ee.Number(numiter).multiply(4).subtract(1), 4\n",
    ").map(lambda x: results.get(x))\n",
    "\n",
    "print_pres_abs_sizes(testing_datasets, numiter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "GvEH7C9SzI0M"
   },
   "outputs": [],
   "source": [
    "def get_acc(hsm, t_data, grain_size):\n",
    "    pr_prob_vals = hsm.sampleRegions(\n",
    "        collection=t_data, properties=[\"PresAbs\"], scale=grain_size\n",
    "    )\n",
    "    seq = ee.List.sequence(start=0, end=1, count=25)  # Divide 0 to 1 into 25 intervals\n",
    "\n",
    "    def calculate_metrics(cutoff):\n",
    "        # Each element of the seq list is passed as cutoff(threshold value)\n",
    "\n",
    "        # Observed present = TP + FN\n",
    "        pres = pr_prob_vals.filterMetadata(\"PresAbs\", \"equals\", 1)\n",
    "\n",
    "        # TP (True Positive)\n",
    "        tp = ee.Number(\n",
    "            pres.filterMetadata(\"classification\", \"greater_than\", cutoff).size()\n",
    "        )\n",
    "\n",
    "        # TPR (True Positive Rate) = Recall = Sensitivity = TP / (TP + FN) = TP / Observed present\n",
    "        tpr = tp.divide(pres.size())\n",
    "\n",
    "        # Observed absent = FP + TN\n",
    "        abs = pr_prob_vals.filterMetadata(\"PresAbs\", \"equals\", 0)\n",
    "\n",
    "        # FN (False Negative)\n",
    "        fn = ee.Number(\n",
    "            pres.filterMetadata(\"classification\", \"less_than\", cutoff).size()\n",
    "        )\n",
    "\n",
    "        # TNR (True Negative Rate) = Specificity = TN  / (FP + TN) = TN / Observed absent\n",
    "        tn = ee.Number(abs.filterMetadata(\"classification\", \"less_than\", cutoff).size())\n",
    "        tnr = tn.divide(abs.size())\n",
    "\n",
    "        # FP (False Positive)\n",
    "        fp = ee.Number(\n",
    "            abs.filterMetadata(\"classification\", \"greater_than\", cutoff).size()\n",
    "        )\n",
    "\n",
    "        # FPR (False Positive Rate) = FP / (FP + TN) = FP / Observed absent\n",
    "        fpr = fp.divide(abs.size())\n",
    "\n",
    "        # Precision = TP / (TP + FP) = TP / Predicted present\n",
    "        precision = tp.divide(tp.add(fp))\n",
    "\n",
    "        # SUMSS = SUM of Sensitivity and Specificity\n",
    "        sumss = tpr.add(tnr)\n",
    "\n",
    "        return ee.Feature(\n",
    "            None,\n",
    "            {\n",
    "                \"cutoff\": cutoff,\n",
    "                \"TP\": tp,\n",
    "                \"TN\": tn,\n",
    "                \"FP\": fp,\n",
    "                \"FN\": fn,\n",
    "                \"TPR\": tpr,\n",
    "                \"TNR\": tnr,\n",
    "                \"FPR\": fpr,\n",
    "                \"Precision\": precision,\n",
    "                \"SUMSS\": sumss,\n",
    "            },\n",
    "        )\n",
    "\n",
    "    return ee.FeatureCollection(seq.map(calculate_metrics))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "O6JACvJbzrn7"
   },
   "outputs": [],
   "source": [
    "def calculate_and_print_auc_metrics(images, testing_datasets, grain_size, numiter):\n",
    "    # Calculate AUC-ROC and AUC-PR\n",
    "    def calculate_auc_metrics(x):\n",
    "        hsm = ee.Image(images.get(x))\n",
    "        t_data = ee.FeatureCollection(testing_datasets.get(x))\n",
    "        acc = get_acc(hsm, t_data, grain_size)\n",
    "\n",
    "        # Calculate AUC-ROC\n",
    "        x = ee.Array(acc.aggregate_array(\"FPR\"))\n",
    "        y = ee.Array(acc.aggregate_array(\"TPR\"))\n",
    "        x1 = x.slice(0, 1).subtract(x.slice(0, 0, -1))\n",
    "        y1 = y.slice(0, 1).add(y.slice(0, 0, -1))\n",
    "        auc_roc = x1.multiply(y1).multiply(0.5).reduce(\"sum\", [0]).abs().toList().get(0)\n",
    "\n",
    "        # Calculate AUC-PR\n",
    "        x = ee.Array(acc.aggregate_array(\"TPR\"))\n",
    "        y = ee.Array(acc.aggregate_array(\"Precision\"))\n",
    "        x1 = x.slice(0, 1).subtract(x.slice(0, 0, -1))\n",
    "        y1 = y.slice(0, 1).add(y.slice(0, 0, -1))\n",
    "        auc_pr = x1.multiply(y1).multiply(0.5).reduce(\"sum\", [0]).abs().toList().get(0)\n",
    "\n",
    "        return (auc_roc, auc_pr)\n",
    "\n",
    "    auc_metrics = (\n",
    "        ee.List.sequence(0, ee.Number(numiter).subtract(1), 1)\n",
    "        .map(calculate_auc_metrics)\n",
    "        .getInfo()\n",
    "    )\n",
    "\n",
    "    # Print AUC-ROC and AUC-PR for each iteration\n",
    "    df = pd.DataFrame(auc_metrics, columns=[\"AUC-ROC\", \"AUC-PR\"])\n",
    "    df.index = [f\"Iteration {i + 1}\" for i in range(len(df))]\n",
    "    df.to_csv(\"auc_metrics.csv\", index_label=\"Iteration\")\n",
    "    print(df)\n",
    "\n",
    "    # Calculate mean and standard deviation of AUC-ROC and AUC-PR\n",
    "    mean_auc_roc, std_auc_roc = df[\"AUC-ROC\"].mean(), df[\"AUC-ROC\"].std()\n",
    "    mean_auc_pr, std_auc_pr = df[\"AUC-PR\"].mean(), df[\"AUC-PR\"].std()\n",
    "    print(f\"Mean AUC-ROC = {mean_auc_roc:.4f} ± {std_auc_roc:.4f}\")\n",
    "    print(f\"Mean AUC-PR = {mean_auc_pr:.4f} ± {std_auc_pr:.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "WVM4DTduz41j"
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "# Calculate AUC-ROC and AUC-PR\n",
    "calculate_and_print_auc_metrics(images, testing_datasets, grain_size, numiter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QDlaa47s1WjC"
   },
   "source": [
    "This tutorial has provided a practical example of using Google Earth Engine (GEE) for Species Distribution Modeling (SDM). An important takeaway is the versatility and flexibility of GEE in the field of SDM. Leveraging Earth Engine's powerful geospatial data processing capabilities opens up endless possibilities for researchers and conservationists to understand and preserve biodiversity on our planet. By applying the knowledge and skills gained from this tutorial, individuals can explore and contribute to this fascinating field of ecological research."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "Species Distribution Modeling",
   "provenance": []
  },
  "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.12.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}