m204-grdiv1-pm25.ipynb

{
  "cells": [
    {
      "cell_type": "raw",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Particulate Matter Across Socioeconomic Strata of Countries\"\n",
        "author: \"Juan F. Martinez\"\n",
        "format: html\n",
        "bibliography: pm25-grdi-references.bib\n",
        "---"
      ],
      "id": "71b5f1ae"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Overview\n",
        "\n",
        "In this lesson, you will use NASA socioeconomic and environmental Earthdata available at NASA [Socioeconomic Data and Applications Center](https://earthdata.nasa.gov/centers/sedac-daac) (SEDAC) to examine air quality by measuring concentrations of **particulate matter (PM)** data in different country-level administrative areas. We will focus on countries not meeting international health guidelines, and we will compare the PM concentration between areas with different levels of socioeconomic deprivation, a proxy for poverty, to determine if there is a correlation between exposure to air pollutants and poverty.\n",
        "\n",
        "This lesson walks through the process of calculating and visualizing zonal statistics for a set of countries using spatial raster data, with a focus on PM2.5 concentrations and **Global Gridded Relative Deprivation Index (GRDI)**, tailored to different socioeconomic areas within each country. We will subset the GRDI and PM2.5 raster data for each country and compute relevant zonal statistics such as mean, median, and various percentiles for PM2.5 concentrations. These statistics are stored as a *GeoDataFrame*, a data table format that can store spatial data, which is later used to create a choropleth map that visualizes specific PM2.5 concentration metrics across countries. The lesson includes a detailed analysis of PM2.5 concentrations within different GRDI quartiles for selected countries. This involves clipping the raster data to each country's geometry, filtering the data based on the GRDI quartiles, and calculating the mean PM2.5 levels for each quartile. The results are then visualized using customized plots to highlight the relationship between air quality and GRDI metrics across the selected countries.\n",
        "\n",
        "## Learning Objectives\n",
        "\n",
        "After completing this lesson, you should be able to:\n",
        "\n",
        "-   Gain a general understanding of what is particulate matter (PM) in the air and how it impacts human health.\n",
        "-   Learn about global socioeconomic dimensions of deprivation and how they are spatially represented.\n",
        "-   Find statistical thresholds in socioeconomic data.\n",
        "-   Perform zonal statistics to summarize spatial data\n",
        "-   Resample spatial data to harmoniza and compare socioeconomic data against environmental data.\n",
        "-   Display data on a map to get a general understanding of the spatial distribution of data.\n",
        "-   Summarize spatial data into table plots to compare how air quality differs in different socioeconomic conditions of international administrative areas.\n",
        "\n",
        "## Introduction\n",
        "\n",
        "As we have learned from previous lessons in this module, air quality is an important factor that if unmitigated can pose significant health risks, particularly to children’s health. In this module, we will use global particulate matter (PM) spatial data and GADM administrative boundaries to determine the average PM2.5 concentrations in each country and determine if they are meeting the [World Health Organization’s (WHO’s)](https://www.who.int/) air quality guidelines. Furthermore we will use spatial population data to normalize the air pollution concentrations per capita. We will subset 10 countries to analyze how areas with different levels of multidimensional deprivation, a proxy for poverty, compare with regards to PM2.5 concentrations.\n",
        "\n",
        "### Background\n",
        "\n",
        "Air pollution is now recognized as the single biggest environmental threat to human health. Air pollution affects different aspects of health and impacts everyone in low- middle- and high-income countries alike. It is estimated that pollution is responsible for 9 million deaths per year [@fuller2022]. The burden of disease attributable to air pollution is now on a par with other major global health risks such as unhealthy diet and tobacco smoking [@who2021; @who2022]. Air pollution also places risks upon several of the United Nations’ 2023 Sustainable Development Goals (SDGs), including SDG 3 (Good health and well-being), SDG 11 (Sustainable Communities), and SDG 15 (Life on Land), among others [@un_sdgs_2024]. Monitoring air quality is crucial for tracking the well-being of people around the world and for the advancement of global sustainable development policies that address these environmental and social risks.\n",
        "\n",
        "Particulate matter (PM₂.₅ and PM₁₀) can penetrate deep into the lungs and bloodstream, leading to respiratory and cardiovascular problems and being associated with conditions such as heart disease and lung cancer. [@who2022]. Benefits from improved air quality include prevention of air pollution-related premature deaths, chronic diseases, damages to ecosystems and crops, as well as the economic benefits for human health from air quality improvement [@calvin2023]. The United Nations Intergovernmental Panel on Climate Change’s (IPCC) latest Climate report acknowledges the largest adaptation gaps exist in projects that manage complex dynamics between air quality and climate risks, and that reductions in greenhouse gas (GHG) emissions would lead to improvements in air quality within a few years [@calvin2023]. Progress has been made to improve air quality, particularly in high-income countries, however; 99% of the world’s population is living in places where the World Health Organization’s (WHO’s) air quality guidelines are not being met [@who2021; @who2022].\n",
        "\n",
        "The WHO Global Air Quality Guidelines, first introduced in 1987, recommends safe levels for six key pollutants—PM₂.₅, PM₁₀, ozone (O₃), nitrogen dioxide (NO₂), sulfur dioxide (SO₂), and carbon monoxide (CO) which represent critical risks to human health globally. Recently updated guidelines respond to new scientific findings demonstrating that even lower concentrations of air pollutants than previously understood can lead to severe health impacts. Table 1 shows the latest Air Quality Guidelines (AQG) for PM2.5.\n",
        "\n",
        "We focus on PM2.5 because we will be working with the SEDAC Global Annual PM2.5 annual data.\n",
        "\n",
        "Table 1. WHO Recommended AQG levels and interim targets\n",
        "\n",
        "<table style=\"border-collapse: collapse; width: 100%;\">\n",
        "  <tr style=\"border: 1px solid black;\">\n",
        "    <th rowspan=\"2\" style=\"border: 1px solid black; padding: 8px;\">Pollutant</th>\n",
        "    <th rowspan=\"2\" style=\"border: 1px solid black; padding: 8px;\">Averaging Time</th>\n",
        "    <th colspan=\"4\" style=\"border: 1px solid black; padding: 8px;\">Interim Target</th>\n",
        "    <th rowspan=\"2\" style=\"border: 1px solid black; padding: 8px;\">AQG Level</th>\n",
        "  </tr>\n",
        "  <tr style=\"border: 1px solid black;\">\n",
        "    <th style=\"padding: 8px;\">1</th>\n",
        "    <th style=\"padding: 8px;\">2</th>\n",
        "    <th style=\"padding: 8px;\">3</th>\n",
        "    <th style=\"padding: 8px;\">4</th>\n",
        "  </tr>\n",
        "  <tr style=\"border: 1px solid black;\">\n",
        "    <td rowspan=\"2\" style=\"border: 1px solid black; padding: 8px;\">PM2.5, µg/m³</td>\n",
        "    <td style=\"border: 1px solid black; padding: 8px;\">Annual</td>\n",
        "    <td style=\"border: 1px solid black; padding: 8px;\">35</td>\n",
        "    <td style=\"border: 1px solid black; padding: 8px;\">25</td>\n",
        "    <td style=\"border: 1px solid black; padding: 8px;\">15</td>\n",
        "    <td style=\"border: 1px solid black; padding: 8px;\">10</td>\n",
        "    <td style=\"border: 1px solid black; padding: 8px;\">5</td>\n",
        "  </tr>\n",
        "  <tr style=\"border: 1px solid black;\">\n",
        "    <td style=\"border: 1px solid black; padding: 8px;\">24-hour*</td>\n",
        "    <td style=\"border: 1px solid black; padding: 8px;\">75</td>\n",
        "    <td style=\"border: 1px solid black; padding: 8px;\">50</td>\n",
        "    <td style=\"border: 1px solid black; padding: 8px;\">37.5</td>\n",
        "    <td style=\"border: 1px solid black; padding: 8px;\">25</td>\n",
        "    <td style=\"border: 1px solid black; padding: 8px;\">15</td>\n",
        "  </tr>\n",
        "</table>\n",
        "\n",
        "\\*99th percentile (i.e. 3-4 exceedance days per year).\n",
        "\n",
        "Inequality, political instability, and increased cost of living are some of the causes preventing significant progress to reduce poverty and meet the SDGs. This continued environmental inequality can also be seen in the significant majority (89%) of premature deaths attributed to air pollution, which occur in low- and middle- income countries (World Health Organization, 2022). In an account of 110 countries, 1.1 billion people (18%) are considered to live in “acute” multidimensional poverty, with half living in Sub-Saharan Africa and a third living in South Asia (UNDP (United Nations Development Programme), 2023). Studies have shown that persons with lower socioeconomic status disproportionately experience higher concentrations of pollutants (Flanagan et al., 2019). The level of poverty experienced by vulnerable communities can compound the risks posed by pollution leading to worsening well being and health outcomes.\n",
        "\n",
        "## Data Collection and Integration\n",
        "\n",
        "To determine whether countries are meeting WHO recommended AQG, we can use data of PM2.5 concentrations from the **Global (GL) Annual PM2.5 Grids from MODIS, MISR and SeaWiFS Aerosol Optical Depth (AOD), v4.03 (1998 – 2019)** dataset, which can can be downloaded from the Socioeconomic Data and Applications Center ([SEDAC](https://doi.org/10.7927/fx80-4n39))[@centerforinternationalearthscienceinformationnetwork-ciesin-columbiauniversity2022]. This data can be summarized using [GADM.org](https://gadm.org/) data, which provides administrative areas of many countries. We can determine the average concentration for the entire country area; however, many areas many not be populated by people, so we will use population data to subset areas where at least one person is living. This will also be used to normalize the data to a per-capita value and provide us with a better representation of the PM2.5 concentrations that people are experiencing.\n",
        "\n",
        "once we determine the top 10 countries with the highest average per capita PM2.5 concentrations, we will use those results to further classify the data into socioeconomic strata in each country. We can compare the PM2.5 per capita concentrations in these different areas to observe the correlation between socioeconomic poverty and air pollution. We employ the **Global Gridded Relative Deprivation Index (GRDI), v1 (2010 – 2020)** dataset can be downloaded from [SEDAC](https://doi.org/10.7927/3xxe-ap97) as well [@centerforinternationalearthscienceinformationnetwork-ciesin-columbiauniversity2022a].\n",
        "\n",
        "### Preparing Computing Environment and Variables\n",
        "\n",
        "Importing python packages required: The python packages are required for the remainder of the lesson. Please review the Python documentation of these packages for detailed information.\n"
      ],
      "id": "473048e9"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "import xarray as xr\n",
        "import rioxarray\n",
        "import rasterstats\n",
        "from rasterio.enums import Resampling\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import geopandas as gpd\n",
        "import pygadm\n",
        "import plotly.graph_objects as go "
      ],
      "id": "fdf78241",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Load the GRDIv1 and PM2.5 data from local sources:\n"
      ],
      "id": "dd6d1704"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Load rasters\n",
        "grdi_path = r\"Z:\\Sedac\\GRDI\\data\\povmap-grdi-v1-geotiff\\final data\\povmap-grdi-v1.tif\"\n",
        "pm25_path = r\"data\\sdei-global-annual-gwr-pm2-5-modis-misr-seawifs-aod-v4-gl-03-2019-geotiff\\sdei-global-annual-gwr-pm2-5-modis-misr-seawifs-aod-v4-gl-03-2019.tif\""
      ],
      "id": "2a168db6",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Using the package rasterio to load the data into memory. This allows us to read the data and use it for processing.\n"
      ],
      "id": "9707d4d9"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Open the input and reference rasters\n",
        "grdi_raster = rioxarray.open_rasterio(grdi_path, mask_and_scale=True)\n",
        "pm25_raster = rioxarray.open_rasterio(pm25_path, mask_and_scale=True)"
      ],
      "id": "482c1b36",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Matching Data Points using Bilinear Resample\n",
        "\n",
        "The GRDI raster and PM2.5 rasters are incompatible in resolution. One method of harmonizing data is by using the `Resampling` tool with a *bilinear* method. In this case, we reduce, or coarsen, the resolution of the GRDI raster to match the PM2.5 raster.\n"
      ],
      "id": "48668330"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Resample the input raster to match the reference raster\n",
        "grdi_raster = grdi_raster.rio.reproject_match(pm25_raster,resampling=Resampling.bilinear)"
      ],
      "id": "01a72493",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Previewing Spatial Data in a Plot\n",
        "\n",
        "Once the data rasters have been resampled to the same resolution, we can quickly plot them to view the data we are working with on a map. "
      ],
      "id": "128ecf38"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Plotting the rasters\n",
        "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 20))\n",
        "\n",
        "# Plot the original GRDI raster in the first subplot\n",
        "im1 = ax1.imshow(grdi_raster.values[0], cmap='viridis', interpolation='nearest')\n",
        "ax1.set_title('Original GRDI Raster')\n",
        "fig.colorbar(im1, ax=ax1, orientation='horizontal', label='GRDI Values')\n",
        "\n",
        "# Plot the PM2.5 raster in the second subplot\n",
        "im2 = ax2.imshow(pm25_raster.values[0], cmap='hot', interpolation='nearest')\n",
        "ax2.set_title('PM2.5 Raster')\n",
        "fig.colorbar(im2, ax=ax2, orientation='horizontal', label='PM2.5 Values')\n",
        "\n",
        "\n",
        "# Show the plots\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "id": "17d62105",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Working with administrative Data\n",
        "\n",
        "`pygadm` is a package that has international administrative units from levels 0 to 2. We can search the available countries by listing the Names.\n"
      ],
      "id": "ce43012e"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "country_table = gpd.GeoDataFrame(pygadm.Names())\n",
        "len(country_table)"
      ],
      "id": "fdc91810",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Some available areas with a unique `GID_0` code share Names; therefore we drop the rows that contain digits.\n"
      ],
      "id": "30a640cc"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "country_table = country_table[~country_table['GID_0'].str.contains('\\d', na=False)]\n",
        "len(country_table)"
      ],
      "id": "d9804460",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Subset Data From a Table\n",
        "\n",
        "Doing Zonal statistics for more than 200 countries may take a while, therefore, we can subset the data randomly with the `.sample()` method.\n"
      ],
      "id": "5cd12c20"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "country_sample = country_table.sample(n=30)\n",
        "country_sample"
      ],
      "id": "ea6f8e1c",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Zonal Statistics for Each Administrative Area\n",
        "\n",
        "`rasterstats` has a funcion `zonal_stats()` which allows you to use vectors to summarize raster data. We summarize GRDIv1 data to calculate the following statistics: count, minimum, mean, max, median, standard deviation, range, and percentiles 20, 40, 60, and 80.\n",
        "\n",
        "\n",
        "We can **define a custom function** that can allow us to use the zonal statistics process multiple times. A custom function can be created using the `def FUNCTION_NAME(PARAMETER1, PARAMETER2):` fuction to define what the fucntion will do.\n"
      ],
      "id": "7817141f"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "def calculate_country_stats(country_sample, grdi_raster, pm25_raster=None):\n",
        "    \"\"\"\n",
        "    Calculate statistics for each country in the sample.\n",
        "\n",
        "    Parameters:\n",
        "    - country_sample: A pandas DataFrame containing country information with 'NAME_0' and 'GID_0' columns, in this case the country_table.\n",
        "    - grdi_raster: A raster object with which to perform the zonal statistics.\n",
        "    - pm25_raster: (Optional) A raster object for PM2.5 data. If provided, statistics will also be calculated for this raster.\n",
        "\n",
        "    Returns:\n",
        "    - stats_results: A GeoDataFrame containing the statistics for each country.\n",
        "    \"\"\"\n",
        "    stats_results = gpd.GeoDataFrame()\n",
        "\n",
        "    for index, row in country_sample.iloc[:].iterrows():\n",
        "        country = row['NAME_0']\n",
        "        country_GID = row['GID_0']\n",
        "        try:\n",
        "            country_poly = pygadm.Items(admin=country_GID, content_level=0)\n",
        "            if isinstance(country_poly, gpd.GeoDataFrame):\n",
        "                country_geometry = country_poly.geometry.iloc[0]  # Ensure single geometry\n",
        "            else:\n",
        "                country_geometry = country_poly\n",
        "        except Exception as e:\n",
        "            print(country, \"skipped due to error:\", e)\n",
        "            continue\n",
        "\n",
        "        # Create a mask for the polygons and perform zonal statistics on GRDI raster\n",
        "        grdi_country_zs = rasterstats.zonal_stats(\n",
        "            country_geometry, grdi_raster.values[0], \n",
        "            affine=grdi_raster.rio.transform(), \n",
        "            stats=\" min mean max percentile_20 percentile_40 percentile_60 percentile_80\",\n",
        "            nodata_value=grdi_raster.rio.nodata\n",
        "        )\n",
        "\n",
        "        # Uncomment and update the following lines if you want to include PM2.5 statistics\n",
        "        if pm25_raster is not None:\n",
        "            pm25_country_zs = rasterstats.zonal_stats(\n",
        "                country_geometry, pm25_raster.values[0], \n",
        "                affine=pm25_raster.rio.transform(), \n",
        "                stats=\"mean\",\n",
        "                nodata_value=pm25_raster.rio.nodata\n",
        "                )\n",
        "\n",
        "        # Extract statistics into a dictionary\n",
        "        country_stats = {\n",
        "            'Country_Name': country,\n",
        "            'Country_GID' : country_GID,\n",
        "            # 'GRDI_Count': grdi_country_zs[0]['count'],\n",
        "            'GRDI_Min': grdi_country_zs[0]['min'],\n",
        "            'GRDI_Mean': grdi_country_zs[0]['mean'],\n",
        "            'GRDI_Max': grdi_country_zs[0]['max'],\n",
        "            # 'GRDI_Median': grdi_country_zs[0]['median'],\n",
        "            # 'GRDI_Std': grdi_country_zs[0]['std'],\n",
        "            # 'GRDI_Range': grdi_country_zs[0]['range'],\n",
        "            'GRDI_P20': grdi_country_zs[0]['percentile_20'],\n",
        "            'GRDI_P40': grdi_country_zs[0]['percentile_40'],\n",
        "            'GRDI_P60': grdi_country_zs[0]['percentile_60'],\n",
        "            'GRDI_P80': grdi_country_zs[0]['percentile_80']#,\n",
        "            # 'geometry' : country_poly['geometry'].iloc[0]\n",
        "        }\n",
        "\n",
        "        # If PM2.5 statistics are calculated, add them to the dictionary\n",
        "        if pm25_raster is not None:\n",
        "            country_stats.update({\n",
        "                'PM25_Mean': pm25_country_zs[0]['mean']\n",
        "            })\n",
        "        \n",
        "        # Filter out None values from `country_stats`\n",
        "        country_stats = {k: v for k, v in country_stats.items() if pd.notnull(v)}\n",
        "\n",
        "        country_stats_gdf = gpd.GeoDataFrame([country_stats], geometry=[country_geometry])\n",
        "\n",
        "        # Append to results if country_stats_gdf has non-empty columns\n",
        "        if not country_stats_gdf.empty and country_stats_gdf.notna().any().any():\n",
        "            stats_results = pd.concat([stats_results, country_stats_gdf], ignore_index=True)\n",
        "\n",
        "    # Drop any rows in `stats_results` that have NaN in any cell\n",
        "    stats_results = stats_results.dropna()\n",
        "    # Reset the index\n",
        "    stats_results = stats_results.reset_index(drop=True)\n",
        "    return stats_results\n"
      ],
      "id": "07c9e489",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's put our defined function to use:\n"
      ],
      "id": "f28e6230"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "stats_results = calculate_country_stats(country_sample, grdi_raster, pm25_raster)"
      ],
      "id": "49425d47",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's use the `.head()` method from Pandas to check the top of our table\n"
      ],
      "id": "a1d4090b"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "stats_results.head()"
      ],
      "id": "618dfc4b",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "So far, we have collected the overall average PM2.5 for the countries. Based on the WHO Air Quality guidelines for PM2.5, we subset the dataframe to keep only countries with the mean PM2.5 value greater than or equal to 5.0:\n"
      ],
      "id": "f589351f"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "stats_results = stats_results[stats_results[\"PM25_Mean\"] >= 5.0]\n",
        "stats_results"
      ],
      "id": "acf6eb8c",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Assuming the column with names is called 'Name'\n",
        "# Define the bins and labels\n",
        "bins = [5, 10, 15, 25, 30, float('inf')]\n",
        "labels = [\"5-10\", \"10-15\", \"15-25\", \"25-30\", \">30\"]\n",
        "\n",
        "# Use .loc to add the new column with binned ranges, avoiding SettingWithCopyWarning\n",
        "stats_results.loc[:, 'PM25_range'] = pd.cut(stats_results[\"PM25_Mean\"], bins=bins, labels=labels, right=False).astype(str)\n",
        "\n",
        "# Group by the 'PM25_range' and list names in each range\n",
        "grouped = stats_results.groupby('PM25_range', observed=True)['Country_Name'].agg(['count', list]).reset_index()\n",
        "grouped.columns = ['PM25_range', 'Count', 'Names']  # Rename columns for clarity\n",
        "\n",
        "# Set display option to prevent truncation of long lists in the output\n",
        "pd.set_option('display.max_colwidth', None)\n",
        "\n",
        "print(grouped)\n",
        "\n",
        "# Reset Pandas to default display options\n",
        "pd.reset_option(\"display.max_colwidth\")\n"
      ],
      "id": "ff5c45d5",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "From the table above, we can see how many and which countries are within the WHO AQG, or are within interim stages of air quality improvement. \n",
        "\n",
        "Below, choose an attribute, or column, to display it in a map plot. In this case, I'm choosing the PM2.5 Mean.\n"
      ],
      "id": "a30036b6"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "column_chosen = 'PM25_Mean' # 'GRDI_Max' #GRDI_Max, GRDI_Min, GRDI_Median\n",
        "# Plotting\n",
        "fig, ax = plt.subplots(1, 1, figsize=(15, 10))\n",
        "stats_results.plot(column=column_chosen, ax=ax, legend=True,\n",
        "    legend_kwds={'label': f\"{column_chosen} per country.\",\n",
        "                      'orientation': \"horizontal\"})\n",
        "ax.set_title(f'Choropleth Map Showing {column_chosen} per country')\n",
        "ax.set_axis_off()  # Turn off the axis numbers and ticks\n",
        "plt.show()"
      ],
      "id": "9badbf26",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Selecting Data by Column\n",
        "\n",
        "Start my creating a list of countries that you are interested in to Subset data from the DataFrame that match the values in the `NAME_0` column. The `.isin()` mehthod checks each element in the DataFrame's column for the item present in the list and returns matching rows.\n"
      ],
      "id": "318f87eb"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# selected_countries = [\"Algeria\", \"Somalia\", \"Colombia\", \"Timor Leste\", \"Finland\", \"Nicaragua\", \"United Kingdom\", \"Mali\"]\n",
        "# selected_countries = [\"Anguilla\", \"Armenia\", \"Angola\", \"Argentina\", \"Albania\", \"United Arab Emirates\", \"American Samoa\", \"Australia\" ]\n",
        "selected_countries = [\"Algeria\", \"Somalia\", \"Colombia\", \"Timor Leste\", \"Finland\", \"Nicaragua\", \"United Kingdom\", \"Mali\", \"Armenia\", \"Argentina\",  \"Albania\", \"United Arab Emirates\", \"Indonesia\", \"Qatar\"]\n",
        "\n",
        "#use the list above to subset the country_table DataFrame by the column NAME_0 \n",
        "selected_countries = country_table[country_table['NAME_0'].isin(selected_countries)]"
      ],
      "id": "439417e2",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Using a Defined Custom Function\n",
        "\n",
        "Recalling the defined fucntion `calculate_country_stats`, we can use our `selected_countries` list, and the GRDI and PM2.5 rasters, to create a new table of zonal statistics.\n"
      ],
      "id": "aca2c1ef"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "stats_results = calculate_country_stats(selected_countries, grdi_raster, pm25_raster)"
      ],
      "id": "67b7f21e",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Show the head of the table again:\n"
      ],
      "id": "efc15ece"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "stats_results.head()"
      ],
      "id": "f8b44cc6",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the map again choosing a column to plot:\n"
      ],
      "id": "deca2841"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "column_chosen = 'PM25_Mean'  #'GRDI_Max' #GRDI_Max, GRDI_Min, GRDI_Median\n",
        "stats_results.plot(column=column_chosen, legend=True)\n",
        "plt.show() "
      ],
      "id": "d400cb35",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Creating a Table with Results\n",
        "\n",
        "We can create a list of **tuples** that we can use to refer to the statistical values, and the name, color, and symbol we want to assign. In this case, we are using the GRDI zonal statistics of each country we selected that include the Mean, Minimum, Maximum, and interquartiles.\n"
      ],
      "id": "8bd78d02"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# List of GRDI values and their corresponding properties\n",
        "#column, value name, color, symbol\n",
        "grdi_data = [\n",
        "    ('GRDI_Mean', 'Mean', 'orange', 'diamond'),\n",
        "    ('GRDI_Min', 'Min', 'gray', '152'),\n",
        "    ('GRDI_Max', 'Max', 'gray', '151'),\n",
        "    ('GRDI_P20', 'Q20', 'blue', '142'),\n",
        "    ('GRDI_P40', 'Q40', 'purple', '142'),\n",
        "    ('GRDI_P60', 'Q60', 'green', '142'),\n",
        "    ('GRDI_P80', 'Q80', 'red', '142')\n",
        "]"
      ],
      "id": "82802c8b",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can create a figure to display the data based on the names colors and symbols we selected.\n"
      ],
      "id": "40286e6a"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Create a figure\n",
        "fig = go.Figure()\n",
        "\n",
        "# Add traces to the figure based on the data\n",
        "for col, name, color, symbol in grdi_data:\n",
        "    fig.add_trace(go.Scatter(\n",
        "        x=stats_results[col],\n",
        "        y=stats_results['Country_Name'],\n",
        "        mode='markers',\n",
        "        name=name,\n",
        "        marker=dict(color=color, size=10, symbol=symbol)\n",
        "    ))\n",
        "\n",
        "# Customize layout\n",
        "fig.update_layout(\n",
        "    title='GRDI Statistics by Country',\n",
        "    xaxis_title='GRDI Values',\n",
        "    yaxis_title='Country Name',\n",
        "    yaxis=dict(tickmode='linear'),\n",
        "    legend_title='Statistics',\n",
        "    yaxis_type='category',\n",
        "    xaxis=dict(tickvals=[0, 20, 40, 60, 80, 100])\n",
        ")\n",
        "\n",
        "# Show plot\n",
        "fig.show()"
      ],
      "id": "f2a418ec",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Summarizing PM2.5 Values by Socioeconomic Deprivation\n",
        "\n",
        "Considering the GRDI quartile values as a level of socieoeconomic deprivation within each country, we can use the `stats_results` GeoDataFrame, the GRDI raster, and the PM2.5 raster to calculate the Mean PM.25 value within each of those areas in each country. This can describe how the air quality for different socioeconomic strata compare within the country, as well as against other countries.\n",
        "\n",
        "The results will be added to the `stats_results` with the corresponting columns.\n"
      ],
      "id": "b38e1c2d"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# iterate through the stats_results table rows\n",
        "for index, row in stats_results.iloc[:].iterrows():\n",
        "    #isolate each country's respective row\n",
        "    row_df = gpd.GeoDataFrame([row], geometry='geometry').reset_index(drop=True)\n",
        "    print(row_df.loc[0,'Country_GID'])\n",
        "    try:\n",
        "        #use rioxarray to clip the GRDI and PM2.5 rasters by the geometry of the respective country.\n",
        "        grdi_country = grdi_raster.rio.clip(row_df.geometry, grdi_raster.rio.crs)\n",
        "        pm25_country = pm25_raster.rio.clip(row_df.geometry, grdi_raster.rio.crs)\n",
        "    except:\n",
        "        print('Error in clip')\n",
        "        continue\n",
        "\n",
        "    #Applying squeeze() to this array removes the singleton dimension, reducing it to a 2D array with dimensions (rows, columns)\n",
        "    grdi_country= grdi_country.squeeze()\n",
        "    pm25_country= pm25_country.squeeze()\n",
        "\n",
        "\n",
        "    # Subset the GRDI raster where values fall between each GRDI quintiles\n",
        "    grdi_countryQ1 = grdi_country.where((grdi_country >= row_df.loc[0, 'GRDI_Min']) & (grdi_country <= row_df.loc[0, 'GRDI_P20']))\n",
        "    grdi_countryQ2 = grdi_country.where((grdi_country >= row_df.loc[0, 'GRDI_P20']) & (grdi_country <= row_df.loc[0, 'GRDI_P40']))\n",
        "    grdi_countryQ3 = grdi_country.where((grdi_country >= row_df.loc[0, 'GRDI_P40']) & (grdi_country <= row_df.loc[0, 'GRDI_P60']))\n",
        "    grdi_countryQ4 = grdi_country.where((grdi_country >= row_df.loc[0, 'GRDI_P60']) & (grdi_country <= row_df.loc[0, 'GRDI_P80']))\n",
        "    grdi_countryQ5 = grdi_country.where((grdi_country >= row_df.loc[0, 'GRDI_P80']) & (grdi_country <= row_df.loc[0, 'GRDI_Max']))\n",
        "\n",
        "\n",
        "    # Mask the PM2.5 raster using the above GRDI quartile rasters, keeping only the cells that intersect\n",
        "    pm25_countryQ1 = pm25_country.where(grdi_countryQ1.notnull())\n",
        "    pm25_countryQ2 = pm25_country.where(grdi_countryQ2.notnull())\n",
        "    pm25_countryQ3 = pm25_country.where(grdi_countryQ3.notnull())\n",
        "    pm25_countryQ4 = pm25_country.where(grdi_countryQ4.notnull())\n",
        "    pm25_countryQ5 = pm25_country.where(grdi_countryQ5.notnull())\n",
        "\n",
        "    #Find the mean value of of the intersected PM2.5 rasters in each quartile\n",
        "    pm25_countryQ1v = pm25_countryQ1.mean().item()\n",
        "    pm25_countryQ2v = pm25_countryQ2.mean().item()\n",
        "    pm25_countryQ3v = pm25_countryQ3.mean().item()\n",
        "    pm25_countryQ4v = pm25_countryQ4.mean().item()\n",
        "    pm25_countryQ5v = pm25_countryQ5.mean().item()\n",
        "\n",
        "    #add the resuts to the stats_results table in the respective column\n",
        "    stats_results.at[index, 'PM25_Q1'] = pm25_countryQ1v\n",
        "    stats_results.at[index, 'PM25_Q2'] = pm25_countryQ2v\n",
        "    stats_results.at[index, 'PM25_Q3'] = pm25_countryQ3v\n",
        "    stats_results.at[index, 'PM25_Q4'] = pm25_countryQ4v\n",
        "    stats_results.at[index, 'PM25_Q5'] = pm25_countryQ5v"
      ],
      "id": "719ad519",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot Results of Mean PM2.5 in Socieceonomic Deprivation Quartiles per country\n",
        "\n",
        "Similarly, we create a list of tuples of how we want to display the data, and create a figure based on the tuples. This plot would show each country in the y axis and the Log of Mean PM2.5 values in each country's GRDI quarties.\n"
      ],
      "id": "79ae6f31"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# List of GRDI values and their corresponding properties\n",
        "#column, value name, color, symbol\n",
        "plot_data =[\n",
        "    ('PM25_Q1', 'Q1', '#440154', '6'),  # Light Blue\n",
        "    ('PM25_Q2', 'Q2', '#31688E', '5'),  # Light Green\n",
        "    ('PM25_Q3', 'Q3', '#35B779', '7'),  # Yellow\n",
        "    ('PM25_Q4', 'Q4', '#FDE725', '8'),  # Orange\n",
        "    ('PM25_Q5', 'Q5', '#FF0000', '1')   # Red\n",
        "]\n",
        "\n",
        "# Create a figure\n",
        "fig = go.Figure()\n",
        "\n",
        "# Add traces to the figure.\n",
        "for col, name, color, symbol in plot_data:\n",
        "    xlog  = np.log(stats_results[col])\n",
        "    fig.add_trace(go.Scatter(\n",
        "        x=stats_results[col], #xlog,\n",
        "        y=stats_results['Country_Name'],\n",
        "        mode='markers+text',  # Add 'text' to mode\n",
        "        text=[f'<b>{name}</b>' for _ in stats_results[col]],  # Repeat name for each point\n",
        "        name=name,\n",
        "        textfont=dict(color=color, size=12),\n",
        "        textposition='top center',  # Position the text above the symbol\n",
        "        marker_color=color,\n",
        "        marker_line_color=\"midnightblue\",\n",
        "        marker_symbol=symbol,\n",
        "        marker_size=14,\n",
        "        marker_line_width=2,\n",
        "        marker_opacity=0.6\n",
        "        ))\n",
        "fig.update_traces(textposition='top center')\n",
        "\n",
        "    # Customize layout\n",
        "fig.update_layout(\n",
        "    title='Mean PM2.5 in each GRDI Quartile by Country',\n",
        "    xaxis_title='PM2.5 Mean Values',\n",
        "    yaxis_title='Country Name',\n",
        "    yaxis=dict(tickmode='linear'),\n",
        "    legend_title='Statistics',\n",
        "    yaxis_type='category',\n",
        "    xaxis=dict(rangemode=\"tozero\"),\n",
        "    \n",
        "    #xaxis=dict(tickvals=[0, 20, 40, 60, 80, 100])\n",
        "    )\n",
        "\n",
        "# Show plot\n",
        "fig.show()"
      ],
      "id": "9074d2d9",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Use the plotly controls to take a closer look at the results.\n",
        "\n",
        "With this shapely plot, We can examine differences between countries and PM2.5 values. The plot displays the coutnries on the Y axis and log values of the average PM2.5 value on the X axis. Each country displays PM2.5 values averaged within the quartile areas based on GRDI values of each country. A higher quartile (Q) implies a higher degree of deprivation, 1 being the lowest and 5 the highest.\n",
        "\n",
        "::: {.callout-tip style=\"color: #9f0b64;\"}\n",
        "#### Knowledge Check\n",
        "1.  Which countries did you identify to be over the WHO AQG values? \n",
        "2.  Which countries have higher contentrations of air pollution in lower socioecomnomic strata? \n",
        "3.  Which are the opposite?\n",
        "\n",
        ":::\n",
        "\n",
        "Congratulations! .... Now you should be able to:\n",
        "\n",
        "-   Understand particulate matter (PM2.5) in the air and its effects on human health.\n",
        "-   Explore the global socioeconomic dimensions of deprivation and their spatial representations.\n",
        "-   Identify statistical thresholds in socioeconomic datasets.\n",
        "-   Use zonal statistics to summarize spatial data.\n",
        "-   Create maps to visualize spatial data distribution.\n",
        "-   Resample spatial data to harmonize the compare socioeconomic data with environmental data.\n",
        "-   Summarize spatial data in tables to analyze how air quality differs across different socioeconomic conditions in administrative areas around the world.\n",
        "\n",
        "\n",
        "## Module 2: Air Quality Home\n",
        "\n",
        "In this lesson, we explored ....\n",
        "\n",
        "[Module 2: Air Quality](https://ciesin-geospatial.github.io/TOPSTSCHOOL-air-quality/index.html){.btn .btn-primary .btn role=\"button\"}"
      ],
      "id": "d7ee025b"
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "language": "python",
      "display_name": "Python 3 (ipykernel)"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}