m203-svi-ejscreen.ipynb

{
  "cells": [
    {
      "cell_type": "raw",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"SVI, TRI, and Health Outcomes\"\n",
        "format: html\n",
        "---"
      ],
      "id": "03529526"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Overview\n",
        "\n",
        "In this lesson, you will use....\n",
        "\n",
        "## Learning Objectives\n",
        "\n",
        "After completing this lesson, you should be able to:\n",
        "\n",
        "-   Manipulating and processing raster data; mainly loading and cropping.\n",
        "-   Exploring SVI/raster by mapping and processing different sub datasets (overall/racial/socioeconomic).\n",
        "-   Connecting to API’s for pollution and public health data.\n",
        "-   Geolocating tabular data with latitude and longitude.\n",
        "-   Joining tabular data to spatial data (probably polygons).\n",
        "-   Connecting to OpenStreetMap API for basemaps.\n",
        "-   Creating regional basemaps.\n",
        "-   Layering spatial data on basemaps.\n",
        "    -   Semi-transparent SVI, TRI points, PLACES boundaries.\n",
        "-   Creating maps with graduated point symbols for point source pollution releases.\n",
        "-   Creating choropleths with census tract polygons for health outcomes data.\n",
        "-   Zonal statistics with census tracts for:\n",
        "    -   Mean SVI values\n",
        "\n",
        "    -   Health outcomes summaries weighted by SVI.\n",
        "-   (Maybe) Interpolate pollution load raster surface from point data.\n",
        "-   Spatial correlation analysis between SVI, point pollution, and health outcomes. .\n",
        "\n",
        "## Introduction\n",
        "\n",
        "Air Quality and Environmental Justice:\n",
        "\n",
        "In many urban areas, air quality issues disproportionately affect low-income communities and communities of color. This disparity is a key focus of environmental justice efforts. In cities like Detroit and Chicago, industrial facilities, highways, and other pollution sources are often concentrated in disadvantaged neighborhoods.\n",
        "\n",
        "NIMBY and Environmental Justice:\n",
        "\n",
        "The NIMBY phenomenon can sometimes conflict with environmental justice goals. While wealthier communities may successfully oppose new polluting facilities in their areas, this can lead to these facilities being placed in less affluent neighborhoods with less political clout, exacerbating environmental injustices.\n",
        "\n",
        "Community-Driven Science:\n",
        "\n",
        "In response to these challenges, community-driven science initiatives have emerged. These efforts involve local residents in collecting data on air quality and other environmental factors, often using low-cost sensors and mobile monitoring techniques. This approach helps fill gaps in official monitoring networks and empowers communities to advocate for themselves.\n",
        "\n",
        "Detroit and Chicago Examples:\n",
        "\n",
        "Detroit: The city has a history of industrial pollution and is working to address air quality issues in areas like Southwest Detroit, where residents face higher exposure to pollutants from nearby industries and heavy truck traffic.\n",
        "\n",
        "Chicago: The city has seen community efforts to address air pollution in areas like the Southeast Side, where residents have fought against polluting industries and advocated for stricter regulations.\n",
        "\n",
        "These topics are interconnected, as community-driven science often arises in response to environmental justice concerns, while NIMBY attitudes can influence the distribution of pollution sources across urban areas.\n",
        "\n",
        "EPA (Environmental Protection Agency):\n",
        "\n",
        "The EPA is the primary federal agency responsible for environmental protection in the United States. It sets and enforces air quality standards, including the National Ambient Air Quality Standards (NAAQS) for six criteria pollutants. The EPA also maintains the AirNow system, which provides real-time air quality data to the public. CDC (Centers for Disease Control and Prevention): While primarily focused on public health, the CDC plays a crucial role in understanding the health impacts of air pollution. It conducts research on the relationships between air quality and various health outcomes, and provides guidance on protecting public health from environmental hazards. Importance of Open Data: Open data is crucial for community-driven solutions in several ways:\n",
        "\n",
        "Transparency: Open data allows communities to verify official reports and hold authorities accountable. Accessibility: When air quality data is freely available, communities can use it to inform local decision-making and advocacy efforts. Innovation: Open data enables researchers, activists, and tech developers to create new tools and analyses that can benefit communities. Collaboration: Open data facilitates collaboration between communities, scientists, and policymakers, leading to more effective solutions.\n",
        "\n",
        "Examples of Open Data in Action:\n",
        "\n",
        "-   The EPA's Air Quality System (AQS) database is publicly accessible, allowing researchers and community groups to analyze historical air quality trends.\n",
        "\n",
        "-   In Chicago, the Array of Things project has installed sensors throughout the city, providing open data on various environmental factors including air quality.\n",
        "\n",
        "-   Detroit's Community Air Monitoring Project uses low-cost sensors to collect and share air quality data in areas underserved by official monitoring stations.\n",
        "\n",
        "-   Community-Driven Solutions: With access to open data, communities can:\n",
        "\n",
        "-   Identify local air quality hotspots that may be missed by sparse official monitoring networks.\n",
        "\n",
        "-   Correlate air quality data with health outcomes to strengthen advocacy efforts.\n",
        "\n",
        "-   Develop targeted interventions, such as promoting indoor air filtration on high-pollution days.\n",
        "\n",
        "-   Create custom alerts and information systems tailored to local needs.\n",
        "\n",
        "Challenges and Opportunities: While open data provides many benefits, challenges remain:\n",
        "\n",
        "-   Ensuring data quality and consistency, especially when integrating data from various sources.\n",
        "\n",
        "-   Addressing the \"digital divide\" to ensure all community members can access and use the data. Balancing the need for detailed local data with privacy concerns.\n",
        "\n",
        "-   Building capacity within communities to effectively use and interpret complex environmental data.\n",
        "\n",
        "The EPA, CDC, and other agencies are increasingly recognizing the value of community-driven science and open data. This has led to initiatives like the EPA's Community-Focused Exposure and Risk Screening Tool (C-FERST) and the CDC's Environmental Public Health Tracking Network, which aim to make environmental health data more accessible to communities.\n",
        "\n",
        "EGLE and Environmental Justice in Detroit\n",
        "\n",
        "-   Certainly. Detroit, Michigan has indeed seen significant environmental justice efforts, particularly in collaboration with EGLE (Michigan Department of Environment, Great Lakes, and Energy). Here are some examples:\n",
        "\n",
        "-   Detroit Environmental Agenda: This community-led initiative works closely with EGLE to address environmental concerns in Detroit. It focuses on issues like air quality, water quality, and waste management.\n",
        "\n",
        "-   48217 Community Air Monitoring Project: Named after the zip code of a heavily industrialized area in Southwest Detroit, this project involves community members working with EGLE to monitor air quality using low-cost sensors.\n",
        "\n",
        "-   Detroit Climate Action Plan: Developed in partnership with EGLE, this plan addresses climate change impacts on the city, with a focus on vulnerable communities.\n",
        "\n",
        "-   Delray Neighborhood Initiatives: EGLE has been involved in efforts to address air quality concerns in the Delray neighborhood, which is impacted by industrial emissions and heavy truck traffic.\n",
        "\n",
        "-   Green Door Initiative: This Detroit-based organization collaborates with EGLE on various environmental justice projects, including lead abatement and air quality improvement efforts.\n",
        "\n",
        "-   Detroit River Sediment Cleanup: EGLE has been involved in efforts to clean up contaminated sediments in the Detroit River, which disproportionately affects nearby low-income communities.\n",
        "\n",
        "-   Asthma Prevention Programs: EGLE supports community-based asthma prevention programs in Detroit, recognizing the link between air quality and asthma rates in disadvantaged neighborhoods.\n",
        "\n",
        "These examples demonstrate the ongoing collaboration between community groups, local government, and EGLE to address environmental justice concerns in Detroit. However, as with many environmental justice efforts, challenges remain and work is ongoing.\n",
        "\n",
        "::: {.callout-tip style=\"color: #5a7a2b;\"}\n",
        "## Data Science Review\n",
        "\n",
        "This lesson uses the [`pandas`](https://pandas.pydata.org/), [`geopandas`](https://geopandas.org/), [`matplotlib`](https://matplotlib.org/), [`numpy`](https://numpy.org/), [`requests`](https://docs.python-requests.org/en/latest/), [`contextily`](https://contextily.readthedocs.io/en/latest/), [`pygris`](https://pygris.readthedocs.io/en/latest/), [`rasterio`](https://rasterio.readthedocs.io/en/latest/), and [`xarray`](https://xarray.pydata.org/) modules. If you'd like to learn more about the functions used in this lesson, you can refer to the documentation on their respective websites.\n",
        "\n",
        "The `pandas` module is essential for data manipulation and analysis, while `geopandas` extends its functionality to handle geospatial data. `matplotlib` is used for creating static, animated, and interactive visualizations. `numpy` provides support for large, multi-dimensional arrays and matrices, along with a collection of mathematical functions to operate on these arrays.\n",
        "\n",
        "The `requests` module allows you to send HTTP requests and interact with APIs easily. `contextily` adds basemaps to your plots, enhancing the visual context of your geospatial data. `pygris` simplifies the process of working with US Census Bureau TIGER/Line shapefiles.\n",
        "\n",
        "`rasterio` and `xarray` are used for working with geospatial raster data. `rasterio` reads and writes geospatial raster datasets, while `xarray` introduces labels in the form of dimensions, coordinates, and attributes on top of raw NumPy-like arrays, making it easier to work with labeled multi-dimensional arrays.\n",
        "\n",
        "Make sure these modules are installed before you begin working with the code in this document. You can install them using pip:\n",
        ":::\n",
        "\n",
        "## Load the Data\n",
        "\n",
        "## SVI and Detroit Metro\n",
        "\n",
        "SVI has a primary overall SVI score, but also provides sublayers. These include minority, socioeconomic, housing, and household data.\n"
      ],
      "id": "0e83552e"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "import xarray as xr\n",
        "import rasterio\n",
        "import rasterio.mask\n",
        "import geopandas as gpd\n",
        "import matplotlib.pyplot as plt\n",
        "import pygris\n",
        "import numpy as np\n",
        "\n",
        "# Fetch Detroit metro area counties using pygris\n",
        "metro_counties = pygris.counties(state=\"MI\", year=2022)\n",
        "detroit_metro = metro_counties[metro_counties['NAME'].isin([\n",
        "    'Wayne', 'Oakland', 'Macomb', 'Livingston', 'St. Clair', 'Lapeer'\n",
        "])]\n",
        "\n",
        "# Dissolve the counties into a single polygon\n",
        "detroit_metro = detroit_metro.dissolve(by='STATEFP')\n",
        "\n",
        "# Convert to GeoDataFrame\n",
        "detroit_metro = gpd.GeoDataFrame(detroit_metro, geometry='geometry', crs='EPSG:4269')\n",
        "\n",
        "# Specify the TIF files\n",
        "tif_files = [\n",
        "    \"data/svi/svi_2020_tract_overall_wgs84.tif\",\n",
        "    \"data/svi/svi_2020_tract_minority_wgs84.tif\",\n",
        "    \"data/svi/svi_2020_tract_socioeconomic_wgs84.tif\",\n",
        "    \"data/svi/svi_2020_tract_housing_wgs84.tif\",\n",
        "    \"data/svi/svi_2020_tract_household_wgs84.tif\"\n",
        "]\n",
        "\n",
        "# Create an empty list to store the individual DataArrays\n",
        "data_arrays = []\n",
        "\n",
        "# Read each TIF file, clip it to Detroit metro's extent, and append it to the list\n",
        "for file in tif_files:\n",
        "    with rasterio.open(file) as src:\n",
        "        # Reproject Detroit metro boundary to match the raster CRS\n",
        "        metro_reprojected = detroit_metro.to_crs(src.crs)\n",
        "        \n",
        "        # Clip the raster to Detroit metro's geometry\n",
        "        out_image, out_transform = rasterio.mask.mask(src, metro_reprojected.geometry, crop=True)\n",
        "        out_meta = src.meta.copy()\n",
        "        \n",
        "        # Update the metadata\n",
        "        out_meta.update({\"driver\": \"GTiff\",\n",
        "                         \"height\": out_image.shape[1],\n",
        "                         \"width\": out_image.shape[2],\n",
        "                         \"transform\": out_transform})\n",
        "        \n",
        "        # Create coordinates\n",
        "        height = out_meta['height']\n",
        "        width = out_meta['width']\n",
        "        cols, rows = np.meshgrid(np.arange(width), np.arange(height))\n",
        "        xs, ys = rasterio.transform.xy(out_transform, rows, cols)\n",
        "        \n",
        "        # Convert lists to numpy arrays\n",
        "        xs = np.array(xs)\n",
        "        ys = np.array(ys)\n",
        "        \n",
        "        # Reshape coordinates to match dimensions of the raster\n",
        "        xs = xs.reshape(height, width)\n",
        "        ys = ys.reshape(height, width)\n",
        "        \n",
        "        # Create a DataArray from the clipped data\n",
        "        da = xr.DataArray(out_image[0],  # Use the first band\n",
        "                          coords={'y': ('y', ys[:, 0]),\n",
        "                                  'x': ('x', xs[0, :])},\n",
        "                          dims=['y', 'x'])\n",
        "        da.attrs['crs'] = str(src.crs)  # Convert CRS to string\n",
        "        da.attrs['transform'] = out_transform\n",
        "        data_arrays.append(da)\n",
        "\n",
        "# Combine all DataArrays into a single DataSet\n",
        "ds = xr.concat(data_arrays, dim='layer')\n",
        "\n",
        "# Rename the layers using the appropriate dimension\n",
        "ds = ds.assign_coords(layer=('layer', ['Overall', 'Minority', 'Socioeconomic', 'Housing', 'Household']))\n",
        "\n",
        "# Define the colorbar limits\n",
        "vmin, vmax = 0, 1\n",
        "\n",
        "# Create a multipanel plot\n",
        "fig, axes = plt.subplots(3, 2, figsize=(15, 20))\n",
        "axes = axes.flatten()\n",
        "\n",
        "# Plot each layer\n",
        "for i, layer in enumerate(ds.layer.values):\n",
        "    # Plot with custom color limits\n",
        "    im = ds.sel(layer=layer).plot(ax=axes[i], add_colorbar=False, vmin=vmin, vmax=vmax, cmap='viridis')\n",
        "    axes[i].set_title(layer)\n",
        "    \n",
        "    # Plot Detroit metro boundary\n",
        "    metro_reprojected.boundary.plot(ax=axes[i], color='red', linewidth=1)\n",
        "\n",
        "# Remove the extra subplot\n",
        "fig.delaxes(axes[5])\n",
        "\n",
        "# Add a single colorbar\n",
        "cbar_ax = fig.add_axes([0.92, 0.15, 0.02, 0.7])\n",
        "cbar = fig.colorbar(im, cax=cbar_ax, label='SVI Score')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "id": "14c63ede",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## TRI API Column Definitions\n",
        "\n",
        "| Column Name                      | Description                                                             |\n",
        "|------------------------------------|------------------------------------|\n",
        "| `FACILITY_NAME`                  | The name of the facility reporting to TRI                               |\n",
        "| `TRI_FACILITY_ID`                | A unique identifier for the facility in the TRI database                |\n",
        "| `STREET_ADDRESS`                 | The street address of the facility                                      |\n",
        "| `CITY_NAME`                      | The city where the facility is located                                  |\n",
        "| `COUNTY_NAME`                    | The county where the facility is located                                |\n",
        "| `STATE_ABBR`                     | The two-letter abbreviation for the state where the facility is located |\n",
        "| `ZIP_CODE`                       | The ZIP code of the facility's location                                 |\n",
        "| `PREF_LATITUDE`                  | The preferred latitude coordinate of the facility                       |\n",
        "| `PREF_LONGITUDE`                 | The preferred longitude coordinate of the facility                      |\n",
        "| `PARENT_CO_NAME`                 | The name of the parent company, if applicable                           |\n",
        "| `INDUSTRY_SECTOR_CODE`           | A code representing the industry sector of the facility                 |\n",
        "| `INDUSTRY_SECTOR`                | A description of the industry \n",
        "\n",
        "**Note**: The availability and exact names of these columns may vary depending on the specific TRI API endpoint and query parameters used. Always refer to the official EPA TRI documentation for the most up-to-date and comprehensive information.\n",
        "\n",
        "### Important Considerations\n",
        "\n",
        "1.  Not all columns may be present in every API response.\n",
        "2.  Column names may have slight variations (e.g., with or without underscores).\n",
        "3.  The EPA occasionally updates their API and data structure.\n",
        "4.  Some columns related to chemical releases and waste management may have additional variations or breakdowns.\n",
        "5.  Numeric values (like release amounts) are typically reported in pounds, but always verify the units.\n",
        "6.  For coordinates (`PREF_LATITUDE` and `PREF_LONGITUDE`), be aware that these are the preferred coordinates, which may have been adjusted for accuracy or privacy reasons.\n",
        "\n",
        "## ENVIROFACTS TRI Facilities w/ API\n"
      ],
      "id": "d30689f3"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "import xarray as xr\n",
        "import rasterio\n",
        "import rasterio.mask\n",
        "import geopandas as gpd\n",
        "import matplotlib.pyplot as plt\n",
        "import pygris\n",
        "import numpy as np\n",
        "from shapely.geometry import box\n",
        "import pandas as pd\n",
        "import requests\n",
        "import contextily as ctx\n",
        "\n",
        "# Fetch Detroit metro area counties using pygris\n",
        "metro_counties = pygris.counties(state=\"MI\", year=2022)\n",
        "detroit_metro = metro_counties[metro_counties['NAME'].isin([\n",
        "    'Wayne', 'Oakland', 'Macomb', 'Livingston', 'St. Clair', 'Lapeer'\n",
        "])]\n",
        "\n",
        "# Dissolve the counties into a single polygon\n",
        "detroit_metro = detroit_metro.dissolve(by='STATEFP')\n",
        "\n",
        "# Get the bounding box\n",
        "bbox = detroit_metro.total_bounds\n",
        "\n",
        "# Print the bounding box coordinates\n",
        "print(\"Bounding Box:\")\n",
        "print(f\"Minimum X (Longitude): {bbox[0]}\")\n",
        "print(f\"Minimum Y (Latitude): {bbox[1]}\")\n",
        "print(f\"Maximum X (Longitude): {bbox[2]}\")\n",
        "print(f\"Maximum Y (Latitude): {bbox[3]}\")\n",
        "\n",
        "# Create the bounding box as a polygon\n",
        "bbox_polygon = gpd.GeoDataFrame(\n",
        "    geometry=[box(*bbox)],\n",
        "    crs=detroit_metro.crs\n",
        ")\n",
        "\n",
        "# Fetch TRI facility data from EPA API for each county\n",
        "counties = ['Wayne', 'Oakland', 'Macomb', 'Livingston', 'St. Clair', 'Lapeer']\n",
        "tri_data = []\n",
        "\n",
        "for county in counties:\n",
        "    api_url = f\"https://data.epa.gov/efservice/tri_facility/state_abbr/MI/county_name/{county}/JSON\"\n",
        "    response = requests.get(api_url)\n",
        "    if response.status_code == 200:\n",
        "        county_data = response.json()\n",
        "        tri_data.extend(county_data)\n",
        "    else:\n",
        "        print(f\"Failed to fetch data for {county} County. Status code: {response.status_code}\")\n",
        "\n",
        "# Convert TRI data to a DataFrame\n",
        "tri_df = pd.DataFrame(tri_data)\n",
        "\n",
        "print(f\"Number of facilities fetched: {len(tri_df)}\")\n",
        "\n",
        "# Create a copy of the dataframe to avoid SettingWithCopyWarning\n",
        "tri_df_clean = tri_df.copy()\n",
        "\n",
        "# Remove facilities with empty latitude or longitude values\n",
        "tri_df_clean = tri_df_clean.dropna(subset=['pref_latitude', 'pref_longitude'])\n",
        "\n",
        "print(f\"Number of facilities after removing empty coordinates: {len(tri_df_clean)}\")\n",
        "\n",
        "# Convert latitude and longitude to numeric type\n",
        "tri_df_clean['pref_latitude'] = pd.to_numeric(tri_df_clean['pref_latitude'], errors='coerce')\n",
        "tri_df_clean['pref_longitude'] = pd.to_numeric(tri_df_clean['pref_longitude'], errors='coerce')\n",
        "\n",
        "# Function to correct longitude\n",
        "def correct_longitude(lon):\n",
        "    if lon > 0:\n",
        "        return -lon\n",
        "    return lon\n",
        "\n",
        "# Apply longitude correction\n",
        "tri_df_clean['pref_longitude'] = tri_df_clean['pref_longitude'].apply(correct_longitude)\n",
        "\n",
        "# Calculate IQR for longitude\n",
        "Q1 = tri_df_clean['pref_longitude'].quantile(0.25)\n",
        "Q3 = tri_df_clean['pref_longitude'].quantile(0.75)\n",
        "IQR = Q3 - Q1\n",
        "\n",
        "# Define bounds for outliers\n",
        "lower_bound = Q1 - 1.5 * IQR\n",
        "upper_bound = Q3 + 1.5 * IQR\n",
        "\n",
        "# Remove outliers\n",
        "tri_df_clean = tri_df_clean[(tri_df_clean['pref_longitude'] >= lower_bound) & \n",
        "                            (tri_df_clean['pref_longitude'] <= upper_bound)]\n",
        "\n",
        "print(f\"Number of facilities after removing longitude outliers: {len(tri_df_clean)}\")\n",
        "\n",
        "# Create a GeoDataFrame from the cleaned TRI data\n",
        "detroit_tri = gpd.GeoDataFrame(\n",
        "    tri_df_clean, \n",
        "    geometry=gpd.points_from_xy(tri_df_clean.pref_longitude, tri_df_clean.pref_latitude),\n",
        "    crs=\"EPSG:4326\"\n",
        ")\n",
        "\n",
        "# Reproject data to Web Mercator for contextily\n",
        "detroit_metro = detroit_metro.to_crs(epsg=3857)\n",
        "bbox_polygon = bbox_polygon.to_crs(epsg=3857)\n",
        "detroit_tri = detroit_tri.to_crs(epsg=3857)\n",
        "\n",
        "# Create the plot\n",
        "fig, ax = plt.subplots(figsize=(15, 15))\n",
        "\n",
        "# Plot the metro area and bounding box\n",
        "detroit_metro.plot(ax=ax, color='lightblue', edgecolor='black', alpha=0.5)\n",
        "bbox_polygon.boundary.plot(ax=ax, color='red', linewidth=2)\n",
        "\n",
        "# Plot TRI facilities\n",
        "detroit_tri.plot(ax=ax, color='red', markersize=50, alpha=0.7)\n",
        "\n",
        "# Add the basemap\n",
        "ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik)\n",
        "\n",
        "# Set the extent of the map to the bounding box\n",
        "ax.set_xlim(bbox_polygon.total_bounds[0], bbox_polygon.total_bounds[2])\n",
        "ax.set_ylim(bbox_polygon.total_bounds[1], bbox_polygon.total_bounds[3])\n",
        "\n",
        "# Remove axes\n",
        "ax.set_axis_off()\n",
        "\n",
        "plt.title(\"Detroit Metro Area TRI Facilities\", fontsize=16)\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n",
        "print(f\"Final number of TRI facilities in the Detroit metro area: {len(detroit_tri)}\")"
      ],
      "id": "19cef76c",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Regulated Facilities w/ ICIS-AIR\n"
      ],
      "id": "143f0a89"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "import pandas as pd\n",
        "import requests\n",
        "import time\n",
        "\n",
        "# Base URL for ECHO ICIS-AIR API\n",
        "base_url = \"https://echodata.epa.gov/echo/air_rest_services\"\n",
        "\n",
        "# Parameters for the initial API call\n",
        "params = {\n",
        "    \"output\": \"JSON\",\n",
        "    \"p_st\": \"MI\"\n",
        "}\n",
        "\n",
        "def get_facilities():\n",
        "    response = requests.get(f\"{base_url}.get_facilities\", params=params)\n",
        "    if response.status_code == 200:\n",
        "        data = response.json()\n",
        "        if 'Results' in data:\n",
        "            qid = data['Results']['QueryID']\n",
        "            print(f\"Query ID: {qid}\")\n",
        "            print(f\"Total Facilities: {data['Results']['QueryRows']}\")\n",
        "            return qid\n",
        "    print(\"Failed to get facilities and QID\")\n",
        "    return None\n",
        "\n",
        "def get_facility_data(qid):\n",
        "    all_facilities = []\n",
        "    page = 1\n",
        "    while True:\n",
        "        params = {\"qid\": qid, \"pageno\": page, \"output\": \"JSON\"}\n",
        "        response = requests.get(f\"{base_url}.get_qid\", params=params)\n",
        "        if response.status_code == 200:\n",
        "            data = response.json()\n",
        "            if 'Results' in data and 'Facilities' in data['Results']:\n",
        "                facilities = data['Results']['Facilities']\n",
        "                if not facilities:  # No more facilities to retrieve\n",
        "                    break\n",
        "                all_facilities.extend(facilities)\n",
        "                print(f\"Retrieved page {page}\")\n",
        "                page += 1\n",
        "            else:\n",
        "                break\n",
        "        else:\n",
        "            print(f\"Failed to retrieve page {page}\")\n",
        "            break\n",
        "    return all_facilities\n",
        "\n",
        "# Step 1: Get the Query ID\n",
        "qid = get_facilities()\n",
        "\n",
        "if qid:\n",
        "    # Step 2: Use get_qid to retrieve all facility data\n",
        "    print(\"Retrieving facility data...\")\n",
        "    facilities = get_facility_data(qid)\n",
        "    \n",
        "    # Convert to DataFrame\n",
        "    df_icis_air = pd.DataFrame(facilities)\n",
        "    \n",
        "    print(f\"\\nSuccessfully retrieved {len(df_icis_air)} ICIS-AIR facilities for Michigan\")\n",
        "    print(\"\\nColumns in the dataset:\")\n",
        "    print(df_icis_air.columns)\n",
        "    \n",
        "    # Display the first few rows\n",
        "    print(\"\\nFirst few rows of the data:\")\n",
        "    print(df_icis_air.head())\n",
        "    \n",
        "    # Save to CSV\n",
        "else:\n",
        "    print(\"Failed to retrieve facility data\")\n",
        "\n",
        "    import pandas as pd\n",
        "\n",
        "# List of Detroit metro counties\n",
        "metro_counties = ['Wayne', 'Oakland', 'Macomb', 'Livingston', 'St. Clair', 'Lapeer']\n",
        "\n",
        "# Subset the dataframe to include only the Detroit metro counties\n",
        "df_detroit_metro = df_icis_air[df_icis_air['AIRCounty'].isin(metro_counties)]\n",
        "\n",
        "# Print information about the subset\n",
        "print(f\"Total ICIS-AIR facilities in Michigan: {len(df_icis_air)}\")\n",
        "print(f\"ICIS-AIR facilities in Detroit metro area: {len(df_detroit_metro)}\")\n",
        "\n",
        "# Display the count of facilities in each metro county\n",
        "print(\"\\nFacilities per county:\")\n",
        "print(df_detroit_metro['AIRCounty'].value_counts())\n",
        "\n",
        "# Display the first few rows of the subset\n",
        "print(\"\\nFirst few rows of the Detroit metro ICIS-AIR facilities:\")\n",
        "print(df_detroit_metro.head())\n",
        "\n",
        "# Additional information: unique values in AIRCounty\n",
        "print(\"\\nUnique values in AIRCounty column:\")\n",
        "print(df_icis_air['AIRCounty'].unique())"
      ],
      "id": "a94d6331",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Create geopandas and plot.\n"
      ],
      "id": "2524853e"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Count records with missing coordinate values\n",
        "missing_coords = df_detroit_metro[(df_detroit_metro['FacLat'].isnull()) | (df_detroit_metro['FacLong'].isnull())]\n",
        "print(f\"Number of ICIS-AIR records with missing coordinates: {len(missing_coords)}\")\n",
        "\n",
        "# Remove records with missing coordinates\n",
        "df_detroit_metro = df_detroit_metro.dropna(subset=['FacLat', 'FacLong'])\n",
        "\n",
        "# Create a GeoDataFrame for ICIS-AIR facilities\n",
        "gdf_icis_air = gpd.GeoDataFrame(\n",
        "    df_detroit_metro, \n",
        "    geometry=gpd.points_from_xy(df_detroit_metro.FacLong, df_detroit_metro.FacLat),\n",
        "    crs=\"EPSG:4326\"\n",
        ")\n",
        "\n",
        "# Reproject ICIS-AIR data to Web Mercator\n",
        "gdf_icis_air = gdf_icis_air.to_crs(epsg=3857)\n",
        "\n",
        "# Create the plot\n",
        "fig, ax = plt.subplots(figsize=(15, 15))\n",
        "\n",
        "# Plot the metro area and bounding box (reusing objects from earlier)\n",
        "detroit_metro.plot(ax=ax, facecolor='none', edgecolor='blue', linewidth=2)\n",
        "bbox_polygon.boundary.plot(ax=ax, color='red', linewidth=2)\n",
        "\n",
        "# Plot ICIS-AIR facilities\n",
        "gdf_icis_air.plot(ax=ax, color='cyan', markersize=50, alpha=0.7, label='ICIS-AIR Facilities')\n",
        "\n",
        "# Plot TRI facilities (reusing the detroit_tri object from earlier)\n",
        "detroit_tri.plot(ax=ax, color='purple', markersize=50, alpha=0.7, label='TRI Facilities')\n",
        "\n",
        "# Add the basemap\n",
        "ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik)\n",
        "\n",
        "# Set the extent of the map to the bounding box\n",
        "ax.set_xlim(bbox_polygon.total_bounds[0], bbox_polygon.total_bounds[2])\n",
        "ax.set_ylim(bbox_polygon.total_bounds[1], bbox_polygon.total_bounds[3])\n",
        "\n",
        "# Remove axes\n",
        "ax.set_axis_off()\n",
        "\n",
        "# Add legend\n",
        "ax.legend()\n",
        "\n",
        "plt.title(\"Detroit Metro Area TRI and ICIS-AIR Facilities\", fontsize=16)\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n",
        "print(f\"Number of TRI facilities plotted: {len(detroit_tri)}\")\n",
        "print(f\"Number of ICIS-AIR facilities plotted: {len(gdf_icis_air)}\")"
      ],
      "id": "a6bd0c86",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Custom TRI Form A Search\n"
      ],
      "id": "1bd04739"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# URL of the CSV file\n",
        "url = \"https://dmap-epa-enviro-prod-export.s3.amazonaws.com/338211710.CSV\"\n",
        "\n",
        "# Read the CSV file directly with pandas\n",
        "df_tri_custom = pd.read_csv(url)\n",
        "print(f\"Successfully read CSV. Number of records: {len(df_tri_custom)}\")\n",
        "\n",
        "# Display information about the dataset\n",
        "print(\"\\nColumns in the dataset:\")\n",
        "print(df_tri_custom.columns)\n",
        "\n",
        "# Assuming the latitude and longitude columns are named 'LATITUDE' and 'LONGITUDE'\n",
        "# Adjust these names if they're different in your CSV\n",
        "lat_col = 'LATITUDE'\n",
        "lon_col = 'LONGITUDE'\n",
        "release_col = 'AIR_TOTAL_RELEASE'  # Adjust this to the actual column name for air releases\n",
        "\n",
        "# Remove records with missing coordinates or air release data\n",
        "df_tri_clean = df_tri_custom.dropna(subset=[lat_col, lon_col, release_col])\n",
        "\n",
        "print(f\"\\nNumber of records after removing missing data: {len(df_tri_clean)}\")\n",
        "\n",
        "# Create a GeoDataFrame\n",
        "gdf_tri_custom = gpd.GeoDataFrame(\n",
        "    df_tri_clean, \n",
        "    geometry=gpd.points_from_xy(df_tri_clean[lon_col], df_tri_clean[lat_col]),\n",
        "    crs=\"EPSG:4326\"\n",
        ")\n",
        "\n",
        "# Reproject to Web Mercator\n",
        "gdf_tri_custom = gdf_tri_custom.to_crs(epsg=3857)\n",
        "\n",
        "# Create the plot\n",
        "fig, ax = plt.subplots(figsize=(15, 15))\n",
        "\n",
        "# Plot the metro area and bounding box (reusing objects from earlier)\n",
        "detroit_metro.plot(ax=ax, facecolor='none', edgecolor='blue', linewidth=2)\n",
        "bbox_polygon.boundary.plot(ax=ax, color='orangered', linewidth=2)\n",
        "\n",
        "# Plot TRI facilities with graduated symbols based on air releases\n",
        "scatter = ax.scatter(gdf_tri_custom.geometry.x, gdf_tri_custom.geometry.y, \n",
        "                     # s=gdf_tri_custom[release_col]/100,  # Adjust the scaling factor as needed\n",
        "                     c='orangered',  # Static fill color\n",
        "                     edgecolor='yellow',  # Outline color\n",
        "                     linewidth=1,  # Adjust the outline width as needed\n",
        "                     alpha=0.7)\n",
        "\n",
        "# Add the basemap\n",
        "ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik)\n",
        "\n",
        "# Set the extent of the map to the bounding box\n",
        "ax.set_xlim(bbox_polygon.total_bounds[0], bbox_polygon.total_bounds[2])\n",
        "ax.set_ylim(bbox_polygon.total_bounds[1], bbox_polygon.total_bounds[3])\n",
        "\n",
        "# Remove axes\n",
        "ax.set_axis_off()\n",
        "\n",
        "# Add a legend for symbol sizes\n",
        "legend_sizes = [1000, 10000, 100000]  # Example sizes, adjust based on your data\n",
        "legend_elements = [plt.scatter([], [], s=size/100, c='orangered', edgecolor='yellow', \n",
        "                               linewidth=1, alpha=1, label=f'{size:,}') \n",
        "                   for size in legend_sizes]\n",
        "ax.legend(handles=legend_elements, title='Total Air Releases (lbs)', \n",
        "          loc='lower right', title_fontsize=12, fontsize=10)\n",
        "\n",
        "plt.title(\"Detroit Metro Area TRI Facilities - Total Air Releases (Custom Data)\", fontsize=16)\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n",
        "print(f\"\\nNumber of TRI facilities plotted: {len(gdf_tri_custom)}\")\n",
        "print(f\"Total air releases: {gdf_tri_custom[release_col].sum():,.2f} lbs\")\n",
        "print(f\"Average air release per facility: {gdf_tri_custom[release_col].mean():,.2f} lbs\")"
      ],
      "id": "5dbcc5a5",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Rasterizing Pollution Sums\n"
      ],
      "id": "14305c89"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "import numpy as np\n",
        "import geopandas as gpd\n",
        "import matplotlib.pyplot as plt\n",
        "import rasterio\n",
        "from rasterio import features\n",
        "from rasterio.transform import from_origin\n",
        "from matplotlib.colors import BoundaryNorm, ListedColormap\n",
        "import xarray as xr\n",
        "import rioxarray  # This imports rioxarray and adds the .rio accessor to xarray objects\n",
        "\n",
        "# Ensure gdf_tri_custom and detroit_metro are in the correct CRS (should be EPSG:3857 for Web Mercator)\n",
        "gdf_tri_custom = gdf_tri_custom.to_crs(epsg=3857)\n",
        "detroit_metro = detroit_metro.to_crs(epsg=3857)\n",
        "\n",
        "# Get the bounds of the Detroit metro area\n",
        "minx, miny, maxx, maxy = detroit_metro.total_bounds\n",
        "\n",
        "# Define the resolution (100m)\n",
        "resolution = 5000\n",
        "\n",
        "# Calculate the number of cells\n",
        "nx = int((maxx - minx) / resolution)\n",
        "ny = int((maxy - miny) / resolution)\n",
        "\n",
        "# Create the transform for the raster\n",
        "transform = from_origin(minx, maxy, resolution, resolution)\n",
        "\n",
        "# Prepare geometries and values for rasterization\n",
        "shapes = ((geom, value) for geom, value in zip(gdf_tri_custom.geometry, gdf_tri_custom.AIR_TOTAL_RELEASE))\n",
        "\n",
        "# Rasterize the point data\n",
        "raster = features.rasterize(shapes=shapes, \n",
        "                            out_shape=(ny, nx), \n",
        "                            transform=transform, \n",
        "                            fill=0, \n",
        "                            all_touched=True, \n",
        "                            merge_alg=rasterio.enums.MergeAlg.add)\n",
        "\n",
        "# Convert the raster to an xarray DataArray\n",
        "# Note: We use ny and nx here to ensure the coordinates match the raster shape\n",
        "raster_da = xr.DataArray(raster, \n",
        "                         coords={'y': np.linspace(maxy, miny, ny),\n",
        "                                 'x': np.linspace(minx, maxx, nx)},\n",
        "                         dims=['y', 'x'])\n",
        "raster_da.rio.write_crs(detroit_metro.crs, inplace=True)\n",
        "\n",
        "# Clip the raster with the Detroit metro boundary\n",
        "clipped_raster = raster_da.rio.clip(detroit_metro.geometry.values, detroit_metro.crs, drop=False, all_touched=True)\n",
        "\n",
        "# Define the breaks for the discrete scale\n",
        "breaks = [0, 1, 10, 100, 1000, 10000, 100000, 250000, 500000]\n",
        "\n",
        "# Create a custom colormap\n",
        "colors = ['#FFFFFF', '#FFFFCC', '#FFEDA0', '#FED976', '#FEB24C', '#FD8D3C', '#FC4E2A', '#E31A1C', '#B10026']\n",
        "cmap = ListedColormap(colors)\n",
        "\n",
        "# Create a normalization based on the breaks\n",
        "norm = BoundaryNorm(breaks, cmap.N)\n",
        "\n",
        "# Create the plot\n",
        "fig, ax = plt.subplots(figsize=(15, 15))\n",
        "\n",
        "# Plot the clipped raster with the custom colormap and norm\n",
        "im = ax.imshow(clipped_raster, extent=[minx, maxx, miny, maxy], origin='upper', \n",
        "               cmap=cmap, norm=norm)\n",
        "\n",
        "# Add colorbar with discrete labels\n",
        "cbar = plt.colorbar(im, ax=ax, extend='max', \n",
        "                    label='Total Air Releases (pounds)', \n",
        "                    ticks=breaks)\n",
        "cbar.ax.set_yticklabels([f'{b:,}' for b in breaks])\n",
        "\n",
        "# Plot the TRI facility points\n",
        "gdf_tri_custom.plot(ax=ax, color='blue', markersize=20, alpha=0.7)\n",
        "\n",
        "# Plot the Detroit metro boundary\n",
        "detroit_metro.boundary.plot(ax=ax, color='black', linewidth=2)\n",
        "\n",
        "# Set the extent to match the Detroit metro area\n",
        "ax.set_xlim(minx, maxx)\n",
        "ax.set_ylim(miny, maxy)\n",
        "\n",
        "# Add title and labels\n",
        "ax.set_title('TRI Air Total Release (100m resolution sum) with Facility Locations', fontsize=16)\n",
        "ax.set_xlabel('X Coordinate')\n",
        "ax.set_ylabel('Y Coordinate')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n",
        "print(f\"Number of TRI facilities plotted: {len(gdf_tri_custom)}\")\n",
        "print(f\"Total air releases: {gdf_tri_custom['AIR_TOTAL_RELEASE'].sum():,.2f}\")\n",
        "print(f\"Maximum cell value in raster: {clipped_raster.max().values:,.2f}\")"
      ],
      "id": "4eeddf9c",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Air Release Vulnerability Index\n"
      ],
      "id": "9cfbe1af"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "import numpy as np\n",
        "import xarray as xr\n",
        "import rioxarray\n",
        "import geopandas as gpd\n",
        "import matplotlib.pyplot as plt\n",
        "from rasterio.enums import Resampling\n",
        "from scipy import stats\n",
        "\n",
        "# Select the 'Overall' layer\n",
        "svi_overall = ds.sel(layer='Overall')\n",
        "\n",
        "# Convert to rioxarray for geospatial operations\n",
        "svi_overall = svi_overall.rio.write_crs(\"EPSG:4326\")\n",
        "\n",
        "# Reproject SVI to match the CRS of the air release raster\n",
        "svi_reprojected = svi_overall.rio.reproject_match(clipped_raster)\n",
        "\n",
        "# Clip SVI raster to the Detroit metro boundary\n",
        "svi_clipped = svi_reprojected.rio.clip(detroit_metro.geometry.values, detroit_metro.crs, drop=True, all_touched=True)\n",
        "\n",
        "# Disaggregate the air release data to match the resolution of the SVI data\n",
        "air_release_disaggregated = clipped_raster.rio.reproject_match(\n",
        "    svi_clipped,\n",
        "    resampling=Resampling.bilinear\n",
        ")\n",
        "\n",
        "# Calculate raster correlation between SVI overall and air release within Detroit metro boundary\n",
        "svi_flat = svi_clipped.values.flatten()\n",
        "air_flat = air_release_disaggregated.values.flatten()\n",
        "\n",
        "# Remove NaN values\n",
        "mask = ~np.isnan(svi_flat) & ~np.isnan(air_flat)\n",
        "svi_flat = svi_flat[mask]\n",
        "air_flat = air_flat[mask]\n",
        "\n",
        "correlation, p_value = stats.pearsonr(svi_flat, air_flat)\n",
        "\n",
        "print(f\"Raster Correlation between SVI Overall and Air Release within Detroit metro: {correlation:.4f}\")\n",
        "print(f\"P-value: {p_value:.4f}\")\n",
        "\n",
        "# New correlation analysis\n",
        "# Ensure gdf_tri_custom is in the same CRS as svi_clipped\n",
        "gdf_tri_custom = gdf_tri_custom.to_crs(svi_clipped.rio.crs)\n",
        "\n",
        "# Extract SVI values at TRI facility locations\n",
        "svi_values = []\n",
        "for point in gdf_tri_custom.geometry:\n",
        "    svi_value = svi_clipped.sel(x=point.x, y=point.y, method=\"nearest\").values\n",
        "    svi_values.append(svi_value)\n",
        "\n",
        "# Add SVI values to gdf_tri_custom\n",
        "gdf_tri_custom['SVI_OVERALL'] = svi_values\n",
        "\n",
        "# Perform correlation analysis\n",
        "correlation_points, p_value_points = stats.pearsonr(gdf_tri_custom['AIR_TOTAL_RELEASE'], gdf_tri_custom['SVI_OVERALL'])\n",
        "\n",
        "print(f\"Point-based Correlation between SVI Overall and AIR_TOTAL_RELEASE: {correlation_points:.4f}\")\n",
        "print(f\"P-value: {p_value_points:.4f}\")\n",
        "\n",
        "# Log1p transform the air release data and scale to 0-1\n",
        "air_release_log = np.log1p(air_release_disaggregated)\n",
        "air_release_scaled = (air_release_log - air_release_log.min()) / (air_release_log.max() - air_release_log.min())\n",
        "\n",
        "# Multiply scaled air release data with SVI data\n",
        "vulnerability_indicator = air_release_scaled * svi_clipped\n",
        "\n",
        "# Create the plots\n",
        "fig, axs = plt.subplots(3, 1, figsize=(15, 45))\n",
        "\n",
        "# Plot SVI Overall\n",
        "im1 = svi_clipped.plot(ax=axs[0], cmap='viridis', vmin=0, vmax=1, add_colorbar=False)\n",
        "plt.colorbar(im1, ax=axs[0], label='SVI Overall')\n",
        "axs[0].set_title('Social Vulnerability Index (Overall)', fontsize=16)\n",
        "detroit_metro.boundary.plot(ax=axs[0], color='black', linewidth=2)\n",
        "\n",
        "# Plot Original Air Release (log-transformed for better visualization)\n",
        "im2 = np.log1p(air_release_disaggregated).plot(ax=axs[1], cmap='YlOrRd', add_colorbar=False)\n",
        "plt.colorbar(im2, ax=axs[1], label='Log(Air Release + 1)')\n",
        "axs[1].set_title('Air Release (Log-transformed)', fontsize=16)\n",
        "detroit_metro.boundary.plot(ax=axs[1], color='black', linewidth=2)\n",
        "\n",
        "# Plot Air Release Vulnerability Indicator\n",
        "im3 = vulnerability_indicator.plot(ax=axs[2], cmap='YlOrRd', vmin=0, vmax=1, add_colorbar=False)\n",
        "plt.colorbar(im3, ax=axs[2], label='Air Release Vulnerability Indicator')\n",
        "axs[2].set_title('Air Release Vulnerability Indicator\\n(Scaled Air Release * SVI)', fontsize=16)\n",
        "detroit_metro.boundary.plot(ax=axs[2], color='black', linewidth=2)\n",
        "\n",
        "for ax in axs:\n",
        "    ax.set_xlabel('Longitude')\n",
        "    ax.set_ylabel('Latitude')\n",
        "    ax.set_xlim(svi_clipped.x.min(), svi_clipped.x.max())\n",
        "    ax.set_ylim(svi_clipped.y.min(), svi_clipped.y.max())\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n",
        "# Print some statistics\n",
        "print(f\"Minimum vulnerability indicator: {vulnerability_indicator.min().values:.4f}\")\n",
        "print(f\"Maximum vulnerability indicator: {vulnerability_indicator.max().values:.4f}\")\n",
        "print(f\"Mean vulnerability indicator: {vulnerability_indicator.mean().values:.4f}\")"
      ],
      "id": "006bdafd",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Convert the vulnerability indicator to a pandas DataFrame\n",
        "vulnerability_df = vulnerability_indicator.to_dataframe(name='index').reset_index()\n",
        "\n",
        "# Sort by index value and get the top 10\n",
        "top_10 = vulnerability_df.sort_values('index', ascending=False).head(10)\n",
        "\n",
        "# Create points from the coordinates\n",
        "top_10['geometry'] = gpd.points_from_xy(top_10.x, top_10.y)\n",
        "top_10_gdf = gpd.GeoDataFrame(top_10, geometry='geometry', crs=vulnerability_indicator.rio.crs)\n",
        "\n",
        "# Create the final map\n",
        "fig, ax = plt.subplots(figsize=(15, 15))\n",
        "\n",
        "# Plot the Detroit metro boundary\n",
        "detroit_metro.boundary.plot(ax=ax, color='black', linewidth=2)\n",
        "\n",
        "# Plot the top 10 points\n",
        "top_10_gdf.plot(ax=ax, color='blue', markersize=100, alpha=0.7)\n",
        "\n",
        "# Add labels to the points\n",
        "for idx, row in top_10_gdf.iterrows():\n",
        "    ax.annotate(f\"#{idx+1}\", (row.geometry.x, row.geometry.y), \n",
        "                xytext=(3, 3), textcoords=\"offset points\", \n",
        "                color='white', fontweight='bold')\n",
        "\n",
        "# Add a basemap\n",
        "ctx.add_basemap(ax, crs=vulnerability_indicator.rio.crs.to_string(), source=ctx.providers.OpenStreetMap.Mapnik)\n",
        "\n",
        "# Set the extent to match the Detroit metro area\n",
        "ax.set_xlim(vulnerability_indicator.x.min(), vulnerability_indicator.x.max())\n",
        "ax.set_ylim(vulnerability_indicator.y.min(), vulnerability_indicator.y.max())\n",
        "\n",
        "ax.set_title('Top 10 Areas with Highest Air Release Vulnerability Index', fontsize=16)\n",
        "ax.set_axis_off()\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n",
        "# Print the coordinates of the top 10 points\n",
        "print(\"Coordinates of the top 10 points:\")\n",
        "for idx, row in top_10_gdf.iterrows():\n",
        "    print(f\"#{idx+1}: ({row.geometry.x}, {row.geometry.y})\")"
      ],
      "id": "4c7d7054",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## PLACES\n",
        "\n",
        "The CDC PLACES (Population Level Analysis and Community Estimates) dataset is a collaboration between the Centers for Disease Control and Prevention (CDC), the Robert Wood Johnson Foundation, and the CDC Foundation. It provides model-based population-level analysis and community estimates of health indicators for all counties, places (incorporated and census designated places), census tracts, and ZIP Code Tabulation Areas (ZCTAs) across the United States.\n",
        "\n",
        "### Key Points\n",
        "\n",
        "1. **Geographic Coverage:** Entire United States, including all 50 states, the District of Columbia, and Puerto Rico.\n",
        "2. **Geographic Granularity:** Multiple levels including counties, cities/towns, census tracts, and ZIP codes.\n",
        "3. **Health Indicators:** Wide range of chronic disease measures related to health outcomes, prevention, and health risk behaviors.\n",
        "4. **Data Sources:** \n",
        "   - Behavioral Risk Factor Surveillance System (BRFSS)\n",
        "   - U.S. Census Bureau's American Community Survey (ACS)\n",
        "5. **Methodology:** Uses small area estimation methods for small geographic areas.\n",
        "6. **Health Measures Include:**\n",
        "   - Chronic diseases: e.g., asthma, COPD, heart disease, diabetes\n",
        "   - Health risk behaviors: e.g., smoking, physical inactivity, binge drinking\n",
        "   - Prevention practices: e.g., health insurance coverage, dental visits, cholesterol screening\n",
        "7. **Socioeconomic Data:** Includes some socioeconomic and demographic variables.\n",
        "8. **Annual Updates:** Providing recent estimates for local areas.\n",
        "\n",
        "This dataset is particularly valuable for:\n",
        "- Public health researchers\n",
        "- Policymakers\n",
        "- Community organizations\n",
        "\n",
        "It provides a standardized way to compare health indicators across different geographic areas and can be used to inform targeted interventions and policy decisions, especially in addressing health disparities at a local level.\n",
        "\n",
        "### Processing\n"
      ],
      "id": "b572e34c"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "import geopandas as gpd\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "import requests\n",
        "import contextily as ctx\n",
        "\n",
        "# Define the GeoJSON API endpoint\n",
        "url = \"https://data.cdc.gov/resource/cwsq-ngmh.geojson\"\n",
        "\n",
        "# Define the Detroit metro area counties\n",
        "detroit_counties = ['Wayne', 'Oakland', 'Macomb', 'Livingston', 'St. Clair', 'Lapeer']\n",
        "\n",
        "# Create the county filter string\n",
        "county_filter = \" OR \".join([f\"countyname = '{county}'\" for county in detroit_counties])\n",
        "\n",
        "# Define the query parameters\n",
        "params = {\n",
        "    \"$where\": f\"stateabbr = 'MI' AND ({county_filter})\",\n",
        "    \"$limit\": 50000  # Adjust if necessary\n",
        "}\n",
        "\n",
        "# Make the API request\n",
        "response = requests.get(url, params=params)\n",
        "\n",
        "if response.status_code == 200:\n",
        "    data = response.json()\n",
        "    print(f\"Successfully retrieved data\")\n",
        "else:\n",
        "    print(f\"Failed to retrieve data. Status code: {response.status_code}\")\n",
        "    print(response.text)\n",
        "\n",
        "# Convert to GeoDataFrame\n",
        "gdf = gpd.read_file(response.text)\n",
        "\n",
        "# Print available health measures\n",
        "print(\"\\nAvailable health measures:\")\n",
        "print(gdf['measure'].unique())\n",
        "\n",
        "# Print the first few rows to see the structure of the data\n",
        "print(\"\\nFirst few rows of the GeoDataFrame:\")\n",
        "print(gdf.head())\n",
        "\n",
        "# Print basic information about the GeoDataFrame\n",
        "print(\"\\nGeoDataFrame Info:\")\n",
        "print(gdf.info())\n",
        "\n",
        "# Create a sample map for one health measure (e.g., Current asthma)\n",
        "fig, ax = plt.subplots(figsize=(15, 15))\n",
        "\n",
        "# Filter for the specific measure and ensure data_value is numeric\n",
        "gdf_asthma = gdf[gdf['measure'] == 'Current asthma among adults'].copy()\n",
        "gdf_asthma['data_value'] = pd.to_numeric(gdf_asthma['data_value'], errors='coerce')\n",
        "\n",
        "# Plot the asthma data\n",
        "gdf_asthma.plot(column='data_value', \n",
        "                ax=ax, \n",
        "                legend=True, \n",
        "                legend_kwds={'label': 'Asthma Prevalence (%)', 'orientation': 'horizontal'},\n",
        "                cmap='YlOrRd',\n",
        "                missing_kwds={'color': 'lightgrey'})\n",
        "\n",
        "# Add county boundaries\n",
        "#gdf_asthma.dissolve(by='countyname').boundary.plot(ax=ax, color='black', linewidth=0.5)\n",
        "\n",
        "# Add basemap\n",
        "ctx.add_basemap(ax, crs=gdf_asthma.crs.to_string(), source=ctx.providers.OpenStreetMap.Mapnik)\n",
        "\n",
        "# Set the extent to match the Detroit metro area\n",
        "ax.set_xlim(gdf_asthma.total_bounds[0], gdf_asthma.total_bounds[2])\n",
        "ax.set_ylim(gdf_asthma.total_bounds[1], gdf_asthma.total_bounds[3])\n",
        "\n",
        "plt.title('Asthma Prevalence in Detroit Metro Area', fontsize=16)\n",
        "ax.axis('off')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n",
        "# Print some statistics for the asthma data\n",
        "print(\"\\nAsthma Statistics:\")\n",
        "print(f\"Average asthma prevalence: {gdf_asthma['data_value'].mean():.2f}%\")\n",
        "print(f\"Minimum asthma prevalence: {gdf_asthma['data_value'].min():.2f}%\")\n",
        "print(f\"Maximum asthma prevalence: {gdf_asthma['data_value'].max():.2f}%\")\n",
        "\n",
        "# Print the number of census tracts per county\n",
        "print(\"\\nNumber of census tracts per county:\")\n",
        "print(gdf_asthma['countyname'].value_counts())"
      ],
      "id": "02afc407",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "import numpy as np\n",
        "import geopandas as gpd\n",
        "import matplotlib.pyplot as plt\n",
        "import contextily as ctx\n",
        "from scipy.interpolate import griddata\n",
        "import rasterio\n",
        "from rasterio.transform import from_origin\n",
        "\n",
        "# Assuming gdf_asthma is already created and contains the asthma data\n",
        "\n",
        "# Ensure gdf_asthma is in EPSG:4326\n",
        "gdf_asthma = gdf_asthma.to_crs(epsg=4326)\n",
        "\n",
        "# Extract coordinates and values\n",
        "X = gdf_asthma.geometry.x.values\n",
        "Y = gdf_asthma.geometry.y.values\n",
        "Z = gdf_asthma['data_value'].values\n",
        "\n",
        "# Remove any NaN values\n",
        "mask = ~np.isnan(Z)\n",
        "X, Y, Z = X[mask], Y[mask], Z[mask]\n",
        "\n",
        "# Create a grid to interpolate over\n",
        "grid_resolution = 0.025  # in degrees\n",
        "x_min, y_min, x_max, y_max = gdf_asthma.total_bounds\n",
        "grid_x = np.arange(x_min, x_max, grid_resolution)\n",
        "grid_y = np.arange(y_min, y_max, grid_resolution)\n",
        "grid_xx, grid_yy = np.meshgrid(grid_x, grid_y)\n",
        "\n",
        "# Perform IDW interpolation\n",
        "points = np.column_stack((X, Y))\n",
        "grid_z = griddata(points, Z, (grid_xx, grid_yy), method='linear')\n",
        "\n",
        "# Create the plot\n",
        "fig, ax = plt.subplots(figsize=(15, 15))\n",
        "\n",
        "# Plot the interpolated data\n",
        "im = ax.imshow(grid_z, extent=[x_min, x_max, y_min, y_max], \n",
        "               origin='lower', cmap='YlOrRd', alpha=0.7)\n",
        "\n",
        "# Add colorbar\n",
        "cbar = plt.colorbar(im, ax=ax, label='Asthma Prevalence')\n",
        "\n",
        "# Add basemap\n",
        "ctx.add_basemap(ax, crs=gdf_asthma.crs.to_string(), source=ctx.providers.OpenStreetMap.Mapnik)\n",
        "\n",
        "# Set the extent to match the Detroit metro area\n",
        "ax.set_xlim(x_min, x_max)\n",
        "ax.set_ylim(y_min, y_max)\n",
        "\n",
        "plt.title('IDW Interpolated Asthma Prevalence in Detroit Metro Area', fontsize=16)\n",
        "ax.axis('off')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n",
        "# Save the interpolated raster\n",
        "transform = from_origin(x_min, y_max, grid_resolution, grid_resolution)\n",
        "with rasterio.open('idw_asthma.tif', 'w', \n",
        "                   driver='GTiff', \n",
        "                   height=grid_z.shape[0], \n",
        "                   width=grid_z.shape[1], \n",
        "                   count=1, \n",
        "                   dtype=grid_z.dtype, \n",
        "                   crs='EPSG:4326',  # Explicitly set the CRS\n",
        "                   transform=transform) as dst:\n",
        "    dst.write(grid_z, 1)\n",
        "\n",
        "# Print some statistics about the interpolated data\n",
        "print(f\"Minimum interpolated value: {np.nanmin(grid_z):.2f}\")\n",
        "print(f\"Maximum interpolated value: {np.nanmax(grid_z):.2f}\")\n",
        "print(f\"Mean interpolated value: {np.nanmean(grid_z):.2f}\")"
      ],
      "id": "6ee59564",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "::: {.callout-tip style=\"color: #5a7a2b;\"}\n",
        "\n",
        "## Data  Review\n",
        "\n",
        "The Toxics Release Inventory (TRI) and the Integrated Compliance Information System for Air (ICIS-AIR) are two important but distinct environmental reporting systems maintained by the U.S. Environmental Protection Agency (EPA). They have several key differences:\n",
        "\n",
        "- **Regulatory Basis**\n",
        "  - TRI: Established under the Emergency Planning and Community Right-to-Know Act (EPCRA) of 1986\n",
        "  - ICIS-AIR: Part of the Clean Air Act (CAA) compliance and enforcement program\n",
        "\n",
        "- **Focus**\n",
        "  - TRI: Tracks the management of certain toxic chemicals that may pose a threat to human health and the environment\n",
        "  - ICIS-AIR: Focuses specifically on air quality and emissions from facilities regulated under the Clean Air Act\n",
        "\n",
        "- **Reported Information**\n",
        "  - TRI: Facilities report on releases, waste management, and pollution prevention activities for specific toxic chemicals\n",
        "  - ICIS-AIR: Tracks emissions data, compliance status, and enforcement actions related to air quality regulations\n",
        "\n",
        "- **Facility Coverage**\n",
        "  - TRI: Covers facilities in specific industries that manufacture, process, or use TRI-listed chemicals above certain thresholds\n",
        "  - ICIS-AIR: Includes a broader range of facilities that emit air pollutants, regardless of the specific chemicals involved\n",
        "\n",
        "- **Reporting Thresholds**\n",
        "  - TRI: Has specific chemical thresholds that trigger reporting requirements\n",
        "  - ICIS-AIR: Generally doesn't have chemical-specific thresholds; requirements are based on overall emissions and facility type\n",
        "\n",
        "- **Public Accessibility**\n",
        "  - TRI: Designed with a strong focus on public right-to-know, with data easily accessible to the public\n",
        "  - ICIS-AIR: While public, it's primarily designed for regulatory and enforcement purposes\n",
        "\n",
        "- **Data Frequency**\n",
        "  - TRI: Annual reporting is required for covered facilities\n",
        "  - ICIS-AIR: May involve more frequent reporting, depending on permit requirements and compliance status\n",
        "\n",
        "- **Scope of Pollutants**\n",
        "  - TRI: Focuses on a specific list of toxic chemicals and chemical categories\n",
        "  - ICIS-AIR: Covers a wider range of air pollutants, including criteria air pollutants and hazardous air pollutants\n",
        "\n",
        "- **Use in Environmental Management**\n",
        "  - TRI: Often used for assessing long-term trends in toxic chemical releases and waste management practices\n",
        "  - ICIS-AIR: More commonly used for day-to-day air quality management and enforcement activities\n",
        "\n",
        "- **Geographic Coverage**\n",
        "  - TRI: Nationwide program with consistent reporting across states\n",
        "  - ICIS-AIR: While national, implementation can vary more by state or local air quality management district\n",
        ":::\n",
        "\n",
        "### Deployment\n",
        "\n",
        "Congratulations! .... Now you should be able to:\n",
        "\n",
        "-   Test test...\n",
        "\n",
        "## Lesson 3\n",
        "\n",
        "In this lesson, we explored ....\n",
        "\n",
        "[Lesson 3](https://ciesin-geospatial.github.io/TOPSTSCHOOL-air-quality/m203-grdiv1-pm25.html){.btn .btn-primary .btn role=\"button\"}"
      ],
      "id": "68890af4"
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "language": "python",
      "display_name": "Python 3 (ipykernel)"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}