eeSDM.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7a8574ef-e8a1-4dbd-8f0b-dcbc0a485de1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee\n",
    "import geemap\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from statsmodels.stats.outliers_influence import variance_inflation_factor\n",
    "import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "14c2eca4-a5b2-48bb-98a6-5f9846667a46",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_data_distribution(gdf, h_size=12):\n",
    "    \n",
    "    # Create a figure with two subplots: Yearly data distribution (left) and Monthly data distribution (right)\n",
    "    plt.figure(figsize=(h_size, h_size-8))\n",
    "    plt.subplot(1, 2, 1)\n",
    "    \n",
    "    # Calculate the counts of data for each year and sort by index\n",
    "    year_counts = gdf['year'].value_counts().sort_index()\n",
    "    \n",
    "    # Create a bar plot for yearly data distribution\n",
    "    plt.bar(year_counts.index, year_counts.values)\n",
    "    plt.xlabel('Year')\n",
    "    plt.ylabel('Count')\n",
    "    plt.title('Yearly Data Distribution')\n",
    "    \n",
    "    # Display the data count within each bar\n",
    "    for i, count in enumerate(year_counts.values):\n",
    "        plt.text(year_counts.index[i], count, str(count), ha='center', va='bottom')\n",
    "    \n",
    "    # Create the second subplot for monthly data distribution\n",
    "    plt.subplot(1, 2, 2)\n",
    "    \n",
    "    # Calculate the counts of data for each month and sort by index\n",
    "    month_counts = gdf['month'].value_counts().sort_index()\n",
    "    \n",
    "    # Create a bar plot for monthly data distribution\n",
    "    plt.bar(month_counts.index, month_counts.values)\n",
    "    plt.xlabel('Month')\n",
    "    plt.ylabel('Count')\n",
    "    plt.title('Monthly Data Distribution')\n",
    "    \n",
    "    # Display the data count within each bar for the second subplot\n",
    "    for i, count in enumerate(month_counts.values):\n",
    "        plt.text(month_counts.index[i], count, str(count), ha='center', va='bottom')\n",
    "\n",
    "    # Set x-axis tick labels to integer format\n",
    "    plt.xticks(month_counts.index, map(int, month_counts.index))\n",
    "    \n",
    "    # Adjust layout for better appearance\n",
    "    plt.tight_layout()\n",
    "    \n",
    "    # Save the plot as an image file\n",
    "    plt.savefig('data_distribution_plot.png')\n",
    "    \n",
    "    # Display the plot\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cb86cace-e7ce-4c8b-9e68-d3a4072f07a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_heatmap(gdf, h_size=8):\n",
    "    # Calculate necessary statistics\n",
    "    statistics = gdf.groupby([\"month\", \"year\"]).size().unstack(fill_value=0)\n",
    "    \n",
    "    # Visualize statistics using a heatmap\n",
    "    plt.figure(figsize=(h_size, h_size-6))\n",
    "    heatmap = plt.imshow(statistics.values, cmap=\"YlOrBr\", origin=\"upper\", aspect=\"auto\")\n",
    "\n",
    "    # Display values on top of each pixel\n",
    "    for i in range(len(statistics.index)):\n",
    "        for j in range(len(statistics.columns)):\n",
    "            plt.text(j, i, statistics.values[i, j], ha=\"center\", va=\"center\", color=\"black\")\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()\n",
    "    print(gdf.groupby([\"month\", \"year\"]).size().unstack(fill_value=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "086a0397-d197-4db6-bf69-bd248abf706d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def remove_duplicates(data_raw, GrainSize):\n",
    "    # Select one presence record per pixel at the chosen spatial resolution (1km) randomly\n",
    "    # Generate a random raster image and reproject it to the specified coordinate system and resolution\n",
    "    random_raster = ee.Image.random().reproject('EPSG:4326', None, GrainSize)\n",
    "    \n",
    "    # Sample presence points with the generated random raster\n",
    "    # Scale parameter is set to 10 for sampling, geometries are included\n",
    "    rand_point_vals = random_raster.sampleRegions(collection=ee.FeatureCollection(data_raw), scale=10, geometries=True)\n",
    "    \n",
    "    # Keep only distinct presence records based on the 'random' property\n",
    "    return rand_point_vals.distinct('random')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fb734ff6-6e0d-4c00-ad35-700f952e050e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_correlation_heatmap(df, h_size=10):\n",
    "    # Calculate Spearman correlation coefficients\n",
    "    correlation_matrix = df.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",
    "    # Display values on the heatmap\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)  # Adjust fontsize\n",
    "\n",
    "    columns = df.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": 6,
   "id": "5e46aaf9-378f-42d7-b1db-9a19a66ee4c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def filter_variables_by_vif(df, threshold=10):\n",
    "    # Store the original column names\n",
    "    original_columns = df.columns.tolist()\n",
    "    \n",
    "    # Create a copy of column names to track remaining variables\n",
    "    remaining_columns = original_columns[:]\n",
    "    \n",
    "    # Perform VIF-based variable selection iteratively\n",
    "    while True:\n",
    "        # Create a subset of the DataFrame using remaining variables\n",
    "        vif_data = df[remaining_columns]\n",
    "        \n",
    "        # Calculate VIF values for each remaining variable\n",
    "        vif_values = [variance_inflation_factor(vif_data.values, i) for i in range(vif_data.shape[1])]\n",
    "        \n",
    "        # Find the index of the variable with the highest VIF\n",
    "        max_vif_index = vif_values.index(max(vif_values))\n",
    "        max_vif = max(vif_values)\n",
    "        \n",
    "        # Check if the highest VIF is below the specified threshold\n",
    "        if max_vif < threshold:\n",
    "            break\n",
    "        \n",
    "        # Print information about the variable being removed\n",
    "        print(f\"Removing '{remaining_columns[max_vif_index]}' with VIF {max_vif:.2f}\")\n",
    "        \n",
    "        # Remove the variable with the highest VIF from the list of remaining variables\n",
    "        del remaining_columns[max_vif_index]\n",
    "    \n",
    "    # Create a new DataFrame with the filtered variables\n",
    "    filtered_data = df[remaining_columns]\n",
    "    bands = filtered_data.columns.tolist()\n",
    "    print('Bands:', bands)\n",
    "    \n",
    "    return filtered_data, bands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b1ebe380-8914-4c18-b36d-b57eccd7e87d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_pa_full_area(Data, GrainSize, AOI):\n",
    "    presence_mask = Data.reduceToImage(\n",
    "        properties=['random'],\n",
    "        reducer=ee.Reducer.first()\n",
    "    ).reproject('EPSG:4326', None, ee.Number(GrainSize)).mask().neq(1).selfMask()\n",
    "    AreaForPA = presence_mask.updateMask(ee.Image(\"USGS/SRTMGL1_003\").gt(0)).clip(AOI)\n",
    "    return AreaForPA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "17d59f18-b34b-40c4-b8ea-f7675aa5abf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_pa_spatial_constraint(Data, GrainSize, AOI, distance=50000, maxError=1000):\n",
    "    presence_mask = Data.reduceToImage(\n",
    "        properties=['random'],\n",
    "        reducer=ee.Reducer.first()\n",
    "    ).reproject('EPSG:4326', None, ee.Number(GrainSize)).mask().neq(1).selfMask()\n",
    "    presence_buffer_mask = Data.geometry().buffer(distance, maxError)\n",
    "    AreaForPA = presence_mask.clip(presence_buffer_mask).updateMask(ee.Image(\"USGS/SRTMGL1_003\").gt(0)).clip(AOI)\n",
    "    return AreaForPA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8aff2805-7130-497b-9e95-9f87408dfb3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_pa_environmental_profiling(Data, GrainSize, AOI, predictors):\n",
    "    presence_mask = Data.reduceToImage(\n",
    "        properties=['random'],\n",
    "        reducer=ee.Reducer.first()\n",
    "    ).reproject('EPSG:4326', None, ee.Number(GrainSize)).mask().neq(1).selfMask()\n",
    "    # Extract predictor variable values from a random subset of the presence data\n",
    "    PixelVals = predictors.sampleRegions(\n",
    "        collection=Data.randomColumn().sort('random').limit(100),\n",
    "        properties=[],\n",
    "        tileScale=16,\n",
    "        scale=GrainSize\n",
    "    )\n",
    "    \n",
    "    # Train a k-means clustering model\n",
    "    clusterer = ee.Clusterer.wekaKMeans(\n",
    "        nClusters=2,\n",
    "        distanceFunction=\"Euclidean\"\n",
    "    ).train(PixelVals)\n",
    "    \n",
    "    # Use the trained clusterer to assign pixels to clusters\n",
    "    Clresult = predictors.cluster(clusterer)\n",
    "    \n",
    "    # Obtain cluster IDs similar to those of the presence data\n",
    "    clustID = Clresult.sampleRegions(\n",
    "        collection=Data.randomColumn().sort('random').limit(200),\n",
    "        properties=[],\n",
    "        tileScale=16,\n",
    "        scale=GrainSize\n",
    "    )\n",
    "    \n",
    "    # Use the opposite cluster, define pseudo-absence allowed area\n",
    "    clustID = ee.FeatureCollection(clustID).reduceColumns(ee.Reducer.mode(), ['cluster'])\n",
    "    clustID = ee.Number(clustID.get('mode')).subtract(1).abs()\n",
    "    cl_mask = Clresult.select(['cluster']).eq(clustID)\n",
    "    AreaForPA = presence_mask.updateMask(cl_mask).clip(AOI)\n",
    "    return AreaForPA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2f2072d3-5eec-4906-b80c-1fa42b7051ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "def createGrid(geometry, scale):\n",
    "    # Create an image with pixel longitude and latitude vals\n",
    "    lonLat = ee.Image.pixelLonLat()\n",
    "    \n",
    "    # Convert longitude and latitude images to integer grids\n",
    "    lonGrid = lonLat.select('longitude').multiply(100000).toInt()\n",
    "    latGrid = lonLat.select('latitude').multiply(100000).toInt()\n",
    "    \n",
    "    # Reduce grids to vectors to create the grid polygons\n",
    "    rawGrid = lonGrid.multiply(latGrid).reduceToVectors(\n",
    "        geometry=geometry.buffer(distance=20000, maxError=1000),\n",
    "        scale=scale,\n",
    "        geometryType='polygon')\n",
    "    \n",
    "    # Apply the grid to a raster image to obtain mean values within each grid cell\n",
    "    Grid = ee.Image(\"USGS/SRTMGL1_003\").gt(0).reduceRegions(\n",
    "        collection=rawGrid,\n",
    "        reducer=ee.Reducer.mean()).filter(ee.Filter.neq('mean', None))\n",
    "    \n",
    "    return Grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "722f8154-7966-4ccf-9826-fe7f10cb24ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "def batchSDM(Grid, Data, AreaForPA, GrainSize, bands, predictors, numiter, split=0.7, seed=None):\n",
    "    # Fit SDM\n",
    "    # Generate a list of random seeds\n",
    "    if seed is not None:\n",
    "        random.seed(seed)\n",
    "    random_seeds = [random.randint(1, 1000) for _ in range(numiter)]\n",
    "\n",
    "    def SDM(x, split):\n",
    "        Seed = ee.Number(x)\n",
    "    \n",
    "        # Random block partition for training and validation\n",
    "        GRID = ee.FeatureCollection(Grid).randomColumn(seed=Seed).sort('random')\n",
    "        TrainingGrid = GRID.filter(ee.Filter.lt('random', split))  # Training grids\n",
    "        TestingGrid = GRID.filter(ee.Filter.gte('random', split))  # Testing grids\n",
    "        \n",
    "        # Presence points\n",
    "        PresencePoints = ee.FeatureCollection(Data)\n",
    "        PresencePoints = PresencePoints.map(lambda feature: feature.set('PresAbs', 1))\n",
    "        TrPresencePoints = PresencePoints.filter(ee.Filter.bounds(TrainingGrid))  # Training presence points\n",
    "        TePresencePoints = PresencePoints.filter(ee.Filter.bounds(TestingGrid))  # Testing presence points\n",
    "        \n",
    "        # Pseudo-absence points\n",
    "        TrPseudoAbsPoints = AreaForPA.sample(region=TrainingGrid,\n",
    "                                             scale=GrainSize,\n",
    "                                             numPixels=TrPresencePoints.size().add(300),\n",
    "                                             seed=Seed,\n",
    "                                             geometries=True)\n",
    "        TrPseudoAbsPoints = TrPseudoAbsPoints.randomColumn().sort('random').limit(ee.Number(TrPresencePoints.size()))\n",
    "        TrPseudoAbsPoints = TrPseudoAbsPoints.map(lambda feature: feature.set('PresAbs', 0))\n",
    "        \n",
    "        TePseudoAbsPoints = AreaForPA.sample(region=TestingGrid,\n",
    "                                             scale=GrainSize,\n",
    "                                             numPixels=TePresencePoints.size().add(100),\n",
    "                                             seed=Seed,\n",
    "                                             geometries=True)\n",
    "        TePseudoAbsPoints = TePseudoAbsPoints.randomColumn().sort('random').limit(ee.Number(TePresencePoints.size()))\n",
    "        TePseudoAbsPoints = TePseudoAbsPoints.map(lambda feature: feature.set('PresAbs', 0))\n",
    "    \n",
    "        # Merge training presence and pseudo-absence points\n",
    "        trainingPartition = TrPresencePoints.merge(TrPseudoAbsPoints)\n",
    "        testingPartition = TePresencePoints.merge(TePseudoAbsPoints)\n",
    "    \n",
    "        # Extract covariate values of predictor images from training partition\n",
    "        trainPixelVals = predictors.sampleRegions(collection=trainingPartition,\n",
    "                                                  properties=['PresAbs'],\n",
    "                                                  scale=GrainSize,\n",
    "                                                  tileScale=16,\n",
    "                                                  geometries=True)\n",
    "    \n",
    "        # Random Forest classifier\n",
    "        Classifier = ee.Classifier.smileRandomForest(\n",
    "            numberOfTrees=500,\n",
    "            variablesPerSplit=None,\n",
    "            minLeafPopulation=10,\n",
    "            bagFraction=0.5,\n",
    "            maxNodes=None,\n",
    "            seed=Seed\n",
    "        )\n",
    "        \n",
    "        # Presence probability classifier\n",
    "        ClassifierPr = Classifier.setOutputMode('PROBABILITY').train(trainPixelVals, 'PresAbs', bands)\n",
    "        ClassifiedImgPr = predictors.select(bands).classify(ClassifierPr)\n",
    "        \n",
    "        # Binary presence/absence classifier\n",
    "        ClassifierBin = Classifier.setOutputMode('CLASSIFICATION').train(trainPixelVals, 'PresAbs', bands)\n",
    "        ClassifiedImgBin = predictors.select(bands).classify(ClassifierBin)\n",
    "    \n",
    "        # Variable importance\n",
    "        varImportance = ClassifierPr.explain().get('importance')\n",
    "    \n",
    "        return [varImportance, ClassifiedImgPr, ClassifiedImgBin, trainingPartition, testingPartition]\n",
    "    \n",
    "    # Run SDM for each seed and concatenate results\n",
    "    results_list = [SDM(x, split) for x in random_seeds]\n",
    "    flattened_results = ee.List(results_list).flatten()\n",
    "    \n",
    "    return flattened_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "fc444f6f-8b5e-4639-bbb9-28e0c47acb35",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_avg_variable_importance(results, numiter, h_size=8):\n",
    "    # Calculate Average Variable Importance\n",
    "    varImportance = ee.List.sequence(0, ee.Number(numiter).multiply(5).subtract(1), 5).map(lambda x: results.get(x)).getInfo()\n",
    "    avg_varImportance = {key: sum(d[key] for d in varImportance) / len(varImportance) for key in varImportance[0]}\n",
    "\n",
    "    # Sort variables by importance in descending order\n",
    "    avg_varImportance = {key: value for key, value in sorted(avg_varImportance.items(), key=lambda item: item[1], reverse=False)}\n",
    "\n",
    "    # Create a horizontal bar plot\n",
    "    plt.figure(figsize=(h_size, h_size-4))  # Adjust the size as desired\n",
    "    plt.barh(list(avg_varImportance.keys()), list(avg_varImportance.values()))\n",
    "    plt.xlabel('Importance')\n",
    "    plt.ylabel('Variable')\n",
    "    plt.title('Average Variable Importance')\n",
    "\n",
    "    # Display values on top of bars (rounded to two decimal places)\n",
    "    for i, value in enumerate(avg_varImportance.values()):\n",
    "        plt.text(value, i, f'{value:.2f}', va='center')\n",
    "\n",
    "    # Extend x-axis range by 5 units beyond the maximum value\n",
    "    max_value = max(avg_varImportance.values())\n",
    "    plt.xlim(0, max_value + 5)\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.savefig('avg_variable_importance_plot.png')\n",
    "    plt.show()\n",
    "    print(avg_varImportance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "156a4b0d-1863-49b9-862b-b76dee549c1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_pres_abs_sizes(TestingDatasets, numiter):\n",
    "    # Checking if there are sufficient presence and pseudo-absence points\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 = ee.List.sequence(0, ee.Number(numiter).subtract(1), 1).map(get_pres_abs_size).getInfo()\n",
    "    \n",
    "    for i, sizes in enumerate(sizes_info):\n",
    "        presence_size = sizes[0]\n",
    "        pseudo_absence_size = sizes[1]\n",
    "        print(f'Iteration {i + 1}: Presence Size = {presence_size}, Pseudo-absence Size = {pseudo_absence_size}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c66c58b1-7235-4b7f-a972-3ca84b140115",
   "metadata": {},
   "outputs": [],
   "source": [
    "def getAcc(HSM, TData, GrainSize):\n",
    "    Pr_Prob_Vals = HSM.sampleRegions(collection=TData, properties=['PresAbs'], scale=GrainSize, tileScale=16)\n",
    "    seq = ee.List.sequence(start=0, end=1, count=25) # Divide 0 to 1 into 25 intervals\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(Pres.filterMetadata('classification', 'greater_than', cutoff).size())\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(Pres.filterMetadata('classification', 'less_than', cutoff).size())\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(Abs.filterMetadata('classification', 'greater_than', cutoff).size())\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": 15,
   "id": "3411ab4c-31c6-49a9-ac73-8b7e6b775f80",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_and_print_auc_metrics(images, TestingDatasets, GrainSize, numiter):\n",
    "    # Calculate AUC-ROC and AUC-PR\n",
    "    def calculate_auc_metrics(x):\n",
    "        HSM = ee.Image(images.get(x))\n",
    "        TData = ee.FeatureCollection(TestingDatasets.get(x))\n",
    "        Acc = getAcc(HSM, TData, GrainSize)\n",
    "\n",
    "        # 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",
    "        # 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 = ee.List.sequence(0, ee.Number(numiter).subtract(1), 1).map(calculate_auc_metrics).getInfo()\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": 16,
   "id": "d9f2ecba-11ae-4edc-b28d-613ca93deff6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_and_print_ss_metrics(images, TestingDatasets, GrainSize, numiter):\n",
    "    # Calculate Sensitivity and Specificity\n",
    "    def calculate_ss_metrics(x):\n",
    "        HSM = ee.Image(images.get(x))\n",
    "        TData = ee.FeatureCollection(TestingDatasets.get(x))\n",
    "        Acc = getAcc(HSM, TData, GrainSize)\n",
    "\n",
    "        # Sensitivity\n",
    "        sensitivity = Acc.sort('SUMSS', False).first().get('TPR')\n",
    "        \n",
    "        # Specificity\n",
    "        specificity = Acc.sort('SUMSS', False).first().get('TNR')\n",
    "        \n",
    "        return (sensitivity, specificity)\n",
    "        \n",
    "    ss_metrics = ee.List.sequence(0, ee.Number(numiter).subtract(1), 1).map(calculate_ss_metrics).getInfo()\n",
    "\n",
    "    # Print Sensitivity and Specificity for each iteration\n",
    "    df = pd.DataFrame(ss_metrics, columns=['Sensitivity', 'Specificity'])\n",
    "    df.index = [f'Iteration {i + 1}' for i in range(len(df))]\n",
    "    df.to_csv('ss_metrics.csv', index_label='Iteration')\n",
    "    print(df)\n",
    "\n",
    "    # Calculate mean and standard deviation of Sensitivity and Specificity\n",
    "    mean_sensitivity, std_sensitivity = df['Sensitivity'].mean(), df['Sensitivity'].std()\n",
    "    mean_specificity, std_specificity = df['Specificity'].mean(), df['Specificity'].std()\n",
    "    print(f'Mean Sensitivity = {mean_sensitivity:.4f} ± {std_sensitivity:.4f}')\n",
    "    print(f'Mean Specificity = {mean_specificity:.4f} ± {std_specificity:.4f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "c63dcfb2-a1bc-45e0-a53f-66e7bda1240b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_roc_pr_curves(images, TestingDatasets, GrainSize, numiter):\n",
    "    # Plot ROC and PR curves\n",
    "    def get_roc_pr_values(x):\n",
    "        HSM = ee.Image(images.get(x))\n",
    "        TData = ee.FeatureCollection(TestingDatasets.get(x))        \n",
    "        Acc = getAcc(HSM, TData, GrainSize)\n",
    "\n",
    "        # Get ROC and PR values\n",
    "        varFPR = Acc.aggregate_array('FPR')\n",
    "        varTPR = Acc.aggregate_array('TPR')\n",
    "        varPrecision = Acc.aggregate_array('Precision')\n",
    "        \n",
    "        return (varFPR, varTPR, varPrecision)\n",
    "\n",
    "    # Retrieve ROC and PR values for each iteration and model\n",
    "    roc_pr_values = np.array(ee.List.sequence(0, ee.Number(numiter).subtract(1), 1).map(get_roc_pr_values).getInfo())\n",
    "\n",
    "    all_model_data = []\n",
    "    for model_data in roc_pr_values:\n",
    "        # Transpose data to match the desired format\n",
    "        transposed_data = np.transpose(model_data)\n",
    "        all_model_data.append(transposed_data)\n",
    "    \n",
    "    # Calculate mean and standard deviation across models\n",
    "    mean_data = np.mean(all_model_data, axis=0)\n",
    "    std_data = np.std(all_model_data, axis=0)\n",
    "    \n",
    "    # Create DataFrames for mean and standard deviation data\n",
    "    df_mean = pd.DataFrame(mean_data, columns=['FPR', 'TPR', 'Precision'])\n",
    "    df_std = pd.DataFrame(std_data, columns=['FPR', 'TPR', 'Precision'])\n",
    "    \n",
    "    # Set font size and style\n",
    "    plt.rcParams.update({'font.size': 14})\n",
    "    \n",
    "    # Create a figure with two subplots\n",
    "    plt.figure(figsize=(16, 6))\n",
    "    \n",
    "    # Plot mean ROC Curve\n",
    "    plt.subplot(1, 2, 1)\n",
    "    plt.plot(df_mean['FPR'], df_mean['TPR'], color='blue', lw=2, label='Mean ROC Curve')\n",
    "    plt.fill_between(df_mean['FPR'], df_mean['TPR'] - df_std['TPR'], df_mean['TPR'] + df_std['TPR'], color='gray', alpha=0.2)\n",
    "    plt.plot([0, 1], [0, 1], linestyle='--', color='gray')\n",
    "    plt.xlim([0, 1])\n",
    "    plt.ylim([0, 1])\n",
    "    plt.xlabel('False Positive Rate (FPR)')\n",
    "    plt.ylabel('True Positive Rate (TPR)')\n",
    "    plt.title('Receiver Operating Characteristic (ROC) curve')\n",
    "    plt.grid(True)\n",
    "    plt.legend(loc='lower right')\n",
    "    \n",
    "    # Plot mean PR Curve\n",
    "    plt.subplot(1, 2, 2)\n",
    "    plt.plot(df_mean['TPR'], df_mean['Precision'], color='blue', lw=2, label='Mean PR Curve')\n",
    "    plt.fill_between(df_mean['TPR'], df_mean['Precision'] - df_std['Precision'], df_mean['Precision'] + df_std['Precision'], color='gray', alpha=0.2)\n",
    "    plt.xlim([0, 1])\n",
    "    plt.ylim([0, 1])\n",
    "    plt.xlabel('Recall')\n",
    "    plt.ylabel('Precision')\n",
    "    plt.title('Precision-Recall (PR) curve')\n",
    "    plt.grid(True)\n",
    "    plt.legend(loc='lower right')\n",
    "    plt.tight_layout()\n",
    "    plt.savefig('roc_pr_curves_plot.png')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "1272c630-73a7-4e28-855b-888b60e9b35b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_DistributionMap2(images, TestingDatasets, GrainSize, numiter, ModelAverage):\n",
    "    # Potential Distribution Map using the optimal threshold\n",
    "    def get_metrics(x):\n",
    "        HSM = ee.Image(images.get(x))\n",
    "        TData = ee.FeatureCollection(TestingDatasets.get(x))\n",
    "        Acc = getAcc(HSM, TData, GrainSize)\n",
    "        return Acc.sort('SUMSS', False).first()\n",
    "    \n",
    "    metrics = ee.List.sequence(0, ee.Number(numiter).subtract(1), 1).map(get_metrics)\n",
    "    MeanThresh = ee.Number(ee.FeatureCollection(metrics).aggregate_array(\"cutoff\").reduce(ee.Reducer.mean()))\n",
    "    print('Mean threshold:', MeanThresh.getInfo())\n",
    "    \n",
    "    DistributionMap2 = ModelAverage.gte(MeanThresh)\n",
    "    \n",
    "    return DistributionMap2"
   ]
  }
 ],
 "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.10.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}