PNG_MPM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bbd73003",
   "metadata": {},
   "source": [
    "# Spatio-temporal mineral prospectivity modelling of New Guinea\n",
    "\n",
    "### Ehsan Farahbakhsh<sup>1</sup>, Sabin Zahirovic<sup>1</sup>, Brent I. A. McInnes<sup>2</sup>, Sara Polanco<sup>1</sup>, Fabian Kohlmann<sup>3</sup>, Maria Seton<sup>1</sup>, R. Dietmar M&uuml;ller<sup>1</sup>\n",
    "\n",
    "<sup>1</sup>*EarthByte Group, School of Geosciences, The University of Sydney, Sydney, Australia*\n",
    "\n",
    "<sup>2</sup>*John de Laeter Centre, Faculty of Science and Engineering, Curtin University, Perth, Australia*\n",
    "\n",
    "<sup>3</sup>*Lithodat Pty. Ltd., Melbourne, Australia*\n",
    "\n",
    "This notebook enables the user to create a spatio-temporal mineral prospectivity model of New Guinea. It comprises two main sections; in the first section, kinematic features are extracted, and in the second section, machine learning algorithms are applied to create a prospectivity model.\n",
    "\n",
    "### Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87623f46",
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup the working environment\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import cartopy.crs as ccrs\n",
    "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n",
    "import cmcrameri.cm as ccm\n",
    "from collections import deque\n",
    "import contextily as cx\n",
    "import cv2\n",
    "import geopandas as gpd\n",
    "import gplately\n",
    "from ipywidgets import interact\n",
    "import math\n",
    "import matplotlib as mpl\n",
    "import matplotlib.cm as cm\n",
    "from matplotlib.gridspec import GridSpec\n",
    "from matplotlib.lines import Line2D\n",
    "import matplotlib.patheffects as pe\n",
    "from matplotlib.patches import Patch\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as mticker\n",
    "from netCDF4 import Dataset\n",
    "import numpy as np\n",
    "import os\n",
    "from osgeo import gdal\n",
    "from osgeo import osr\n",
    "import pandas as pd\n",
    "import pickle\n",
    "import plotly.express as px\n",
    "from ptt.subduction_convergence import subduction_convergence_over_time\n",
    "from pulearn import BaggingPuClassifier\n",
    "from gplately import pygplates\n",
    "from rasterio.plot import show\n",
    "import rioxarray as rxr\n",
    "from scipy import ndimage, stats\n",
    "from scipy.cluster import hierarchy\n",
    "from scipy.interpolate import griddata, make_interp_spline\n",
    "from scipy.ndimage import gaussian_filter1d\n",
    "import scipy.spatial\n",
    "import seaborn as sns\n",
    "import shapefile\n",
    "from shapely.geometry import LineString, Point\n",
    "import statistics\n",
    "\n",
    "# machine learning\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.experimental import enable_iterative_imputer\n",
    "from sklearn.impute import IterativeImputer\n",
    "from sklearn.metrics import accuracy_score, auc, confusion_matrix, f1_score, precision_score, recall_score, roc_auc_score, roc_curve\n",
    "from sklearn.metrics import ConfusionMatrixDisplay\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import MinMaxScaler, StandardScaler\n",
    "from skopt import BayesSearchCV\n",
    "from skopt.space import Categorical, Integer, Real\n",
    "\n",
    "# from lib_stamp_muller2016 import *\n",
    "from lib_stamp_muller2019 import *\n",
    "\n",
    "# load parameters\n",
    "# from parameters_muller2016 import parameters\n",
    "from parameters_muller2019 import parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c864753a",
   "metadata": {},
   "source": [
    "### Extract Convergence Kinematic Features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "221235cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "plate_motion_model = parameters['plate_motion_model']\n",
    "\n",
    "# get start time, end time, and time step from parameters.py\n",
    "time_period = parameters['time']\n",
    "start_time = time_period['start']\n",
    "end_time = time_period['end']\n",
    "time_step = time_period['step']\n",
    "time_steps = list(range(start_time, end_time+1, time_step))\n",
    "\n",
    "conv_dir = parameters['convergence_data_dir']\n",
    "conv_prefix = parameters['convergence_data_filename_prefix']\n",
    "conv_ext = parameters['convergence_data_filename_ext']\n",
    "\n",
    "trench_points_features(\n",
    "    start_time,\n",
    "    end_time,\n",
    "    time_step,\n",
    "    conv_dir,\n",
    "    conv_prefix,\n",
    "    conv_ext,\n",
    "    plate_motion_model,\n",
    "    random_state=42\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3c90744",
   "metadata": {},
   "source": [
    "### Create the Plate Reconstruction Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1d7b53f",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plate_motion_model == 'muller2016':\n",
    "    rotation_files = parameters['rotation_file']\n",
    "    topology_files = parameters['topology_file']\n",
    "elif plate_motion_model == 'muller2019':\n",
    "    rotation_files = [os.path.join(dirpath, f) for (dirpath, dirnames, filenames) in os.walk(parameters['rotation_dir']) for f in filenames]\n",
    "    topology_files = [os.path.join(dirpath, f) for (dirpath, dirnames, filenames) in os.walk(parameters['topology_dir']) for f in filenames]\n",
    "\n",
    "coastlines = parameters['coastlines_file']\n",
    "static_polygons = parameters['static_polygons_file']\n",
    "continents = parameters['coastlines_file']\n",
    "cob = parameters['cob_file']\n",
    "\n",
    "rotation_model = pygplates.RotationModel(rotation_files)\n",
    "\n",
    "topology_features = pygplates.FeatureCollection()\n",
    "for topology_file in topology_files:\n",
    "    topology_features.add(pygplates.FeatureCollection(topology_file))\n",
    "\n",
    "# use the PlateReconstruction object to create a plate motion model\n",
    "model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74917245",
   "metadata": {},
   "source": [
    "### Plot Kinematic Features in a Global Scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf3ec89b",
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_features = parameters['selected_features']\n",
    "selected_features_plot = selected_features.copy()\n",
    "selected_features_plot.remove('distance_deg')\n",
    "\n",
    "agegrid_dir = parameters['agegrid_dir']\n",
    "\n",
    "extent_globe = [-180, 180, -90, 90]\n",
    "\n",
    "@interact\n",
    "def show_map(time=time_steps, feature=selected_features_plot):\n",
    "    # call the PlotTopologies object\n",
    "    gplot = gplately.PlotTopologies(model, coastlines, continents, cob, time=time)\n",
    "    \n",
    "    if plate_motion_model == 'muller2016':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2016_AREPS_v1.17_AgeGrid-{time}.nc'\n",
    "    elif plate_motion_model == 'muller2019':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2019_Tectonics_v2.0_AgeGrid-{time}.nc'\n",
    "    \n",
    "    agegrid = gplately.grids.read_netcdf_grid(agegrid_file)\n",
    "    \n",
    "    features_t = pd.read_csv(f'{conv_dir}/{conv_prefix}_{time}.00.{conv_ext}', index_col=False)\n",
    "\n",
    "    # dual colour bars\n",
    "    fig = plt.figure(figsize=(16, 12))\n",
    "    gs = GridSpec(2, 2, hspace=-0.4, wspace=0.4, height_ratios=[1, 0.02])\n",
    "    ax = fig.add_subplot(gs[0, :], projection=ccrs.Mollweide(central_longitude=150))\n",
    "    \n",
    "    ax.set_global()\n",
    "\n",
    "    im = gplot.plot_grid(ax, agegrid.data, cmap=ccm.lapaz_r, vmin=0, vmax=230, alpha=0.5, zorder=1) # cmap: Blues, viridis, winter\n",
    "\n",
    "    gplot.plot_continents(ax, edgecolor='none', facecolor='tan', zorder=2) # facecolor: tan\n",
    "    gplot.plot_ridges(ax, color='red', alpha=0.5, zorder=3)\n",
    "    gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=True, alpha=0.1, zorder=4)\n",
    "\n",
    "    sc = ax.scatter(features_t['trench_lon'], features_t['trench_lat'], 50, marker='.',\n",
    "                    c=features_t[feature], cmap=ccm.hawaii_r, transform=ccrs.PlateCarree(), zorder=5) # cmap: Spectral_r, YlOrRd\n",
    "\n",
    "    gplot.plot_trenches(ax, color='k', alpha=0.3, zorder=6)\n",
    "    gplot.plot_subduction_teeth(ax, spacing=0.05, color='k', alpha=0.3, zorder=7)\n",
    "    \n",
    "    ax.gridlines(linestyle=':')\n",
    "    \n",
    "    cax1 = fig.add_subplot(gs[1, 0])\n",
    "    cax2 = fig.add_subplot(gs[1, 1])\n",
    "    \n",
    "    fig.colorbar(sc, cax=cax2, orientation='horizontal', label=feature, extend='both')\n",
    "    fig.colorbar(im, cax=cax1, orientation='horizontal', label='Seafloor age (Ma)', extend='max')\n",
    "    \n",
    "    # Define custom legend handles\n",
    "    custom_handles = [\n",
    "        Patch(facecolor='tan', edgecolor='none', label='Continental Crust'),  # Custom handle for the filled polygon\n",
    "        Line2D([0], [0], color='red', lw=2, label='Mid-Ocean Ridge')  # Custom handle for the line (ridge)\n",
    "    ]\n",
    "\n",
    "    # Add the custom legend to the plot\n",
    "    ax.legend(handles=custom_handles, loc='lower left')\n",
    "    \n",
    "    ax.set_title(f'Subduction Zones {time} Ma')\n",
    "        \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8a1ecfd",
   "metadata": {},
   "source": [
    "### Clip Trench Points based on the Target Extent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37a851f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "target_extent_anim_file = parameters['target_extent_animation']\n",
    "target_extent_anim_gdf = gpd.read_file(target_extent_anim_file)\n",
    "target_extent_anim_bounds = target_extent_anim_gdf.bounds\n",
    "target_extent_anim = [target_extent_anim_bounds.loc[0]['minx'], target_extent_anim_bounds.loc[0]['maxx'],\n",
    "                      target_extent_anim_bounds.loc[0]['miny'], target_extent_anim_bounds.loc[0]['maxy']]\n",
    "\n",
    "if target_extent_anim[0] < -180:\n",
    "    target_extent_anim[0] = -180\n",
    "if target_extent_anim[1] > 180:\n",
    "    target_extent_anim[1] = 180\n",
    "if target_extent_anim[2] < -90:\n",
    "    target_extent_anim[2] = -90\n",
    "if target_extent_anim[3] > 90:\n",
    "    target_extent_anim[3] = 90\n",
    "\n",
    "features_target_extent_anim_files_lst = []\n",
    "\n",
    "for time in time_steps:\n",
    "    features_target_extent_anim_files_lst.append(f'{conv_dir}/{conv_prefix}_target_extent_animation_{time}.00.{conv_ext}')\n",
    "\n",
    "for features_target_extent_anim_file in features_target_extent_anim_files_lst:\n",
    "    if not os.path.isfile(features_target_extent_anim_file):\n",
    "        time = features_target_extent_anim_files_lst.index(features_target_extent_anim_file)\n",
    "        features_t = pd.read_csv(f'{conv_dir}/{conv_prefix}_{time}.00.{conv_ext}', index_col=False)\n",
    "        features_target_extent_anim_lst = []\n",
    "\n",
    "        for i in range(features_t.shape[0]):\n",
    "            x = features_t.iloc[i]['trench_lon']\n",
    "            y = features_t.iloc[i]['trench_lat']\n",
    "            p = Point((x, y))\n",
    "            if p.within(target_extent_anim_gdf.geometry[0]):\n",
    "                features_target_extent_anim_lst.append(features_t.iloc[i].values)\n",
    "\n",
    "        features_target_extent_anim = pd.DataFrame(np.row_stack(features_target_extent_anim_lst), columns=features_t.columns)\n",
    "        features_target_extent_anim.to_csv(features_target_extent_anim_file, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6aeeb548",
   "metadata": {},
   "outputs": [],
   "source": [
    "target_extent_map_file = parameters['target_extent_map']\n",
    "target_extent_map_gdf = gpd.read_file(target_extent_map_file)\n",
    "target_extent_map_bounds = target_extent_map_gdf.bounds\n",
    "target_extent_map = [target_extent_map_bounds.loc[0]['minx'], target_extent_map_bounds.loc[0]['maxx'],\n",
    "                     target_extent_map_bounds.loc[0]['miny'], target_extent_map_bounds.loc[0]['maxy']]\n",
    "\n",
    "if target_extent_map[0] < -180:\n",
    "    target_extent_map[0] = -180\n",
    "if target_extent_map[1] > 180:\n",
    "    target_extent_map[1] = 180\n",
    "if target_extent_map[2] < -90:\n",
    "    target_extent_map[2] = -90\n",
    "if target_extent_map[3] > 90:\n",
    "    target_extent_map[3] = 90\n",
    "\n",
    "features_target_extent_map_files_lst = []\n",
    "\n",
    "for time in time_steps:\n",
    "    features_target_extent_map_files_lst.append(f'{conv_dir}/{conv_prefix}_target_extent_map_{time}.00.{conv_ext}')\n",
    "\n",
    "for features_target_extent_map_file in features_target_extent_map_files_lst:\n",
    "    if not os.path.isfile(features_target_extent_map_file):\n",
    "        time = features_target_extent_map_files_lst.index(features_target_extent_map_file)\n",
    "        features_t = pd.read_csv(f'{conv_dir}/{conv_prefix}_{time}.00.{conv_ext}', index_col=False)\n",
    "        features_target_extent_map_lst = []\n",
    "\n",
    "        for i in range(features_t.shape[0]):\n",
    "            x = features_t.iloc[i]['trench_lon']\n",
    "            y = features_t.iloc[i]['trench_lat']\n",
    "            p = Point((x, y))\n",
    "            if p.within(target_extent_map_gdf.geometry[0]):\n",
    "                features_target_extent_map_lst.append(features_t.iloc[i].values)\n",
    "\n",
    "        features_target_extent_map = pd.DataFrame(np.row_stack(features_target_extent_map_lst), columns=features_t.columns)\n",
    "        features_target_extent_map.to_csv(features_target_extent_map_file, index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6da1411f",
   "metadata": {},
   "source": [
    "### Plot Kinematic Features based on the Target Extent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d761705",
   "metadata": {},
   "outputs": [],
   "source": [
    "def set_ax(ax, extent, interval_x, interval_y, font_size=None, stock_img=True, order=None):\n",
    "    if stock_img:\n",
    "        ax.stock_img()\n",
    "\n",
    "    ax.set_extent(extent)\n",
    "\n",
    "    gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,\n",
    "                      linewidth=1, color='gray', alpha=0.5, linestyle='--', \n",
    "                      xlocs=np.arange(-180, 180, interval_x), ylocs=np.arange(-90, 90, interval_y), zorder=order)\n",
    "    gl.xlabels_top = False\n",
    "    gl.ylabels_right = False\n",
    "    gl.xformatter = LONGITUDE_FORMATTER\n",
    "    gl.yformatter = LATITUDE_FORMATTER\n",
    "    gl.xlabel_style = {'color': 'gray', 'weight': 'bold', 'fontsize': font_size}\n",
    "    gl.ylabel_style = {'color': 'gray', 'weight': 'bold', 'fontsize': font_size}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e766e299",
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.LambertAzimuthalEqualArea(150, 0)\n",
    "\n",
    "@interact\n",
    "def show_map(time=time_steps, feature=selected_features_plot):\n",
    "    # call the PlotTopologies object\n",
    "    gplot = gplately.PlotTopologies(model, coastlines, continents, cob, time=time)\n",
    "    \n",
    "    if plate_motion_model == 'muller2016':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2016_AREPS_v1.17_AgeGrid-{time}.nc'\n",
    "    elif plate_motion_model == 'muller2019':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2019_Tectonics_v2.0_AgeGrid-{time}.nc'\n",
    "    \n",
    "    agegrid = gplately.grids.read_netcdf_grid(agegrid_file)\n",
    "    \n",
    "    features_t = pd.read_csv(f'{conv_dir}/{conv_prefix}_target_extent_map_{time}.00.{conv_ext}', index_col=False)\n",
    "\n",
    "    # dual colour bars\n",
    "    fig = plt.figure(figsize=(6, 8))\n",
    "    gs = GridSpec(2, 2, hspace=-0.75, wspace=0.1, height_ratios=[1, 0.01])\n",
    "    ax = fig.add_subplot(gs[0, :], projection=proj)\n",
    "    \n",
    "    set_ax(ax, target_extent_map, 10, 5, stock_img=False, order=8)\n",
    "\n",
    "    im = gplot.plot_grid(ax, agegrid.data, cmap=ccm.lapaz_r, vmin=0, vmax=153, alpha=0.5, zorder=1)\n",
    "\n",
    "    gplot.plot_continents(ax, edgecolor='gray', facecolor='tan', zorder=2)\n",
    "    gplot.plot_ridges(ax, color='red', alpha=0.5, zorder=3)\n",
    "    gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=False, regrid_shape=10, alpha=0.2, zorder=4)\n",
    "\n",
    "    sc = ax.scatter(features_t['trench_lon'], features_t['trench_lat'], 50, marker='.',\n",
    "                    c=features_t[feature], cmap=ccm.hawaii_r, transform=ccrs.PlateCarree(), zorder=5)\n",
    "\n",
    "    gplot.plot_trenches(ax, color='k', alpha=0.3, zorder=6)\n",
    "    gplot.plot_subduction_teeth(ax, spacing=0.03, color='k', alpha=0.3, zorder=7)\n",
    "        \n",
    "    cax1 = fig.add_subplot(gs[1, 0])\n",
    "    cax2 = fig.add_subplot(gs[1, 1])\n",
    "    \n",
    "    fig.colorbar(sc, cax=cax2, orientation='horizontal', label=feature, extend='both')\n",
    "    fig.colorbar(im, cax=cax1, orientation='horizontal', label='Seafloor age (Ma)', extend='max')\n",
    "    \n",
    "    # Define custom legend handles\n",
    "    custom_handles = [\n",
    "        Patch(facecolor='tan', edgecolor='gray', label='Continental Crust'),  # Custom handle for the filled polygon\n",
    "        Line2D([0], [0], color='red', lw=2, label='Mid-Ocean Ridge')  # Custom handle for the line (ridge)\n",
    "    ]\n",
    "\n",
    "    # Add the custom legend to the plot\n",
    "    legend = ax.legend(handles=custom_handles, loc='upper right', bbox_to_anchor=(1.4, 1), borderaxespad=0.)\n",
    "    \n",
    "#     plt.savefig(\n",
    "#     f'./figures/muller2019/features/conv_angle_deg_edited_.png',\n",
    "#     bbox_inches='tight',\n",
    "#     pad_inches=0.1,\n",
    "#     dpi=150\n",
    "#     )\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "724a596b",
   "metadata": {},
   "source": [
    "### Reconstruct Mineral Occurrences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70e7d703",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "min_occ_file = parameters['min_occ_file']\n",
    "coreg_input_dir = parameters['coreg_input_dir']\n",
    "coreg_input_files = parameters['coreg_input_files']\n",
    "min_occ_data_file = coreg_input_dir + coreg_input_files[0] # returns 'mineral_occurrences'\n",
    "\n",
    "if os.path.isfile(min_occ_data_file):\n",
    "    min_occ_data = pd.read_csv(min_occ_data_file, index_col=False)\n",
    "else:\n",
    "    # id, lon, lat, age, and plate id of mineral occurrences\n",
    "    min_occ_data = process_real_deposits(min_occ_file, start_time, end_time, time_step, plate_motion_model)\n",
    "    # save the attributes of mineral occurrences\n",
    "    min_occ_data.to_csv(min_occ_data_file, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a12bccf4",
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.LambertAzimuthalEqualArea(150, 0)\n",
    "\n",
    "@interact\n",
    "def show_map(time=time_steps, feature=selected_features_plot):\n",
    "    lons_lats_recon = []\n",
    "    \n",
    "    for min_occ in min_occ_data.iterrows():\n",
    "        if time == 0:\n",
    "            lons_lats_recon.append((min_occ[1]['lon'], min_occ[1]['lat'], min_occ[1]['weight']))\n",
    "        elif int(min_occ[1]['age']) < time:\n",
    "            lons_lats_recon.append((np.nan, np.nan, np.nan))\n",
    "        elif int(min_occ[1]['age']) == time:\n",
    "            lons_lats_recon.append((min_occ[1]['lon_recon'], min_occ[1]['lat_recon'], min_occ[1]['weight']))\n",
    "        else:\n",
    "            lat_lon_recon = get_recon_ccords([min_occ[1]['lon']],\n",
    "                                             [min_occ[1]['lat']],\n",
    "                                             plate_motion_model='muller2019',\n",
    "                                             time=time)[0]\n",
    "            lons_lats_recon.append(tuple((lat_lon_recon[1], lat_lon_recon[0], min_occ[1]['weight'])))\n",
    "\n",
    "    # call the PlotTopologies object\n",
    "    gplot = gplately.PlotTopologies(model, coastlines, continents, cob, time=time)\n",
    "    \n",
    "    if plate_motion_model == 'muller2016':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2016_AREPS_v1.17_AgeGrid-{time}.nc'\n",
    "    elif plate_motion_model == 'muller2019':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2019_Tectonics_v2.0_AgeGrid-{time}.nc'\n",
    "    \n",
    "    agegrid = gplately.grids.read_netcdf_grid(agegrid_file)\n",
    "    \n",
    "    features_t = pd.read_csv(f'{conv_dir}/{conv_prefix}_{time}.00.{conv_ext}', index_col=False)\n",
    "\n",
    "    # dual colour bars\n",
    "    fig = plt.figure(figsize=(6, 8))\n",
    "    gs = GridSpec(2, 2, hspace=-0.75, wspace=0.1, height_ratios=[1, 0.01])\n",
    "    ax = fig.add_subplot(gs[0, :], projection=proj)\n",
    "    \n",
    "    set_ax(ax, target_extent_map, 10, 5, stock_img=False, order=9)\n",
    "\n",
    "    im = gplot.plot_grid(ax, agegrid.data, cmap=ccm.lapaz_r, vmin=0, vmax=153, alpha=0.5, zorder=1)\n",
    "\n",
    "    gplot.plot_continents(ax, edgecolor='gray', facecolor='tan', zorder=2)\n",
    "    gplot.plot_ridges(ax, color='red', alpha=0.5, zorder=3)\n",
    "    gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=False, regrid_shape=10, alpha=0.2, zorder=4)\n",
    "\n",
    "    sc0 = ax.scatter(features_t['trench_lon'], features_t['trench_lat'], 50, marker='.',\n",
    "                    c=features_t[feature], cmap=ccm.hawaii_r, transform=ccrs.PlateCarree(), zorder=5)\n",
    "    \n",
    "    sc1 = ax.scatter(\n",
    "        [coords[0] for coords in lons_lats_recon],\n",
    "        [coords[1] for coords in lons_lats_recon],\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        marker='o',\n",
    "        facecolor='yellow',\n",
    "        edgecolor='black',\n",
    "        s=[x * 20 for x in [coords[2] for coords in lons_lats_recon]],\n",
    "        alpha=0.7,\n",
    "        zorder=6\n",
    "    )\n",
    "\n",
    "    gplot.plot_trenches(ax, color='k', alpha=0.3, zorder=7)\n",
    "    gplot.plot_subduction_teeth(ax, spacing=0.03, color='k', alpha=0.3, zorder=8)\n",
    "    \n",
    "    cax1 = fig.add_subplot(gs[1, 0])\n",
    "    cax2 = fig.add_subplot(gs[1, 1])\n",
    "    \n",
    "    fig.colorbar(sc0, cax=cax2, orientation='horizontal', label=feature, extend='both')\n",
    "    fig.colorbar(im, cax=cax1, orientation='horizontal', label='Seafloor age (Ma)', extend='max')\n",
    "    \n",
    "    # Define custom legend handles\n",
    "    custom_handles = [\n",
    "        Patch(facecolor='tan', edgecolor='gray', label='Continental Crust'),  # Custom handle for the filled polygon\n",
    "        Line2D([0], [0], color='red', lw=2, label='Mid-Ocean Ridge'),  # Custom handle for the line (ridge)\n",
    "        Line2D([0], [0], marker='o', markerfacecolor='yellow', markeredgecolor='black', markersize=5, linestyle='None', label='Mineral Occurrence')  # Custom handle for mineral occurrences\n",
    "    ]\n",
    "\n",
    "    # Add the custom legend to the plot\n",
    "    legend = ax.legend(handles=custom_handles, loc='upper right', bbox_to_anchor=(1.42, 1), borderaxespad=0.)\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c1cf702",
   "metadata": {},
   "source": [
    "### Create Buffer Zones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c57e7a72",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if plate_motion_model == 'muller2016':\n",
    "    rotation_files = parameters['rotation_files']\n",
    "    topology_files = parameters['topology_files']\n",
    "elif plate_motion_model == 'muller2019':\n",
    "    rotation_files = [os.path.join(dirpath, f) for (dirpath, dirnames, filenames) in os.walk(parameters['rotation_dir']) for f in filenames]\n",
    "    topology_files = [os.path.join(dirpath, f) for (dirpath, dirnames, filenames) in os.walk(parameters['topology_dir']) for f in filenames]\n",
    "\n",
    "buffer_zones_files_lst = []\n",
    "\n",
    "for time in time_steps:\n",
    "    buffer_zones_files_lst.append(f'{coreg_input_dir}buffer_zones/buffer_zone_{time}_Ma.shp')\n",
    "\n",
    "for buffer_zone_file in buffer_zones_files_lst:\n",
    "    if not os.path.isfile(buffer_zone_file):\n",
    "        index = buffer_zones_files_lst.index(buffer_zone_file)\n",
    "        \n",
    "        resolved_topologies = []\n",
    "        shared_boundary_sections = []\n",
    "        # use pygplates to resolve the topologies\n",
    "        pygplates.resolve_topologies(topology_files, rotation_files, resolved_topologies, time_steps[index], shared_boundary_sections)\n",
    "\n",
    "        # subduction zones\n",
    "        subduction_geoms = []\n",
    "        get_subduction_geometries(subduction_geoms, shared_boundary_sections)\n",
    "\n",
    "        _, buffer_zone = generate_buffer_zones(subduction_geoms, width=3)\n",
    "        buffer_zone.to_file(buffer_zone_file)\n",
    "        print(f'Buffer zones saved to {buffer_zone_file}')\n",
    "\n",
    "buffer_zones_lst = []\n",
    "buffer_zones_clipped_lst = []\n",
    "\n",
    "for buffer_zone_file in buffer_zones_files_lst:\n",
    "    buffer_zone = gpd.read_file(buffer_zone_file)\n",
    "    buffer_zone_clipped = buffer_zone.clip(target_extent_anim_gdf)\n",
    "    buffer_zones_lst.append(buffer_zone)\n",
    "    buffer_zones_clipped_lst.append(buffer_zone_clipped)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbe2289d",
   "metadata": {},
   "source": [
    "### Plot Buffer Zones in a Global Scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f545b43",
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.Mollweide(central_longitude=150)\n",
    "\n",
    "@interact\n",
    "def show_map(time=time_steps):\n",
    "    # call the PlotTopologies object\n",
    "    gplot = gplately.PlotTopologies(model, coastlines, continents, cob, time=time)\n",
    "    \n",
    "    if plate_motion_model == 'muller2016':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2016_AREPS_v1.17_AgeGrid-{time}.nc'\n",
    "    elif plate_motion_model == 'muller2019':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2019_Tectonics_v2.0_AgeGrid-{time}.nc'\n",
    "    \n",
    "    agegrid = gplately.grids.read_netcdf_grid(agegrid_file)\n",
    "    \n",
    "    # single colour bar\n",
    "    fig = plt.figure(figsize=(16, 12))\n",
    "    ax = plt.axes(projection=proj)\n",
    "    ax.set_global()\n",
    "\n",
    "    im = gplot.plot_grid(ax, agegrid.data, cmap=ccm.lapaz_r, vmin=0, vmax=230, alpha=0.5, zorder=1)\n",
    "\n",
    "    gplot.plot_continents(ax, edgecolor='none', facecolor='tan', zorder=2)\n",
    "    gplot.plot_ridges(ax, color='red', alpha=0.5, zorder=3)\n",
    "    gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=True, alpha=0.1, zorder=4)\n",
    "\n",
    "    buffer_zones_lst[time_steps.index(time)].plot(\n",
    "        ax=ax,\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        edgecolor='none',\n",
    "        facecolor='gray',\n",
    "        linewidth=1,\n",
    "        alpha=0.7,\n",
    "        zorder=5,\n",
    "    )\n",
    "    \n",
    "    gplot.plot_trenches(ax, color='k', zorder=6)\n",
    "    gplot.plot_subduction_teeth(ax, spacing=0.05, color='k', zorder=7)\n",
    "    \n",
    "    ax.gridlines(linestyle=':')\n",
    "    \n",
    "    fig.colorbar(im, orientation='horizontal', shrink=0.4, pad=0.05, label='Seafloor Age (Ma)', extend='max')\n",
    "    \n",
    "    # Define custom legend handles\n",
    "    custom_handles = [\n",
    "        Patch(facecolor='tan', edgecolor='none', label='Continental Crust'),  # Custom handle for the filled polygon\n",
    "        Patch(facecolor='gray', edgecolor='none', label='Target Areas in\\nBack-Arc Basins'),  # Custom handle for the buffer zones\n",
    "        Line2D([0], [0], color='red', lw=2, label='Mid-Ocean Ridge')  # Custom handle for the line (ridge)\n",
    "    ]\n",
    "\n",
    "    # Add the custom legend to the plot\n",
    "    ax.legend(handles=custom_handles, loc='lower left')\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "477e6148",
   "metadata": {},
   "source": [
    "### Plot Buffer Zones based on the Target Extent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e48a6a9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.LambertAzimuthalEqualArea(150, 0)\n",
    "\n",
    "@interact\n",
    "def show_map(time=time_steps):\n",
    "    # call the PlotTopologies object\n",
    "    gplot = gplately.PlotTopologies(model, coastlines, continents, cob, time=time)\n",
    "    \n",
    "    if plate_motion_model == 'muller2016':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2016_AREPS_v1.17_AgeGrid-{time}.nc'\n",
    "    elif plate_motion_model == 'muller2019':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2019_Tectonics_v2.0_AgeGrid-{time}.nc'\n",
    "    \n",
    "    agegrid = gplately.grids.read_netcdf_grid(agegrid_file)\n",
    "    \n",
    "    # single colour bar\n",
    "    fig = plt.figure(figsize=(6, 8))\n",
    "    ax = plt.axes(projection=proj)\n",
    "    \n",
    "    set_ax(ax, target_extent_map, 10, 5, stock_img=False, order=8)\n",
    "\n",
    "    im = gplot.plot_grid(ax, agegrid.data, cmap=ccm.lapaz_r, vmin=0, vmax=153, alpha=0.5, zorder=1)\n",
    "\n",
    "    gplot.plot_continents(ax, edgecolor='gray', facecolor='tan', zorder=2)\n",
    "    gplot.plot_ridges(ax, color='red', alpha=0.5, zorder=3)\n",
    "    gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=False, regrid_shape=10, alpha=0.2, zorder=4)\n",
    "    \n",
    "    buffer_zones_clipped_lst[time_steps.index(time)].plot(\n",
    "        ax=ax,\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        edgecolor='none',\n",
    "        facecolor='gray',\n",
    "        linewidth=1,\n",
    "        alpha=0.7,\n",
    "        zorder=5\n",
    "    )\n",
    "    \n",
    "    gplot.plot_trenches(ax, color='k', alpha=0.3, zorder=6)\n",
    "    gplot.plot_subduction_teeth(ax, spacing=0.03, color='k', alpha=0.3, zorder=7)\n",
    "    \n",
    "    fig.colorbar(im, orientation='horizontal', shrink=0.4, pad=0.05, label='Seafloor Age (Ma)', extend='max')\n",
    "    \n",
    "    # Define custom legend handles\n",
    "    custom_handles = [\n",
    "        Patch(facecolor='tan', edgecolor='gray', label='Continental Crust'),  # Custom handle for the filled polygon\n",
    "        Patch(facecolor='gray', edgecolor='none', label='Target Areas in\\nBack-Arc Basins'),  # Custom handle for the buffer zones\n",
    "        Line2D([0], [0], color='red', lw=2, label='Mid-Ocean Ridge')  # Custom handle for the line (ridge)\n",
    "    ]\n",
    "\n",
    "    # Add the custom legend to the plot\n",
    "    legend = ax.legend(handles=custom_handles, loc='upper right', bbox_to_anchor=(1.4, 1), borderaxespad=0.)\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c83d3e3",
   "metadata": {},
   "source": [
    "### Generate Random Samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3004219c",
   "metadata": {},
   "outputs": [],
   "source": [
    "random_data_file = coreg_input_dir + coreg_input_files[1]\n",
    "num_features = len(selected_features)\n",
    "\n",
    "if os.path.isfile(random_data_file):\n",
    "    random_data = pd.read_csv(random_data_file, index_col=False)\n",
    "    time_steps_random = random_data['age'].tolist()\n",
    "else:\n",
    "    time_steps_random, random_data = generate_random_samples(buffer_zones_clipped_lst,\n",
    "                                          start_time=start_time,\n",
    "                                          end_time=end_time,\n",
    "                                          time_step=time_step,\n",
    "                                          num_features=num_features,\n",
    "                                          num_features_factor=5,\n",
    "                                          rand_factor=20,\n",
    "                                          plate_motion_model='muller2019',\n",
    "                                          random_state=42\n",
    "                                         )\n",
    "    random_data.to_csv(random_data_file, index=False, float_format='%.4f')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fffa0e2",
   "metadata": {},
   "source": [
    "### Plot Random Samples based on the Target Extent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5915c649",
   "metadata": {},
   "outputs": [],
   "source": [
    "time_steps_random_sorted = sorted(time_steps_random)\n",
    "time_steps_random_sorted = [*set(time_steps_random_sorted)]\n",
    "\n",
    "proj = ccrs.LambertAzimuthalEqualArea(150, 0)\n",
    "\n",
    "@interact\n",
    "def show_map(time=time_steps_random_sorted):\n",
    "    # call the PlotTopologies object\n",
    "    gplot = gplately.PlotTopologies(model, coastlines, continents, cob, time=time)\n",
    "    \n",
    "    if plate_motion_model == 'muller2016':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2016_AREPS_v1.17_AgeGrid-{time}.nc'\n",
    "    elif plate_motion_model == 'muller2019':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2019_Tectonics_v2.0_AgeGrid-{time}.nc'\n",
    "    \n",
    "    agegrid = gplately.grids.read_netcdf_grid(agegrid_file)\n",
    "    \n",
    "    # single colour bar\n",
    "    fig = plt.figure(figsize=(6, 8))\n",
    "    ax = plt.axes(projection=proj)\n",
    "    \n",
    "    set_ax(ax, target_extent_map, 10, 5, stock_img=False, order=9)\n",
    "\n",
    "    im = gplot.plot_grid(ax, agegrid.data, cmap=ccm.lapaz_r, vmin=0, vmax=153, alpha=0.5, zorder=1)\n",
    "\n",
    "    gplot.plot_continents(ax, edgecolor='gray', facecolor='tan', zorder=2)\n",
    "    gplot.plot_ridges(ax, color='red', alpha=0.5, zorder=3)\n",
    "    gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=False, regrid_shape=10, alpha=0.2, zorder=4)\n",
    "    \n",
    "    buffer_zones_clipped_lst[time_steps.index(time)].plot(\n",
    "        ax=ax,\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        edgecolor='none',\n",
    "        facecolor='gray',\n",
    "        linewidth=1,\n",
    "        alpha=0.7,\n",
    "        zorder=5\n",
    "    )\n",
    "    \n",
    "    random_samples = random_data.loc[random_data['age'] == time]\n",
    "    ax.scatter(\n",
    "        random_samples['lon'],\n",
    "        random_samples['lat'],\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        marker='X',\n",
    "        edgecolor='black',\n",
    "        facecolor='cyan',\n",
    "        s=50,\n",
    "        zorder=8\n",
    "    )\n",
    "    \n",
    "    gplot.plot_trenches(ax, color='k', alpha=0.3, zorder=6)\n",
    "    gplot.plot_subduction_teeth(ax, spacing=0.03, color='k', alpha=0.3, zorder=7)\n",
    "    \n",
    "    fig.colorbar(im, orientation='horizontal', shrink=0.4, pad=0.05, label='Seafloor Age (Ma)', extend='max')\n",
    "    \n",
    "    # Define custom legend handles\n",
    "    custom_handles = [\n",
    "        Patch(facecolor='tan', edgecolor='gray', label='Continental Crust'),  # Custom handle for the filled polygon\n",
    "        Patch(facecolor='gray', edgecolor='none', label='Target Areas in\\nBack-Arc Basins'),  # Custom handle for the buffer zones\n",
    "        Line2D([0], [0], color='red', lw=2, label='Mid-Ocean Ridge')  # Custom handle for the line (ridge)\n",
    "    ]\n",
    "\n",
    "    # Add the custom legend to the plot\n",
    "    legend = ax.legend(handles=custom_handles, loc='upper right', bbox_to_anchor=(1.4, 1), borderaxespad=0.)\n",
    "        \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2be7a737",
   "metadata": {},
   "source": [
    "### Generate Target Points in Back-Arc Basins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5c20a13",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "target_points_coreg_in_files_lst = []\n",
    "mask_coords_files_lst = []\n",
    "\n",
    "for time in time_steps:\n",
    "    target_points_coreg_in_files_lst.append(coreg_input_dir + coreg_input_files[2] + f'_{time}_Ma.csv')\n",
    "    mask_coords_files_lst.append(coreg_input_dir + f'mask_{time}_Ma.csv')\n",
    "\n",
    "for target_points_file, mask_coords_file in zip(target_points_coreg_in_files_lst, mask_coords_files_lst):\n",
    "    if not(os.path.isfile(target_points_file) and os.path.isfile(mask_coords_file)):\n",
    "        index = target_points_coreg_in_files_lst.index(target_points_file)\n",
    "\n",
    "        # generate target points\n",
    "        target_points, mask_coords, nx, ny = generate_samples(buffer_zones_clipped_lst[index], 0.2, 0.2, # dist_x and dist_y\n",
    "                                                              time_steps[index], plate_motion_model='muller2019')\n",
    "        # save the attributes of target points\n",
    "        target_points.to_csv(target_points_file, index=False, float_format='%.4f')\n",
    "        # save the mask\n",
    "        mask_coords.to_csv(mask_coords_file, index=False, float_format='%.4f')\n",
    "        print(f'Target points saved to {target_points_file}')\n",
    "\n",
    "target_points_coreg_in_lst = []\n",
    "mask_coords_lst = []\n",
    "\n",
    "for target_points_file, mask_coords_file in zip(target_points_coreg_in_files_lst, mask_coords_files_lst):\n",
    "    target_points_coreg_in_lst.append(pd.read_csv(target_points_file, index_col=False))\n",
    "    mask_coords_lst.append(pd.read_csv(mask_coords_file, index_col=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c88c5884",
   "metadata": {},
   "source": [
    "### Plot Target Points based on the Target Extent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec5e9718",
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.LambertAzimuthalEqualArea(150, 0)\n",
    "\n",
    "@interact\n",
    "def show_map(time=time_steps):\n",
    "    # call the PlotTopologies object\n",
    "    gplot = gplately.PlotTopologies(model, coastlines, continents, cob, time=time)\n",
    "    \n",
    "    if plate_motion_model == 'muller2016':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2016_AREPS_v1.17_AgeGrid-{time}.nc'\n",
    "    elif plate_motion_model == 'muller2019':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2019_Tectonics_v2.0_AgeGrid-{time}.nc'\n",
    "\n",
    "    agegrid = gplately.grids.read_netcdf_grid(agegrid_file)\n",
    "    \n",
    "    # single colour bar\n",
    "    fig = plt.figure(figsize=(6, 8))\n",
    "    ax = plt.axes(projection=proj)\n",
    "    \n",
    "    set_ax(ax, target_extent_map, 10, 5, stock_img=False, order=9)\n",
    "\n",
    "    im = gplot.plot_grid(ax, agegrid.data, cmap=ccm.lapaz_r, vmin=0, vmax=153, alpha=0.5, zorder=1)\n",
    "\n",
    "    gplot.plot_continents(ax, edgecolor='gray', facecolor='tan', zorder=2)\n",
    "    gplot.plot_ridges(ax, color='red', alpha=0.5, zorder=3)\n",
    "    gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=False, regrid_shape=10, alpha=0.2, zorder=4)\n",
    "    \n",
    "    buffer_zones_clipped_lst[time_steps.index(time)].plot(\n",
    "        ax=ax,\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        edgecolor='none',\n",
    "        facecolor='gray',\n",
    "        linewidth=1,\n",
    "        alpha=0.7,\n",
    "        zorder=5\n",
    "    )\n",
    "    \n",
    "    ax.scatter(\n",
    "        target_points_coreg_in_lst[time_steps.index(time)]['lon'],\n",
    "        target_points_coreg_in_lst[time_steps.index(time)]['lat'],\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        marker='.',\n",
    "        c='red',\n",
    "        s=1,\n",
    "        zorder=8\n",
    "    )\n",
    "        \n",
    "    gplot.plot_trenches(ax, color='k', alpha=0.3, zorder=6)\n",
    "    gplot.plot_subduction_teeth(ax, spacing=0.03, color='k', alpha=0.3, zorder=7)\n",
    "    \n",
    "    fig.colorbar(im, orientation='horizontal', shrink=0.4, pad=0.05, label='Seafloor Age (Ma)', extend='max')\n",
    "    \n",
    "    # Define custom legend handles\n",
    "    custom_handles = [\n",
    "        Patch(facecolor='tan', edgecolor='gray', label='Continental Crust'),  # Custom handle for the filled polygon\n",
    "        Patch(facecolor='gray', edgecolor='none', label='Target Areas in\\nBack-Arc Basins'),  # Custom handle for the buffer zones\n",
    "        Line2D([0], [0], color='red', lw=2, label='Mid-Ocean Ridge'),  # Custom handle for the line (ridge)\n",
    "        Line2D([0], [0], marker='.', markerfacecolor='red', markeredgecolor='none', markersize=10, linestyle='None', label='Target Points')\n",
    "    ]\n",
    "\n",
    "    # Add the custom legend to the plot\n",
    "    legend = ax.legend(handles=custom_handles, loc='upper right', bbox_to_anchor=(1.4, 1), borderaxespad=0.)\n",
    "            \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6df7421",
   "metadata": {},
   "source": [
    "### Generate Target Points Inside New Guinea Island"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e1aa703",
   "metadata": {},
   "outputs": [],
   "source": [
    "target_polygon = parameters['target_polygon']\n",
    "target_polygon_gdf = gpd.read_file(target_polygon)\n",
    "\n",
    "target_points_ng_0_file = coreg_input_dir + coreg_input_files[3] + '_0_Ma.csv'\n",
    "mask_coords_ng_0_file = coreg_input_dir + 'mask_ng_0_Ma.csv'\n",
    "\n",
    "if os.path.isfile(target_points_ng_0_file) and os.path.isfile(mask_coords_ng_0_file):\n",
    "    target_points_ng_0 = pd.read_csv(target_points_ng_0_file, index_col=False)\n",
    "    mask_coords_ng_0 = pd.read_csv(mask_coords_ng_0_file, index_col=False)\n",
    "else:\n",
    "    target_points_ng_0, mask_coords_ng_0, nx_ng_0, ny_ng_0 = generate_samples(target_polygon_gdf, 0.2, 0.2, # dist_x and dist_y\n",
    "                                                                              time=0, plate_motion_model='muller2019')\n",
    "    target_points_ng_0.to_csv(target_points_ng_0_file, index=False, float_format='%.4f')\n",
    "    mask_coords_ng_0.to_csv(mask_coords_ng_0_file, index=False, float_format='%.4f')\n",
    "\n",
    "proj = ccrs.LambertAzimuthalEqualArea(150, 0)\n",
    "\n",
    "# call the PlotTopologies object\n",
    "gplot = gplately.PlotTopologies(model, coastlines, continents, cob, time=0)\n",
    "\n",
    "if plate_motion_model == 'muller2016':\n",
    "    agegrid_file = agegrid_dir + f'Muller_etal_2016_AREPS_v1.17_AgeGrid-0.nc'\n",
    "elif plate_motion_model == 'muller2019':\n",
    "    agegrid_file = agegrid_dir + f'Muller_etal_2019_Tectonics_v2.0_AgeGrid-0.nc'\n",
    "\n",
    "agegrid = gplately.grids.read_netcdf_grid(agegrid_file)\n",
    "\n",
    "ng_extent_bounds = target_polygon_gdf.bounds\n",
    "ng_extent = [ng_extent_bounds.loc[0]['minx'], ng_extent_bounds.loc[0]['maxx'],\n",
    "             ng_extent_bounds.loc[0]['miny'], ng_extent_bounds.loc[0]['maxy']]\n",
    "\n",
    "# single colour bar\n",
    "fig = plt.figure(figsize=(6, 8))\n",
    "ax = plt.axes(projection=proj)\n",
    "\n",
    "set_ax(ax, ng_extent, 10, 5, stock_img=False, order=9)\n",
    "\n",
    "im = gplot.plot_grid(ax, agegrid.data, cmap=ccm.lapaz_r, vmin=0, vmax=153, alpha=0.5, zorder=1)\n",
    "\n",
    "gplot.plot_continents(ax, edgecolor='gray', facecolor='tan', zorder=2)\n",
    "gplot.plot_ridges(ax, color='red', alpha=0.5, zorder=3)\n",
    "gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=False, regrid_shape=10, alpha=0.2, zorder=4)\n",
    "\n",
    "ax.scatter(\n",
    "    target_points_ng_0['lon'],\n",
    "    target_points_ng_0['lat'],\n",
    "    transform=ccrs.PlateCarree(),\n",
    "    marker='.',\n",
    "    c='green',\n",
    "    s=1,\n",
    "    zorder=8\n",
    ")\n",
    "\n",
    "gplot.plot_trenches(ax, color='k', alpha=0.3, zorder=6)\n",
    "gplot.plot_subduction_teeth(ax, spacing=0.03, color='k', alpha=0.3, zorder=7)\n",
    "\n",
    "fig.colorbar(im, orientation='horizontal', shrink=0.4, pad=0.05, label='Seafloor Age (Ma)', extend='max')\n",
    "\n",
    "# Define custom legend handles\n",
    "custom_handles = [\n",
    "    Patch(facecolor='tan', edgecolor='gray', label='Continental Crust'),  # Custom handle for the filled polygon\n",
    "    Line2D([0], [0], color='red', lw=2, label='Mid-Ocean Ridge'),  # Custom handle for the line (ridge)\n",
    "    Line2D([0], [0], marker='.', markerfacecolor='green', markeredgecolor='none', markersize=10, linestyle='None', label='Target Points')\n",
    "    \n",
    "]\n",
    "\n",
    "# Add the custom legend to the plot\n",
    "legend = ax.legend(handles=custom_handles, loc='upper right', bbox_to_anchor=(1.4, 1), borderaxespad=0.)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f615218-db0e-4272-a4cd-a060a9580346",
   "metadata": {},
   "outputs": [],
   "source": [
    "target_points_ng_coreg_in_files_lst = []\n",
    "\n",
    "time_steps_ng = time_steps.copy()\n",
    "time_steps_ng.remove(0)\n",
    "\n",
    "for time in time_steps_ng:\n",
    "    target_points_ng_coreg_in_files_lst.append(coreg_input_dir + coreg_input_files[3] + f'_{time}_Ma.csv')\n",
    "\n",
    "for target_points_ng_file in target_points_ng_coreg_in_files_lst:\n",
    "    if not os.path.isfile(target_points_ng_file):\n",
    "        index = target_points_ng_coreg_in_files_lst.index(target_points_ng_file)\n",
    "        # generate target points\n",
    "        target_points_ng = generate_samples_polygon(target_points_ng_0_file, time_steps_ng[index], plate_motion_model='muller2019')\n",
    "        # save the attributes of target points\n",
    "        target_points_ng.to_csv(target_points_ng_file, index=False, float_format='%.4f')\n",
    "        print(f'Target points saved to {target_points_ng_file}')\n",
    "\n",
    "target_points_ng_coreg_in_lst = []\n",
    "\n",
    "for target_points_ng_file in target_points_ng_coreg_in_files_lst:\n",
    "    target_points_ng_coreg_in_lst.append(pd.read_csv(target_points_ng_file, index_col=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f96b67a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.LambertAzimuthalEqualArea(150, 0)\n",
    "\n",
    "@interact\n",
    "def show_map(time=time_steps_ng):\n",
    "    # call the PlotTopologies object\n",
    "    gplot = gplately.PlotTopologies(model, coastlines, continents, cob, time=time)\n",
    "    \n",
    "    if plate_motion_model == 'muller2016':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2016_AREPS_v1.17_AgeGrid-{time}.nc'\n",
    "    elif plate_motion_model == 'muller2019':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2019_Tectonics_v2.0_AgeGrid-{time}.nc'\n",
    "\n",
    "    agegrid = gplately.grids.read_netcdf_grid(agegrid_file)\n",
    "    \n",
    "    # single colour bar\n",
    "    fig = plt.figure(figsize=(6, 8))\n",
    "    ax = plt.axes(projection=proj)\n",
    "    \n",
    "    set_ax(ax, target_extent_map, 10, 5, stock_img=False, order=9)\n",
    "\n",
    "    im = gplot.plot_grid(ax, agegrid.data, cmap=ccm.lapaz_r, vmin=0, vmax=153, alpha=0.5, zorder=1)\n",
    "\n",
    "    gplot.plot_continents(ax, edgecolor='gray', facecolor='tan', zorder=2)\n",
    "    gplot.plot_ridges(ax, color='red', alpha=0.5, zorder=3)\n",
    "    gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=False, regrid_shape=10, alpha=0.2, zorder=4)\n",
    "    \n",
    "    ax.scatter(\n",
    "        target_points_ng_coreg_in_lst[time_steps_ng.index(time)]['lon'],\n",
    "        target_points_ng_coreg_in_lst[time_steps_ng.index(time)]['lat'],\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        marker='.',\n",
    "        c='green',\n",
    "        s=0.1,\n",
    "        zorder=8\n",
    "    )\n",
    "        \n",
    "    gplot.plot_trenches(ax, color='k', alpha=0.3, zorder=6)\n",
    "    gplot.plot_subduction_teeth(ax, spacing=0.03, color='k', alpha=0.3, zorder=7)\n",
    "    \n",
    "    fig.colorbar(im, orientation='horizontal', shrink=0.4, pad=0.05, label='Seafloor Age (Ma)', extend='max')\n",
    "    \n",
    "    # Define custom legend handles\n",
    "    custom_handles = [\n",
    "        Patch(facecolor='tan', edgecolor='gray', label='Continental Crust'),  # Custom handle for the filled polygon\n",
    "        Line2D([0], [0], color='red', lw=2, label='Mid-Ocean Ridge'),  # Custom handle for the line (ridge)\n",
    "        Line2D([0], [0], marker='.', markerfacecolor='green', markeredgecolor='none', markersize=10, linestyle='None', label='Target Points')\n",
    "    ]\n",
    "\n",
    "    # Add the custom legend to the plot\n",
    "    legend = ax.legend(handles=custom_handles, loc='upper right', bbox_to_anchor=(1.4, 1), borderaxespad=0.)\n",
    "            \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "783efad6",
   "metadata": {},
   "source": [
    "### Coregistration and Data Wrangling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01944d11",
   "metadata": {},
   "outputs": [],
   "source": [
    "coreg_output_dir = parameters['coreg_output_dir']\n",
    "positive_data_file = coreg_output_dir + coreg_input_files[0]\n",
    "unlabelled_data_file = coreg_output_dir + coreg_input_files[1]\n",
    "target_points_coreg_out_files_lst = []\n",
    "target_points_ng_coreg_out_files_lst = []\n",
    "\n",
    "for time in time_steps:\n",
    "    target_points_coreg_out_files_lst.append(coreg_output_dir + coreg_input_files[2] + f'_{time}_Ma.csv')\n",
    "    \n",
    "for time in time_steps:\n",
    "    target_points_ng_coreg_out_files_lst.append(coreg_output_dir + coreg_input_files[3] + f'_{time}_Ma.csv')\n",
    "\n",
    "coregistration(\n",
    "    coreg_input_dir,\n",
    "    coreg_output_dir,\n",
    "    coreg_input_files,\n",
    "    conv_dir,\n",
    "    conv_prefix,\n",
    "    conv_ext,\n",
    "    time_steps=time_steps,\n",
    "    search_radius=3\n",
    ")\n",
    "\n",
    "positive_data = pd.read_csv(positive_data_file, index_col=False)\n",
    "unlabelled_data = pd.read_csv(unlabelled_data_file, index_col=False)\n",
    "\n",
    "target_points_coreg_out_lst = []\n",
    "for file_name in target_points_coreg_out_files_lst:\n",
    "    target_points_coreg_out_lst.append(pd.read_csv(file_name, index_col=False))\n",
    "\n",
    "target_points_ng_coreg_out_lst = []\n",
    "for file_name in target_points_ng_coreg_out_files_lst:\n",
    "    target_points_ng_coreg_out_lst.append(pd.read_csv(file_name, index_col=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9216f5b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ml_input_dir = parameters['ml_input_dir']\n",
    "\n",
    "positive_data['label'] = 1\n",
    "unlabelled_data['label'] = 0\n",
    "\n",
    "positive_features = positive_data[selected_features]\n",
    "unlabelled_features = unlabelled_data[selected_features]\n",
    "features_all = pd.concat([positive_features, unlabelled_features]).reset_index(drop=True)\n",
    "\n",
    "# save correlation csv file\n",
    "corr_file = ml_input_dir + 'correlation.csv'\n",
    "\n",
    "if os.path.isfile(corr_file):\n",
    "    corr = pd.read_csv(corr_file, index_col=0)\n",
    "else:\n",
    "    corr = features_all.corr(method='spearman').round(3)\n",
    "    corr.to_csv(corr_file, index=True)\n",
    "\n",
    "corr.style.background_gradient(cmap='coolwarm', axis=None).format('{:.3}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22b120db",
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_features_names = parameters['selected_features_names_nounit_01']\n",
    "corr.columns = selected_features_names\n",
    "corr.index = selected_features_names\n",
    "\n",
    "f = plt.figure(figsize=(20, 15))\n",
    "plt.matshow(corr, fignum=f.number, cmap='coolwarm', vmin=-1, vmax=1)\n",
    "plt.xticks(range(corr.select_dtypes(['number']).shape[1]), corr.select_dtypes(['number']).columns, fontsize=14, rotation=90)\n",
    "plt.yticks(range(corr.select_dtypes(['number']).shape[1]), corr.select_dtypes(['number']).columns, fontsize=14)\n",
    "cb = plt.colorbar(aspect=50)\n",
    "cb.ax.tick_params(labelsize=14)\n",
    "\n",
    "# plt.savefig(\n",
    "#     f'./figures/muller2019/correlation.png',\n",
    "#     bbox_inches='tight',\n",
    "#     pad_inches=0.1,\n",
    "#     dpi=150\n",
    "#     )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7961141e-b6c4-4131-8e2a-a2c633a4c8e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def analyze_correlations(file_path, threshold=0.7):\n",
    "    # Read the correlation matrix\n",
    "    corr_matrix = pd.read_csv(file_path, index_col=0)\n",
    "    \n",
    "    # Dictionary to store correlations\n",
    "    correlations = {}\n",
    "    \n",
    "    for column in corr_matrix.columns:\n",
    "        positive_corr = []\n",
    "        negative_corr = []\n",
    "        feature = corr_matrix[column]\n",
    "        \n",
    "        for i in range(feature.shape[0]):\n",
    "            if abs(feature[i]) >= threshold and feature.index[i] != column:\n",
    "                if feature[i] > 0:\n",
    "                    positive_corr.append((feature.index[i], feature[i]))\n",
    "                else:\n",
    "                    negative_corr.append((feature.index[i], feature[i]))\n",
    "        \n",
    "        if positive_corr or negative_corr:\n",
    "            correlations[column] = {\n",
    "                'positive': sorted(positive_corr, key=lambda x: x[1], reverse=True),\n",
    "                'negative': sorted(negative_corr, key=lambda x: x[1])\n",
    "            }\n",
    "    \n",
    "    return correlations\n",
    "\n",
    "def generate_report(correlations, threshold):\n",
    "    print(f\"Correlation Analysis Report (Threshold: {threshold})\")\n",
    "    print(\"=\" * 50)\n",
    "    \n",
    "    for feature, corr in correlations.items():\n",
    "        print(f\"\\nFeature: {feature}\")\n",
    "        print(\"-\" * 30)\n",
    "        \n",
    "        if corr['positive']:\n",
    "            print(\"Positive Correlations:\")\n",
    "            for c, value in corr['positive']:\n",
    "                print(f\"  {c}: {value:.3f}\")\n",
    "        \n",
    "        if corr['negative']:\n",
    "            print(\"Negative Correlations:\")\n",
    "            for c, value in corr['negative']:\n",
    "                print(f\"  {c}: {value:.3f}\")\n",
    "    \n",
    "    print(\"\\nTotal features with strong correlations:\", len(correlations))\n",
    "\n",
    "# Main execution\n",
    "file_path = ml_input_dir + 'correlation.csv'\n",
    "threshold = 0.7\n",
    "\n",
    "correlations = analyze_correlations(file_path, threshold)\n",
    "generate_report(correlations, threshold)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff7138da",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert correlations to distances\n",
    "distance_matrix = 1 - np.abs(corr)\n",
    "\n",
    "# Perform hierarchical clustering\n",
    "linkage_matrix = hierarchy.linkage(\n",
    "    hierarchy.distance.squareform(distance_matrix),\n",
    "    method='complete'\n",
    ")\n",
    "\n",
    "# Create figure\n",
    "plt.figure(figsize=(12, 8))\n",
    "\n",
    "# Create dendrogram\n",
    "dendrogram = hierarchy.dendrogram(\n",
    "    linkage_matrix,\n",
    "    labels=corr.columns,\n",
    "    orientation='right',\n",
    "    leaf_font_size=12,\n",
    "    leaf_rotation=0\n",
    ")\n",
    "\n",
    "# Get the axes object\n",
    "ax = plt.gca()\n",
    "\n",
    "# Get y-axis tick labels\n",
    "labels = ax.get_yticklabels()\n",
    "\n",
    "# Get leaf colors from dendrogram\n",
    "# The dendrogram returns a dictionary containing the color information\n",
    "leaf_colors = {}\n",
    "for i, leaf_color in enumerate(dendrogram['leaves_color_list']):\n",
    "    leaf_colors[dendrogram['ivl'][i]] = leaf_color\n",
    "\n",
    "# Set each label's color to match its corresponding line\n",
    "for label in labels:\n",
    "    label.set_color(leaf_colors[label.get_text()])\n",
    "\n",
    "# Customize the plot\n",
    "plt.title('Hierarchical Clustering Dendrogram (Spearman Correlation)', pad=20)\n",
    "plt.xlabel('|Spearman Correlation|')\n",
    "\n",
    "# Get current x-axis limits\n",
    "current_xlim = plt.xlim()\n",
    "\n",
    "# Create new tick positions (keep the same scale but show different labels)\n",
    "tick_positions = np.linspace(0, current_xlim[1], 6)\n",
    "tick_labels = [f'{x:.1f}' for x in np.linspace(1.0, 0.0, 6)]\n",
    "\n",
    "# Set new ticks and labels\n",
    "plt.xticks(tick_positions, tick_labels)\n",
    "\n",
    "# Add vertical line at x=0.7\n",
    "# Need to convert from our display scale (1.0 to 0.0) to the actual scale (0 to max)\n",
    "line_position = current_xlim[1] * (1 - 0.7)  # Convert from display scale to actual scale\n",
    "plt.axvline(x=line_position, color='r', linestyle='--', alpha=0.5)\n",
    "plt.annotate(f'0.7', xy=(0.315, -0.06), xycoords='data',\n",
    "             xytext=(0, -20), textcoords='offset points',\n",
    "             ha='center', va='top', color='r',\n",
    "             arrowprops=dict(arrowstyle='->', color='r'))\n",
    "\n",
    "# Adjust layout\n",
    "plt.tight_layout()\n",
    "\n",
    "# plt.savefig(\n",
    "#     f'./figures/muller2019/dendrogram.png',\n",
    "#     bbox_inches='tight',\n",
    "#     pad_inches=0.1,\n",
    "#     dpi=150\n",
    "# )\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e3776bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xy_train_file = ml_input_dir + 'Xy_train.csv'\n",
    "Xy_pos_test_file = ml_input_dir + 'Xy_test.csv'\n",
    "\n",
    "if os.path.isfile(Xy_train_file) and os.path.isfile(Xy_pos_test_file):\n",
    "    Xy_train = pd.read_csv(Xy_train_file, index_col=False)\n",
    "    features_list = Xy_train.columns.tolist()\n",
    "    features_list = [e for e in features_list if e not in ('label', 'weight')]\n",
    "    Xy_pos_test = pd.read_csv(Xy_pos_test_file, index_col=False)\n",
    "    print('Training data file already exists!')\n",
    "else:\n",
    "    positive_labels = positive_data[positive_data.columns[-1]]\n",
    "    unlabelled_labels = unlabelled_data[unlabelled_data.columns[-1]]\n",
    "\n",
    "    positive_weights = positive_data['weight']\n",
    "    unlabelled_weights = unlabelled_data['weight']\n",
    "    \n",
    "    labels = pd.concat([positive_labels, unlabelled_labels]).reset_index(drop=True)\n",
    "    weights = pd.concat([positive_weights, unlabelled_weights]).reset_index(drop=True)\n",
    "    \n",
    "    features_labels_all = pd.concat([features_all, weights, labels], axis=1).reset_index(drop=True)\n",
    "\n",
    "    # drop highly correlated features\n",
    "    # create a correlation matrix\n",
    "    corr_matrix = features_all.corr(method='spearman').abs()\n",
    "    # select the upper triangle of the correlation matrix\n",
    "    corr_upper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(bool))\n",
    "    # find features with the correlation greater than 0.7\n",
    "    corr_drop = [column for column in corr_upper.columns if any(corr_upper[column] > 0.7)]\n",
    "    print('List of the features removed due to high correlation with other features:', corr_drop)\n",
    "    # drop features\n",
    "    features = features_all.drop(corr_drop, axis=1)\n",
    "    features_list = features.columns.tolist()\n",
    "    features_labels = pd.concat([features, weights, labels], axis=1).reset_index(drop=True)\n",
    "    features_labels_list = features_labels.columns.tolist()\n",
    "\n",
    "    positive_data = features_labels[features_labels['label']==1]\n",
    "    unlabelled_data = features_labels[features_labels['label']==0]\n",
    "\n",
    "    positive_features = positive_data[positive_data.columns[:-1]]\n",
    "    unlabelled_features = unlabelled_data[unlabelled_data.columns[:-1]]\n",
    "\n",
    "    positive_labels = positive_data[positive_data.columns[-1]]\n",
    "    unlabelled_labels = unlabelled_data[unlabelled_data.columns[-1]]\n",
    "\n",
    "    # split positive samples into training and test datasets\n",
    "    X_pos_train, X_pos_test, y_pos_train, y_pos_test = train_test_split(positive_features, positive_labels, train_size=0.8, random_state=42)\n",
    "    X_train = np.vstack((X_pos_train, unlabelled_features))\n",
    "    y_train = np.vstack((y_pos_train.values.reshape(-1, 1), unlabelled_labels.values.reshape(-1, 1)))\n",
    "    \n",
    "    Xy_train_original = np.hstack((X_train, y_train))\n",
    "    Xy_train_original = pd.DataFrame(Xy_train_original, columns=features_labels_list)\n",
    "    Xy_train_original.to_csv(ml_input_dir + 'Xy_train_original.csv', index=False)\n",
    "\n",
    "    Xy_train_features = Xy_train_original[features_list]\n",
    "    Xy_train_labels = Xy_train_original['label']\n",
    "    Xy_train_weights = Xy_train_original['weight']\n",
    "\n",
    "    st_scaler = StandardScaler()\n",
    "    X_train = st_scaler.fit_transform(Xy_train_features)\n",
    "    Xy_train = np.hstack((X_train, Xy_train_weights.values.reshape(-1, 1), Xy_train_labels.values.reshape(-1, 1)))\n",
    "    Xy_train = pd.DataFrame(Xy_train, columns=features_labels_list)\n",
    "\n",
    "    Xy_pos_test_original = np.hstack((X_pos_test, y_pos_test.values.reshape(-1, 1)))\n",
    "    Xy_pos_test_original = pd.DataFrame(Xy_pos_test_original, columns=features_labels_list)\n",
    "    Xy_pos_test_original.to_csv(ml_input_dir + 'Xy_pos_test_original.csv', index=False)\n",
    "    \n",
    "    Xy_pos_test_features = Xy_pos_test_original[features_list]\n",
    "    Xy_pos_test_labels = Xy_pos_test_original['label']\n",
    "    Xy_pos_test_weights = Xy_pos_test_original['weight']\n",
    "    \n",
    "    X_pos_test = st_scaler.transform(Xy_pos_test_features)\n",
    "    Xy_pos_test = np.hstack((X_pos_test, Xy_pos_test_weights.values.reshape(-1, 1), Xy_pos_test_labels.values.reshape(-1, 1)))\n",
    "    Xy_pos_test = pd.DataFrame(Xy_pos_test, columns=features_labels_list)\n",
    "\n",
    "    Xy_train.to_csv(ml_input_dir + 'Xy_train.csv', index=False)\n",
    "    Xy_pos_test.to_csv(ml_input_dir + 'Xy_pos_test.csv', index=False)\n",
    "\n",
    "    # save the standard scaler model\n",
    "    with open(ml_input_dir + 'st_scaler.pkl', 'wb') as f:\n",
    "        pickle.dump(st_scaler, f)\n",
    "\n",
    "    print('\\nNumber of features reduced from', num_features, 'to', len(features_list))\n",
    "    print('Number of positive samples:', positive_features.shape[0])\n",
    "\n",
    "    print('\\nNumber of training samples:', X_train.shape[0]), print('Number of training labels:', y_train.shape[0])\n",
    "    print('Number of positive testing samples:', X_pos_test.shape[0]), print('Number of testing labels:', y_pos_test.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a834e425",
   "metadata": {},
   "outputs": [],
   "source": [
    "target_points_ml_in_files_lst = []\n",
    "for time in time_steps:\n",
    "    target_points_ml_in_files_lst.append(ml_input_dir + coreg_input_files[2] + f'_{time}_Ma.csv')\n",
    "    \n",
    "target_points_ml_in_lst = []\n",
    "for target_points_ml_in_file in target_points_ml_in_files_lst:\n",
    "    if os.path.isfile(target_points_ml_in_file):\n",
    "        target_points_ml_in_lst.append(pd.read_csv(target_points_ml_in_file, index_col=False))\n",
    "    if not os.path.isfile(target_points_ml_in_file):\n",
    "        index = target_points_ml_in_files_lst.index(target_points_ml_in_file)\n",
    "        target_points_ml_in = target_points_coreg_out_lst[index][selected_features]\n",
    "        target_points_ml_in = target_points_ml_in[target_points_ml_in.columns.intersection(features_list)]\n",
    "        \n",
    "        try:\n",
    "            target_points_ml_in = st_scaler.transform(target_points_ml_in)\n",
    "        except:\n",
    "            # load the model\n",
    "            with open(ml_input_dir + 'st_scaler.pkl', 'rb') as f:\n",
    "                st_scaler = pickle.load(f)\n",
    "            target_points_ml_in = st_scaler.transform(target_points_ml_in)\n",
    "            \n",
    "        target_points_ml_in = pd.DataFrame(target_points_ml_in, columns=features_list)\n",
    "        target_points_ml_in_lst.append(target_points_ml_in)\n",
    "        target_points_ml_in.to_csv(target_points_ml_in_file, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb2c6512",
   "metadata": {},
   "outputs": [],
   "source": [
    "target_points_ng_ml_in_files_lst = []\n",
    "for time in time_steps:\n",
    "    target_points_ng_ml_in_files_lst.append(ml_input_dir + coreg_input_files[3] + f'_{time}_Ma.csv')\n",
    "    \n",
    "target_points_ng_ml_in_lst = []\n",
    "for target_points_ng_ml_in_file in target_points_ng_ml_in_files_lst:\n",
    "    if os.path.isfile(target_points_ng_ml_in_file):\n",
    "        target_points_ng_ml_in_lst.append(pd.read_csv(target_points_ng_ml_in_file, index_col=False))\n",
    "    if not os.path.isfile(target_points_ng_ml_in_file):\n",
    "        index = target_points_ng_ml_in_files_lst.index(target_points_ng_ml_in_file)\n",
    "        target_points_ng_ml_in = target_points_ng_coreg_out_lst[index][selected_features]\n",
    "        target_points_ng_ml_in = target_points_ng_ml_in[target_points_ng_ml_in.columns.intersection(features_list)]\n",
    "        \n",
    "        try:\n",
    "            target_points_ng_ml_in = st_scaler.transform(target_points_ng_ml_in)\n",
    "        except:\n",
    "            # load the model\n",
    "            with open(ml_input_dir + 'st_scaler.pkl', 'rb') as f:\n",
    "                st_scaler = pickle.load(f)\n",
    "            target_points_ng_ml_in = st_scaler.transform(target_points_ng_ml_in)\n",
    "            \n",
    "        target_points_ng_ml_in = pd.DataFrame(target_points_ng_ml_in, columns=features_list)\n",
    "        target_points_ng_ml_in_lst.append(target_points_ng_ml_in)\n",
    "        target_points_ng_ml_in.to_csv(target_points_ng_ml_in_file, index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9925cbed",
   "metadata": {},
   "source": [
    "### Machine Learning\n",
    "#### Positive and Unlabelled Bagging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09d1711e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ml_output_dir = parameters['ml_output_dir']\n",
    "model_pub_file = ml_output_dir + 'model_pub.pkl'\n",
    "\n",
    "if os.path.isfile(model_pub_file):\n",
    "    print('A PUB model file exists. Attempting to load...')\n",
    "    try:\n",
    "        with open(model_pub_file, 'rb') as f:\n",
    "            model_pub = pickle.load(f)\n",
    "        print('Model loaded successfully.')\n",
    "    except (ModuleNotFoundError, ImportError) as e:\n",
    "        print(f\"Error loading the model: {e}\")\n",
    "        print(\"This may be due to a version mismatch. Proceeding to train a new model.\")\n",
    "        model_pub = None\n",
    "else:\n",
    "    print('No existing model found. Proceeding to train a new model.')\n",
    "    model_pub = None\n",
    "\n",
    "if model_pub is None:\n",
    "    # Random Forest model structure\n",
    "    rf_pub = RandomForestClassifier(n_jobs=-1, random_state=42)\n",
    "    pub = BaggingPuClassifier(rf_pub, n_jobs=-1, random_state=42)\n",
    "\n",
    "    n_fold = 10\n",
    "\n",
    "    features = Xy_train[Xy_train.columns[:-1]]\n",
    "    labels = Xy_train[Xy_train.columns[-1]]\n",
    "    X_train, X_test, y_train, y_test = train_test_split(features, labels, train_size=0.8, random_state=42)\n",
    "    \n",
    "    Xy_pub_train = np.hstack((X_train, y_train.values.reshape(-1, 1)))\n",
    "    Xy_pub_train = pd.DataFrame(Xy_pub_train, columns=Xy_train.columns)\n",
    "    Xy_pub_train.to_csv(ml_input_dir + 'Xy_pub_train.csv', index=False)\n",
    "    \n",
    "    Xy_pub_test = np.hstack((X_test, y_test.values.reshape(-1, 1)))\n",
    "    Xy_pub_test = pd.DataFrame(Xy_pub_test, columns=Xy_train.columns)\n",
    "    Xy_pub_test.to_csv(ml_input_dir + 'Xy_pub_test.csv', index=False)\n",
    "\n",
    "    weights_train = X_train[X_train.columns[-1]]\n",
    "    X_train = X_train[X_train.columns[:-1]]\n",
    "    weights_test = X_test[X_test.columns[-1]]\n",
    "    X_test = X_test[X_test.columns[:-1]]\n",
    "\n",
    "    search_space = {\n",
    "    'estimator__bootstrap': Categorical([True, False]), # values for boostrap can be either True or False\n",
    "    'estimator__max_depth': Integer(5, 20), # values of max_depth are integers\n",
    "    'estimator__max_features': Categorical([None, 'sqrt','log2']), \n",
    "    'estimator__min_samples_leaf': Integer(2, 20),\n",
    "    'estimator__min_samples_split': Integer(2, 30),\n",
    "    'estimator__n_estimators': Integer(10, 200),\n",
    "    'max_samples': Integer(int(0.5*(len(y_train)-sum(y_train))), int(0.9*(len(y_train)-sum(y_train))))\n",
    "    }\n",
    "\n",
    "    pub_bayes_search = BayesSearchCV(pub,\n",
    "                                     search_space,\n",
    "                                     n_iter=100, # specify how many iterations\n",
    "                                     scoring='f1',\n",
    "                                     n_jobs=4,\n",
    "                                     cv=n_fold,\n",
    "                                     verbose=1,\n",
    "                                     random_state=42)\n",
    "    pub_bayes_search.fit(X_train, y_train, sample_weight=weights_train)\n",
    "    \n",
    "    # Extract the optimization results\n",
    "    optimization_results = pub_bayes_search.cv_results_['mean_test_score']\n",
    "    \n",
    "    model_pub = pub_bayes_search.best_estimator_\n",
    "    model_pub_acc = pub_bayes_search.best_score_    \n",
    "    print('The highest F1-score during cross validation:', model_pub_acc)\n",
    "    \n",
    "    # save the model\n",
    "    with open(model_pub_file, 'wb') as f:\n",
    "        pickle.dump(model_pub, f)\n",
    "    \n",
    "    # Plot the Bayesian optimization progress\n",
    "    plt.figure(figsize=(12, 6))\n",
    "    plt.plot(range(1, len(optimization_results) + 1), optimization_results, marker='o', color='black', markerfacecolor='red')\n",
    "    plt.xlim(0, 51)\n",
    "    plt.xlabel('Bayesian Optimization Iteration')\n",
    "    plt.ylabel('Mean Test Precision')\n",
    "    plt.title('Bayesian Optimization Progress')\n",
    "    plt.grid(True, linestyle=':')\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(ml_output_dir + 'bayesian_optimization_pub_progress.png')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c1227e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(model_pub)\n",
    "\n",
    "Xy_pub_test_file = ml_input_dir + 'Xy_pub_test.csv'\n",
    "Xy_pub_test = pd.read_csv(Xy_pub_test_file, index_col=False)\n",
    "\n",
    "features = Xy_pub_test[Xy_pub_test.columns[:-2]]\n",
    "labels = Xy_pub_test[Xy_pub_test.columns[-1]]\n",
    "weights = Xy_pub_test[Xy_pub_test.columns[-2]]\n",
    "\n",
    "X_pred = model_pub.predict(features)\n",
    "\n",
    "# assuming that all zeros in labels are true negative samples\n",
    "cMatrix = confusion_matrix(labels, X_pred)\n",
    "X_pred_acc = accuracy_score(labels, X_pred, sample_weight=weights)\n",
    "X_pred_pre = precision_score(labels, X_pred, sample_weight=weights)\n",
    "X_pred_rec = recall_score(labels, X_pred, sample_weight=weights)\n",
    "X_pred_f1 = f1_score(labels, X_pred, sample_weight=weights)\n",
    "\n",
    "print('Confusion matrix:\\n', cMatrix)\n",
    "print('Accuracy:', X_pred_acc)\n",
    "print('Precision:', X_pred_pre)\n",
    "print('Recall:', X_pred_rec)\n",
    "print('F1-Score:', X_pred_f1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f9bf957",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xy_pos_test = pd.read_csv(ml_input_dir + 'Xy_pos_test.csv', index_col=False)\n",
    "X_pos_test = Xy_pos_test[Xy_pos_test.columns[:-2]]\n",
    "y_pos_test = Xy_pos_test[Xy_pos_test.columns[-1]]\n",
    "weights_test = Xy_pos_test[Xy_pos_test.columns[-2]]\n",
    "X_pos_pred = model_pub.predict(X_pos_test)\n",
    "X_pos_pred_acc = accuracy_score(y_pos_test, X_pos_pred, sample_weight=weights_test)\n",
    "print('Accuracy:', X_pos_pred_acc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f07e9879",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xy_train = pd.read_csv(ml_input_dir + 'Xy_train.csv', index_col=False)\n",
    "Xy_train_features = Xy_train[Xy_train.columns[:-2]]\n",
    "Xy_train_labels = Xy_train[Xy_train.columns[-1]]\n",
    "Xy_train_weights = Xy_train[Xy_train.columns[-2]]\n",
    "Xy_train_labels_new = model_pub.predict(Xy_train_features)\n",
    "\n",
    "Xy_train_labels_new_acc = accuracy_score(Xy_train_labels, Xy_train_labels_new, sample_weight=Xy_train_weights)\n",
    "Xy_train_labels_new_pre = precision_score(Xy_train_labels, Xy_train_labels_new, sample_weight=Xy_train_weights)\n",
    "Xy_train_labels_new_rec = recall_score(Xy_train_labels, Xy_train_labels_new, sample_weight=Xy_train_weights)\n",
    "Xy_train_labels_new_f1 = f1_score(Xy_train_labels, Xy_train_labels_new, sample_weight=Xy_train_weights)\n",
    "\n",
    "# assuming that all zeros in Xy_train_labels are true negative samples\n",
    "print('Accuracy:', Xy_train_labels_new_acc)\n",
    "print('Precision:', Xy_train_labels_new_pre)\n",
    "print('Recall:', Xy_train_labels_new_rec)\n",
    "print('F1-Score:', Xy_train_labels_new_f1)\n",
    "\n",
    "mask = Xy_train_labels == 1\n",
    "Xy_train_labels_new[mask] = 1\n",
    "Xy_train['label'] = Xy_train_labels_new\n",
    "Xy_train.to_csv(ml_output_dir + 'Xy_train_new.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1bd8e96",
   "metadata": {},
   "source": [
    "#### Random Forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37751b0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_rf_file = ml_output_dir + 'model_rf.pkl'\n",
    "\n",
    "if os.path.isfile(model_rf_file):\n",
    "    print('A model file exists. Attempting to load...')\n",
    "    try:\n",
    "        with open(model_rf_file, 'rb') as f:\n",
    "            model_rf = pickle.load(f)\n",
    "        print('Model loaded successfully.')\n",
    "    except (ModuleNotFoundError, ImportError) as e:\n",
    "        print(f\"Error loading the model: {e}\")\n",
    "        print(\"This may be due to a version mismatch. Proceeding to train a new model.\")\n",
    "        model_rf = None\n",
    "else:\n",
    "    print('No existing model found. Proceeding to train a new model.')\n",
    "    model_rf = None\n",
    "\n",
    "if model_rf is None:\n",
    "    # Random Forest model structure\n",
    "    rf = RandomForestClassifier(n_jobs=-1, random_state=42)\n",
    "\n",
    "    n_fold = 10\n",
    "\n",
    "    Xy_train_file = ml_output_dir + 'Xy_train_new.csv'\n",
    "    Xy_train = pd.read_csv(Xy_train_file, index_col=False)\n",
    "    features = Xy_train[Xy_train.columns[:-1]]\n",
    "    labels = Xy_train[Xy_train.columns[-1]]\n",
    "    X_train, X_test, y_train, y_test = train_test_split(features, labels, train_size=0.8, random_state=42)\n",
    "\n",
    "    Xy_rf_train = np.hstack((X_train, y_train.values.reshape(-1, 1)))\n",
    "    Xy_rf_train = pd.DataFrame(Xy_rf_train, columns=Xy_train.columns)\n",
    "    Xy_rf_train.to_csv(ml_input_dir + 'Xy_rf_train.csv', index=False)\n",
    "    \n",
    "    Xy_rf_test = np.hstack((X_test, y_test.values.reshape(-1, 1)))\n",
    "    Xy_rf_test = pd.DataFrame(Xy_rf_test, columns=Xy_train.columns)\n",
    "    Xy_rf_test.to_csv(ml_input_dir + 'Xy_rf_test.csv', index=False)\n",
    "\n",
    "    weights_train = X_train[X_train.columns[-1]]\n",
    "    X_train = X_train[X_train.columns[:-1]]\n",
    "    weights_test = X_test[X_test.columns[-1]]\n",
    "    X_test = X_test[X_test.columns[:-1]]\n",
    "\n",
    "    search_space = {\n",
    "    'bootstrap': Categorical([True, False]), # values for boostrap can be either True or False\n",
    "    'max_depth': Integer(5, 20), # values of max_depth are integers\n",
    "    'max_features': Categorical([None, 'sqrt','log2']), \n",
    "    'min_samples_leaf': Integer(2, 20),\n",
    "    'min_samples_split': Integer(2, 30),\n",
    "    'n_estimators': Integer(10, 200)\n",
    "    }\n",
    "\n",
    "    rf_bayes_search = BayesSearchCV(rf,\n",
    "                                    search_space,\n",
    "                                    n_iter=100, # specify how many iterations\n",
    "                                    scoring='f1',\n",
    "                                    n_jobs=4,\n",
    "                                    cv=n_fold,\n",
    "                                    verbose=1,\n",
    "                                    random_state=42)\n",
    "    rf_bayes_search.fit(X_train, y_train, sample_weight=weights_train)\n",
    "    \n",
    "    # Extract the optimization results\n",
    "    optimization_results = rf_bayes_search.cv_results_['mean_test_score']\n",
    "    \n",
    "    model_rf = rf_bayes_search.best_estimator_\n",
    "    model_rf_acc = rf_bayes_search.best_score_ \n",
    "    print('The highest F1-score during cross validation:', model_rf_acc)\n",
    "    \n",
    "    # save the model\n",
    "    with open(ml_output_dir + 'model_rf.pkl', 'wb') as f:\n",
    "        pickle.dump(model_rf, f)\n",
    "    \n",
    "    importances = []\n",
    "    estimators = model_rf.estimators_\n",
    "    importances = [estimators[j].feature_importances_.reshape(-1, 1) for j in range(len(estimators))]\n",
    "    importances = np.hstack(importances)\n",
    "\n",
    "    # Plot the Bayesian optimization progress\n",
    "    plt.figure(figsize=(12, 6))\n",
    "    plt.plot(range(1, len(optimization_results) + 1), optimization_results, marker='o', color='black', markerfacecolor='red')\n",
    "    plt.xlim(0, 51)\n",
    "    plt.xlabel('Bayesian Optimization Iteration')\n",
    "    plt.ylabel('Mean Test Accuracy')\n",
    "    plt.title('Bayesian Optimization Progress')\n",
    "    plt.grid(True, linestyle=':')\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(ml_output_dir + 'bayesian_optimization_rf_progress.png')\n",
    "    plt.show()\n",
    "else:\n",
    "    print(\"Using the loaded model. No new training performed.\")\n",
    "    print(\"To retrain the model, delete or rename the existing model file and run the script again.\")\n",
    "\n",
    "print(\"Script execution completed.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "486fa8e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(model_rf)\n",
    "\n",
    "Xy_rf_test_file = ml_input_dir + 'Xy_rf_test.csv'\n",
    "Xy_rf_test = pd.read_csv(Xy_rf_test_file, index_col=False)\n",
    "\n",
    "features = Xy_rf_test[Xy_rf_test.columns[:-2]]\n",
    "labels = Xy_rf_test[Xy_rf_test.columns[-1]]\n",
    "weights = Xy_rf_test[Xy_rf_test.columns[-2]]\n",
    "\n",
    "X_pred = model_rf.predict(features)\n",
    "\n",
    "cMatrix = confusion_matrix(labels, X_pred)\n",
    "X_pred_acc = accuracy_score(labels, X_pred, sample_weight=weights)\n",
    "X_pred_pre = precision_score(labels, X_pred, sample_weight=weights)\n",
    "X_pred_rec = recall_score(labels, X_pred, sample_weight=weights)\n",
    "X_pred_f1 = f1_score(labels, X_pred, sample_weight=weights)\n",
    "\n",
    "print('Confusion matrix:\\n', cMatrix)\n",
    "print('Accuracy:', X_pred_acc)\n",
    "print('Precision:', X_pred_pre)\n",
    "print('Recall:', X_pred_rec)\n",
    "print('F1-Score:', X_pred_f1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98e2669f",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xy_pos_test = pd.read_csv(ml_input_dir + 'Xy_pos_test.csv', index_col=False)\n",
    "X_pos_test = Xy_pos_test[Xy_pos_test.columns[:-2]]\n",
    "y_pos_test = Xy_pos_test[Xy_pos_test.columns[-1]]\n",
    "weights_test = Xy_pos_test[Xy_pos_test.columns[-2]]\n",
    "X_pos_pred = model_rf.predict(X_pos_test)\n",
    "X_pos_pred_acc = accuracy_score(y_pos_test, X_pos_pred, sample_weight=weights_test)\n",
    "print('Accuracy:', X_pos_pred_acc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8a8a7b2",
   "metadata": {},
   "source": [
    "#### ROC Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e88aba2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def roc_plot(y_test, z_test, n_classes, labels_name, average='macro'):\n",
    "    fpr = {}\n",
    "    tpr = {}\n",
    "    roc_auc = {}\n",
    "\n",
    "    y_test_dummies = pd.get_dummies(y_test).values\n",
    "    \n",
    "    for i in range(n_classes):\n",
    "        fpr[i], tpr[i], _ = roc_curve(y_test_dummies[:, i], z_test[:, i])\n",
    "        roc_auc[i] = auc(fpr[i], tpr[i])\n",
    "\n",
    "    # roc for each class\n",
    "    fig, ax = plt.subplots(figsize=(6, 6))\n",
    "    ax.plot([0, 1], [0, 1], 'k--')\n",
    "    ax.set_xlim([0.0, 1.0])\n",
    "    ax.set_ylim([0.0, 1.05])\n",
    "    ax.set_xlabel('False Positive Rate')\n",
    "    ax.set_ylabel('True Positive Rate')\n",
    "    ax.set_title('Receiver Operating Characteristic')\n",
    "    \n",
    "    for i in range(n_classes):\n",
    "        ax.plot(fpr[i], tpr[i], label='{}, AUC = {}'.format(labels_name[i], '{0:.4f}'.format(roc_auc[i])))\n",
    "    \n",
    "    ax.legend(loc='best')\n",
    "    ax.grid(alpha=0.5)\n",
    "    sns.despine()\n",
    "    plt.show()\n",
    "    print('ROC AUC score:', roc_auc_score(y_test_dummies, z_test, average=average))\n",
    "\n",
    "# predict for the test dataset\n",
    "z_test = model_rf.predict_proba(features)\n",
    "\n",
    "labels_name = ['Non-mineralised', 'Mineralised']\n",
    "roc_plot(labels, z_test, 2, labels_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2604e9b8",
   "metadata": {},
   "source": [
    "#### Feature Importance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78326980",
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_importances_file = ml_output_dir + 'feature_importances_muller2019.csv'\n",
    "\n",
    "if os.path.isfile(feature_importances_file):\n",
    "    feature_importances = pd.read_csv(feature_importances_file, index_col=False).to_numpy().tolist()\n",
    "else:\n",
    "    output_features_index = [selected_features.index(feature) for feature in features_list]\n",
    "    selected_features_names = parameters['selected_features_names_nounit_01']\n",
    "    selected_features_names = [selected_features_names[i] for i in output_features_index]\n",
    "\n",
    "    importances_mean = importances.mean(axis=1)\n",
    "    importances_var = importances.var(axis=1)\n",
    "\n",
    "    feature_importances = [(feature, round(importance, 4)) for feature, importance in zip(selected_features_names, importances_mean)]\n",
    "    feature_importances = sorted(feature_importances, key=lambda x:x[1], reverse=True)\n",
    "    feature_importances_df = pd.DataFrame(feature_importances, columns=['Feature', 'Importance'])\n",
    "    feature_importances_df['Variance'] = importances_var\n",
    "    \n",
    "    feature_importances_df.to_csv(feature_importances_file, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6956dc91",
   "metadata": {},
   "outputs": [],
   "source": [
    "# list of features sorted from most to least important\n",
    "sorted_importances = [importance[1] for importance in feature_importances]\n",
    "# cumulative importance\n",
    "cumulative_importances = np.cumsum(sorted_importances)\n",
    "\n",
    "x_values = list(range(len(feature_importances)))\n",
    "x_values = [x+1 for x in x_values]\n",
    "\n",
    "fig = plt.figure(figsize=(4, 4))\n",
    "ax2 = fig.add_subplot(111)\n",
    "ax1 = ax2.twinx()\n",
    "ax2.set_facecolor('whitesmoke')\n",
    "\n",
    "ax2.bar(x_values, sorted_importances, edgecolor='gray', facecolor='LightSalmon', width=1, alpha=0.5)\n",
    "ax1.plot(x_values, cumulative_importances, 'k--')\n",
    "\n",
    "plt.xlim(0.5, len(cumulative_importances)+0.5)\n",
    "ax1.set_ylim(0, 1.05)\n",
    "\n",
    "ax2.set_ylabel('Feature Importance')\n",
    "ax1.set_ylabel('Cumulative Importance')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a285014",
   "metadata": {},
   "outputs": [],
   "source": [
    "# print significant features above some threshold\n",
    "feature_importances.sort(key=lambda x:x[1])\n",
    "ft_imps = [x[1] for x in feature_importances]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_facecolor('whitesmoke')\n",
    "bar = ax.barh(range(len(ft_imps)), ft_imps)\n",
    "\n",
    "def gradientbars(bars, data):\n",
    "    ax = bars[0].axes\n",
    "    lim = ax.get_xlim()+ax.get_ylim()\n",
    "    for bar in bars:\n",
    "        bar.set_zorder(1)\n",
    "        bar.set_facecolor('none')\n",
    "        bar.set_edgecolor('black')\n",
    "        x, y = bar.get_xy()\n",
    "        w, h = bar.get_width(), bar.get_height()\n",
    "        cmap = plt.get_cmap('coolwarm')\n",
    "        grad = np.atleast_2d(np.linspace(0, 1*w/max(data), 256))\n",
    "        ax.imshow(grad, extent=[x, x+w, y, y+h], aspect='auto', zorder=0, norm=mpl.colors.NoNorm(vmin=0, vmax=1), cmap=cmap, alpha=0.8)\n",
    "        manual_labels = [x[0] for x in feature_importances]\n",
    "        ax.set_yticks(np.arange(0, len(data), 1).tolist())\n",
    "        ax.set_yticklabels(manual_labels, minor=False)\n",
    "    ax.axis(lim)\n",
    "    ax.set_xlabel('Feature Importance')\n",
    "\n",
    "gradientbars(bar, ft_imps)\n",
    "plt.gca().yaxis.grid(False)\n",
    "\n",
    "# plt.savefig(\n",
    "#     f'./figures/muller2019/importances.png',\n",
    "#     bbox_inches='tight',\n",
    "#     pad_inches=0.1,\n",
    "#     dpi=150\n",
    "#     )\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70b09a24",
   "metadata": {},
   "source": [
    "### Histogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2d852f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xy_train_new = pd.read_csv(ml_output_dir + 'Xy_train_new.csv', index_col=False)\n",
    "Xy_train_original = pd.read_csv(ml_input_dir + 'Xy_train_original.csv', index_col=False)\n",
    "\n",
    "@interact\n",
    "def show_map(feature=Xy_train_new.columns):\n",
    "    fig = plt.figure(figsize=(6, 6))\n",
    "    ax1 = fig.add_subplot(111)\n",
    "    ax1.set_facecolor('whitesmoke')\n",
    "    ax2 = ax1.twiny()\n",
    "    ax3 = ax2.twinx()\n",
    "    \n",
    "    # Calculate bin edges based on the entire dataset\n",
    "    min_val = Xy_train_new[feature].min()\n",
    "    max_val = Xy_train_new[feature].max()\n",
    "    bin_edges = np.linspace(min_val, max_val, 26)  # 25 bins, 26 edges\n",
    "    \n",
    "    # Histogram for original data (transparent)\n",
    "    ax1.hist(Xy_train_original[feature], bins=25, alpha=0.0)\n",
    "    \n",
    "    # Histograms for standardized data with black edges and consistent bin widths\n",
    "    h1 = ax2.hist(Xy_train_new.loc[Xy_train_new['label']==0][feature], bins=bin_edges, \n",
    "                  color='LightSalmon', label='Negative', \n",
    "                  edgecolor='black', linewidth=1)\n",
    "    h2 = ax2.hist(Xy_train_new.loc[Xy_train_new['label']==1][feature], bins=bin_edges, \n",
    "                  color='DarkSeaGreen', label='Positive', alpha=0.8, \n",
    "                  edgecolor='black', linewidth=1)\n",
    "    \n",
    "    kde_x = np.linspace(min_val, max_val, 100)\n",
    "    \n",
    "    kde1 = stats.gaussian_kde(Xy_train_new.loc[Xy_train_new['label']==0][feature])\n",
    "    kde2 = stats.gaussian_kde(Xy_train_new.loc[Xy_train_new['label']==1][feature])\n",
    "    \n",
    "    k1 = ax3.plot(kde_x, kde1(kde_x), color='LightSalmon')\n",
    "    k2 = ax3.plot(kde_x, kde2(kde_x), color='DarkSeaGreen')\n",
    "    \n",
    "    # Calculate maximum values for y-axes\n",
    "    max_freq = max(np.max(h1[0]), np.max(h2[0]))\n",
    "    max_density = max(np.max(kde1(kde_x)), np.max(kde2(kde_x)))\n",
    "    \n",
    "    # Set y-axis limits with some padding\n",
    "    ax2.set_ylim(0, max_freq * 1.1)\n",
    "    ax3.set_ylim(0, max_density * 1.1)\n",
    "    \n",
    "    # Adjust tick locations for both y-axes\n",
    "    ax2.yaxis.set_major_locator(plt.MaxNLocator(5))\n",
    "    ax3.yaxis.set_major_locator(plt.MaxNLocator(5))\n",
    "    \n",
    "    ax2.legend(loc='upper right')\n",
    "    \n",
    "    ax1.set_xlabel(feature + '\\n(Actual)')\n",
    "    ax2.set_xlabel(feature + '\\n(Standardised)')\n",
    "    ax1.set_ylabel('Frequency')\n",
    "    ax3.set_ylabel('Probability density')\n",
    "\n",
    "#     plt.savefig(\n",
    "#         f'./figures/muller2019/histograms/conv_angle_deg_edited.png',\n",
    "#         bbox_inches='tight',\n",
    "#         pad_inches=0.1,\n",
    "#         dpi=150\n",
    "#     )\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c5cfe6e",
   "metadata": {},
   "source": [
    "### Scatter Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df5cf140",
   "metadata": {},
   "outputs": [],
   "source": [
    "features_1 = Xy_train_new.columns.drop('label')\n",
    "features_2 = deque(features_1)\n",
    "features_2.rotate()\n",
    "features_3 = deque(features_2)\n",
    "features_3.rotate()\n",
    "\n",
    "@interact\n",
    "def show_map(feature_1=features_1, feature_2=features_2, feature_3=features_3):\n",
    "\n",
    "    fig = plt.figure(figsize=(6, 6))\n",
    "    gs = GridSpec(2, 2, hspace=0.4, wspace=0.2, height_ratios=[1, 0.03])\n",
    "    ax1 = fig.add_subplot(gs[0, :])\n",
    "    ax2 = ax1.twinx()\n",
    "    ax3 = ax2.twiny()\n",
    "\n",
    "    min1 = Xy_train_original[feature_1].min()\n",
    "    max1 = Xy_train_original[feature_1].max()\n",
    "    min2 = Xy_train_original[feature_2].min()\n",
    "    max2 = Xy_train_original[feature_2].max()\n",
    "\n",
    "    x_range = np.linspace(min1, max1+0.1*max1, num=100)\n",
    "    y_range = np.linspace(min2, max2+0.1*max2, num=100)\n",
    "    grid_x, grid_y = np.meshgrid(x_range, y_range)\n",
    "    \n",
    "    grid_data = griddata(list(zip(Xy_train_original[feature_1], Xy_train_original[feature_2])), Xy_train_original[feature_3],\n",
    "                         (grid_x, grid_y), method='nearest', fill_value=0)\n",
    "    grid_data = ndimage.gaussian_filter(grid_data, sigma=3)\n",
    "    cb1 = ax1.imshow(grid_data.T, extent=(min1, max1, min2, max2), origin='lower', aspect='auto', cmap=plt.cm.Spectral_r, alpha=0)\n",
    "\n",
    "    sc1 = ax1.scatter(Xy_train_original.loc[Xy_train_original['label']==0, feature_1],\n",
    "                      Xy_train_original.loc[Xy_train_original['label']==0, feature_2], 40, marker='.', c='blue', alpha=0)\n",
    "    sc2 = ax1.scatter(Xy_train_original.loc[Xy_train_original['label']==1, feature_1],\n",
    "                      Xy_train_original.loc[Xy_train_original['label']==1, feature_2], 40, marker='.', c='orange', alpha=0)\n",
    "\n",
    "    min1 = Xy_train_new[feature_1].min()\n",
    "    max1 = Xy_train_new[feature_1].max()\n",
    "    min2 = Xy_train_new[feature_2].min()\n",
    "    max2 = Xy_train_new[feature_2].max()\n",
    "    \n",
    "    x_range = np.linspace(min1, max1+0.1*max1, num=100)\n",
    "    y_range = np.linspace(min2, max2+0.1*max2, num=100)\n",
    "    grid_x, grid_y = np.meshgrid(x_range, y_range)\n",
    "    \n",
    "    grid_data = griddata(list(zip(Xy_train_new[feature_1], Xy_train_new[feature_2])), Xy_train_new[feature_3],\n",
    "                         (grid_x, grid_y), method='nearest', fill_value=0)\n",
    "    grid_data = ndimage.gaussian_filter(grid_data, sigma=3)\n",
    "    cb2 = ax3.imshow(grid_data.T, extent=(min1, max1, min2, max2), origin='lower', aspect='auto', cmap=plt.cm.Spectral_r, alpha=0.7)\n",
    "    \n",
    "    sc3 = ax3.scatter(Xy_train_new.loc[Xy_train_new['label']==1, feature_1],\n",
    "                      Xy_train_new.loc[Xy_train_new['label']==1, feature_2],\n",
    "                      100, marker='.', facecolor='mediumseagreen', edgecolor='black')\n",
    "    sc4 = ax3.scatter(Xy_train_new.loc[Xy_train_new['label']==0, feature_1],\n",
    "                      Xy_train_new.loc[Xy_train_new['label']==0, feature_2],\n",
    "                      100, marker='.', facecolor='orangered', edgecolor='black', alpha=0.7)\n",
    "        \n",
    "    ax3.legend([sc3, sc4], ['Positive', 'Negative'], loc='best',  borderaxespad=0.1, fontsize=8) # numpoints=1\n",
    "\n",
    "    ax1.set_xlabel(feature_1 + '\\n(Actual)')\n",
    "    ax1.set_ylabel(feature_2 + '\\n(Actual)')\n",
    "    ax2.set_ylabel(feature_2 + '\\n(Standardised)')\n",
    "    ax3.set_xlabel(feature_1 + '\\n(Standardised)')\n",
    "        \n",
    "    cax1 = fig.add_subplot(gs[1, 1])\n",
    "    cax2 = fig.add_subplot(gs[1, 0])\n",
    "        \n",
    "    fig.colorbar(cb2, cax=cax2, orientation='horizontal', label=feature_3 + '\\n(Standardised)', extend='both')\n",
    "    fig.colorbar(cb1, cax=cax1, orientation='horizontal', label=feature_3 + '\\n(Actual)', extend='both')\n",
    "    \n",
    "#     plt.savefig(\n",
    "#         f'./figures/muller2019/features_three/dist_nearest_edge_deg.png',\n",
    "#         bbox_inches='tight',\n",
    "#         pad_inches=0.1,\n",
    "#         dpi=150\n",
    "#     )\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d9c8c4c",
   "metadata": {},
   "source": [
    "### Boxplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82e07e30",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xy_train_pivot = Xy_train_new.pivot(columns=['label'])\n",
    "Xy_train_original_pivot = Xy_train_original.pivot(columns=['label'])\n",
    "nb_groups1 = Xy_train_new['label'].nunique()\n",
    "nb_groups2 = Xy_train_original['label'].nunique()\n",
    "\n",
    "@interact\n",
    "def show_map(feature=Xy_train_new.columns):\n",
    "    bplot1 = [Xy_train_pivot[feature][var].dropna() for var in Xy_train_pivot[feature]]\n",
    "    bplot2 = [Xy_train_original_pivot[feature][var].dropna() for var in Xy_train_original_pivot[feature]]\n",
    "    fig, ax1 = plt.subplots(figsize=(6, 6))\n",
    "    ax1.set_facecolor('whitesmoke')\n",
    "    \n",
    "    # Define colors for negative and positive samples\n",
    "    colors = ['LightSalmon', 'DarkSeaGreen']\n",
    "    \n",
    "    # Create box plots for standardized data with different colors\n",
    "    bp1 = ax1.boxplot(bplot1, positions=np.arange(nb_groups1), patch_artist=True,\n",
    "                      whis=(5, 95), widths=0.2,\n",
    "                      flierprops=dict(marker='.', markeredgecolor='black', fillstyle=None),\n",
    "                      medianprops=dict(color='black'))\n",
    "    \n",
    "    # Color the boxes\n",
    "    for patch, color in zip(bp1['boxes'], colors):\n",
    "        patch.set_facecolor(color)\n",
    "    \n",
    "    # Create transparent box plots for original data\n",
    "    box_param2 = dict(whis=(5, 95), widths=0, patch_artist=True,\n",
    "                      flierprops=dict(marker='.', markeredgecolor='none', fillstyle=None),\n",
    "                      medianprops=dict(color='none'), whiskerprops=dict(color='none'),\n",
    "                      boxprops=dict(facecolor='none', edgecolor='none'))\n",
    "    ax2 = ax1.twinx()\n",
    "    ax2.boxplot(bplot2, positions=np.arange(nb_groups2), **box_param2)\n",
    "    \n",
    "    # Format x ticks\n",
    "    labelsize = 12\n",
    "    ax1.set_xticks(np.arange(nb_groups1))\n",
    "    ax1.set_xticklabels(['Negative', 'Positive'])\n",
    "    ax1.tick_params(axis='x', labelsize=labelsize)\n",
    "    \n",
    "    # Format y ticks\n",
    "    ax1.tick_params(axis='y', labelsize=labelsize)\n",
    "    ax2.tick_params(axis='y', labelsize=labelsize)\n",
    "    \n",
    "    # Format axes labels\n",
    "    label_fmt = dict(size=12, labelpad=15)\n",
    "    ax1.set_xlabel(feature, **label_fmt)\n",
    "    ax1.set_ylabel(feature + '\\n(Standardised)', **label_fmt)\n",
    "    ax2.set_ylabel(feature + '\\n(Actual)', **label_fmt)\n",
    "        \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24584512",
   "metadata": {},
   "source": [
    "### Violin Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ec90b60",
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to calculate whiskers (for outliers)\n",
    "def calculate_whiskers(data):\n",
    "    q1, q3 = np.percentile(data, [25, 75])\n",
    "    iqr = q3 - q1\n",
    "    whisker_low = q1 - 1.5 * iqr\n",
    "    whisker_high = q3 + 1.5 * iqr\n",
    "    return whisker_low, whisker_high\n",
    "\n",
    "# pivot dataframes\n",
    "Xy_train_new_pivot = Xy_train_new.pivot(columns=['label'])\n",
    "Xy_train_original_pivot = Xy_train_original.pivot(columns=['label'])\n",
    "\n",
    "# calculate the number of unique groups\n",
    "nb_groups1 = Xy_train_new['label'].nunique()\n",
    "nb_groups2 = Xy_train_original['label'].nunique()\n",
    "\n",
    "colors = ['LightSalmon', 'DarkSeaGreen']\n",
    "\n",
    "@interact\n",
    "def show_map(feature=Xy_train_new.columns):\n",
    "    vplot1_data = [Xy_train_new_pivot[feature][var].dropna() for var in Xy_train_new_pivot[feature]]\n",
    "    vplot2_data = [Xy_train_original_pivot[feature][var].dropna() for var in Xy_train_original_pivot[feature]]\n",
    "\n",
    "    fig, ax1 = plt.subplots(figsize=(6, 6))\n",
    "    ax1.set_facecolor('whitesmoke')\n",
    "\n",
    "    # create violin plots\n",
    "    vplot1_parts = ax1.violinplot(vplot1_data, positions=np.arange(nb_groups1))\n",
    "    \n",
    "    for i, part in enumerate(vplot1_parts['bodies']):\n",
    "        part.set_facecolor(colors[i])\n",
    "        part.set_edgecolor(colors[i])\n",
    "        \n",
    "    vplot1_parts['cbars'].set_edgecolor('black')\n",
    "    vplot1_parts['cmins'].set_edgecolor('black')\n",
    "    vplot1_parts['cmaxes'].set_edgecolor('black')\n",
    "    \n",
    "    ax2 = ax1.twinx()\n",
    "    vplot2_parts = ax2.violinplot(vplot2_data, positions=np.arange(nb_groups2))\n",
    "    \n",
    "    for part in vplot2_parts['bodies']:\n",
    "        part.set_facecolor('none')\n",
    "        part.set_edgecolor('none')\n",
    "        \n",
    "    vplot2_parts['cbars'].set_edgecolor('none')\n",
    "    vplot2_parts['cmins'].set_edgecolor('none')\n",
    "    vplot2_parts['cmaxes'].set_edgecolor('none')\n",
    "\n",
    "    # plotting outliers for the first violin plot\n",
    "    for i, data in enumerate(vplot1_data):\n",
    "        low, high = calculate_whiskers(data)\n",
    "        outliers = data[(data > high) | (data < low)]\n",
    "        ax1.scatter([i]*len(outliers), outliers, facecolor='red', edgecolor='black', s=20)\n",
    "\n",
    "    # format x ticks\n",
    "    labelsize = 12\n",
    "    ax1.set_xticks(np.arange(nb_groups1))\n",
    "    ax1.set_xticklabels(['Negative', 'Positive'])\n",
    "\n",
    "    # format axes labels\n",
    "    ax1.set_ylabel(feature + '\\n(Standardised)')\n",
    "    ax2.set_ylabel(feature + '\\n(Actual)')\n",
    "        \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "684b22d7",
   "metadata": {},
   "source": [
    "### Box-Violin Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69616fd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# pivot dataframes\n",
    "Xy_train_pivot = Xy_train_new.pivot(columns=['label'])\n",
    "Xy_train_original_pivot = Xy_train_original.pivot(columns=['label'])\n",
    "nb_groups1 = Xy_train_new['label'].nunique()\n",
    "nb_groups2 = Xy_train_original['label'].nunique()\n",
    "\n",
    "@interact\n",
    "def show_combined_plot(feature=Xy_train_new.columns):\n",
    "    # Prepare data\n",
    "    plot_data1 = [Xy_train_pivot[feature][var].dropna() for var in Xy_train_pivot[feature]]\n",
    "    plot_data2 = [Xy_train_original_pivot[feature][var].dropna() for var in Xy_train_original_pivot[feature]]\n",
    "    \n",
    "    # Create figure and axes\n",
    "    fig, ax1 = plt.subplots(figsize=(8, 6))\n",
    "    ax1.set_facecolor('whitesmoke')\n",
    "    \n",
    "    # Define colors\n",
    "    colors = ['LightSalmon', 'DarkSeaGreen']\n",
    "    \n",
    "    # Create violin plots first (semi-transparent)\n",
    "    vplot1_parts = ax1.violinplot(plot_data1, positions=np.arange(nb_groups1))\n",
    "    \n",
    "    # Style violin plots\n",
    "    for i, body in enumerate(vplot1_parts['bodies']):\n",
    "        body.set_facecolor(colors[i])\n",
    "        body.set_edgecolor(colors[i])\n",
    "        \n",
    "    vplot1_parts['cbars'].set_edgecolor('black')\n",
    "    vplot1_parts['cmins'].set_edgecolor('black')\n",
    "    vplot1_parts['cmaxes'].set_edgecolor('black')\n",
    "    \n",
    "    # Create box plots on top with gray whiskers and no caps\n",
    "    bp1 = ax1.boxplot(plot_data1, positions=np.arange(nb_groups1), patch_artist=True,\n",
    "                      whis=(5, 95), widths=0.2,\n",
    "                      flierprops=dict(marker='.', markersize=8, markerfacecolor='red', markeredgecolor='black', fillstyle=None),\n",
    "                      medianprops=dict(color='black', linewidth=1.5),\n",
    "                      whiskerprops=dict(color='black', linestyle='-'),  # Gray whiskers\n",
    "                      capprops=dict(visible=False))  # Remove whisker caps\n",
    "    \n",
    "    # Color the boxes\n",
    "    for patch, color in zip(bp1['boxes'], colors):\n",
    "        patch.set_facecolor(color)\n",
    "    \n",
    "    # Create transparent box plots for original data\n",
    "    ax2 = ax1.twinx()\n",
    "    box_param2 = dict(whis=(5, 95), widths=0, patch_artist=True,\n",
    "                     flierprops=dict(marker='.', markeredgecolor='none', fillstyle=None),\n",
    "                     medianprops=dict(color='none'), whiskerprops=dict(color='none'),\n",
    "                     boxprops=dict(facecolor='none', edgecolor='none'))\n",
    "    ax2.boxplot(plot_data2, positions=np.arange(nb_groups2), **box_param2)\n",
    "        \n",
    "    # Format axes\n",
    "    ax1.set_xticks(np.arange(nb_groups1))\n",
    "    ax1.set_xticklabels(['Negative', 'Positive'])\n",
    "    ax1.tick_params(axis='x')\n",
    "    ax1.tick_params(axis='y')\n",
    "    ax2.tick_params(axis='y')\n",
    "    \n",
    "    # Format labels\n",
    "    ax1.set_ylabel(feature + '\\n(Standardised)')\n",
    "    ax2.set_ylabel(feature + '\\n(Actual)')\n",
    "    \n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10afdeb6",
   "metadata": {},
   "source": [
    "### Probability in Back-Arc Basins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1a4b5c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "target_points_prob_files_lst = []\n",
    "target_points_prob_lst = []\n",
    "\n",
    "for time in time_steps:\n",
    "    target_points_prob_files_lst.append(ml_output_dir + f'target_points_prob_{time}_Ma.csv')\n",
    "\n",
    "for i, target_points_prob_file in enumerate(target_points_prob_files_lst):\n",
    "    if not(os.path.isfile(target_points_prob_file)):\n",
    "        df = target_points_ml_in_lst[target_points_prob_files_lst.index(target_points_prob_file)]\n",
    "        probs = model_rf.predict_proba(df)[:, 1].reshape(-1, 1)\n",
    "\n",
    "        mm_scaler1 = MinMaxScaler()\n",
    "        probs_scaled = mm_scaler1.fit_transform(probs)\n",
    "\n",
    "        df_xy = df.copy()\n",
    "        df_xy['lon'] = target_points_coreg_out_lst[i]['lon'].to_numpy()\n",
    "        df_xy['lat'] = target_points_coreg_out_lst[i]['lat'].to_numpy()\n",
    "        df_xy['prob'] = probs_scaled\n",
    "        df_xy.to_csv(target_points_prob_file, index=False)\n",
    "        target_points_prob_lst.append(df_xy)\n",
    "        \n",
    "for target_points_prob_file in target_points_prob_files_lst:\n",
    "    target_points_prob_lst.append(pd.read_csv(target_points_prob_file, index_col=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee829648",
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.LambertAzimuthalEqualArea(150, 0)\n",
    "\n",
    "@interact\n",
    "def show_map(time=time_steps):\n",
    "    lons_lats_recon = []\n",
    "    \n",
    "    for min_occ in min_occ_data.iterrows():\n",
    "        if time == 0:\n",
    "            lons_lats_recon.append((min_occ[1]['lon'], min_occ[1]['lat'], min_occ[1]['weight']))\n",
    "        elif int(min_occ[1]['age']) < time:\n",
    "            lons_lats_recon.append((np.nan, np.nan, np.nan))\n",
    "        elif int(min_occ[1]['age']) == time:\n",
    "            lons_lats_recon.append((min_occ[1]['lon_recon'], min_occ[1]['lat_recon'], min_occ[1]['weight']))\n",
    "        else:\n",
    "            lat_lon_recon = get_recon_ccords([min_occ[1]['lon']],\n",
    "                                             [min_occ[1]['lat']],\n",
    "                                             plate_motion_model='muller2019',\n",
    "                                             time=time)[0]\n",
    "            lons_lats_recon.append(tuple((lat_lon_recon[1], lat_lon_recon[0], min_occ[1]['weight'])))\n",
    "            \n",
    "    # call the PlotTopologies object\n",
    "    gplot = gplately.PlotTopologies(model, coastlines, continents, cob, time=time)\n",
    "    \n",
    "    if plate_motion_model == 'muller2016':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2016_AREPS_v1.17_AgeGrid-{time}.nc'\n",
    "    elif plate_motion_model == 'muller2019':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2019_Tectonics_v2.0_AgeGrid-{time}.nc'\n",
    "\n",
    "    agegrid = gplately.grids.read_netcdf_grid(agegrid_file)\n",
    "    \n",
    "    plot_x = target_points_prob_lst[time_steps.index(time)]['lon']\n",
    "    plot_y = target_points_prob_lst[time_steps.index(time)]['lat']\n",
    "    \n",
    "    # dual colour bars\n",
    "    fig = plt.figure(figsize=(6, 8))\n",
    "    gs = GridSpec(2, 2, hspace=-0.75, wspace=0.1, height_ratios=[1, 0.01])\n",
    "    ax = fig.add_subplot(gs[0, :], projection=proj)\n",
    "    \n",
    "    set_ax(ax, target_extent_map, 10, 5, stock_img=False, order=9)\n",
    "\n",
    "    im = gplot.plot_grid(ax, agegrid.data, cmap=ccm.lapaz_r, vmin=0, vmax=153, alpha=0.5, zorder=1)\n",
    "\n",
    "    gplot.plot_continents(ax, edgecolor='gray', facecolor='tan', zorder=2)\n",
    "    gplot.plot_ridges(ax, color='red', alpha=0.5, zorder=3)\n",
    "    gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=False, regrid_shape=10, alpha=0.2, zorder=4)\n",
    "    \n",
    "    sc0 = ax.scatter(\n",
    "        plot_x,\n",
    "        plot_y,\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        marker='.',\n",
    "        c=target_points_prob_lst[time_steps.index(time)]['prob'],\n",
    "        s=30,\n",
    "        cmap=ccm.hawaii_r,\n",
    "        zorder=5\n",
    "    )\n",
    "    \n",
    "    sc1 = ax.scatter(\n",
    "        [coords[0] for coords in lons_lats_recon],\n",
    "        [coords[1] for coords in lons_lats_recon],\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        marker='o',\n",
    "        facecolor='yellow',\n",
    "        edgecolor='black',\n",
    "        s=[x * 20 for x in [coords[2] for coords in lons_lats_recon]],\n",
    "        alpha=0.7,\n",
    "        zorder=6\n",
    "    )\n",
    "\n",
    "    gplot.plot_trenches(ax, color='k', alpha=0.3, zorder=7)\n",
    "    gplot.plot_subduction_teeth(ax, spacing=0.03, color='k', alpha=0.3, zorder=8)\n",
    "\n",
    "    cax1 = fig.add_subplot(gs[1, 0])\n",
    "    cax2 = fig.add_subplot(gs[1, 1])\n",
    "    \n",
    "    fig.colorbar(sc0, cax=cax2, orientation='horizontal', label='Mineralisation probability')\n",
    "    fig.colorbar(im, cax=cax1, orientation='horizontal', label='Seafloor age (Ma)', extend='max')\n",
    "    \n",
    "    # Define custom legend handles\n",
    "    custom_handles = [\n",
    "        Patch(facecolor='tan', edgecolor='gray', label='Continental Crust'),  # Custom handle for the filled polygon\n",
    "        Line2D([0], [0], color='red', lw=2, label='Mid-Ocean Ridge'),  # Custom handle for the line (ridge)\n",
    "        Line2D([0], [0], marker='o', markerfacecolor='yellow', markeredgecolor='black', markersize=5, linestyle='None', label='Mineral Occurrence')  # Custom handle for mineral occurrences\n",
    "    ]\n",
    "\n",
    "    # Add the custom legend to the plot\n",
    "    legend = ax.legend(handles=custom_handles, loc='upper right', bbox_to_anchor=(1.42, 1), borderaxespad=0.)\n",
    "    \n",
    "    ax.set_title(f'Porphyry Mineralisation Probability {time} Ma')\n",
    "        \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ea26650",
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_coords_lst = []\n",
    "for mask_coords_file in mask_coords_files_lst:\n",
    "    mask_coords_lst.append(pd.read_csv(mask_coords_file, index_col=False))\n",
    "\n",
    "proj = ccrs.LambertAzimuthalEqualArea(150, 0)\n",
    "\n",
    "@interact\n",
    "def show_map(time=time_steps):\n",
    "    lons_lats_recon = []\n",
    "    \n",
    "    for min_occ in min_occ_data.iterrows():\n",
    "        if time == 0:\n",
    "            lons_lats_recon.append((min_occ[1]['lon'], min_occ[1]['lat'], min_occ[1]['weight']))\n",
    "        elif int(min_occ[1]['age']) < time:\n",
    "            lons_lats_recon.append((np.nan, np.nan, np.nan))\n",
    "        elif int(min_occ[1]['age']) == time:\n",
    "            lons_lats_recon.append((min_occ[1]['lon_recon'], min_occ[1]['lat_recon'], min_occ[1]['weight']))\n",
    "        else:\n",
    "            lat_lon_recon = get_recon_ccords([min_occ[1]['lon']],\n",
    "                                             [min_occ[1]['lat']],\n",
    "                                             plate_motion_model='muller2019',\n",
    "                                             time=time)[0]\n",
    "            lons_lats_recon.append(tuple((lat_lon_recon[1], lat_lon_recon[0], min_occ[1]['weight'])))\n",
    "    \n",
    "    # call the PlotTopologies object\n",
    "    gplot = gplately.PlotTopologies(model, coastlines, continents, cob, time=time)\n",
    "    \n",
    "    if plate_motion_model == 'muller2016':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2016_AREPS_v1.17_AgeGrid-{time}.nc'\n",
    "    elif plate_motion_model == 'muller2019':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2019_Tectonics_v2.0_AgeGrid-{time}.nc'\n",
    "\n",
    "    agegrid = gplately.grids.read_netcdf_grid(agegrid_file)\n",
    "    \n",
    "    mask_coords = mask_coords_lst[time_steps.index(time)]\n",
    "    \n",
    "    probabilities = []\n",
    "    count = 0\n",
    "            \n",
    "    for mask in mask_coords['include']:\n",
    "        if mask:\n",
    "            probabilities.append(target_points_prob_lst[time_steps.index(time)]['prob'][count])\n",
    "            count += 1\n",
    "        else:\n",
    "            probabilities.append(np.nan)\n",
    "    \n",
    "    nx = mask_coords_lst[time_steps.index(time)]['lon'].nunique()\n",
    "    ny = mask_coords_lst[time_steps.index(time)]['lat'].nunique()\n",
    "    \n",
    "    x_min = mask_coords_lst[time_steps.index(time)]['lon'].min()\n",
    "    x_max = mask_coords_lst[time_steps.index(time)]['lon'].max()\n",
    "    y_min = mask_coords_lst[time_steps.index(time)]['lat'].min()\n",
    "    y_max = mask_coords_lst[time_steps.index(time)]['lat'].max()\n",
    "    \n",
    "    probabilities_2d = np.reshape(probabilities, (ny, nx))\n",
    "    probabilities_2d_ud = np.flipud(np.reshape(probabilities, (ny, nx)))\n",
    "    \n",
    "    # dual colour bars\n",
    "    fig = plt.figure(figsize=(6, 8))\n",
    "    gs = GridSpec(2, 2, hspace=-0.75, wspace=0.1, height_ratios=[1, 0.01])\n",
    "    ax = fig.add_subplot(gs[0, :], projection=proj)\n",
    "    \n",
    "    set_ax(ax, target_extent_map, 10, 5, stock_img=False, order=9)\n",
    "\n",
    "    im0 = gplot.plot_grid(ax, agegrid.data, cmap=ccm.lapaz_r, vmin=0, vmax=153, alpha=0.5, zorder=1)\n",
    "\n",
    "    gplot.plot_continents(ax, edgecolor='gray', facecolor='tan', zorder=2)\n",
    "    gplot.plot_ridges(ax, color='red', alpha=0.5, zorder=3)\n",
    "    gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=False, regrid_shape=10, alpha=0.2, zorder=4)\n",
    "    \n",
    "    im1 = plt.imshow(\n",
    "        probabilities_2d,\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        origin='lower',\n",
    "        cmap=ccm.hawaii_r,\n",
    "        interpolation='nearest',\n",
    "        extent=(x_min, x_max, y_min, y_max),\n",
    "        zorder=5\n",
    "    )\n",
    "    \n",
    "    sc1 = ax.scatter(\n",
    "        [coords[0] for coords in lons_lats_recon],\n",
    "        [coords[1] for coords in lons_lats_recon],\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        marker='o',\n",
    "        facecolor='yellow',\n",
    "        edgecolor='black',\n",
    "        s=[x * 20 for x in [coords[2] for coords in lons_lats_recon]],\n",
    "        alpha=0.7,\n",
    "        zorder=6\n",
    "    )\n",
    "\n",
    "    gplot.plot_trenches(ax, color='k', alpha=0.3, zorder=7)\n",
    "    gplot.plot_subduction_teeth(ax, spacing=0.03, color='k', alpha=0.3, zorder=8)\n",
    "    \n",
    "    cax1 = fig.add_subplot(gs[1, 0])\n",
    "    cax2 = fig.add_subplot(gs[1, 1])\n",
    "    \n",
    "    fig.colorbar(im1, cax=cax2, orientation='horizontal', label='Mineralisation probability')\n",
    "    fig.colorbar(im0, cax=cax1, orientation='horizontal', label='Seafloor age (Ma)', extend='max')\n",
    "    \n",
    "    # Define custom legend handles\n",
    "    custom_handles = [\n",
    "        Patch(facecolor='tan', edgecolor='gray', label='Continental Crust'),  # Custom handle for the filled polygon\n",
    "        Line2D([0], [0], color='red', lw=2, label='Mid-Ocean Ridge'),  # Custom handle for the line (ridge)\n",
    "        Line2D([0], [0], marker='o', markerfacecolor='yellow', markeredgecolor='black', markersize=5, linestyle='None', label='Mineral Occurrence')  # Custom handle for mineral occurrences\n",
    "    ]\n",
    "\n",
    "    # Add the custom legend to the plot\n",
    "    legend = ax.legend(handles=custom_handles, loc='upper right', bbox_to_anchor=(1.42, 1), borderaxespad=0.)\n",
    "    \n",
    "    ax.set_title(f'Porphyry Mineralisation Probability {time} Ma')\n",
    "        \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "072bb55c",
   "metadata": {},
   "source": [
    "### Probability in New Guinea"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7083ca0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "target_points_ng_prob_files_lst = []\n",
    "target_points_ng_prob_lst = []\n",
    "\n",
    "for time in time_steps:\n",
    "    target_points_ng_prob_files_lst.append(ml_output_dir + f'target_points_ng_prob_{time}_Ma.csv')\n",
    "\n",
    "for i, target_points_ng_prob_file in enumerate(target_points_ng_prob_files_lst):\n",
    "    if not(os.path.isfile(target_points_ng_prob_file)):\n",
    "        df = target_points_ng_ml_in_lst[target_points_ng_prob_files_lst.index(target_points_ng_prob_file)]\n",
    "        probs = model_rf.predict_proba(df)[:, 1].reshape(-1, 1)\n",
    "\n",
    "        mm_scaler1 = MinMaxScaler()\n",
    "        probs_scaled = mm_scaler1.fit_transform(probs)\n",
    "\n",
    "        df_xy = df.copy()\n",
    "        df_xy['lon'] = target_points_ng_coreg_out_lst[i]['lon'].to_numpy()\n",
    "        df_xy['lat'] = target_points_ng_coreg_out_lst[i]['lat'].to_numpy()\n",
    "        df_xy['prob'] = probs_scaled\n",
    "        df_xy.to_csv(target_points_ng_prob_file, index=False)\n",
    "        target_points_ng_prob_lst.append(df_xy)\n",
    "        \n",
    "for target_points_ng_prob_file in target_points_ng_prob_files_lst:\n",
    "    target_points_ng_prob_lst.append(pd.read_csv(target_points_ng_prob_file, index_col=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d325ffa1",
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.LambertAzimuthalEqualArea(150, 0)\n",
    "\n",
    "@interact\n",
    "def show_map(time=time_steps):\n",
    "    lons_lats_recon = []\n",
    "    \n",
    "    for min_occ in min_occ_data.iterrows():\n",
    "        if time == 0:\n",
    "            lons_lats_recon.append((min_occ[1]['lon'], min_occ[1]['lat'], min_occ[1]['weight']))\n",
    "        elif int(min_occ[1]['age']) < time:\n",
    "            lons_lats_recon.append((np.nan, np.nan, np.nan))\n",
    "        elif int(min_occ[1]['age']) == time:\n",
    "            lons_lats_recon.append((min_occ[1]['lon_recon'], min_occ[1]['lat_recon'], min_occ[1]['weight']))\n",
    "        else:\n",
    "            lat_lon_recon = get_recon_ccords([min_occ[1]['lon']],\n",
    "                                             [min_occ[1]['lat']],\n",
    "                                             plate_motion_model='muller2019',\n",
    "                                             time=time)[0]\n",
    "            lons_lats_recon.append(tuple((lat_lon_recon[1], lat_lon_recon[0], min_occ[1]['weight'])))\n",
    "            \n",
    "    # call the PlotTopologies object\n",
    "    gplot = gplately.PlotTopologies(model, coastlines, continents, cob, time=time)\n",
    "    \n",
    "    if plate_motion_model == 'muller2016':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2016_AREPS_v1.17_AgeGrid-{time}.nc'\n",
    "    elif plate_motion_model == 'muller2019':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2019_Tectonics_v2.0_AgeGrid-{time}.nc'\n",
    "\n",
    "    agegrid = gplately.grids.read_netcdf_grid(agegrid_file)\n",
    "    \n",
    "    plot_x = target_points_ng_prob_lst[time_steps.index(time)]['lon']\n",
    "    plot_y = target_points_ng_prob_lst[time_steps.index(time)]['lat']\n",
    "    \n",
    "    # dual colour bars\n",
    "    fig = plt.figure(figsize=(6, 8))\n",
    "    gs = GridSpec(2, 2, hspace=-0.75, wspace=0.1, height_ratios=[1, 0.01])\n",
    "    ax = fig.add_subplot(gs[0, :], projection=proj)\n",
    "    \n",
    "    set_ax(ax, target_extent_map, 10, 5, stock_img=False, order=9)\n",
    "\n",
    "    im = gplot.plot_grid(ax, agegrid.data, cmap=ccm.lapaz_r, vmin=0, vmax=230, alpha=0.5, zorder=1)\n",
    "\n",
    "    gplot.plot_continents(ax, edgecolor='gray', facecolor='tan', zorder=2)\n",
    "    gplot.plot_ridges(ax, color='red', alpha=0.5, zorder=3)\n",
    "    gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=False, regrid_shape=10, alpha=0.2, zorder=4)\n",
    "    \n",
    "    sc0 = ax.scatter(\n",
    "        plot_x,\n",
    "        plot_y,\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        marker='.',\n",
    "        c=target_points_ng_prob_lst[time_steps.index(time)]['prob'],\n",
    "        s=30,\n",
    "        cmap=ccm.hawaii_r,\n",
    "        zorder=5\n",
    "    )\n",
    "    \n",
    "    sc1 = ax.scatter(\n",
    "        [coords[0] for coords in lons_lats_recon],\n",
    "        [coords[1] for coords in lons_lats_recon],\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        marker='o',\n",
    "        facecolor='yellow',\n",
    "        edgecolor='black',\n",
    "        s=[x * 20 for x in [coords[2] for coords in lons_lats_recon]],\n",
    "        alpha=0.7,\n",
    "        zorder=6\n",
    "    )\n",
    "\n",
    "    gplot.plot_trenches(ax, color='k', alpha=0.3, zorder=7)\n",
    "    gplot.plot_subduction_teeth(ax, spacing=0.03, color='k', alpha=0.3, zorder=8)\n",
    "\n",
    "    cax1 = fig.add_subplot(gs[1, 0])\n",
    "    cax2 = fig.add_subplot(gs[1, 1])\n",
    "    \n",
    "    fig.colorbar(sc0, cax=cax2, orientation='horizontal', label='Mineralisation probability')\n",
    "    fig.colorbar(im, cax=cax1, orientation='horizontal', label='Seafloor age (Ma)', extend='max')\n",
    "    \n",
    "    # Define custom legend handles\n",
    "    custom_handles = [\n",
    "        Patch(facecolor='tan', edgecolor='gray', label='Continental Crust'),  # Custom handle for the filled polygon\n",
    "        Line2D([0], [0], color='red', lw=2, label='Mid-Ocean Ridge'),  # Custom handle for the line (ridge)\n",
    "        Line2D([0], [0], marker='o', markerfacecolor='yellow', markeredgecolor='black', markersize=5, linestyle='None', label='Mineral Occurrence')  # Custom handle for mineral occurrences\n",
    "    ]\n",
    "\n",
    "    # Add the custom legend to the plot\n",
    "    legend = ax.legend(handles=custom_handles, loc='upper right', bbox_to_anchor=(1.42, 1), borderaxespad=0.)\n",
    "    \n",
    "    ax.set_title(f'Porphyry Mineralisation Probability {time} Ma')\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9551659",
   "metadata": {},
   "outputs": [],
   "source": [
    "target_points_ng_probs_file = ml_output_dir + 'target_points_ng_prob.csv'\n",
    "\n",
    "if not(os.path.isfile(target_points_ng_probs_file)):\n",
    "    prob_lst = []\n",
    "\n",
    "    for row in target_points_ng_0.iterrows():\n",
    "        index_initial = row[1]['index']\n",
    "        index_lst = []\n",
    "        probs = []\n",
    "\n",
    "        for time in time_steps:\n",
    "            target_points_ng_coreg_out = target_points_ng_coreg_out_lst[time]\n",
    "            target_points_ng_coreg_out_index = target_points_ng_coreg_out['index'].tolist()\n",
    "            try:\n",
    "                index_lst.append(target_points_ng_coreg_out_index.index(index_initial))\n",
    "            except ValueError:\n",
    "                index_lst.append(None)\n",
    "\n",
    "        for time in time_steps:\n",
    "            target_points_ng_prob = target_points_ng_prob_lst[time]\n",
    "            if index_lst[time] != None:\n",
    "                target_points_ng_prob_val = target_points_ng_prob.iloc[index_lst[time]]['prob']\n",
    "                probs.append(target_points_ng_prob_val)\n",
    "            else:\n",
    "                probs.append(np.nan)\n",
    "\n",
    "        prob_lst.append(probs)\n",
    "\n",
    "    prob_arr = np.array(prob_lst)\n",
    "    ma_list = [f\"{i}_ma\" for i in range(31)]\n",
    "    prob_df = pd.DataFrame(prob_arr, columns=ma_list)\n",
    "\n",
    "    target_points_ng_probs = pd.DataFrame()\n",
    "    target_points_ng_probs['lon'] = target_points_ng_0['lon']\n",
    "    target_points_ng_probs['lat'] = target_points_ng_0['lat']\n",
    "    target_points_ng_probs = pd.concat([target_points_ng_probs, prob_df], axis=1)\n",
    "\n",
    "    probs = []\n",
    "\n",
    "    for row in target_points_ng_probs.iterrows():\n",
    "        prob = row[1][2:]\n",
    "        prob = prob.dropna()\n",
    "        prob_filtered = prob[prob > 0.5]\n",
    "        if not prob_filtered.empty:\n",
    "            prob_result = prob_filtered.median()\n",
    "        else:\n",
    "            prob_result = prob.median()\n",
    "        probs.append(prob_result)\n",
    "\n",
    "    probs = [0 if pd.isna(x) else x for x in probs]\n",
    "    target_points_ng_probs['prob'] = probs\n",
    "    target_points_ng_probs.to_csv(ml_output_dir + 'target_points_ng_prob.csv', index=False)\n",
    "else:\n",
    "    target_points_ng_probs = pd.read_csv(target_points_ng_probs_file, index_col=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fcb30fd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.LambertAzimuthalEqualArea(150, 0)\n",
    "\n",
    "# call the PlotTopologies object\n",
    "gplot = gplately.PlotTopologies(model, coastlines, continents, cob, time=0)\n",
    "\n",
    "if plate_motion_model == 'muller2016':\n",
    "    agegrid_file = agegrid_dir + f'Muller_etal_2016_AREPS_v1.17_AgeGrid-0.nc'\n",
    "elif plate_motion_model == 'muller2019':\n",
    "    agegrid_file = agegrid_dir + f'Muller_etal_2019_Tectonics_v2.0_AgeGrid-0.nc'\n",
    "\n",
    "agegrid = gplately.grids.read_netcdf_grid(agegrid_file)\n",
    "\n",
    "plot_x = target_points_ng_probs['lon']\n",
    "plot_y = target_points_ng_probs['lat']\n",
    "\n",
    "ng_extent_bounds = target_polygon_gdf.bounds\n",
    "ng_extent = [ng_extent_bounds.loc[0]['minx'], ng_extent_bounds.loc[0]['maxx'],\n",
    "             ng_extent_bounds.loc[0]['miny'], ng_extent_bounds.loc[0]['maxy']]\n",
    "\n",
    "# single colour bar\n",
    "fig = plt.figure(figsize=(9, 12))\n",
    "gs = GridSpec(2, 2, hspace=-0.68, wspace=0.1, height_ratios=[1, 0.01])\n",
    "ax = fig.add_subplot(gs[0, :], projection=proj)\n",
    "\n",
    "set_ax(ax, ng_extent, 10, 5, stock_img=False, order=9)\n",
    "\n",
    "im = gplot.plot_grid(ax, agegrid.data, cmap=ccm.lapaz_r, vmin=0, vmax=153, alpha=0.5, zorder=1)\n",
    "\n",
    "gplot.plot_continents(ax, edgecolor='gray', facecolor='tan', zorder=2)\n",
    "gplot.plot_ridges(ax, color='red', alpha=0.5, zorder=3)\n",
    "gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=False, regrid_shape=10, alpha=0.2, zorder=4)\n",
    "\n",
    "sc0 = ax.scatter(\n",
    "    plot_x,\n",
    "    plot_y,\n",
    "    transform=ccrs.PlateCarree(),\n",
    "    marker='.',\n",
    "    c=target_points_ng_probs['prob'],\n",
    "    s=10,\n",
    "    cmap=ccm.hawaii_r,\n",
    "    zorder=5\n",
    ")\n",
    "\n",
    "sc1 = ax.scatter(\n",
    "    min_occ_data['lon'],\n",
    "    min_occ_data['lat'],\n",
    "    transform=ccrs.PlateCarree(),\n",
    "    marker='o',\n",
    "    facecolor='yellow',\n",
    "    edgecolor='black',\n",
    "    s=20 * min_occ_data['weight'],\n",
    "    alpha=0.7,\n",
    "    zorder=6\n",
    ")\n",
    "\n",
    "gplot.plot_trenches(ax, color='k', alpha=0.3, zorder=6)\n",
    "gplot.plot_subduction_teeth(ax, spacing=0.03, color='k', alpha=0.3, zorder=7)\n",
    "\n",
    "cax1 = fig.add_subplot(gs[1, 0])\n",
    "cax2 = fig.add_subplot(gs[1, 1])\n",
    "\n",
    "fig.colorbar(sc0, cax=cax2, orientation='horizontal', label='Mineralisation Probability')\n",
    "fig.colorbar(im, cax=cax1, orientation='horizontal', label='Seafloor age (Ma)', extend='max')\n",
    "\n",
    "# Define custom legend handles\n",
    "custom_handles = [\n",
    "    Patch(facecolor='tan', edgecolor='gray', label='Continental Crust'),  # Custom handle for the filled polygon\n",
    "    Line2D([0], [0], color='red', lw=2, label='Mid-Ocean Ridge'),  # Custom handle for the line (ridge)\n",
    "    Line2D([0], [0], marker='o', markerfacecolor='yellow', markeredgecolor='black', markersize=5, linestyle='None', label='Mineral Occurrence')  # Custom handle for mineral occurrences\n",
    "]\n",
    "\n",
    "# Add the custom legend to the plot\n",
    "legend = ax.legend(handles=custom_handles, loc='lower left')\n",
    "\n",
    "ax.set_title(f'Porphyry Mineralisation Probability')\n",
    "\n",
    "# plt.savefig(\n",
    "# f'./figures/muller2019/target_points.png',\n",
    "# bbox_inches='tight',\n",
    "# pad_inches=0.1,\n",
    "# dpi=150\n",
    "# )\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "078a5c55",
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.LambertAzimuthalEqualArea(150, 0)\n",
    "\n",
    "# call the PlotTopologies object\n",
    "gplot = gplately.PlotTopologies(model, coastlines, continents, cob, time=0)\n",
    "\n",
    "if plate_motion_model == 'muller2016':\n",
    "    agegrid_file = agegrid_dir + 'Muller_etal_2016_AREPS_v1.17_AgeGrid-0.nc'\n",
    "elif plate_motion_model == 'muller2019':\n",
    "    agegrid_file = agegrid_dir + 'Muller_etal_2019_Tectonics_v2.0_AgeGrid-0.nc'\n",
    "\n",
    "agegrid = gplately.grids.read_netcdf_grid(agegrid_file)\n",
    "\n",
    "probabilities = []\n",
    "count = 0\n",
    "\n",
    "for mask in mask_coords_ng_0['include']:\n",
    "    if mask:\n",
    "        probabilities.append(target_points_ng_probs['prob'][count])\n",
    "        count += 1\n",
    "    else:\n",
    "        probabilities.append(np.nan)\n",
    "\n",
    "nx = mask_coords_ng_0['lon'].nunique()\n",
    "ny = mask_coords_ng_0['lat'].nunique()\n",
    "\n",
    "x_min = mask_coords_ng_0['lon'].min()\n",
    "x_max = mask_coords_ng_0['lon'].max()\n",
    "y_min = mask_coords_ng_0['lat'].min()\n",
    "y_max = mask_coords_ng_0['lat'].max()\n",
    "\n",
    "probabilities_2d = np.reshape(probabilities, (ny, nx))\n",
    "probabilities_2d_ud = np.flipud(np.reshape(probabilities, (ny, nx)))\n",
    "\n",
    "# dual colour bars\n",
    "fig = plt.figure(figsize=(9, 12))\n",
    "gs = GridSpec(2, 2, hspace=-0.68, wspace=0.1, height_ratios=[1, 0.01])\n",
    "ax = fig.add_subplot(gs[0, :], projection=proj)\n",
    "\n",
    "set_ax(ax, ng_extent, 5, 3, stock_img=False, order=9)\n",
    "\n",
    "im0 = gplot.plot_grid(ax, agegrid.data, cmap=ccm.lapaz_r, vmin=0, vmax=153, alpha=0.5, zorder=1)\n",
    "\n",
    "gplot.plot_continents(ax, edgecolor='gray', facecolor='tan', zorder=2)\n",
    "gplot.plot_ridges(ax, color='red', alpha=0.5, zorder=3)\n",
    "gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=False, regrid_shape=10, alpha=0.2, zorder=4)\n",
    "\n",
    "im1 = plt.imshow(\n",
    "    probabilities_2d,\n",
    "    transform=ccrs.PlateCarree(),\n",
    "    origin='lower',\n",
    "    cmap=ccm.hawaii_r,\n",
    "    interpolation='nearest',\n",
    "    extent=(x_min, x_max, y_min, y_max),\n",
    "    zorder=5\n",
    ")\n",
    "\n",
    "sc1 = ax.scatter(\n",
    "    min_occ_data['lon'],\n",
    "    min_occ_data['lat'],\n",
    "    transform=ccrs.PlateCarree(),\n",
    "    marker='o',\n",
    "    facecolor='yellow',\n",
    "    edgecolor='black',\n",
    "    s=20 * min_occ_data['weight'],\n",
    "    alpha=0.7,\n",
    "    zorder=6\n",
    ")\n",
    "\n",
    "gplot.plot_trenches(ax, color='k', alpha=0.5, zorder=7)\n",
    "gplot.plot_subduction_teeth(ax, spacing=0.01, color='k', alpha=0.5, zorder=8)\n",
    "\n",
    "cax1 = fig.add_subplot(gs[1, 0])\n",
    "cax2 = fig.add_subplot(gs[1, 1])\n",
    "\n",
    "fig.colorbar(im1, cax=cax2, orientation='horizontal', label='Mineralisation probability')\n",
    "fig.colorbar(im0, cax=cax1, orientation='horizontal', label='Seafloor age (Ma)', extend='max')\n",
    "\n",
    "# Define custom legend handles\n",
    "custom_handles = [\n",
    "    Patch(facecolor='tan', edgecolor='gray', label='Continental Crust'),  # Custom handle for the filled polygon\n",
    "    Line2D([0], [0], color='red', lw=2, label='Mid-Ocean Ridge'),  # Custom handle for the line (ridge)\n",
    "    Line2D([0], [0], marker='o', markerfacecolor='yellow', markeredgecolor='black', markersize=5, linestyle='None', label='Mineral Occurrence')  # Custom handle for mineral occurrences\n",
    "]\n",
    "\n",
    "# Add the custom legend to the plot\n",
    "legend = ax.legend(handles=custom_handles, loc='lower left')\n",
    "\n",
    "ax.set_title(f'Porphyry Mineralisation Probability')\n",
    "\n",
    "# plt.savefig(\n",
    "# f'./figures/muller2019/prob-pa/prob_ng.png',\n",
    "# bbox_inches='tight',\n",
    "# pad_inches=0.1,\n",
    "# dpi=150\n",
    "# )\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "654f21e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "target_points_ng_prob_map_file = ml_output_dir + 'target_points_ng_prob.tif'\n",
    "\n",
    "if not(os.path.isfile(target_points_ng_prob_map_file)):\n",
    "    # export the map to a GeoTIFF file\n",
    "    xmin, ymin, xmax, ymax = [x_min, y_min, x_max, y_max]\n",
    "    geotransform = (xmin, 0.2, 0, ymax, 0, -0.2)\n",
    "    driver = gdal.GetDriverByName('GTiff')\n",
    "    dataset = driver.Create(target_points_ng_prob_map_file, nx, ny, 1, gdal.GDT_Float32)\n",
    "    dataset.SetGeoTransform(geotransform)\n",
    "    srs = osr.SpatialReference()\n",
    "    srs.ImportFromEPSG(4326)\n",
    "    dataset.SetProjection(srs.ExportToWkt())\n",
    "    dataset.GetRasterBand(1).WriteArray(probabilities_2d_ud)\n",
    "    dataset.FlushCache()\n",
    "    dataset = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f87a32dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the geotiff file (mineral prospectivity map) using rioxarray\n",
    "prospectivity_data = rxr.open_rasterio(ml_output_dir + 'target_points_ng_prob.tif', masked=True).squeeze()\n",
    "\n",
    "# Load the shapefile (mineral occurrences)\n",
    "mineral_occurrences = gpd.read_file('../GIS/min_occ_ng.shp')\n",
    "\n",
    "# Extract probability values at known mineral occurrence locations\n",
    "occurrence_probabilities = []\n",
    "for geom in mineral_occurrences.geometry:\n",
    "    x, y = geom.x, geom.y\n",
    "    # Extract probability value at this point\n",
    "    prob_value = prospectivity_data.sel(x=x, y=y, method='nearest').item()\n",
    "    occurrence_probabilities.append(prob_value)\n",
    "\n",
    "# Convert to numpy array for easier manipulation\n",
    "occurrence_probabilities = np.array(occurrence_probabilities)\n",
    "occurrence_probabilities = occurrence_probabilities[~np.isnan(occurrence_probabilities)]\n",
    "\n",
    "# Calculate total study area and the number of cells\n",
    "total_area = np.count_nonzero(~np.isnan(prospectivity_data))\n",
    "total_occurrences = len(occurrence_probabilities)\n",
    "\n",
    "# Define the probability thresholds\n",
    "thresholds = np.linspace(0, 1, 100)\n",
    "\n",
    "# Initialize lists to store results\n",
    "occurrence_percentages = []\n",
    "area_percentages = []\n",
    "\n",
    "# Calculate percentages for each threshold\n",
    "for threshold in thresholds:\n",
    "    occurrence_percentage = np.sum(occurrence_probabilities >= threshold) / total_occurrences * 100\n",
    "    area_percentage = np.sum(prospectivity_data >= threshold) / total_area * 100\n",
    "    occurrence_percentages.append(occurrence_percentage)\n",
    "    area_percentages.append(area_percentage)\n",
    "\n",
    "# Convert lists to numpy arrays for plotting\n",
    "occurrence_percentages = np.array(occurrence_percentages)\n",
    "area_percentages = np.array(area_percentages)\n",
    "\n",
    "# Set the font sizes\n",
    "SMALL_SIZE = 12\n",
    "MEDIUM_SIZE = 14\n",
    "LARGE_SIZE = 16\n",
    "\n",
    "# Create the P-A plot\n",
    "fig, ax1 = plt.subplots(figsize=(10, 6))\n",
    "\n",
    "# Plot the occurrence percentages\n",
    "ax1.plot(thresholds, occurrence_percentages, 'b-', linewidth=2)\n",
    "ax1.set_xlabel('Probability', fontsize=MEDIUM_SIZE)\n",
    "ax1.set_ylabel('Percentage of Known Mineral Occurrences', color='b', fontsize=MEDIUM_SIZE)\n",
    "ax1.set_ylim(0, 100)\n",
    "ax1.set_xlim(0, 1)\n",
    "ax1.grid(True)\n",
    "\n",
    "# Increase tick label sizes\n",
    "ax1.tick_params(axis='both', labelsize=SMALL_SIZE)\n",
    "\n",
    "# Create a second y-axis to plot the area percentages\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(thresholds, area_percentages, 'r-', linewidth=2)\n",
    "ax2.set_ylabel('Percentage of Study Area', color='r', fontsize=MEDIUM_SIZE)\n",
    "ax2.set_ylim(100, 0)  # Inverted y-axis\n",
    "ax2.tick_params(axis='y', labelsize=SMALL_SIZE)\n",
    "\n",
    "# Add horizontal and vertical lines\n",
    "intersection_x = 0.79\n",
    "intersection_y1 = 60\n",
    "intersection_y2 = 40\n",
    "h_line = ax1.hlines(y=intersection_y1, xmin=0, xmax=1, color='g', linestyles=':')\n",
    "v_line = ax1.vlines(x=intersection_x, ymin=0, ymax=100, color='g', linestyles=':')\n",
    "\n",
    "# Add a text box to explain the intersection point with larger font\n",
    "ax1.annotate(f'At probability {intersection_x}:\\n{intersection_y1}% of occurrences\\n{intersection_y2}% of area',\n",
    "             xy=(0.72, 0.82), xycoords='axes fraction',\n",
    "             ha='center', va='top',\n",
    "             fontsize=SMALL_SIZE,\n",
    "             bbox=dict(boxstyle='round', fc='white', ec='gray', alpha=0.8))\n",
    "\n",
    "plt.tight_layout()  # Adjust layout to prevent clipping of labels\n",
    "\n",
    "# Adjust the plot margins to make room for annotations\n",
    "plt.subplots_adjust(left=0.15, right=0.85)\n",
    "\n",
    "# plt.savefig(\n",
    "# f'./figures/muller2019/prob-pa/pa_ng.png',\n",
    "# bbox_inches='tight',\n",
    "# pad_inches=0.1,\n",
    "# dpi=150\n",
    "# )\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5cd4ab9",
   "metadata": {},
   "source": [
    "### Coregistration of Mineral Occurrences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "290e4b7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "min_occ_prob_dir = parameters['min_occ_prob_dir']\n",
    "min_occ_prob_files_lst = []\n",
    "min_occ_prob_tran_files_lst = []\n",
    "\n",
    "for index in min_occ_data['index']:\n",
    "    min_occ_prob_files_lst.append(min_occ_prob_dir + f'min_occ_features_{index}.csv')\n",
    "    min_occ_prob_tran_files_lst.append(min_occ_prob_dir + f'min_occ_features_tran_{index}.csv')\n",
    "\n",
    "coregistration_point(\n",
    "    min_occ_data,\n",
    "    conv_dir,\n",
    "    conv_prefix,\n",
    "    conv_ext,\n",
    "    min_occ_prob_dir,\n",
    "    file_prefix='min_occ_features',\n",
    "    time_steps=time_steps,\n",
    "    search_radius=3,\n",
    "    plate_motion_model='muller2019'\n",
    ")\n",
    "\n",
    "for min_occ_prob_file, min_occ_prob_tran_file in zip(min_occ_prob_files_lst, min_occ_prob_tran_files_lst):\n",
    "    if not os.path.isfile(min_occ_prob_tran_file):\n",
    "        min_occ_prob = pd.read_csv(min_occ_prob_file, index_col=False)\n",
    "        min_occ_prob_tran = min_occ_prob.copy()\n",
    "        probs = []\n",
    "\n",
    "        for i, row in min_occ_prob.iterrows():\n",
    "            row_features = row[features_list]\n",
    "                \n",
    "            if row_features.isnull().values.any():\n",
    "                probs.append(np.nan)\n",
    "            else:\n",
    "                try:\n",
    "                    row_features = st_scaler.transform(row_features.values.reshape(1, -1))\n",
    "                except:\n",
    "                    # load the model\n",
    "                    with open(ml_input_dir + 'st_scaler.pkl', 'rb') as f:\n",
    "                        st_scaler = pickle.load(f)\n",
    "                    row_features = st_scaler.transform(row_features.values.reshape(1, -1))\n",
    "\n",
    "                min_occ_prob_tran.loc[min_occ_prob_tran['age'] == i, features_list] = row_features[0].tolist()\n",
    "                prob = model_rf.predict_proba(row_features)[0, 1]\n",
    "                probs.append(prob)\n",
    "                \n",
    "        mm_scaler2 = MinMaxScaler()\n",
    "        probs_scaled = mm_scaler2.fit_transform(np.array(probs).reshape(-1, 1))\n",
    "        min_occ_prob['prob'] = probs_scaled\n",
    "        min_occ_prob_tran['prob'] = probs_scaled\n",
    "        min_occ_prob.to_csv(min_occ_prob_file, index=False)\n",
    "        min_occ_prob_tran.to_csv(min_occ_prob_tran_file, index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ca01929",
   "metadata": {},
   "source": [
    "### Probability Changes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8aea302",
   "metadata": {},
   "outputs": [],
   "source": [
    "@interact\n",
    "def show_map(file1=min_occ_prob_files_lst, file2=min_occ_prob_tran_files_lst, feature=features_list):\n",
    "    df1 = pd.read_csv(file1, index_col=False)\n",
    "    df2 = pd.read_csv(file2, index_col=False)\n",
    "\n",
    "    fig = plt.figure(figsize=(12, 4))\n",
    "    ax2 = fig.add_subplot(121, xlim=[df1['age'].max(), 0])\n",
    "    ax1 = ax2.twinx()\n",
    "\n",
    "    ax2.plot(df1['age'], df1['prob'], c='red')\n",
    "    \n",
    "    index1 = min_occ_prob_files_lst.index(file1)\n",
    "    age1 = min_occ_data.iloc[index1]['age']\n",
    "    ax2.vlines(x=age1, ymin=0, ymax=1, color='k', linestyles=':')\n",
    "    \n",
    "    ax2.set_ylim(0, 1)\n",
    "    ax1.plot(df1['age'], df1[feature], c='blue')\n",
    "    \n",
    "    ax4 = fig.add_subplot(122, xlim=[df2['age'].max(), 0])\n",
    "    ax3 = ax4.twinx()\n",
    "\n",
    "    ax4.plot(df2['age'], df2['prob'], c='red')\n",
    "    \n",
    "    index2 = min_occ_prob_tran_files_lst.index(file2)\n",
    "    age2 = min_occ_data.iloc[index2]['age']\n",
    "    ax4.vlines(x=age2, ymin=0, ymax=1, color='k', linestyles=':')\n",
    "    \n",
    "    ax4.set_ylim(0, 1)\n",
    "    ax3.plot(df2['age'], df2[feature], c='blue')\n",
    "\n",
    "    ax2.set_ylabel('Probability')\n",
    "    ax1.set_ylabel(feature + ' (Actual)')\n",
    "    ax4.set_ylabel('Probability')\n",
    "    ax3.set_ylabel(feature + ' (Standardised)')\n",
    "    \n",
    "    fig.tight_layout()\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6639d647",
   "metadata": {},
   "source": [
    "#### Smoothed Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4b26913",
   "metadata": {},
   "outputs": [],
   "source": [
    "# smoothing function for the plots\n",
    "def smooth_data(x, y, points=200):\n",
    "    x_new = np.linspace(x.min(), x.max(), points)\n",
    "    spl = make_interp_spline(x, y, k=2)  # b-spline\n",
    "    y_smooth = spl(x_new)\n",
    "    return x_new, y_smooth\n",
    "\n",
    "@interact\n",
    "def show_map(file1=min_occ_prob_files_lst, feature=features_list):\n",
    "    df1 = pd.read_csv(file1, index_col=False)\n",
    "    df1 = df1.replace([np.inf, -np.inf], np.nan).dropna()\n",
    "\n",
    "    fig = plt.figure(figsize=(6, 4))\n",
    "    ax1 = fig.add_subplot(111, xlim=[df1['age'].max(), 0])\n",
    "    ax2 = ax1.twinx()\n",
    "    \n",
    "    ax1.set_facecolor('whitesmoke')\n",
    "\n",
    "    # smooth and plot probability\n",
    "    age_smooth, prob_smooth = smooth_data(df1['age'], df1['prob'])\n",
    "    ln1 = ax1.plot(age_smooth, prob_smooth, c='orangered')\n",
    "\n",
    "    index1 = min_occ_prob_files_lst.index(file1)\n",
    "    age1 = min_occ_data.iloc[index1]['age']\n",
    "    ln2 = ax1.vlines(x=age1, ymin=0, ymax=1, color='k', linestyles=':', label='Age of Formation')\n",
    "\n",
    "    ax1.set_ylim(0, 1)\n",
    "\n",
    "    # smooth and plot feature\n",
    "    age_smooth, feature_smooth = smooth_data(df1['age'], df1[feature])\n",
    "    ln3 = ax2.plot(age_smooth, feature_smooth, c='royalblue')\n",
    "    \n",
    "    # Define custom legend handles\n",
    "    custom_handles = [\n",
    "        Line2D([0], [0], color='k', linestyle=':', label='Age of formation')  # Custom handle for the line (age of formation)\n",
    "    ]\n",
    "\n",
    "    # Add the custom legend to the plot\n",
    "    legend = ax2.legend(handles=custom_handles, loc='lower left')\n",
    "\n",
    "    ax1.set_xlabel('Age (Ma)')\n",
    "    ax1.set_ylabel('Mineralisation probability', color='orangered')\n",
    "    ax2.set_ylabel('Orthogonal component of the\\nrelative motion vector (cm/yr)', color='royalblue')\n",
    "        \n",
    "    fig.tight_layout()\n",
    "    \n",
    "#     plt.savefig(\n",
    "#         f'./figures/muller2019/panguna/conv_ortho_cm_yr.png',\n",
    "#         bbox_inches='tight',\n",
    "#         pad_inches=0.1,\n",
    "#         dpi=150\n",
    "#     )\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d8a7bcd",
   "metadata": {},
   "source": [
    "### Traceplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "731c496d",
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.LambertAzimuthalEqualArea(150, 0)\n",
    "\n",
    "@interact\n",
    "def show_map(file=min_occ_prob_files_lst, time=time_steps):\n",
    "    df = pd.read_csv(file, index_col=False)\n",
    "    lons = df.loc[df['age'] >= time]['lon'].tolist()\n",
    "    lats = df.loc[df['age'] >= time]['lat'].tolist()\n",
    "    bm = df.loc[df['age'] >= time]['before_mineralisation'].tolist()\n",
    "    val = df.loc[df['age'] >= time]['valid'].tolist()\n",
    "    \n",
    "    lons_inval = []\n",
    "    lats_inval = []\n",
    "    lons_bm = []\n",
    "    lats_bm = []\n",
    "    lons_am = []\n",
    "    lats_am = []\n",
    "    \n",
    "    # colour of the last point\n",
    "    if not val[0]:\n",
    "        last_point = 'invalid'\n",
    "    elif bm[0]:\n",
    "        last_point = 'before_mineralisation'\n",
    "    else:\n",
    "        last_point = 'after_mineralisation'\n",
    "    \n",
    "    for index in range(1, len(lons)):\n",
    "        if not val[index]:\n",
    "            lons_inval.append(lons[index])\n",
    "            lats_inval.append(lats[index])\n",
    "        elif bm[index]:\n",
    "            lons_bm.append(lons[index])\n",
    "            lats_bm.append(lats[index])\n",
    "        else:\n",
    "            lons_am.append(lons[index])\n",
    "            lats_am.append(lats[index])\n",
    "        \n",
    "    # call the PlotTopologies object\n",
    "    gplot = gplately.PlotTopologies(model, coastlines, continents, cob, time=time)\n",
    "\n",
    "    if plate_motion_model == 'muller2016':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2016_AREPS_v1.17_AgeGrid-{time}.nc'\n",
    "    elif plate_motion_model == 'muller2019':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2019_Tectonics_v2.0_AgeGrid-{time}.nc'\n",
    "        \n",
    "    agegrid = gplately.grids.read_netcdf_grid(agegrid_file)\n",
    "    \n",
    "    mask_coords = mask_coords_lst[time_steps.index(time)]\n",
    "    \n",
    "    probabilities = []\n",
    "    count = 0\n",
    "            \n",
    "    for mask in mask_coords['include']:\n",
    "        if mask:\n",
    "            probabilities.append(target_points_prob_lst[time_steps.index(time)]['prob'][count])\n",
    "            count += 1\n",
    "        else:\n",
    "            probabilities.append(np.nan)\n",
    "    \n",
    "    nx = mask_coords_lst[time_steps.index(time)]['lon'].nunique()\n",
    "    ny = mask_coords_lst[time_steps.index(time)]['lat'].nunique()\n",
    "    \n",
    "    x_min = mask_coords_lst[time_steps.index(time)]['lon'].min()\n",
    "    x_max = mask_coords_lst[time_steps.index(time)]['lon'].max()\n",
    "    y_min = mask_coords_lst[time_steps.index(time)]['lat'].min()\n",
    "    y_max = mask_coords_lst[time_steps.index(time)]['lat'].max()\n",
    "    \n",
    "    probabilities_2d = np.reshape(probabilities, (ny, nx))\n",
    "    probabilities_2d_ud = np.flipud(np.reshape(probabilities, (ny, nx)))\n",
    "    \n",
    "    # single colour bar\n",
    "    fig = plt.figure(figsize=(6, 8))\n",
    "    gs = GridSpec(2, 2, hspace=-0.4, wspace=0.1, height_ratios=[1, 0.02])\n",
    "    ax = fig.add_subplot(gs[0, :], projection=proj)\n",
    "    \n",
    "    set_ax(ax, target_extent_anim, 15, 15, stock_img=False, order=10)\n",
    "\n",
    "    im0 = gplot.plot_grid(ax, agegrid.data, cmap=ccm.lapaz_r, vmin=0, vmax=153, alpha=0.5, zorder=1)\n",
    "    \n",
    "    gplot.plot_continents(ax, edgecolor='gray', facecolor='tan', zorder=2)\n",
    "    gplot.plot_ridges_and_transforms(ax, color='red', alpha=0.5, zorder=3)\n",
    "    gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=False, regrid_shape=10, alpha=0.2, zorder=4)\n",
    "    \n",
    "    im1 = plt.imshow(\n",
    "        probabilities_2d,\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        origin='lower',\n",
    "        cmap=ccm.hawaii_r,\n",
    "        interpolation='nearest',\n",
    "        extent=(x_min, x_max, y_min, y_max),\n",
    "        alpha=0.7,\n",
    "        zorder=5\n",
    "    )\n",
    "    \n",
    "    if last_point == 'invalid':\n",
    "        sc = ax.scatter(lons[0], lats[0], transform=ccrs.PlateCarree(), marker='*', facecolor='gray', s=20, zorder=6)\n",
    "    elif last_point == 'before_mineralisation':\n",
    "        sc = ax.scatter(lons[0], lats[0], transform=ccrs.PlateCarree(), marker='*', facecolor='royalblue', s=20, zorder=6)\n",
    "    else:\n",
    "        sc = ax.scatter(lons[0], lats[0], transform=ccrs.PlateCarree(), marker='*', facecolor='lime', s=20, zorder=6)\n",
    "    \n",
    "    sc = ax.scatter(lons_bm, lats_bm, transform=ccrs.PlateCarree(), marker='.', facecolor='royalblue', s=20, zorder=7)\n",
    "    sc = ax.scatter(lons_am, lats_am, transform=ccrs.PlateCarree(), marker='.', facecolor='lime', s=20, zorder=7)\n",
    "    sc = ax.scatter(lons_inval, lats_inval, transform=ccrs.PlateCarree(), marker='.', facecolor='gray', s=20, zorder=7)\n",
    "    \n",
    "    gplot.plot_trenches(ax, color='k', alpha=0.3, zorder=8)\n",
    "    gplot.plot_subduction_teeth(ax, spacing=0.03, color='k', alpha=0.3, zorder=9)\n",
    "    \n",
    "    cax1 = fig.add_subplot(gs[1, 0])\n",
    "    cax2 = fig.add_subplot(gs[1, 1])\n",
    "    \n",
    "    fig.colorbar(im1, cax=cax2, orientation='horizontal', label='Mineralisation Probability')\n",
    "    fig.colorbar(im0, cax=cax1, orientation='horizontal', label='Seafloor age (Ma)', extend='max')\n",
    "\n",
    "    # Define custom legend handles\n",
    "    custom_handles = [\n",
    "        Patch(facecolor='tan', edgecolor='gray', label='Continental Crust'),  # Custom handle for the filled polygon\n",
    "        Line2D([0], [0], color='red', lw=2, label='Mid-Ocean Ridge'),  # Custom handle for the line (ridge)\n",
    "        Line2D([0], [0], marker='.', markerfacecolor='royalblue', markeredgecolor='none', markersize=10, linestyle='None', label='Pre-Mineralisation'),\n",
    "        Line2D([0], [0], marker='.', markerfacecolor='lime', markeredgecolor='none', markersize=10, linestyle='None', label='Post-Mineralisation')\n",
    "    ]\n",
    "\n",
    "    # Add the custom legend to the plot\n",
    "    legend = ax.legend(handles=custom_handles, loc='lower left')\n",
    "    legend.set_zorder(11)\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "986bb5bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.LambertAzimuthalEqualArea(150, 0)\n",
    "\n",
    "@interact\n",
    "def show_map(file=min_occ_prob_files_lst, time=time_steps):\n",
    "    df = pd.read_csv(file, index_col=False)\n",
    "    lons = df.loc[df['age'] >= time]['lon'].tolist()\n",
    "    lats = df.loc[df['age'] >= time]['lat'].tolist()\n",
    "    bm = df.loc[df['age'] >= time]['before_mineralisation'].tolist()\n",
    "    val = df.loc[df['age'] >= time]['valid'].tolist()\n",
    "        \n",
    "    # colour of the last point\n",
    "    if not val[0]:\n",
    "        last_point = 'invalid'\n",
    "    elif bm[0]:\n",
    "        last_point = 'before_mineralisation'\n",
    "    else:\n",
    "        last_point = 'after_mineralisation'\n",
    "\n",
    "    # list of invalid lines\n",
    "    lons_inval = []\n",
    "    lats_inval = []\n",
    "    lons_inval_temp = []\n",
    "    lats_inval_temp = []\n",
    "\n",
    "    for index2 in range(len(lons)-1, -1, -1):\n",
    "        if not val[index2]:\n",
    "            lons_inval_temp.append(lons[index2])\n",
    "            lats_inval_temp.append(lats[index2])\n",
    "            if index2 == 0:\n",
    "                lons_inval.append(lons_inval_temp)\n",
    "                lats_inval.append(lats_inval_temp)\n",
    "        else:\n",
    "            if len(lons_inval_temp) != 0:\n",
    "                lons_inval_temp.append(lons[index2])\n",
    "                lats_inval_temp.append(lats[index2])\n",
    "                lons_inval.append(lons_inval_temp)\n",
    "                lats_inval.append(lats_inval_temp)\n",
    "                lons_inval_temp = []\n",
    "                lats_inval_temp = []\n",
    "            else:\n",
    "                continue\n",
    "\n",
    "    # list of lines created before mineralisation\n",
    "    lons_bm = []\n",
    "    lats_bm = []\n",
    "    lons_bm_temp = []\n",
    "    lats_bm_temp = []\n",
    "\n",
    "    for index2 in range(len(lons)-1, -1, -1):\n",
    "        if bm[index2] and val[index2]:\n",
    "            lons_bm_temp.append(lons[index2])\n",
    "            lats_bm_temp.append(lats[index2])\n",
    "            if index2 == 0:\n",
    "                lons_bm.append(lons_bm_temp)\n",
    "                lats_bm.append(lats_bm_temp)\n",
    "        else:\n",
    "            if len(lons_bm_temp) != 0:\n",
    "                lons_bm_temp.append(lons[index2])\n",
    "                lats_bm_temp.append(lats[index2])\n",
    "                lons_bm.append(lons_bm_temp)\n",
    "                lats_bm.append(lats_bm_temp)\n",
    "                lons_bm_temp = []\n",
    "                lats_bm_temp = []\n",
    "            else:\n",
    "                continue\n",
    "\n",
    "    # list of lines created after mineralisation\n",
    "    lons_am = []\n",
    "    lats_am = []\n",
    "    lons_am_temp = []\n",
    "    lats_am_temp = []\n",
    "\n",
    "    for index2 in range(len(lons)-1, -1, -1):\n",
    "        if not bm[index2] and val[index2]:\n",
    "            lons_am_temp.append(lons[index2])\n",
    "            lats_am_temp.append(lats[index2])\n",
    "            if index2 == 0:\n",
    "                lons_am.append(lons_am_temp)\n",
    "                lats_am.append(lats_am_temp)\n",
    "        else:\n",
    "            if len(lons_am_temp) != 0:\n",
    "                lons_am_temp.append(lons[index2])\n",
    "                lats_am_temp.append(lats[index2])\n",
    "                lons_am.append(lons_am_temp)\n",
    "                lats_am.append(lats_am_temp)\n",
    "                lons_am_temp = []\n",
    "                lats_am_temp = []\n",
    "            else:\n",
    "                continue\n",
    "\n",
    "    # call the PlotTopologies object\n",
    "    gplot = gplately.PlotTopologies(model, coastlines, continents, cob, time=time)\n",
    "    \n",
    "    if plate_motion_model == 'muller2016':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2016_AREPS_v1.17_AgeGrid-{time}.nc'\n",
    "    elif plate_motion_model == 'muller2019':\n",
    "        agegrid_file = agegrid_dir + f'Muller_etal_2019_Tectonics_v2.0_AgeGrid-{time}.nc'\n",
    "\n",
    "    agegrid = gplately.grids.read_netcdf_grid(agegrid_file)\n",
    "    \n",
    "    mask_coords = mask_coords_lst[time_steps.index(time)]\n",
    "    \n",
    "    probabilities = []\n",
    "    count = 0\n",
    "            \n",
    "    for mask in mask_coords['include']:\n",
    "        if mask:\n",
    "            probabilities.append(target_points_prob_lst[time_steps.index(time)]['prob'][count])\n",
    "            count += 1\n",
    "        else:\n",
    "            probabilities.append(np.nan)\n",
    "    \n",
    "    nx = mask_coords_lst[time_steps.index(time)]['lon'].nunique()\n",
    "    ny = mask_coords_lst[time_steps.index(time)]['lat'].nunique()\n",
    "    \n",
    "    x_min = mask_coords_lst[time_steps.index(time)]['lon'].min()\n",
    "    x_max = mask_coords_lst[time_steps.index(time)]['lon'].max()\n",
    "    y_min = mask_coords_lst[time_steps.index(time)]['lat'].min()\n",
    "    y_max = mask_coords_lst[time_steps.index(time)]['lat'].max()\n",
    "    \n",
    "    probabilities_2d = np.reshape(probabilities, (ny, nx))\n",
    "    probabilities_2d_ud = np.flipud(np.reshape(probabilities, (ny, nx)))\n",
    "    \n",
    "    # single colour bar\n",
    "    fig = plt.figure(figsize=(6, 8))\n",
    "    gs = GridSpec(2, 2, hspace=-0.4, wspace=0.1, height_ratios=[1, 0.02])\n",
    "    ax = fig.add_subplot(gs[0, :], projection=proj)\n",
    "    \n",
    "    set_ax(ax, target_extent_anim, 15, 15, stock_img=False, order=10)\n",
    "\n",
    "    im0 = gplot.plot_grid(ax, agegrid.data, cmap=ccm.lapaz_r, vmin=0, vmax=153, alpha=0.5, zorder=1)\n",
    "    \n",
    "    gplot.plot_continents(ax, edgecolor='gray', facecolor='tan', zorder=2)\n",
    "    gplot.plot_ridges_and_transforms(ax, color='red', alpha=0.5, zorder=3)\n",
    "    gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=False, regrid_shape=10, alpha=0.2, zorder=4)\n",
    "    \n",
    "    im1 = plt.imshow(\n",
    "        probabilities_2d,\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        origin='lower',\n",
    "        cmap=ccm.hawaii_r,\n",
    "        interpolation='nearest',\n",
    "        extent=(x_min, x_max, y_min, y_max),\n",
    "        alpha=0.7,\n",
    "        zorder=5\n",
    "    )\n",
    "        \n",
    "    for lons_, lats_ in zip(lons_bm, lats_bm):\n",
    "        sc = ax.plot(lons_, lats_, transform=ccrs.PlateCarree(), color='royalblue', zorder=6)\n",
    "        \n",
    "    for lons_, lats_ in zip(lons_am, lats_am):\n",
    "        sc = ax.plot(lons_, lats_, transform=ccrs.PlateCarree(), color='lime', zorder=6)\n",
    "        \n",
    "    for lons_, lats_ in zip(lons_inval, lats_inval):\n",
    "        sc = ax.plot(lons_, lats_, transform=ccrs.PlateCarree(), color='gray', zorder=6)\n",
    "        \n",
    "    if last_point == 'invalid':\n",
    "        sc = ax.scatter(lons[0], lats[0], transform=ccrs.PlateCarree(), marker='*', facecolor='gray', s=20, zorder=7)\n",
    "    elif last_point == 'before_mineralisation':\n",
    "        sc = ax.scatter(lons[0], lats[0], transform=ccrs.PlateCarree(), marker='*', facecolor='royalblue', s=20, zorder=7)\n",
    "    else:\n",
    "        sc = ax.scatter(lons[0], lats[0], transform=ccrs.PlateCarree(), marker='*', facecolor='lime', s=20, zorder=7)\n",
    "    \n",
    "    gplot.plot_trenches(ax, color='k', alpha=0.3, zorder=8)\n",
    "    gplot.plot_subduction_teeth(ax, spacing=0.03, color='k', alpha=0.3, zorder=9)\n",
    "    \n",
    "    cax1 = fig.add_subplot(gs[1, 0])\n",
    "    cax2 = fig.add_subplot(gs[1, 1])\n",
    "    \n",
    "    fig.colorbar(im1, cax=cax2, orientation='horizontal', label='Mineralisation Probability')\n",
    "    fig.colorbar(im0, cax=cax1, orientation='horizontal', label='Seafloor age (Ma)', extend='max')\n",
    "    \n",
    "    # Define custom legend handles\n",
    "    custom_handles = [\n",
    "        Patch(facecolor='tan', edgecolor='gray', label='Continental Crust'),  # Custom handle for the filled polygon\n",
    "        Line2D([0], [0], color='red', lw=2, label='Mid-Ocean Ridge'),  # Custom handle for the line (ridge)\n",
    "        Line2D([0], [0], color='royalblue', lw=2, label='Pre-Mineralisation'),\n",
    "        Line2D([0], [0], color='lime', lw=2, label='Post-Mineralisation')\n",
    "    ]\n",
    "\n",
    "    # Add the custom legend to the plot\n",
    "    legend = ax.legend(handles=custom_handles, loc='lower left')\n",
    "    legend.set_zorder(11)\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae456ba5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "png_jun23",
   "language": "python",
   "name": "png_jun23"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}