04 Hotspots to classes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This notebook classifies each HUC basin by *Gi** value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import geopandas as gpd\n",
    "import glob\n",
    "from tqdm import tqdm\n",
    "\n",
    "%matplotlib inline\n",
    "import statsmodels as sm\n",
    "from statsmodels.imputation import mice\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "element_list = [\n",
    "    \"ag\",\n",
    "    \"al\",\n",
    "    \"as\",\n",
    "    \"au\",\n",
    "    \"b\",\n",
    "    \"ba\",\n",
    "    \"be\",\n",
    "    \"bi\",\n",
    "    \"ca\",\n",
    "    \"cd\",\n",
    "    \"ce\",\n",
    "    \"co\",\n",
    "    \"cr\",\n",
    "    \"cs\",\n",
    "    \"cu\",\n",
    "    \"dy\",\n",
    "    \"er\",\n",
    "    \"eu\",\n",
    "    \"fe\",\n",
    "    \"ga\",\n",
    "    \"gd\",\n",
    "    \"ge\",\n",
    "    \"hf\",\n",
    "    \"hg\",\n",
    "    \"ho\",\n",
    "    \"k\",\n",
    "    \"la\",\n",
    "    \"li\",\n",
    "    \"lu\",\n",
    "    \"mg\",\n",
    "    \"mn\",\n",
    "    \"mo\",\n",
    "    \"na\",\n",
    "    \"nb\",\n",
    "    \"nd\",\n",
    "    \"ni\",\n",
    "    \"os\",\n",
    "    \"p\",\n",
    "    \"pb\",\n",
    "    \"pd\",\n",
    "    \"pr\",\n",
    "    \"pt\",\n",
    "    \"rb\",\n",
    "    \"re\",\n",
    "    \"rh\",\n",
    "    \"ru\",\n",
    "    \"sb\",\n",
    "    \"sc\",\n",
    "    \"se\",\n",
    "    \"sm\",\n",
    "    \"sn\",\n",
    "    \"sr\",\n",
    "    \"ta\",\n",
    "    \"tb\",\n",
    "    \"te\",\n",
    "    \"th\",\n",
    "    \"ti\",\n",
    "    \"tl\",\n",
    "    \"tm\",\n",
    "    \"u\",\n",
    "    \"v\",\n",
    "    \"w\",\n",
    "    \"y\",\n",
    "    \"yb\",\n",
    "    \"zn\",\n",
    "    \"zr\",\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All rock *Gi** values should be calculated in ArcMap and exported to an `xls` file for each element, the next cell reads in the excel files for rock hotspots, and reads in the sediment hotspots that were calculated in `03 Sediment hotspots.ipynb`. It combines the two for each element and classifies as `known`, `missed`, `potential`, `below`, and `background`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this cell does all elements for rock at once in a batch\n",
    "for element in element_list:\n",
    "    try:\n",
    "        huc12 = gpd.read_file(\n",
    "            r\"J:\\Geology\\WSGS\\GIS\\huc_12.shp\"\n",
    "        )  # reads in the hydraulic basin units and adds an area column\n",
    "        huc12 = huc12.to_crs({\"init\": \"epsg:3732\"})\n",
    "\n",
    "        frame = pd.read_excel(\n",
    "            \"J:/Geology/WSGS/Projects/Critical Minerals/hotspot analysis/rock hotspot joined with hucs/\"\n",
    "            + element\n",
    "            + \".xls\",\n",
    "            index_col=[0],\n",
    "            sheet=0,\n",
    "        )\n",
    "        frame = frame.drop(\n",
    "            columns=[\n",
    "                \"Join_Count\",\n",
    "                \"TARGET_FID\",\n",
    "                \"AREA\",\n",
    "                \"PERIMETER\",\n",
    "                \"WY_HU12_\",\n",
    "                \"WY_HU12_ID\",\n",
    "                \"HUC_8\",\n",
    "                \"HUC_10\",\n",
    "                \"ACRES\",\n",
    "                \"STATES\",\n",
    "                \"NCONTRB_A\",\n",
    "                \"HU_10_DS\",\n",
    "                \"HU_10_NAME\",\n",
    "                \"HU_10_MOD\",\n",
    "                \"HU_10_TYPE\",\n",
    "                \"HU_12_DS\",\n",
    "                \"HU_12_NAME\",\n",
    "                \"HU_12_MOD\",\n",
    "                \"HU_12_TYPE\",\n",
    "                \"HU_2_NAME\",\n",
    "                \"HU_4_NAME\",\n",
    "                \"HU_6_NAME\",\n",
    "                \"HU_8_NAME\",\n",
    "                \"Shape_Leng\",\n",
    "                \"SOURCE_ID\",\n",
    "                \"GiPValue\",\n",
    "                \"NNeighbors\",\n",
    "                \"Gi_Bin\",\n",
    "                \"Shape_Length\",\n",
    "                \"Shape_Area\",\n",
    "            ]\n",
    "        )\n",
    "        summary = frame.groupby(\"HUC_12\", as_index=False).sum()\n",
    "        middle = huc12.merge(summary.astype(object), on=\"HUC_12\", how=\"right\")\n",
    "        combined = huc12.merge(\n",
    "            middle[[\"HUC_12\", \"GiZScore\"]].astype(str), on=\"HUC_12\", how=\"left\"\n",
    "        )\n",
    "        combined.dropna(inplace=True)\n",
    "        plt.figure(figsize=(18, 12))\n",
    "\n",
    "        combined[\"zscore\"] = combined.GiZScore.astype(float).values\n",
    "        combined.to_file(\n",
    "            driver=\"ESRI Shapefile\",\n",
    "            filename=\"/results/\" + element + \"_rock_heatmap.shp\",\n",
    "        )\n",
    "\n",
    "        frame = pd.read_csv(\n",
    "            \"/sed hotspot joined with hucs/\" + element + \"_ppm.csv\"\n",
    "        )\n",
    "        summary = frame.groupby(\"HUC_12\", as_index=False).sum()\n",
    "        middle = huc12.merge(summary.astype(object), on=\"HUC_12\", how=\"right\")\n",
    "        sed_comb = huc12.merge(\n",
    "            middle[[\"HUC_12\", \"z_score\"]].astype(str), on=\"HUC_12\", how=\"left\"\n",
    "        )\n",
    "        crs = {\"init\": \"epsg:3732\"}\n",
    "        sed_comb.dropna(inplace=True)\n",
    "        sed_comb = sed_comb.to_crs(crs)\n",
    "        sed_comb[\"zscore\"] = sed_comb.z_score.astype(float).values\n",
    "\n",
    "        sed_comb.to_file(\n",
    "            driver=\"ESRI Shapefile\",\n",
    "            filename=\"/results/\" + element + \"_sed_heatmap.shp\",\n",
    "        )\n",
    "\n",
    "        missed = sed_comb[\n",
    "            (sed_comb[\"zscore\"] >= 2)\n",
    "            & sed_comb[\"HUC_12\"].isin(combined[combined[\"zscore\"] < 2].HUC_12)\n",
    "        ]\n",
    "        missed.drop([\"z_score\"], axis=1, inplace=True)\n",
    "        missed[\"class\"] = \"missed\"\n",
    "\n",
    "        known = combined[combined[\"zscore\"] >= 2]\n",
    "        known.drop([\"GiZScore\"], axis=1, inplace=True)\n",
    "        known[\"class\"] = \"known\"\n",
    "\n",
    "        potential = sed_comb[\n",
    "            (sed_comb[\"zscore\"] >= 2)\n",
    "            & (~sed_comb[\"HUC_12\"].isin(missed.HUC_12))\n",
    "            & (~sed_comb[\"HUC_12\"].isin(known.HUC_12))\n",
    "        ]\n",
    "        potential.drop([\"z_score\"], axis=1, inplace=True)\n",
    "        potential[\"class\"] = \"potential\"\n",
    "\n",
    "        below = sed_comb[\n",
    "            (sed_comb[\"zscore\"] <= -2)\n",
    "            & (\n",
    "                sed_comb[\"HUC_12\"].isin(\n",
    "                    combined[combined[\"zscore\"] <= -2].HUC_12\n",
    "                )\n",
    "            )\n",
    "        ]\n",
    "        below.drop([\"z_score\"], axis=1, inplace=True)\n",
    "        below[\"class\"] = \"below\"\n",
    "\n",
    "        bg = huc12[\n",
    "            (~huc12[\"HUC_12\"].isin(below.HUC_12))\n",
    "            & (~huc12[\"HUC_12\"].isin(potential.HUC_12))\n",
    "            & (~huc12[\"HUC_12\"].isin(known.HUC_12))\n",
    "            & (~huc12[\"HUC_12\"].isin(missed.HUC_12))\n",
    "        ]\n",
    "        bg[\"class\"] = \"background\"\n",
    "\n",
    "        classes = pd.concat([missed, known, potential, bg, below], sort=True)\n",
    "        classes.to_file(\n",
    "            driver=\"ESRI Shapefile\",\n",
    "            filename=\"/results/\" + element + \"_classes.shp\",\n",
    "        )\n",
    "    except:\n",
    "        print(element + \" was not successfully written to shapefile\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "element = \"co\"\n",
    "# this cell creates individual rock heatmaps\n",
    "huc12 = gpd.read_file(\n",
    "    r\"\\huc_12.shp\"\n",
    ")  # reads in the hydraulic basin units and adds an area column\n",
    "huc12 = huc12.to_crs({\"init\": \"epsg:3732\"})\n",
    "\n",
    "frame = pd.read_excel(\n",
    "    \"/rock hotspot joined with hucs/\"\n",
    "    + element\n",
    "    + \"_SpatialJoin_TableToExcel.xls\",\n",
    "    index_col=[0],\n",
    "    sheet_name=element + \"_SpatialJoin_TableToExcel\",\n",
    ")\n",
    "frame = frame.drop(\n",
    "    columns=[\n",
    "        \"Join_Count\",\n",
    "        \"TARGET_FID\",\n",
    "        \"AREA\",\n",
    "        \"PERIMETER\",\n",
    "        \"WY_HU12_\",\n",
    "        \"WY_HU12_ID\",\n",
    "        \"HUC_8\",\n",
    "        \"HUC_10\",\n",
    "        \"ACRES\",\n",
    "        \"STATES\",\n",
    "        \"NCONTRB_A\",\n",
    "        \"HU_10_DS\",\n",
    "        \"HU_10_NAME\",\n",
    "        \"HU_10_MOD\",\n",
    "        \"HU_10_TYPE\",\n",
    "        \"HU_12_DS\",\n",
    "        \"HU_12_NAME\",\n",
    "        \"HU_12_MOD\",\n",
    "        \"HU_12_TYPE\",\n",
    "        \"HU_2_NAME\",\n",
    "        \"HU_4_NAME\",\n",
    "        \"HU_6_NAME\",\n",
    "        \"HU_8_NAME\",\n",
    "        \"Shape_Leng\",\n",
    "        \"SOURCE_ID\",\n",
    "        \"GiPValue\",\n",
    "        \"NNeighbors\",\n",
    "        \"Gi_Bin\",\n",
    "        \"Shape_Length\",\n",
    "        \"Shape_Area\",\n",
    "    ]\n",
    ")\n",
    "summary = frame.groupby(\"HUC_12\", as_index=False).sum()\n",
    "middle = huc12.merge(summary.astype(object), on=\"HUC_12\", how=\"right\")\n",
    "combined = huc12.merge(\n",
    "    middle[[\"HUC_12\", \"GiZScore\"]].astype(str), on=\"HUC_12\", how=\"left\"\n",
    ")\n",
    "combined.dropna(inplace=True)\n",
    "plt.figure(figsize=(18, 12))\n",
    "combined.plot(\n",
    "    column=pd.to_numeric(combined[\"GiZScore\"]),\n",
    "    cmap=\"Greens\",\n",
    "    legend=True,\n",
    "    figsize=(20, 10),\n",
    ")\n",
    "plt.title(\"Rock Sample Gi* Scores\")\n",
    "\n",
    "combined[\"zscore\"] = combined.GiZScore.astype(float).values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this does the same as above, but for sediment samples\n",
    "frame = pd.read_csv(\"/sed hotspot joined with hucs/\" + element + \"_ppm.csv\")\n",
    "summary = frame.groupby(\"HUC_12\", as_index=False).sum()\n",
    "middle = huc12.merge(summary.astype(object), on=\"HUC_12\", how=\"right\")\n",
    "sed_comb = huc12.merge(\n",
    "    middle[[\"HUC_12\", \"z_score\"]].astype(str), on=\"HUC_12\", how=\"left\"\n",
    ")\n",
    "crs = {\"init\": \"epsg:3732\"}\n",
    "sed_comb.dropna(inplace=True)\n",
    "sed_comb = sed_comb.to_crs(crs)\n",
    "sed_comb.plot(\n",
    "    column=pd.to_numeric(sed_comb[\"z_score\"]),\n",
    "    cmap=\"Oranges\",\n",
    "    legend=True,\n",
    "    figsize=(20, 10),\n",
    ")\n",
    "plt.title(\"Sediment Sample Gi* Scores\")\n",
    "sed_comb[\"zscore\"] = sed_comb.z_score.astype(float).values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then create a map of each class and assign values to each HUC based on cumulative *Gi** values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "missed = sed_comb[(sed_comb['zscore']>=2) & sed_comb['HUC_12'].isin(combined[combined['zscore']<2].HUC_12)]\n",
    "missed[]\n",
    "missed.drop(['z_score'], axis=1, inplace=True)\n",
    "missed['class'] = 'missed'\n",
    "\n",
    "known = combined[combined['zscore']>=2]\n",
    "known.drop(['GiZScore'], axis=1, inplace=True)\n",
    "known['class'] = 'known'\n",
    "\n",
    "potential = sed_comb[(sed_comb['zscore']>=2) & (~sed_comb['HUC_12'].isin(missed.HUC_12)) & (~sed_comb['HUC_12'].isin(known.HUC_12))]\n",
    "potential.drop(['z_score'], axis=1, inplace=True)\n",
    "potential['class'] = 'potential'\n",
    "\n",
    "below = sed_comb[(sed_comb['zscore']<=-2) & (sed_comb['HUC_12'].isin(combined[combined['zscore']<=-2].HUC_12))]\n",
    "below.drop(['z_score'], axis=1, inplace=True)\n",
    "below['class'] = 'below'\n",
    "\n",
    "bg = huc12[(~huc12['HUC_12'].isin(below.HUC_12)) & (~huc12['HUC_12'].isin(potential.HUC_12)) & (~huc12['HUC_12'].isin(known.HUC_12)) & (~huc12['HUC_12'].isin(missed.HUC_12))]\n",
    "bg['class'] = 'background'\n",
    "classes = pd.concat([missed, known, potential, bg, below], sort=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "classes[classes[\"class\"] == \"missed\"].plot(\n",
    "    column=\"zscore\", cmap=\"viridis\", legend=True, figsize=(20, 10)\n",
    ")\n",
    "plt.title(\"Missed Areas Gi* Scores\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}