{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "# Rivulet: U. S. EPA Air Quality System\n", "_by Michelle H Wilkerson, Lucas Coletti, and Adelmo Eloy_" ] }, { "cell_type": "markdown", "id": "1", "metadata": {}, "source": [ "### Purpose of this Notebook\n", "\n", "This notebook is focused on **Air Quality** as a phenomenon. While most students understand that poor Air Quality can impact health, they may not know that there are many different kinds of air pollution, each caused by different processes and chemicals. These are reflected by different patterns over the course of a day, week, major event, or year.\n", "\n", "This data tool allows users to connect to the U. S. Environmental Protection Agency's Air Quality System (AQS) API, which provides air quality data for dates up to six months before the present date (the time lag allows for data to be scientifically validated). You are able to search for air quality data streams in an area of interest, identify a date range of interest within that area, and then access a variety of combinations of data that allow for explorations of different pollutants, comparison of AQI across different regions, and more. The datasets this notebook is constucted to fetch can serve as a launching point for examining what air quality is, and what are its underlying mechanistic and compositional complexities.\n", "\n", "
\n", " Click here for more information about APIs\n", "\n", "API stands for **Application Programming Interface**. Think of it as a language to communicate with data centers so you can search and get the data you need. For students, using APIs is like having a direct line to the most advanced scientific sensors on the planet. It allows us to work with the same **live, real-world data** that professional scientists use to track what's going on in the world around us.\n", "\n", "You are welcome to modify and adapt this script. You may find the AQS API documentation [here](https://aqs.epa.gov/aqsweb/documents/data_api.html) and the `pyaqsapi` documentation [here](https://usepa.github.io/pyaqsapi/) helpful.\n", "\n", "This notebook was developed as part of NSF Grant 2445609 to support accessing and processing public datasets for middle and high school classroom activities. It's written to be relatively accessible to beginners, but if you have not interacted with computational notebooks or Python before, you may find navigating this tool difficult. (Check out the Show Your Work project for a gentle introduction to computational notebooks for educators!)\n", "\n", "Our project is focused on supporting data analysis and mechanistic reasoning in science education. In other words, we want students to learn how data provides information about _how scientific mechanisms work_ and how understanding scientific mechanisms can help them to _explain and interpret patterns in data_. This builds on a long history of research on complex systems and agent-based modeling, and more closely connects that work to current expansions of data analysis across subjects." ] }, { "cell_type": "markdown", "id": "2", "metadata": {}, "source": [ "# Part I: Setup\n", "\n", "First, you need to connect to the API and specify what region and time period you are interested in. This section will help you with that." ] }, { "cell_type": "markdown", "id": "3", "metadata": {}, "source": [ "### Connecting with AQS\n", "\n", "Before you get started, you will need an AQS API Key. To get one, use the url `https://aqs.epa.gov/data/api/signup?email=myemail@example.com` (replace myemail@example.com with your email) and paste it into a browser. You will receive a cute sounding API key via email to the address you provided in the URL. \n", "\n", "Copy your API key and change EMAIL and API_KEY in the cell block below to your email and key from the service. \n", "\n", "You can run the notebook using the test email and key that are already provided, however, the test account has a limited number of uses per day and may not work. Register for an account as soon as you know you'd like to use the service. If you lose your key and need a new one, use the same url with the same email address." ] }, { "cell_type": "code", "execution_count": null, "id": "4", "metadata": {}, "outputs": [], "source": [ "# EDIT BELOW: Depending on what tool you are using to run this \n", "# notebook, you may need to replace the \"%\"\" below with \"!\"\n", "!pip install pyaqsapi # This installs the AQS API package from the EPA\n", "\n", "EMAIL = \"test@aqs.api\" # EDIT HERE: with your registered email\n", "API_KEY = \"test\" # EDIT HERE: with your API key" ] }, { "cell_type": "markdown", "id": "5", "metadata": {}, "source": [ "The `pyaqsapi` package provides our requested data as pandas dataframes. Below, we run a query to see if the service is available. We've found this service in particular to be finicky at times. If you don't get a successful response on a query, take a break and try again later. Remember that you are more likely to get a successful response if you use your own API key rather than the test key provided here.\n", "\n", "#### **Note for Google CoLab users:**\n", "If you are running this notebook in Google CoLab, uncomment the first line of the code below. This will make it easier for you to download your data files. There will also be lines to uncomment in each \"Data Fetch\" section, once you are ready to get your data." ] }, { "cell_type": "code", "execution_count": null, "id": "6", "metadata": {}, "outputs": [], "source": [ "# from google.colab import files ## UNCOMMENT THIS LINE IF USING COLAB\n", "import pyaqsapi as aqs # library for getting data\n", "import pandas as pd # library for working with data\n", "\n", "aqs.aqs_credentials(username=EMAIL, key=API_KEY)\n", "\n", "aqs.aqs_is_available()" ] }, { "cell_type": "markdown", "id": "7", "metadata": {}, "source": [ "### Customizing Your Location\n", "\n", "This section allows you to specify a location and a date for which you would like to collect data. You'll need to know the approximate longitude and latitude of the region you are interested in. One easy way to do this is by asking Google, \"What is the longitude and latitude of [area]?\" You will use the code in this section to use a minimum and maximum latitude and longitude to draw a bounding box around your location of interest. This will filter your future queries to focus only on air quality sensors found within that box. If you don't find any sensors, select a different location or increase the bounding box size." ] }, { "cell_type": "markdown", "id": "8", "metadata": {}, "source": [ "Now, identify the specific location you want to explore. AQS can fetch all the air quality data within a bounding box. Let's not get crazy - start with a relatively small bounding box (try one degree for an urban area) and get bigger if you need to. \n", "\n", "(Tip: If you click on a location in Google Maps, you'll see the lat and long for that point in the URL.)" ] }, { "cell_type": "code", "execution_count": null, "id": "9", "metadata": {}, "outputs": [], "source": [ "# Define a bounding box around your target region. \n", "# If it is densely populated, we suggest you start with \n", "# a bounding box that is only one degree in area. \n", "\n", "min_lat = 37 # EDIT HERE\n", "max_lat = 38 # EDIT HERE\n", "min_long = -122.5 # EDIT HERE\n", "max_long = -121.5 # EDIT HERE" ] }, { "cell_type": "markdown", "id": "10", "metadata": {}, "source": [ "Below, we'll create a map with the box you defined, to make sure you're capturing what you want." ] }, { "cell_type": "code", "execution_count": null, "id": "11", "metadata": {}, "outputs": [], "source": [ "%pip install folium #install the mapping library\n", "import folium\n", "\n", "bbox = [[min_lat, min_long], [max_lat, max_long]]\n", "\n", "# Calculate the center of the box to position the map\n", "map_center = [(bbox[0][0] + bbox[1][0]) / 2, (bbox[0][1] + bbox[1][1]) / 2]\n", "\n", "# Create a Folium map object\n", "m = folium.Map(location=map_center, zoom_start=8)\n", "\n", "# Add a rectangle for the bounding box to the map\n", "folium.Rectangle(\n", " bounds=bbox,\n", " color=\"#ff0000\", # Red border\n", " fill=True,\n", " fill_color=\"#ff7800\", # Orange fill\n", " fill_opacity=0.2\n", ").add_to(m)\n", "\n", "m" ] }, { "cell_type": "markdown", "id": "12", "metadata": {}, "source": [ "### Customizing Your Date\n", "\n", "Next, specify a \"target date.\" Different parts of the notebook will collect different time periods of data (usually about a week's worth to a month's worth of data). This part makes sure that the date you choose is included in the dataset. You may want to consider when a major air quality event (like a wildfire, fireworks, weather event, etc.) happened, or you can choose what you understand to be a \"typical\" day to get baseline datasets for your region." ] }, { "cell_type": "code", "execution_count": null, "id": "13", "metadata": {}, "outputs": [], "source": [ "from datetime import datetime, timedelta\n", "\n", "target_date = \"04-10-2020\" # EDIT HERE\n", "target_datetime = datetime.strptime(target_date, \"%m-%d-%Y\")" ] }, { "cell_type": "markdown", "id": "14", "metadata": {}, "source": [ "### Defining Pollutants" ] }, { "cell_type": "markdown", "id": "15", "metadata": {}, "source": [ "AQS has settings for different collections of parameters that reflect different \"classes\" of interest. For example, one parameter is SCHOOL AIR TOXICS, which lists 125 industrial air pollutants that have been tracked near some K-12 schools, or HAZARDOUS AIR POLLUTANTS with a total of 407 pollutants. \n", "\n", "We want to start simple here, and focus on the pollutants that are most well-known and most likely to be monitored across the country. This also helps us make sure we're not taxing the system with big queries when we don't need to. AQI POLLUTANTS includes codes for the \"big five\" pollutants that are used to calculate the Air Quality Index: " ] }, { "cell_type": "code", "execution_count": null, "id": "16", "metadata": {}, "outputs": [], "source": [ "try:\n", " parameter_list = aqs.aqs_parameters_by_class(\"AQI POLLUTANTS\")\n", "except Exception as e:\n", " print(f\"Something didn't work: {e}\")\n", "\n", "parameter_list" ] }, { "cell_type": "markdown", "id": "17", "metadata": {}, "source": [ "You can use this notebook to explore other pollutants by replacing the \"codes\" you see above with codes for the pollutant you're interested in. \n", "\n", "In each section below, we will suggest different combinations of pollutants to explore. You can always come back here to get the code of other pollutant(s) you might want to explore." ] }, { "cell_type": "markdown", "id": "18", "metadata": {}, "source": [ "# Part II: Fetch Data" ] }, { "cell_type": "markdown", "id": "19", "metadata": {}, "source": [ "Now, we are ready to begin fetching data for your investigation. Each section below allows you to download a _different kind of dataset_, designed to support exploration of _different patterns and issues related to air quality_. You can scroll through to learn more about each one, or click a link to go directly to a section.\n", "\n", "[Exploring Pollutant Behaviors: Ozone and PM2.5](#pollutants)\n", "[Reconstructing AQI](#reconstructing)\n", "[Comparing Neighbors](#comparing)" ] }, { "cell_type": "markdown", "id": "20", "metadata": {}, "source": [ "## Exploring Pollutant Behaviors: Ozone and PM2.5\n", " Data can highlight how different air pollutants have different sources and behaviors. For example, Ozone often peaks in the heat of the afternoon (a photochemical reaction), while PM2.5 (smoke and dust) might stay high all night or spike during morning traffic. \n", "\n", "Here, we guide you through downloading a dataset that you can use to compare the behavior of PM2.5 and Ozone, or other pollutants that you might be interested in." ] }, { "cell_type": "code", "execution_count": null, "id": "21", "metadata": {}, "outputs": [], "source": [ "PARAMETERS = ('44201', #ozone\n", " '88101') #pm2.5" ] }, { "cell_type": "markdown", "id": "22", "metadata": {}, "source": [ "### Choose Your Site and Length of Time" ] }, { "cell_type": "markdown", "id": "23", "metadata": {}, "source": [ "First, let's look for the monitoring stations in your custom region that were actually monitoring the pollutants defined above during the time you are interested in." ] }, { "cell_type": "code", "execution_count": null, "id": "24", "metadata": {}, "outputs": [], "source": [ "bdate = target_datetime\n", "edate = target_datetime + timedelta(days=1)\n", "\n", "monitors = aqs.bybox.monitors(\n", " parameter= ','.join(PARAMETERS), # all six pollutants are listed\n", " bdate=bdate, \n", " edate=edate,\n", " minlat=min_lat,\n", " maxlat=max_lat,\n", " minlon=min_long,\n", " maxlon=max_long,\n", ")\n", "\n", "# Filter the monitors so we are showing only the ones that have all requested parameters\n", "relevant_monitors = monitors.groupby(\n", " ['state_code','county_code','site_number']).filter( # make groups of each monitor\n", " lambda group: set(group['parameter_code']) == set(PARAMETERS) # include only groups with both monitors\n", ")\n", "\n", "\n", "relevant_monitors[['state_code','county_code','site_number','address','city_name','county_name','state_name']].drop_duplicates()" ] }, { "cell_type": "markdown", "id": "25", "metadata": {}, "source": [ "Hopefully, at least one monitor in your designated area appeared that has all the pollutants you are looking for. If not, you may want change the pollutants you are looking for (some pollutants, like SO or PM10, and rarely monitored compared to PM2.5 or Ozone).\n", "\n", "We can also adjust how many days of data you will get. This might be helpful depending on the kind of event you want to explore (for example, you may only need a few days to explore the impacts of fireworks on air quality, versus a few weeks if you are tracking the impacts of a longer-lasting wildfire). In the code block below, we specify how many weeks of data we'd like to fetch around our target date. By default, we have set this to fetch two weeks of data, with your target date right in the middle of the time period." ] }, { "cell_type": "code", "execution_count": null, "id": "26", "metadata": {}, "outputs": [], "source": [ "beginningdate = target_datetime - timedelta(weeks=1)\n", "enddate = target_datetime + timedelta(weeks=1)" ] }, { "cell_type": "markdown", "id": "27", "metadata": {}, "source": [ "Now, use list of sites above to find the monitoring station you want to use to get your data. Enter the corresponding stateFIPS (state_code), countycode, and sitenum (site_number) to get your data.\n", "\n", "By default, we have set the code below to get data from a monitoring station in Concord, CA that we know has been monitoring all five AQI pollutants." ] }, { "cell_type": "markdown", "id": "28", "metadata": {}, "source": [ "### Fetch Your Data" ] }, { "cell_type": "code", "execution_count": null, "id": "29", "metadata": {}, "outputs": [], "source": [ "# CUSTOMIZE THIS: Enter the stateFIPS, county_code, and sitenum of your preferred \n", "# monitoring station from the list above. \n", "\n", "aqdata = aqs.bysite.sampledata(parameter=','.join(PARAMETERS),\n", " bdate=beginningdate,\n", " edate=enddate,\n", " stateFIPS=\"06\", # EDIT HERE enter your state_code \n", " countycode=\"001\", # EDIT HERE enter your county_code \n", " sitenum=\"0013\") # EDIT HERE enter your site_number \n", "\n", "aqdata" ] }, { "cell_type": "markdown", "id": "30", "metadata": {}, "source": [ "The data we get from the AQS has dates and times, but not a combined datetime field. This would make it hard for us to analyze our data from hour to hour (for example, there would not be an easy way for graphing programs to know that 11pm Dec 31, 2020 comes before 1am on Jan 1, 2021), so let's combine the information so dates and times are all organized correctly." ] }, { "cell_type": "code", "execution_count": null, "id": "31", "metadata": {}, "outputs": [], "source": [ "aqdata['datetime_local'] = aqdata.apply(lambda row: datetime.strptime(f\"{row['date_local']} {row['time_local']}\", \"%Y-%m-%d %H:%M\"), axis=1)" ] }, { "cell_type": "markdown", "id": "32", "metadata": {}, "source": [ "Now, we can graph the data to see if it's demonstrating patterns that are worth exploring further with your students." ] }, { "cell_type": "markdown", "id": "33", "metadata": {}, "source": [ "### Plot to Explore the Patterns" ] }, { "cell_type": "code", "execution_count": null, "id": "34", "metadata": {}, "outputs": [], "source": [ "import seaborn as sns # a library for graphing\n", "\n", "g = sns.FacetGrid(aqdata, \n", " row=\"parameter_code\", \n", " aspect=4, # this makes the graphs 4x wider than they are tall\n", " sharey=False) # this treats the scale for each pollutant separately\n", "\n", "g.map(sns.lineplot, \"datetime_local\", \"sample_measurement\")" ] }, { "cell_type": "markdown", "id": "35", "metadata": {}, "source": [ "### Export for Further Investigation" ] }, { "cell_type": "markdown", "id": "36", "metadata": {}, "source": [ "If this set of data looks interesting enough to keep working with, use the code below to export to csv. We first pivot the data so that it's in a format that's easier to work with using student-friendly tools like Excel or CODAP." ] }, { "cell_type": "code", "execution_count": null, "id": "37", "metadata": {}, "outputs": [], "source": [ "aqdata_pivoted = aqdata.pivot_table(\n", " index='datetime_local', # Pivot so each row is a datetime\n", " columns='parameter_code',\n", " values='sample_measurement'\n", ").reset_index()\n", "\n", "aqdata_pivoted\n", "\n", "aqdata_pivoted.to_csv(\"pollutantlevels.csv\") \n", "\n", "# files.download(\"pollutantlevels.csv\") ## UNCOMMENT THIS LINE IF USING COLAB" ] }, { "cell_type": "markdown", "id": "38", "metadata": {}, "source": [ "## Reconstructing AQI from Pollutant Levels\n", " The Air Quality Index (AQI) is a **composite index** used to communicate how clean or polluted the air is. Instead of tracking just one chemical, it looks at a set of five main pollutants: Ozone, Carbon Monoxide, Sulfur Dioxide, Nitrogen Dioxide, and particulate matter (PM2.5 and PM10). \n", "\n", "Generating a dataset that includes all of these pollutants is pedagogically powerful because different pollutants have different **\"signature patterns.\"** By reconstructing the AQI ourselves, we can see exactly which pollutant is \"driving\" the risk level at any given time.\n", "\n", "\"alt" ] }, { "cell_type": "markdown", "id": "39", "metadata": {}, "source": [ "### Select all AQI Pollutants\n", "In the block below, we define the `PARAMETERS` we want to fetch. By default, we include all six major pollutants that are tracked by the AQI measure. " ] }, { "cell_type": "code", "execution_count": null, "id": "40", "metadata": {}, "outputs": [], "source": [ "# CUSTOMIZE THIS: Identify the parameters that you want to include in your dataset.\n", "\n", "PARAMETERS = (\n", " '44201', # Ozone\n", " '88101', # PM2.5\n", " '42101', # Carbon monoxide\n", " '42401', # Sulfur dioxide\n", " '42602', # Nitrogen dioxide (NO2)\n", " '81102' # PM10\n", ")" ] }, { "cell_type": "markdown", "id": "41", "metadata": {}, "source": [ "### Find Monitoring Stations that Measure All Pollutants\n", "To generate that most accurate AQI calculation, we need to find a station that reports **all six** the pollutants that are tracked. The code below searches the local region you defined above for these \"complete\" monitoring stations. \n", "\n", "Unfortunately, not many monitoring stations track every pollutant. It's pretty tricky to find a monitoring station that has all the AQI pollutants in practice. If you can't find one for your area, try reducing the parameters above. Usually, the pollutants that are most important for setting the overall AQI are PM2.5, Ozone, and NO2. Try adjusting the required parameter list in the code above to find a \"good enough\" monitoring station for your region." ] }, { "cell_type": "code", "execution_count": null, "id": "42", "metadata": {}, "outputs": [], "source": [ "# RUN THIS: This code figures out what monitoring stations have the information you need.\n", "# If there are no monitoring stations in your area that have all the pollutants, use the\n", "# box above to reduce the number of pollutants you require.\n", "\n", "monitors = aqs.bybox.monitors(\n", " parameter= ','.join(PARAMETERS), # all six pollutants are listed\n", " bdate=target_datetime, \n", " edate=target_datetime + timedelta(days=1),\n", " minlat=min_lat,\n", " maxlat=max_lat,\n", " minlon=min_long,\n", " maxlon=max_long,\n", ")\n", "\n", "# Filter the monitors so we are showing only the ones that have all requested parameters\n", "relevant_monitors = monitors.groupby(\n", " ['state_code','county_code','site_number']).filter( # make groups of each monitor\n", " lambda group: set(group['parameter_code']) == set(PARAMETERS) # include only groups with both monitors\n", ")\n", "\n", "\n", "relevant_monitors[['state_code','county_code','site_number','address','city_name','county_name','state_name']].drop_duplicates()" ] }, { "cell_type": "markdown", "id": "43", "metadata": {}, "source": [ "### Fetch the Time-Series Data\n", "Now we zoom in on a specific station. We are fetching a **two-week window** of data around our target date. This allows students to see several day/night cycles (diurnal patterns), which helps them distinguish between consistent biological/physical rhythms and unique weather events. You can adjust the time period using the next code block." ] }, { "cell_type": "code", "execution_count": null, "id": "44", "metadata": {}, "outputs": [], "source": [ "# CUSTOMIZE THIS: Define the period of time for which you would like to fetch \n", "# data. By default, the code below specifies a beginning date one week before\n", "# the target date, and one week after.\n", "\n", "beginningdate = target_datetime - timedelta(weeks=1)\n", "enddate = target_datetime + timedelta(weeks=1)" ] }, { "cell_type": "markdown", "id": "45", "metadata": {}, "source": [ "In the following block of code, we fetch data from the monitoring station of your choice. Adjust the `stateFIPS` (state_code), `countycode`, and `sitenum` to match a monitoring station from the list above that has all the pollutants you are looking for." ] }, { "cell_type": "code", "execution_count": null, "id": "46", "metadata": {}, "outputs": [], "source": [ "# CUSTOMIZE THIS: Enter the stateFIPS, county_code, and sitenum of your preferred \n", "# monitoring station from the list above. \n", "\n", "aqdata = aqs.bysite.sampledata(parameter=','.join(PARAMETERS),\n", " bdate=beginningdate,\n", " edate=enddate,\n", " stateFIPS=\"06\", # enter the state_code here\n", " countycode=\"013\", # enter the county_code here\n", " sitenum=\"0002\") # enter the site_number here\n", "\n", "aqdata.head()" ] }, { "cell_type": "markdown", "id": "47", "metadata": {}, "source": [ "### Translate Pollutant Levels to the Composite AQI\n", "The next code blocks contain the EPA's official **\"breakpoints\"** and the functions needed to translate raw concentrations (like ppm or µg/m³) into the 0-500 AQI scale. This illustrates how scientists translate complex measurements into simpler measurements to communicate about physical phenomena. To learn more about this conversion, you can check out the EPA's technical documentation [here](https://document.airnow.gov/technical-assistance-document-for-the-reporting-of-daily-air-quailty.pdf)." ] }, { "cell_type": "code", "execution_count": null, "id": "48", "metadata": {}, "outputs": [], "source": [ "# RUN THIS: Define the EPA Air Quality Index (AQI) Breakpoints\n", "# Format: (AQI_Low, AQI_High, Conc_Low, Conc_High)\n", "\n", "AQI_BREAKPOINTS = {\n", " # PM2.5 (24-hr, in µg/m³)\n", " '88101': [\n", " (0, 50, 0.0, 12.0),\n", " (51, 100, 12.1, 35.4),\n", " (101, 150, 35.5, 55.4),\n", " (151, 200, 55.5, 150.4),\n", " (201, 300, 150.5, 250.4),\n", " (301, 400, 250.5, 350.4),\n", " (401, 500, 350.5, 500.4),\n", " ],\n", " # PM10 (24-hr, in µg/m³)\n", " '81102': [\n", " (0, 50, 0, 54),\n", " (51, 100, 55, 154),\n", " (101, 150, 155, 254),\n", " (151, 200, 255, 354),\n", " (201, 300, 355, 424),\n", " (301, 400, 425, 504),\n", " (401, 500, 505, 604),\n", " ],\n", " # O3 (8-hr, in ppm)\n", " '44201': [\n", " (0, 50, 0.000, 0.054),\n", " (51, 100, 0.055, 0.070),\n", " (101, 150, 0.071, 0.085),\n", " (151, 200, 0.086, 0.105),\n", " (201, 300, 0.106, 0.200),\n", " # Note: AQI > 200 for 8-hr O3 is calculated using 1-hr O3.\n", " # This implementation assumes you will provide the 8-hr value\n", " # and will cap at the 201-300 range.\n", " ],\n", " # O3 (1-hr, in ppm) - Used when 8-hr values are high TODO\n", " 'O3_1hr': [\n", " (101, 150, 0.125, 0.164),\n", " (151, 200, 0.165, 0.204),\n", " (201, 300, 0.205, 0.404),\n", " (301, 400, 0.405, 0.504),\n", " (401, 500, 0.505, 0.604),\n", " ],\n", " # CO (8-hr, in ppm)\n", " '42101': [\n", " (0, 50, 0.0, 4.4),\n", " (51, 100, 4.5, 9.4),\n", " (101, 150, 9.5, 12.4),\n", " (151, 200, 12.5, 15.4),\n", " (201, 300, 15.5, 30.4),\n", " (301, 400, 30.5, 40.4),\n", " (401, 500, 40.5, 50.4),\n", " ],\n", " # SO2 (1-hr, in ppb)\n", " '42401': [\n", " (0, 50, 0, 35),\n", " (51, 100, 36, 75),\n", " (101, 150, 76, 185),\n", " (151, 200, 186, 304),\n", " (201, 300, 305, 604),\n", " (301, 400, 605, 804),\n", " (401, 500, 805, 1004),\n", " ],\n", " # NO2 (1-hr, in ppb)\n", " '42602': [\n", " (0, 50, 0, 53),\n", " (51, 100, 54, 100),\n", " (101, 150, 101, 360),\n", " (151, 200, 361, 649),\n", " (201, 300, 650, 1249),\n", " (301, 400, 1250, 1649),\n", " (401, 500, 1650, 2049),\n", " ],\n", "}" ] }, { "cell_type": "markdown", "id": "49", "metadata": {}, "source": [ "The code below takes all the thresholds that are defined above, and creates functions that will use them to calculate the AQI for the data we have fetched. These functions will become the machinery we apply to the data to figure out the overall composite AQI." ] }, { "cell_type": "code", "execution_count": null, "id": "50", "metadata": {}, "outputs": [], "source": [ "import math # we need this library for some of the calculations below\n", "\n", "def calculate_individual_aqi(pollutant_code, concentration): \n", " if concentration is None or concentration < 0 or \\\n", " not math.isfinite(concentration):\n", " return None\n", "\n", " C_p = concentration\n", "\n", " # Find the correct breakpoint category\n", " table = AQI_BREAKPOINTS[pollutant_code]\n", " \n", " for (I_lo, I_hi, C_lo, C_hi) in table:\n", " if C_lo <= C_p <= C_hi:\n", " # Avoid division by zero if C_hi == C_lo\n", " if (C_hi - C_lo) == 0:\n", " aqi = I_lo\n", " else:\n", " # Apply the linear interpolation formula\n", " aqi = ((I_hi - I_lo) / (C_hi - C_lo)) * (C_p - C_lo) + I_lo\n", " \n", " return round(aqi)\n", " \n", " # If concentration is beyond the highest breakpoint,\n", " # use the last category for calculation\n", " (I_lo, I_hi, C_lo, C_hi) = table[-1]\n", " if C_p > C_hi:\n", " if (C_hi - C_lo) == 0:\n", " aqi = I_lo\n", " else:\n", " aqi = ((I_hi - I_lo) / (C_hi - C_lo)) * (C_p - C_lo) + I_lo\n", " return round(aqi)\n", "\n", "def calculate_composite_aqi(concentration_data):\n", " individual_aqis = []\n", " \n", " for pollutant_code, concentration in concentration_data.items():\n", " aqi = calculate_individual_aqi(pollutant_code, concentration)\n", " \n", " if aqi is not None:\n", " individual_aqis.append(aqi)\n", " \n", " # The composite AQI is the highest of the available individual AQIs\n", " if not individual_aqis:\n", " return None\n", " \n", " return max(individual_aqis)\n", "\n", "# Helper to build pollutant -> concentration dict and compute composite AQI\n", "def _row_to_composite(row):\n", " conc = {str(code): row[code] for code in row.index if pd.notnull(row[code])}\n", " return calculate_composite_aqi(conc)" ] }, { "cell_type": "markdown", "id": "51", "metadata": {}, "source": [ "### Calculate the Composite AQI\n", "Now we apply our \"engine\" to the raw data. We first calculate the AQI for each individual pollutant at every time point. " ] }, { "cell_type": "code", "execution_count": null, "id": "52", "metadata": {}, "outputs": [], "source": [ "aqdata['individual_aqi'] = aqdata.apply(\n", " lambda row: calculate_individual_aqi(row['parameter_code'], row['sample_measurement']),\n", " axis=1\n", ")" ] }, { "cell_type": "markdown", "id": "53", "metadata": {}, "source": [ "The data we get from the AQS has dates and times, but not a combined datetime field. This would make it hard for us to analyze our data from hour to hour (for example, there would not be an easy way for graphing programs to know that 11pm Dec 31, 2020 comes before 1am on Jan 1, 2021), so let's combine the information so dates and times are all organized correctly." ] }, { "cell_type": "code", "execution_count": null, "id": "54", "metadata": {}, "outputs": [], "source": [ "combined_series = aqdata['date_local'] + ' ' + aqdata['time_local']\n", "aqdata['datetime_local'] = pd.to_datetime(combined_series, format=\"%Y-%m-%d %H:%M\")" ] }, { "cell_type": "markdown", "id": "55", "metadata": {}, "source": [ "Now, we are going to pivot the dataset so that it is easier to compare pollutants that were measured at the same time. This is what we need to do in order to figure out the composite AQI for each moment of time when pollutants were measured." ] }, { "cell_type": "code", "execution_count": null, "id": "56", "metadata": {}, "outputs": [], "source": [ "# Pivot so each row is a datetime and columns are pollutant codes with their measurements\n", "measurements = aqdata.pivot_table(\n", " index='datetime_local',\n", " columns='parameter_code',\n", " values='sample_measurement',\n", ")" ] }, { "cell_type": "markdown", "id": "57", "metadata": {}, "source": [ "Now, we find the \"highest\" value across all pollutants to determine the **Composite AQI**. The code below also does some cleanup to make sure dates are always formatted in the same way and everything is labeled properly." ] }, { "cell_type": "code", "execution_count": null, "id": "58", "metadata": {}, "outputs": [], "source": [ "# Compute composite AQI for each datetime\n", "measurements['composite_aqi'] = measurements.apply(_row_to_composite, axis=1)\n", "\n", "# Make a tidy timeseries dataframe\n", "aqi_timeseries = measurements[['composite_aqi']].reset_index().sort_values('datetime_local')\n", "\n", "new_rows = pd.DataFrame({\n", " 'datetime_local': aqi_timeseries['datetime_local'],\n", " 'sample_measurement': aqi_timeseries['composite_aqi'],\n", " 'parameter_code': 'COMPOSITE_AQI', # tag so these rows are identifiable\n", " 'parameter': 'Composite AQI',\n", " 'individual_aqi': aqi_timeseries['composite_aqi']\n", "})\n", "\n", "# Concatenate and keep a consistent datetime dtype; sort by datetime if desired.\n", "aqdata = pd.concat([aqdata, new_rows], ignore_index=True, sort=False)\n", "aqdata['datetime_local'] = pd.to_datetime(aqdata['datetime_local'])\n", "aqdata = aqdata.sort_values('datetime_local').reset_index(drop=True)" ] }, { "cell_type": "markdown", "id": "59", "metadata": {}, "source": [ "### Visualize Patterns\n", "Now, we can plot all our pollutants and the composite AQI on the same time scale. \n", "\n", "This is a great way to see which pollutant is the \"limiting factor\" or the primary risk driver at different times of the day. For example, you might see the composite score being driven by Ozone in the afternoon and PM2.5 at night.\n", "\n", "These visual patterns can raise interesting questions about AQI: Does a score of 50 always mean the same thing? Why or why not? What are the relationships between the different pollutants? " ] }, { "cell_type": "code", "execution_count": null, "id": "60", "metadata": {}, "outputs": [], "source": [ "# RUN THIS: Plot a grid of reported levels for each pollutant, plus composite AQI\n", "\n", "import seaborn as sns\n", "\n", "# All that dataframe action messes us the indices. \n", "# We need to reset the index to visualize everything together.\n", "aqdata = aqdata.reset_index()\n", "\n", "# Let's also sort the aqdata so that COMPOSITE_AQI is at the top of the grid.\n", "aqdata = aqdata.sort_values(\"parameter_code\", ascending=False)\n", "\n", "g = sns.FacetGrid(aqdata, \n", " row=\"parameter_code\", \n", " aspect=4, # this makes the graphs 4x wider than they are tall\n", " sharey=False) # this treats the scale for each pollutant separately\n", "\n", "g.map(sns.lineplot, \"datetime_local\", \"sample_measurement\")" ] }, { "cell_type": "markdown", "id": "61", "metadata": {}, "source": [ "### Export for Further Investigation\n", "If the patterns look interesting and you'd like to use this dataset with your students, you can export the data to a `.csv` file. This file can then be uploaded to tools like CODAP or Excel for further analysis." ] }, { "cell_type": "code", "execution_count": null, "id": "62", "metadata": {}, "outputs": [], "source": [ "# Let's pivot so all measures are collected in one row with the datetime \n", "aqdata = aqdata.reset_index().pivot_table(index='datetime_local', columns='parameter_code', values='sample_measurement')\n", "\n", "aqdata.to_csv(\"aqipluspollutantlevels.csv\")\n", "\n", "# files.download(\"aqipluspollutantlevels.csv\") ## UNCOMMENT THIS LINE IF USING COLAB" ] }, { "cell_type": "markdown", "id": "63", "metadata": {}, "source": [ "## Comparing Air Quality Between Neighbors\n", " Another interesting activity to do with air quality data is to explore differences in air quality between different locations in the same region. This may be useful to think about the impacts of particular natural (e.g. a coastal breeze or mountains) or man-made features (e.g. the presence of a freeway or factory) on local air quality patterns over longer periods of time.\n", "\n", "This section helps you conduct a search within the region you defined in Section 1, specifically for finding monitoring stations with the largest differences in mean concentrations between pollutants. " ] }, { "cell_type": "markdown", "id": "64", "metadata": {}, "source": [ "### Specify Your Pollutant and Time Interval" ] }, { "cell_type": "markdown", "id": "65", "metadata": {}, "source": [ "First, let's identify the specific pollutant you want to compare across sites. Here, we look at PM2.5." ] }, { "cell_type": "code", "execution_count": null, "id": "66", "metadata": {}, "outputs": [], "source": [ "pollutants = \"88101\" # EDIT HERE to identify the pollutant you want to compare." ] }, { "cell_type": "markdown", "id": "67", "metadata": {}, "source": [ "Now, let's define a time period to compare. Below, our default is to look at data for a one-month period around the target date you specified in the setup." ] }, { "cell_type": "code", "execution_count": null, "id": "68", "metadata": {}, "outputs": [], "source": [ "# define the month around the target datetime\n", "bdate = target_datetime - timedelta(weeks=2)\n", "edate = target_datetime + timedelta(weeks=2)" ] }, { "cell_type": "markdown", "id": "69", "metadata": {}, "source": [ "### Find Neighboring Monitoring Stations" ] }, { "cell_type": "markdown", "id": "70", "metadata": {}, "source": [ "Let's find the monitoring sites within the specified region that have the most dramatic differences in mean pollutant concentrations.\n", "\n", "The code below takes daily summaries for each site in the region, for the period of time specified above." ] }, { "cell_type": "code", "execution_count": null, "id": "71", "metadata": {}, "outputs": [], "source": [ "# RUN THIS: For all sites within the bounding box\n", "# get the aggregated stats by site of the requested pollutant(s) \n", "\n", "aqsummary = aqs.bybox.dailysummary(parameter=pollutants,\n", " bdate=bdate,\n", " edate=edate,\n", " minlat=min_lat,\n", " maxlat=max_lat,\n", " minlon=min_long,\n", " maxlon=max_long)\n", "\n", "aqsummary.head()" ] }, { "cell_type": "markdown", "id": "72", "metadata": {}, "source": [ "### Compare and Map Means Between Neighbors" ] }, { "cell_type": "markdown", "id": "73", "metadata": {}, "source": [ "Now, we'll compute overall arithmetic mean of those daily summaries above for each site. We'll then rank the sites by mean to see where we might find some extreme differences. This isn't the only information we need, but it's a good first step for now." ] }, { "cell_type": "code", "execution_count": null, "id": "74", "metadata": {}, "outputs": [], "source": [ "# sort the sites from highest to lowest mean pollutant concentration.\n", "meanaq = aqsummary.groupby(\n", " ['state_code', 'county_code', 'site_number']\n", ").mean(numeric_only=True).sort_values(\n", " by=\"arithmetic_mean\", \n", " ascending=False\n", ")\n", "\n", "meanaq[['arithmetic_mean']]" ] }, { "cell_type": "markdown", "id": "75", "metadata": {}, "source": [ "Lets map these sites and color each site by its relative mean pollutant concentration. This might help you find out more about which specific regional differences for this pollutant might be most interesting or compelling for students." ] }, { "cell_type": "code", "execution_count": null, "id": "76", "metadata": {}, "outputs": [], "source": [ "import branca.colormap as cm # for coloring map markers\n", "\n", "linear_scale = cm.LinearColormap(\n", " [\"green\", \"yellow\", \"red\"],\n", " vmin=min(meanaq['arithmetic_mean']), vmax=max(meanaq['arithmetic_mean'])\n", ")\n", "\n", "map_center = [meanaq['latitude'].mean(), meanaq['longitude'].mean()]\n", "m = folium.Map(location=map_center, zoom_start=8)\n", "\n", "# 3. Loop through the DataFrame and add markers\n", "for index, row in meanaq.iterrows():\n", " folium.CircleMarker(\n", " location=[row['latitude'], row['longitude']],\n", " radius=10,\n", " tooltip= index, # Show the name on hover\n", " fill_color=linear_scale(row['arithmetic_mean']),\n", " fill_opacity=1\n", " ).add_to(m)\n", "\n", "m" ] }, { "cell_type": "markdown", "id": "77", "metadata": {}, "source": [ "### Choose Focal Neighbors\n", "\n", "Okay, now that you've got your map and table, it's time to pick your favorites! Look for two or three sites that are close to each other but might have different surroundings (like one near a freeway and another in a park). \n", "\n", "Copy the `state_code`, `county_code`, and `site_number` for your chosen sites and plop them into the list below. We'll pull the full month of data for these neighbors so we can see how they compare day-to-day. This dataset is great for comparing distributions or time-series patterns." ] }, { "cell_type": "code", "execution_count": null, "id": "78", "metadata": {}, "outputs": [], "source": [ "# EDIT HERE: Choose 2-3 sites from the table above to compare.\n", "# Format: [state_code, county_code, site_number]\n", "my_sites = [\n", " ['06', '087', '1005'], \n", " ['06', '013', '1004']\n", "]\n", "\n", "# make a list of all the columns that will be returned\n", "aqs_columns = [\n", " \"state_code\", \"county_code\", \"site_number\", \"parameter_code\", \"poc\",\n", " \"datum\", \"parameter_name\", \"sample_duration\", \"pollutant_standard\",\n", " \"date_local\", \"time_local\", \"date_gmt\", \"time_gmt\", \"sample_measurement\",\n", " \"units_of_measure\", \"sample_frequency\", \"detection_limit\", \"uncertainty\",\n", " \"qualifier\", \"method_type\", \"method_code\", \"method_description\",\n", " \"cbsa_name\", \"state_name\", \"county_name\", \"site_address\", \"local_site_name\",\n", " \"date_of_last_change\", \"latitude\", \"longitude\"\n", "]\n", "\n", "neighbor_comparison = pd.DataFrame(columns=aqs_columns)\n", "\n", "for site in my_sites:\n", " neighbor_data = aqs.bysite.sampledata(parameter=pollutants,\n", " bdate=bdate,\n", " edate=edate,\n", " stateFIPS=site[0], # enter the state_code here\n", " countycode=site[1], # enter the county_code here\n", " sitenum=site[2]) # enter the site_number here\n", "\n", " neighbor_comparison = pd.concat([neighbor_comparison, neighbor_data], ignore_index=True)\n", "\n", "neighbor_comparison.head()" ] }, { "cell_type": "markdown", "id": "79", "metadata": {}, "source": [ "### Visualize Distributions of Pollutants per Neighbor" ] }, { "cell_type": "markdown", "id": "80", "metadata": {}, "source": [ "Let's take a look at what we have! " ] }, { "cell_type": "code", "execution_count": null, "id": "81", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "neighbor_comparison['Site_ID'] = (neighbor_comparison['state_code'].astype(str) + \"-\" + \n", " neighbor_comparison['county_code'].astype(str) + \"-\" + \n", " neighbor_comparison['site_number'].astype(str))\n", "\n", "plt.figure(figsize=(12, 6))\n", "\n", "# Boxplot\n", "sns.boxplot(data=neighbor_comparison, x='Site_ID', y='sample_measurement', color='white')\n", "\n", "# Overlay Dots (Stripplot)\n", "sns.stripplot(data=neighbor_comparison, x='Site_ID', y='sample_measurement', color='blue', alpha=0.4, jitter=True)\n", "\n", "plt.xticks(rotation=45) # Rotate labels so they don't overlap\n", "plt.title('Distribution by Hierarchical Site ID')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "82", "metadata": {}, "source": [ "### Download the Dataset" ] }, { "cell_type": "markdown", "id": "83", "metadata": {}, "source": [ "If you think the distributions you see above are worth more investigation -- to see what's happening over time, to see if the differences are statistically significant, or even to do some further work to try and understand why these differences exist, go ahead and download the dataset below!" ] }, { "cell_type": "code", "execution_count": null, "id": "84", "metadata": {}, "outputs": [], "source": [ "# You might want to check out what's going on over time, so let's datetime it\n", "combined_series = neighbor_comparison['date_local'] + ' ' + neighbor_comparison['time_local']\n", "neighbor_comparison['datetime_local'] = pd.to_datetime(combined_series, format=\"%Y-%m-%d %H:%M\")\n", "\n", "# Let's pivot so all measures are collected in one row with the datetime \n", "neighbor_comparison = neighbor_comparison.reset_index().pivot_table(index='datetime_local', columns=['Site_ID'], values='sample_measurement')\n", "\n", "neighbor_comparison.to_csv(\"neighborcomparison.csv\")\n", "\n", "# files.download(\"neighborcomparison.csv\") ## UNCOMMENT THIS LINE IF USING COLAB" ] }, { "cell_type": "markdown", "id": "85", "metadata": {}, "source": [ "# Credits\n", "\n", "This notebook was developed as part of NSF Grant 2445609 to support accessing and processing public datasets for middle and high school classroom activities.\n", "\n", "Data provided by the U.S. Environmental Protection Agency Air Quality System (AQS). You can learn more about AQS data [here](https://aqs.epa.gov/aqsweb/documents/about_aqs_data.html).\n", "\n", "Special thanks to the authors of the `pyaqsapi` package for making this data so accessible to the Python community.\n", "\n", "Developed by Michelle H Wilkerson, Lucas Coletti, and Adelmo Eloy as part of the Rivulet Project." ] } ], "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.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }