{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load and Visualize EFRU Contaminant data\n",
    "1. Load a CSV containing testing results\n",
    "1. Print out basic metrics about data\n",
    "1. Generate piechart for lead and ash/soot/char detection\n",
    "1. Generate boxplot for lead\n",
    "1. Generate boxplot for all contaminants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "\n",
    "# Initialize logger\n",
    "logging.basicConfig(level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s')\n",
    "logger = logging.getLogger(__name__)\n",
    "\t\n",
    "# UPDATE PATH TO TEST DATA AND DESIRED METALS\n",
    "\n",
    "csv_fpath = '../EFRU_Test_Results_20250823.csv'\n",
    "metals_to_analyze = ['Max Pb', 'As', 'Sb', 'Ba', 'Be', 'Cd', 'Cr', 'Co', 'Cu', 'Hg', 'Mo', 'Ni', 'Se', 'Ag', 'Tl', 'V', 'Zn']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data\n",
    "df = pd.read_csv(csv_fpath)\n",
    "\n",
    "# Verify required columns exist\n",
    "required_cols = ['Latitude', 'Longitude']\n",
    "metal_cols = [f\"{metal}_numeric\" for metal in metals_to_analyze]\n",
    "\n",
    "missing_cols = [col for col in required_cols + metal_cols if col not in df.columns]\n",
    "if missing_cols:\n",
    "    logger.warning(f\"Missing columns in data: {', '.join(missing_cols)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total samples:  201\n",
      "Metal concentration statistics (µg/ft^2):\n",
      "Max Pb: mean =  286.0, median =   94.9, std =  663.8, n_valid = 180\n",
      "As    : mean =   12.0, median =    4.5, std =   18.2, n_valid =  75\n",
      "Sb    : mean =   26.2, median =    5.1, std =   66.3, n_valid =  53\n",
      "Ba    : mean =  515.0, median =  261.1, std =  634.6, n_valid =  57\n",
      "Be    : mean =    1.3, median =    0.0, std =    2.4, n_valid =  55\n",
      "Cd    : mean =    2.5, median =    0.3, std =    6.0, n_valid =  65\n",
      "Cr    : mean =   50.6, median =   20.4, std =   67.3, n_valid =  65\n",
      "Co    : mean =   22.2, median =   10.0, std =   31.7, n_valid =  59\n",
      "Cu    : mean =  325.8, median =  148.6, std =  450.4, n_valid =  59\n",
      "Hg    : mean =   10.6, median =    0.1, std =   68.9, n_valid =  46\n",
      "Mo    : mean =    5.4, median =    0.0, std =   20.2, n_valid =  53\n",
      "Ni    : mean =   51.7, median =   26.9, std =   73.3, n_valid =  61\n",
      "Se    : mean =    6.5, median =    0.0, std =   16.8, n_valid =  54\n",
      "Ag    : mean =    0.8, median =    0.0, std =    2.3, n_valid =  59\n",
      "Tl    : mean =    0.5, median =    0.0, std =    3.2, n_valid =  52\n",
      "V     : mean =   84.5, median =   38.6, std =  123.3, n_valid =  58\n",
      "Zn    : mean = 5606.4, median =  923.5, std = 26767.7, n_valid =  59\n"
     ]
    }
   ],
   "source": [
    "# Print out statistics for each metal\n",
    "\n",
    "print(\"Total samples: \", len(df))\n",
    "print(\"Metal concentration statistics (µg/ft^2):\")\n",
    "for metal in metals_to_analyze:\n",
    "    metal_col = f\"{metal}_numeric\"\n",
    "    if metal_col in df.columns:\n",
    "        # Overall statistics\n",
    "        mean_val = df[metal_col].mean()\n",
    "        median_val = df[metal_col].median()\n",
    "        std_val = df[metal_col].std()\n",
    "        n_valid_samples = df[df[metal_col].notna()].shape[0]\n",
    "\n",
    "        print(f\"{metal:6}: mean = {mean_val:6.1f}, median = {median_val:6.1f}, std = {std_val:6.1f}, n_valid = {n_valid_samples:3d}\")\n",
    "    else:\n",
    "        logger.warning(f\"Metal {metal} not found in data\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot Pie Charts for Wildfire Debris and Lead"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_wildfire_debris_pie(df):\n",
    "    \"\"\"\n",
    "    Create a pie chart showing the detection of Wildfire Debris.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    df : pandas.DataFrame\n",
    "        DataFrame containing contaminant data with 'Wildfire Debris' column\n",
    "    \"\"\"\n",
    "    logger.info(\"Creating pie chart for Wildfire Debris categories\")\n",
    "    \n",
    "    # Get valid categories from config\n",
    "    valid_categories = ['Detected', 'Not Detected', 'Not Tested']\n",
    "    category_colors = ['tomato', 'lightskyblue', 'lightgray']\n",
    "    debris_counts = df['Wildfire Debris'].value_counts()\n",
    "    \n",
    "    # Ensure all valid categories are present in the counts, even if zero\n",
    "    for category in valid_categories:\n",
    "        if category not in debris_counts:\n",
    "            debris_counts[category] = 0\n",
    "    \n",
    "    # Sort by category names to ensure consistent ordering\n",
    "    debris_counts = debris_counts.sort_index()\n",
    "    \n",
    "    # Create figure\n",
    "    fig, ax = plt.subplots(figsize=(8, 8))\n",
    "    \n",
    "    # Create pie chart\n",
    "    wedges, texts, autotexts = ax.pie(\n",
    "        debris_counts,\n",
    "        labels=debris_counts.index,\n",
    "        autopct=lambda p: f'{p:.1f}%' if p > 0 else '',  # Only show percentage for non-zero slices\n",
    "        startangle=90,\n",
    "        shadow=False,\n",
    "        colors=category_colors\n",
    "    )\n",
    "\n",
    "    # Enhance the appearance of text elements\n",
    "    for text in texts:\n",
    "        text.set_fontsize(12)\n",
    "    for autotext in autotexts:\n",
    "        autotext.set_fontsize(12)\n",
    "        autotext.set_weight('bold')\n",
    "    \n",
    "    # Add title\n",
    "    ax.set_title('Distribution of Wildfire Debris Categories', fontsize=14)\n",
    "    \n",
    "    # Add legend showing actual counts, including zeros\n",
    "    legend_labels = [f\"{idx} ({count})\" for idx, count in zip(debris_counts.index, debris_counts)]\n",
    "    ax.legend(wedges, legend_labels, title=\"Categories (Count)\", \n",
    "              loc=\"lower right\", bbox_to_anchor=(1, 0, 0.5, 1), fontsize=10)\n",
    "    \n",
    "    # Equal aspect ratio ensures the pie chart is circular\n",
    "    ax.set_aspect('equal')\n",
    "    \n",
    "    # Adjust layout\n",
    "    plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2025-08-23 22:17:44,811 - INFO - Creating pie chart for Wildfire Debris categories\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 49352 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_wildfire_debris_pie(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_lead_pie(df, bin_edges, log_scale=True):\n",
    "    \"\"\"\n",
    "    Create a pie chart for lead concentration data using bin categories.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    df : pandas.DataFrame\n",
    "        DataFrame containing lead concentration data\n",
    "    bin_edges : array-like\n",
    "        The bin edges for categorizing lead concentrations\n",
    "    log_scale : bool, default=True\n",
    "        Whether to display log-scale information in labels\n",
    "    \"\"\"\n",
    "    logger.info(\"Creating pie chart for lead concentration categories\")\n",
    "    \n",
    "    # Create figure\n",
    "    fig, ax = plt.subplots(figsize=(7.5, 6))\n",
    "    \n",
    "    # Get lead data and create categories\n",
    "    lead_data = df['Max Pb_numeric']\n",
    "    \n",
    "    # Create category labels based on bin edges\n",
    "    category_labels = []\n",
    "    for i in range(len(bin_edges)-1):\n",
    "        if log_scale:\n",
    "            category_labels.append(f\"{bin_edges[i]}-{bin_edges[i+1]} µg/ft²\")\n",
    "        else:\n",
    "            category_labels.append(f\"{bin_edges[i]}-{bin_edges[i+1]} µg/ft²\")\n",
    "    \n",
    "    # Count samples in each category\n",
    "    category_counts = []\n",
    "    for i in range(len(bin_edges)-1):\n",
    "        if i == 0:\n",
    "            # For the first bin, include the lower bound\n",
    "            count = ((lead_data >= bin_edges[i]) & (lead_data < bin_edges[i+1])).sum()\n",
    "        else:\n",
    "            # For other bins, exclude the lower bound\n",
    "            count = ((lead_data >= bin_edges[i]) & (lead_data < bin_edges[i+1])).sum()\n",
    "        category_counts.append(count)\n",
    "    \n",
    "    # Get colors from magma colormap for the number of categories we have\n",
    "    n_categories = len(category_counts)\n",
    "    colors = [mpl.cm.Reds(i/n_categories) for i in range(n_categories)] + ['lightgray']\n",
    "\n",
    "    # Count NaN values for \"Not Tested\" category\n",
    "    not_tested_count = lead_data.isna().sum()\n",
    "    category_counts.append(not_tested_count)\n",
    "\n",
    "    # Add \"Not Tested\" category\n",
    "    category_labels.append(\"Not Tested\")\n",
    "    \n",
    "    if sum(category_counts) != len(lead_data):\n",
    "        logger.warning(f'Category counts sum {sum(category_counts)} does not line up with length of data {len(lead_data)}).')\n",
    "    \n",
    "    # Create pie chart\n",
    "    wedges, texts, autotexts = ax.pie(\n",
    "        category_counts,\n",
    "        labels=category_labels,\n",
    "        autopct=lambda p: f'{p:.1f}%' if p > 0 else '',  # Only show percentage for non-zero slices\n",
    "        startangle=90,\n",
    "        shadow=False,\n",
    "        colors=colors\n",
    "    )\n",
    "    \n",
    "    # Enhance the appearance of text elements\n",
    "    for text in texts:\n",
    "        text.set_fontsize(12)\n",
    "    for autotext in autotexts:\n",
    "        autotext.set_fontsize(12)\n",
    "        autotext.set_weight('bold')\n",
    "    \n",
    "    # Add title\n",
    "    ax.set_title('Lead Concentration Distribution by Category', fontsize=14)\n",
    "    \n",
    "    # Add legend showing actual counts\n",
    "    legend_labels = [f\"{label} ({count})\" for label, count in zip(category_labels, category_counts)]\n",
    "    ax.legend(wedges, legend_labels, title=\"Categories (Count)\", \n",
    "              loc=\"lower right\", bbox_to_anchor=(1, 0, 0.5, 1), fontsize=10)\n",
    "    \n",
    "    # Equal aspect ratio ensures the pie chart is circular\n",
    "    ax.set_aspect('equal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2025-08-23 22:17:44,900 - INFO - Creating pie chart for lead concentration categories\n",
      "2025-08-23 22:17:44,905 - WARNING - Category counts sum 199 does not line up with length of data 201).\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 80468 characters]",
      "text/plain": [
       "<Figure size 750x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_lead_pie(df, (0, 5, 40, 100, 250, 4000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot Boxplots for Lead and Other Contaminants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prep_for_boxplot(df, metals_to_analyze):\n",
    "    \"\"\"\n",
    "    Prepare the data for boxplot visualization by generating plot-ready dataframes.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    df : pandas.DataFrame\n",
    "        DataFrame containing contaminant data with in_boundary column\n",
    "    metals_to_analyze : list\n",
    "        List of metals to analyze\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    pandas.DataFrame\n",
    "        Processed DataFrame in long format for boxplot visualization\n",
    "    list\n",
    "        List of available metal columns\n",
    "    \"\"\"\n",
    "    # Get metal columns\n",
    "    metal_cols = [f\"{metal}_numeric\" for metal in metals_to_analyze]\n",
    "    missing_cols = [col for col in metal_cols if col not in df.columns]\n",
    "    \n",
    "    if missing_cols:\n",
    "        logger.warning(f\"Missing metal columns: {', '.join(missing_cols)}\")\n",
    "        metal_cols = [col for col in metal_cols if col in df.columns]\n",
    "    \n",
    "    # Prepare data for seaborn in long format\n",
    "    plot_data = []\n",
    "    \n",
    "    for metal, col in zip(metals_to_analyze, metal_cols):\n",
    "        if col in df.columns:\n",
    "            temp_df = pd.DataFrame({'Contaminant': metal, 'Concentration': df[col]})\n",
    "            plot_data.append(temp_df)\n",
    "    \n",
    "    plot_data = pd.concat(plot_data, ignore_index=True)\n",
    "    \n",
    "    return plot_data, metal_cols\n",
    "\n",
    "\n",
    "def plot_metals_boxplot(df,  metals_to_analyze, log_scale=True, horizontal_lines_labels=None):\n",
    "    \"\"\"\n",
    "    Create a boxplot with selected metals and stripplot overlay. Optionally add horizontal lines and labels (e.g., for EPA limits).\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    df : pandas.DataFrame\n",
    "        DataFrame containing loaded contaminant data \n",
    "    metals_to_analyze : list\n",
    "        List of metals to include in the plot\n",
    "    log_scale : bool\n",
    "        Whether to use a log scale for the y-axis\n",
    "    horizontal_lines_labels : list\n",
    "        List of tuples containing (line, label) for horizontal lines and labels \n",
    "    \"\"\"\n",
    "    # Prepare data\n",
    "    plot_data, metal_cols = prep_for_boxplot(df, metals_to_analyze)\n",
    "    logger.info(f\"Creating boxplot with metals: {', '.join(metals_to_analyze)}\")\n",
    "    \n",
    "    fig, ax = plt.subplots(figsize=(3 + 0.5*len(metal_cols), 8))\n",
    "    \n",
    "    # Create boxplot\n",
    "    sns.boxplot(\n",
    "        data=plot_data,\n",
    "        x='Contaminant',\n",
    "        y='Concentration',\n",
    "        color='0.85',\n",
    "        ax=ax,\n",
    "        showfliers=False\n",
    "    )\n",
    "    \n",
    "    # Add stripplot to show individual samples\n",
    "    sns.stripplot(\n",
    "        data=plot_data,\n",
    "        x='Contaminant',\n",
    "        y='Concentration',\n",
    "        size=3,\n",
    "        jitter=True,\n",
    "        ax=ax\n",
    "    )\n",
    "    \n",
    "    # Add horizontal lines and labels if provided\n",
    "    if horizontal_lines_labels:\n",
    "        for line, label in horizontal_lines_labels:\n",
    "            ax.axhline(y=line, color='red', linestyle='--', linewidth=1, alpha=0.7)\n",
    "            ax.text(0.6, line, label, color='red', fontsize=10, va='center', ha='left')\n",
    "    \n",
    "    # Add labels and formatting\n",
    "    ax.set_title('Indoor Contamination\\nfor All Positive Detections', fontsize=14)\n",
    "    ax.set_xlabel('Contaminant', fontsize=12)\n",
    "    ax.set_ylabel('Concentration\\n(µg/ft²)')\n",
    "\n",
    "    if log_scale:\n",
    "        ax.set_yscale('log')\n",
    "        ax.set_ylabel('Concentration', fontsize=12)\n",
    "\n",
    "    ax.tick_params(axis='x', labelsize=None, bottom=False)\n",
    "    ax.tick_params(axis='y', labelsize=None)\n",
    "    ax.grid(True, axis='y', linestyle='--', alpha=0.5)\n",
    "    \n",
    "    # Remove top and right spines\n",
    "    sns.despine(ax=ax, top=True, right=True, bottom=True)\n",
    "    \n",
    "    plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2025-08-23 22:17:44,977 - INFO - Creating boxplot with metals: Max Pb\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 40632 characters]",
      "text/plain": [
       "<Figure size 350x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_metals_boxplot(df, ['Max Pb'], horizontal_lines_labels=[(5, 'EPA Lead\\nLimit Floor'), (40, 'EPA Lead\\nLimit Window\\nSills')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2025-08-23 22:17:45,164 - INFO - Creating boxplot with metals: Max Pb, As, Sb, Ba, Be, Cd, Cr, Co, Cu, Hg, Mo, Ni, Se, Ag, Tl, V, Zn\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 86848 characters]",
      "text/plain": [
       "<Figure size 1150x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_metals_boxplot(df, metals_to_analyze)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "firespec",
   "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.11.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
