import xarray as xr
import rioxarray
import rasterstats
from rasterio.enums import Resampling
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import geopandas as gpd
import pygadm
import plotly.graph_objects as go
Particulate Matter Across Socioeconomic Strata of Countries
Overview
In this lesson, you will use NASA socioeconomic and environmental Earthdata available at NASA SEDAC to compare relationships between levels of socioeconomic deprivation agaisnts air quality data of particulate matter (PM) in different international administrative areas.
This lesson walks through the process of calculating and visualizing zonal statistics for a set of countries using raster data, focusing on GRDI country quintiles and PM2.5 concentration levels within these quintile areas. It begins by subsetting data by country and iterating over each country to extract relevant zonal statistics like mean, median, and various percentiles for each quintile. These statistics are stored in a GeoDataFrame, which is later used to create a choropleth map that visualizes specific GRDI 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.
Learning Objectives
After completing this lesson, you should be able to:
- Gain a general understanding of what is particulate matter (PM) in the air and how it impacts human health.
- Learn about global socioeconomic dimensions of deprivation and how they are spatially represented.
- Find statistical thresholds in socioeconomic data.
- Perform zonal statistics to summarize spatial data
- Resample spatial data to harmoniza and compare socioeconomic data against environmental data.
- Display data on a maps to get a general understanding of the spatial distribution of data.
- Summarize spatial data into table plots to compare how air quality differs in different socioeconomic conditions of international administrative areas.
Introduction
Data Collection and Integration
The Global (GL) Annual PM2.5 Grids from MODIS, MISR and SeaWiFS Aerosol Optical Depth (AOD), v4.03 (1998 – 2019) can can be downloaded from the Socioeconomic Data and Applications Center ([SEDAC](https://sedac.ciesin.columbia.edu/)) (Center For International Earth Science Information Network-CIESIN-Columbia University 2022a).
The Global Gridded Relative Deprivation Index (GRDI), v1 (2010 – 2020) dataset can be downloaded from SEDAC as well (Center For International Earth Science Information Network-CIESIN-Columbia University 2022b).
Gather comprehensive datasets from reliable sources such as the US EPA’s EJSCREEN tool and other public health databases like John Hopkins and County Health Rankings. Integrate relevant environmental data (e.g., air pollutant concentrations, pollution source proximity) with health outcomes data (e.g., COVID-19 prevalence, NSCLC incidence). Ensure data compatibility and quality through data cleaning and validation procedures.
Preparing Environment and Variables
Importing python packages required:
Load the GRDIv1 and PM2.5 data from local sources:
# Load rasters
= r"Z:\Sedac\GRDI\data\povmap-grdi-v1-geotiff\final data\povmap-grdi-v1.tif"
grdi_path = r"F:\TOPSSCHOOL\git\TOPSTSCHOOL-air-quality\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" pm25_path
Using the package rasterio to load the data into memory. This allows us to read the data and use it for processing.
# Open the input and reference rasters
= rioxarray.open_rasterio(grdi_path, mask_and_scale=True)
grdi_raster = rioxarray.open_rasterio(pm25_path, mask_and_scale=True) pm25_raster
Matching Data Points using Bilinear Resample
The GRDI raster and PM2.5 rasters are incompatible in resolution. One method of harmonizing data is by using the Resampling
bethod with a bilinear method. In this case, we reduce, or coarsen, the resolution of the GRDI raster to match the PM2.5 raster.
# Resample the input raster to match the reference raster
= grdi_raster.rio.reproject_match(pm25_raster,resampling=Resampling.bilinear) grdi_raster
Previewing Spatial Data in a Plot
# Plotting the rasters
= plt.subplots(2, 1, figsize=(10, 20))
fig, (ax1, ax2)
# Plot the original GRDI raster in the first subplot
= ax1.imshow(grdi_raster.values[0], cmap='viridis', interpolation='nearest')
im1 'Original GRDI Raster')
ax1.set_title(=ax1, orientation='horizontal', label='GRDI Values')
fig.colorbar(im1, ax
# Plot the PM2.5 raster in the second subplot
= ax2.imshow(pm25_raster.values[0], cmap='hot', interpolation='nearest')
im2 'PM2.5 Raster')
ax2.set_title(=ax2, orientation='horizontal', label='PM2.5 Values')
fig.colorbar(im2, ax
# Show the plots
plt.tight_layout() plt.show()
Working with administrative Data
pygadm
is a package that has international administrative units from levels 0 to 2. We can search the available countries by listing the Names.
= gpd.GeoDataFrame(pygadm.Names())
country_table len(country_table)
263
Some available areas with a unique GID_0
code share Names; therefore we drop the rows that contain digits.
= country_table[~country_table['GID_0'].str.contains('\d', na=False)]
country_table len(country_table)
254
Subset Data From a Table
Doing Zonal statistics for more than 200 countries may take a while, therefore, we can subset the data randomly with the .sample()
method.
= country_table.sample(n=15)
country_sample country_sample
NAME_0 | GID_0 | |
---|---|---|
127 | Kiribati | KIR |
80 | Micronesia | FSM |
177 | Oman | OMN |
262 | Zimbabwe | ZWE |
63 | Djibouti | DJI |
9 | Armenia | ARM |
249 | Vanuatu | VUT |
37 | Botswana | BWA |
196 | Russia | RUS |
164 | Mayotte | MYT |
261 | Northern Cyprus | ZNC |
14 | Australia | AUS |
202 | South Georgia and the South Sand | SGS |
221 | Syria | SYR |
129 | South Korea | KOR |
Zonal Statistics for Each Administrative Area
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.
= gpd.GeoDataFrame()
stats_results
for index, row in country_sample.iloc[:].iterrows():
= row['NAME_0']
country = row['GID_0']
country_GID try:
= pygadm.Items(admin=country_GID, content_level=0)
country_poly except:
print(country, " skipped.")
continue
# Create a mask for the polygons
= rasterstats.zonal_stats(country_poly, grdi_raster.values[0], affine=grdi_raster.rio.transform(), stats="count min mean max median std median range percentile_20 percentile_40 percentile_60 percentile_80")
grdi_country_zs# # pm25_country_zs= rasterstats.zonal_stats(country_poly, pm25_arr, affine=pm25_transform, stats="count min mean max median std median range percentile_20 percentile_40 percentile_60 percentile_80", nodata=pm25_raster.nodata)
# # Extract statistics into a dictionary
= {
country_stats 'Country_Name': country,
'Country_GID' : country_GID,
'GRDI_Count': grdi_country_zs[0]['count'],
'GRDI_Min': grdi_country_zs[0]['min'],
'GRDI_Mean': grdi_country_zs[0]['mean'],
'GRDI_Max': grdi_country_zs[0]['max'],
'GRDI_Median': grdi_country_zs[0]['median'],
'GRDI_Std': grdi_country_zs[0]['std'],
'GRDI_Range': grdi_country_zs[0]['range'],
'GRDI_P20': grdi_country_zs[0]['percentile_20'],
'GRDI_P40': grdi_country_zs[0]['percentile_40'],
'GRDI_P60': grdi_country_zs[0]['percentile_60'],
'GRDI_P80': grdi_country_zs[0]['percentile_80'],
# # 'PM25_Count': pm25_country_zs[0]['count'],
# # 'PM25_Min': pm25_country_zs[0]['min'],
# 'PM25_Mean': pm25_country_zs[0]['mean'],
# # 'PM25_Max': pm25_country_zs[0]['max'],
# # 'PM25_Median': pm25_country_zs[0]['median'],
# # 'PM25_Std': pm25_country_zs[0]['std'],
# # 'PM25_Range': pm25_country_zs[0]['range'],
# # 'PM25_P20': pm25_country_zs[0]['percentile_20'],
# # 'PM25_P40': pm25_country_zs[0]['percentile_40'],
# # 'PM25_P60': pm25_country_zs[0]['percentile_60'],
# # 'PM25_P80': pm25_country_zs[0]['percentile_80'],
'geometry' : country_poly['geometry'].iloc[0]
}= gpd.GeoDataFrame([country_stats], geometry='geometry')
country_stats_gdf # stats_results.append(country_stats_gdf)
= pd.concat([stats_results, country_stats_gdf], ignore_index=True) stats_results
Let’s use the .head()
method from Pandas to check the top of our table
stats_results.head()
Country_Name | Country_GID | GRDI_Count | GRDI_Min | GRDI_Mean | GRDI_Max | GRDI_Median | GRDI_Std | GRDI_Range | GRDI_P20 | GRDI_P40 | GRDI_P60 | GRDI_P80 | geometry | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | Kiribati | KIR | 116 | 55.782413 | 80.730006 | 84.099190 | 83.011345 | 5.609290 | 28.316776 | 79.892601 | 82.550774 | 83.338074 | 83.656990 | MULTIPOLYGON (((-143.0905 -17.8498, -143.0858 ... |
1 | Micronesia | FSM | 345 | 32.953396 | 74.996014 | 85.425018 | 77.590477 | 7.864845 | 52.471622 | 71.955696 | 76.526154 | 78.361168 | 79.964798 | MULTIPOLYGON (((154.7811 1.0278, 154.7797 1.02... |
2 | Oman | OMN | 38957 | 5.814715 | 58.435711 | 75.041962 | 64.611206 | 14.829892 | 69.227246 | 48.413792 | 61.939667 | 66.912010 | 68.876038 | MULTIPOLYGON (((54.806 16.9438, 54.8076 16.942... |
3 | Zimbabwe | ZWE | 180198 | 24.202816 | 89.134846 | 97.393387 | 89.752472 | 3.847163 | 73.190571 | 87.505157 | 88.919373 | 90.362289 | 91.690750 | MULTIPOLYGON (((32.7043 -18.9602, 32.7081 -18.... |
4 | Djibouti | DJI | 781 | 26.355202 | 77.480739 | 94.670731 | 85.478348 | 15.927292 | 68.315529 | 70.784187 | 84.738541 | 85.735367 | 85.902008 | MULTIPOLYGON (((43.0834 11.182, 43.0013 11.051... |
Defining a Funtion
We can create 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.
def calculate_country_stats(country_sample, grdi_raster, pm25_raster=None):
"""
Calculate statistics for each country in the sample.
Parameters:
- country_sample: A pandas DataFrame containing country information with 'NAME_0' and 'GID_0' columns, in this case the country_table.
- grdi_raster: A raster object with which to perform the zonal statistics.
- pm25_raster: (Optional) A raster object for PM2.5 data. If provided, statistics will also be calculated for this raster.
Returns:
- stats_results: A GeoDataFrame containing the statistics for each country.
"""
= gpd.GeoDataFrame()
stats_results
for index, row in country_sample.iloc[:].iterrows():
= row['NAME_0']
country = row['GID_0']
country_GID try:
= pygadm.Items(admin=country_GID, content_level=0)
country_poly except Exception as e:
print(country, "skipped due to error:", e)
continue
# Create a mask for the polygons and perform zonal statistics on GRDI raster
= rasterstats.zonal_stats(
grdi_country_zs 0],
country_poly, grdi_raster.values[=grdi_raster.rio.transform(),
affine="count min mean max median std range percentile_20 percentile_40 percentile_60 percentile_80"
stats
)
# Uncomment and update the following lines if you want to include PM2.5 statistics
# if pm25_raster is not None:
# pm25_country_zs = rasterstats.zonal_stats(
# country_poly, pm25_raster.values[0],
# affine=pm25_raster.rio.transform(),
# stats="count min mean max median std range percentile_20 percentile_40 percentile_60 percentile_80",
# nodata=pm25_raster.nodata
# )
# Extract statistics into a dictionary
= {
country_stats 'Country_Name': country,
'Country_GID' : country_GID,
'GRDI_Count': grdi_country_zs[0]['count'],
'GRDI_Min': grdi_country_zs[0]['min'],
'GRDI_Mean': grdi_country_zs[0]['mean'],
'GRDI_Max': grdi_country_zs[0]['max'],
'GRDI_Median': grdi_country_zs[0]['median'],
'GRDI_Std': grdi_country_zs[0]['std'],
'GRDI_Range': grdi_country_zs[0]['range'],
'GRDI_P20': grdi_country_zs[0]['percentile_20'],
'GRDI_P40': grdi_country_zs[0]['percentile_40'],
'GRDI_P60': grdi_country_zs[0]['percentile_60'],
'GRDI_P80': grdi_country_zs[0]['percentile_80'],
'geometry' : country_poly['geometry'].iloc[0]
}
# If PM2.5 statistics are calculated, add them to the dictionary
# if pm25_raster is not None:
# country_stats.update({
# 'PM25_Count': pm25_country_zs[0]['count'],
# 'PM25_Min': pm25_country_zs[0]['min'],
# 'PM25_Mean': pm25_country_zs[0]['mean'],
# 'PM25_Max': pm25_country_zs[0]['max'],
# 'PM25_Median': pm25_country_zs[0]['median'],
# 'PM25_Std': pm25_country_zs[0]['std'],
# 'PM25_Range': pm25_country_zs[0]['range'],
# 'PM25_P20': pm25_country_zs[0]['percentile_20'],
# 'PM25_P40': pm25_country_zs[0]['percentile_40'],
# 'PM25_P60': pm25_country_zs[0]['percentile_60'],
# 'PM25_P80': pm25_country_zs[0]['percentile_80'],
# })
= gpd.GeoDataFrame([country_stats], geometry='geometry')
country_stats_gdf = pd.concat([stats_results, country_stats_gdf], ignore_index=True)
stats_results
return stats_results
From the table above, we can choose an attribute, or column, to display it in a map plot. In this case, I’m choosing the GRDI Max
= 'GRDI_Max' #GRDI_Max, GRDI_Min, GRDI_Median
column_chosen # Plotting
= plt.subplots(1, 1, figsize=(15, 10))
fig, ax =column_chosen, ax=ax, legend=True,
stats_results.plot(column={'label': f"{column_chosen} per country.",
legend_kwds'orientation': "horizontal"})
'Choropleth Map Showing GRDI Mean per country')
ax.set_title(# Turn off the axis numbers and ticks
ax.set_axis_off() plt.show()
Selecting Data by Column
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.
# selected_countries = ["Algeria", "Somalia", "Colombia", "Timor Leste", "Finland", "Nicaragua", "United Kingdom", "Mali"]
# selected_countries = ["Anguilla", "Armenia", "Angola", "Argentina", "Albania", "United Arab Emirates", "American Samoa", "Australia" ]
= ["Algeria", "Somalia", "Colombia", "Timor Leste", "Finland", "Nicaragua", "United Kingdom", "Mali", "Armenia", "Argentina", "Albania", "United Arab Emirates", "Indonesia", "Qatar"]
selected_countries
#use the list above to subset the country_table DataFrame by the column NAME_0
= country_table[country_table['NAME_0'].isin(selected_countries)] selected_countries
Using a Defined Custom Function
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.
= calculate_country_stats(selected_countries, grdi_raster) stats_results
Show the head of the table again:
stats_results.head()
Country_Name | Country_GID | GRDI_Count | GRDI_Min | GRDI_Mean | GRDI_Max | GRDI_Median | GRDI_Std | GRDI_Range | GRDI_P20 | GRDI_P40 | GRDI_P60 | GRDI_P80 | geometry | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | Albania | ALB | 19076 | 8.272310 | 61.513866 | 75.395561 | 66.220139 | 11.644754 | 67.123251 | 55.861691 | 64.445587 | 67.104332 | 68.589424 | MULTIPOLYGON (((20.0541 39.6917, 20.0389 39.69... |
1 | United Arab Emirates | ARE | 18229 | 5.732072 | 42.347647 | 67.470955 | 45.034630 | 17.949307 | 61.738883 | 24.688160 | 38.623692 | 50.710217 | 61.909931 | MULTIPOLYGON (((54.1541 22.7548, 53.3313 22.85... |
2 | Argentina | ARG | 474297 | 7.572341 | 66.341158 | 81.701645 | 68.789696 | 10.416561 | 74.129304 | 66.774643 | 68.278297 | 69.335510 | 71.098465 | MULTIPOLYGON (((-66.5458 -55.061, -66.5486 -55... |
3 | Armenia | ARM | 9108 | 6.944193 | 59.871617 | 73.707901 | 64.272858 | 12.315193 | 66.763708 | 52.692970 | 61.288189 | 66.743027 | 69.004845 | MULTIPOLYGON (((45.8319 39.8311, 45.8448 39.82... |
4 | Colombia | COL | 223557 | 11.956444 | 71.523674 | 84.922409 | 73.350586 | 8.762370 | 72.965965 | 69.966934 | 72.190056 | 74.405403 | 76.502792 | MULTIPOLYGON (((-77.491 4.1451, -77.4985 4.140... |
Plot the map again choosing a column to plot:
= 'GRDI_Max' #GRDI_Max, GRDI_Min, GRDI_Median
column_chosen =column_chosen, legend=True)
stats_results.plot(column plt.show()
Creating a Table with Results
We can create a list of tuples that we can use to refer to the GRDI 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.
# List of GRDI values and their corresponding properties
#column, value name, color, symbol
= [
grdi_data 'GRDI_Mean', 'Mean', 'orange', 'diamond'),
('GRDI_Min', 'Min', 'gray', '152'),
('GRDI_Max', 'Max', 'gray', '151'),
('GRDI_P20', 'Q20', 'blue', '142'),
('GRDI_P40', 'Q40', 'purple', '142'),
('GRDI_P60', 'Q60', 'green', '142'),
('GRDI_P80', 'Q80', 'red', '142')
( ]
We can create a figure to display the data based on the names colors and symbols we selected.
# Create a figure
= go.Figure()
fig
# Add traces to the figure based on the data
for col, name, color, symbol in grdi_data:
fig.add_trace(go.Scatter(=stats_results[col],
x=stats_results['Country_Name'],
y='markers',
mode=name,
name=dict(color=color, size=10, symbol=symbol)
marker
))
# Customize layout
fig.update_layout(='GRDI Statistics by Country',
title='GRDI Values',
xaxis_title='Country Name',
yaxis_title=dict(tickmode='linear'),
yaxis='Statistics',
legend_title='category',
yaxis_type=dict(tickvals=[0, 20, 40, 60, 80, 100])
xaxis
)
# Show plot
fig.show()
Summarizing PM2.5 Values by Socioeconomic Deprivation
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.
The results will be added to the stats_results
with the corresponting columns.
# iterate through the stats_results table rows
for index, row in stats_results.iloc[:].iterrows():
#isolate each country's respective row
= gpd.GeoDataFrame([row], geometry='geometry').reset_index(drop=True)
row_df print(row_df.loc[0,'Country_GID'])
try:
#use rioxarray to clip the GRDI and PM2.5 rasters by the geometry of the respective country.
= grdi_raster.rio.clip(row_df.geometry, grdi_raster.rio.crs)
grdi_country = pm25_raster.rio.clip(row_df.geometry, grdi_raster.rio.crs)
pm25_country except:
print('Error in clip')
continue
#Applying squeeze() to this array removes the singleton dimension, reducing it to a 2D array with dimensions (rows, columns)
= grdi_country.squeeze()
grdi_country= pm25_country.squeeze()
pm25_country
# Subset the GRDI raster where values fall between each GRDI quintiles
= grdi_country.where((grdi_country >= row_df.loc[0, 'GRDI_Min']) & (grdi_country <= row_df.loc[0, 'GRDI_P20']))
grdi_countryQ1 = grdi_country.where((grdi_country >= row_df.loc[0, 'GRDI_P20']) & (grdi_country <= row_df.loc[0, 'GRDI_P40']))
grdi_countryQ2 = grdi_country.where((grdi_country >= row_df.loc[0, 'GRDI_P40']) & (grdi_country <= row_df.loc[0, 'GRDI_P60']))
grdi_countryQ3 = grdi_country.where((grdi_country >= row_df.loc[0, 'GRDI_P60']) & (grdi_country <= row_df.loc[0, 'GRDI_P80']))
grdi_countryQ4 = grdi_country.where((grdi_country >= row_df.loc[0, 'GRDI_P80']) & (grdi_country <= row_df.loc[0, 'GRDI_Max']))
grdi_countryQ5
# Mask the PM2.5 raster using the above GRDI quartile rasters, keeping only the cells that intersect
= pm25_country.where(grdi_countryQ1.notnull())
pm25_countryQ1 = pm25_country.where(grdi_countryQ2.notnull())
pm25_countryQ2 = pm25_country.where(grdi_countryQ3.notnull())
pm25_countryQ3 = pm25_country.where(grdi_countryQ4.notnull())
pm25_countryQ4 = pm25_country.where(grdi_countryQ5.notnull())
pm25_countryQ5
#Find the mean value of of the intersected PM2.5 rasters in each quartile
= pm25_countryQ1.mean().item()
pm25_countryQ1v = pm25_countryQ2.mean().item()
pm25_countryQ2v = pm25_countryQ3.mean().item()
pm25_countryQ3v = pm25_countryQ4.mean().item()
pm25_countryQ4v = pm25_countryQ5.mean().item()
pm25_countryQ5v
#add the resuts to the stats_results table in the respective column
'PM25_Q1'] = pm25_countryQ1v
stats_results.at[index, 'PM25_Q2'] = pm25_countryQ2v
stats_results.at[index, 'PM25_Q3'] = pm25_countryQ3v
stats_results.at[index, 'PM25_Q4'] = pm25_countryQ4v
stats_results.at[index, 'PM25_Q5'] = pm25_countryQ5v stats_results.at[index,
ALB
ARE
ARG
ARM
COL
DZA
FIN
GBR
IDN
MLI
NIC
QAT
SOM
stats_results.head()
Country_Name | Country_GID | GRDI_Count | GRDI_Min | GRDI_Mean | GRDI_Max | GRDI_Median | GRDI_Std | GRDI_Range | GRDI_P20 | GRDI_P40 | GRDI_P60 | GRDI_P80 | geometry | PM25_Q1 | PM25_Q2 | PM25_Q3 | PM25_Q4 | PM25_Q5 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | Albania | ALB | 19076 | 8.272310 | 61.513866 | 75.395561 | 66.220139 | 11.644754 | 67.123251 | 55.861691 | 64.445587 | 67.104332 | 68.589424 | MULTIPOLYGON (((20.0541 39.6917, 20.0389 39.69... | 15.438293 | 15.031855 | 14.699892 | 14.612475 | 15.709009 |
1 | United Arab Emirates | ARE | 18229 | 5.732072 | 42.347647 | 67.470955 | 45.034630 | 17.949307 | 61.738883 | 24.688160 | 38.623692 | 50.710217 | 61.909931 | MULTIPOLYGON (((54.1541 22.7548, 53.3313 22.85... | 47.710175 | 47.893559 | 47.822292 | 48.220722 | 49.729328 |
2 | Argentina | ARG | 474297 | 7.572341 | 66.341158 | 81.701645 | 68.789696 | 10.416561 | 74.129304 | 66.774643 | 68.278297 | 69.335510 | 71.098465 | MULTIPOLYGON (((-66.5458 -55.061, -66.5486 -55... | 7.519970 | 5.924110 | 6.083218 | 7.562082 | 8.672854 |
3 | Armenia | ARM | 9108 | 6.944193 | 59.871617 | 73.707901 | 64.272858 | 12.315193 | 66.763708 | 52.692970 | 61.288189 | 66.743027 | 69.004845 | MULTIPOLYGON (((45.8319 39.8311, 45.8448 39.82... | 19.292490 | 16.728804 | 16.334721 | 16.337336 | 15.342737 |
4 | Colombia | COL | 223557 | 11.956444 | 71.523674 | 84.922409 | 73.350586 | 8.762370 | 72.965965 | 69.966934 | 72.190056 | 74.405403 | 76.502792 | MULTIPOLYGON (((-77.491 4.1451, -77.4985 4.140... | 27.425064 | 29.047859 | 24.132767 | 22.309097 | 20.367249 |
Plot Results of Mean PM2.5 in Socieceonomic Deprivation Quartiles per country
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.
# List of GRDI values and their corresponding properties
#column, value name, color, symbol
=[
plot_data 'PM25_Q1', 'Q1', '#440154', '6'), # Light Blue
('PM25_Q2', 'Q2', '#31688E', '5'), # Light Green
('PM25_Q3', 'Q3', '#35B779', '7'), # Yellow
('PM25_Q4', 'Q4', '#FDE725', '8'), # Orange
('PM25_Q5', 'Q5', '#FF0000', '1') # Red
(
]
# Create a figure
= go.Figure()
fig
# Add traces to the figure.
for col, name, color, symbol in plot_data:
= np.log(stats_results[col])
xlog
fig.add_trace(go.Scatter(=xlog,
x=stats_results['Country_Name'],
y='markers+text', # Add 'text' to mode
mode=[f'<b>{name}</b>' for _ in stats_results[col]], # Repeat name for each point
text=name,
name=dict(color=color, size=12),
textfont='top center', # Position the text above the symbol
textposition=color,
marker_color="midnightblue",
marker_line_color=symbol,
marker_symbol=14,
marker_size=2,
marker_line_width=0.6
marker_opacity
))='top center')
fig.update_traces(textposition
# Customize layout
fig.update_layout(='Mean PM2.5 in each GRDI Quartile by Country',
title='Log of PM2.5 Mean Values',
xaxis_title='Country Name',
yaxis_title=dict(tickmode='linear'),
yaxis='Statistics',
legend_title='category',
yaxis_type=dict(rangemode="tozero"),
xaxis
#xaxis=dict(tickvals=[0, 20, 40, 60, 80, 100])
)
# Show plot
fig.show()
Use the plotly controls to take a closer look at the results.
With this shapely plot, We can examine differences between countires 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.
Congratulations! …. Now you should be able to:
- Test test…
Module 2: Air Quality Home
In this lesson, we explored ….