1. Introduction
Smoggy weather can be frequently observed in many parts of the world nowadays. Fine solid particles suspended in the air play a major role in causing smoggy conditions by absorbing or scattering sunlight. These particles are called aerosols. Sources of aerosols include pollution from factories, smoke from fires, dust from storms, volcanic ash etc (Source: NASA). Aerosols can be detrimental to human health. They can also interfere with temperature, cloud formation and atmospheric circulation patterns. It is, therefore, important to monitor its level and keep it within limits.
Moderate Resolution Imaging Spectroradiometer (MODIS) sensors aboard NASA’s Terra satellite track the Aerosol Optical Depth (AOD) at daily frequency. It is a dimensionless number that is related to the amount of aerosol in the vertical column of atmosphere over the observation location. AOD tells us how much direct sunlight is prevented from reaching the ground by these aerosol particles. A value of 0.01 corresponds to an extremely clean atmosphere, and a value of 0.4 corresponds to hazy air condition. AOD value above 1 represents high concentration of aerosols in the atmosphere during severe pollution. Some natural phenomena like wildfires and sand storms result in high AOD values between 1 and 3. The maximum AOD which is reported by MODIS sensors is 4.
In this blog piece, we will extract daily AOD data (2000–2020) using Google Earth Engine (GEE) and find out the severity of aerosol induced pollution in the city of Delhi.
2. Workflow
The workflow can be divided into three broad sections
- Collecting AOD data using GEE.
- Importing the data as pandas dataframe and plotting.
- Counting the number of days with severe pollution in each year.
3. Collecting AOD data using GEE
We will begin by importing the relevant libraries and authenticating the GEE.
# Import relevant libraries
import ee
# Authenticate GEE and initialize
try:
ee.Initialize()
except Exception as e:
ee.Authenticate()
ee.Initialize()
Next, we import the vector shapefile for the region of interest (RoI) and the raster images for AOD. Daily AOD values are obtained from different images of the collection which are tagged by the date on which the location was observed by the satellite. We define a function reduce_del
to reduce the pixel values contained within the RoI to a single statistic by averaging them. This function is then mapped over all the images in the collection to derive the average AOD for the region for the selected duration. The resulting time series is then exported to Google Drive as a CSV file.
# Import the boundaries of the region of interest
delhi = ee.FeatureCollection("FAO/GAUL/2015/level1")
.filter(ee.Filter.eq('ADM1_NAME','Delhi'))
# Import the raster data on aerosol concentration and apply regional filters.
# It vastly improves run time of the code by reducing the number of images in the collection.
asol = ee.ImageCollection("MODIS/006/MCD19A2_GRANULES")
.filterBounds(delhi)
.filterDate('2000-01-01', '2020-12-31');
# Reduce the image data into a statistic for Delhi
def reduce_del(img):
stat = img.reduceRegion(reducer=ee.Reducer.mean(), geometry=delhi, scale=1000).get('Optical_Depth_047')
date = img.date().format()
mydict = {'date': date, 'stat': stat}
return ee.Feature(None).set(mydict)
# Convert ImageCollection into FeatureCollection and filter the FeatureCollection keeping only notNull values
reduced_del = ee.FeatureCollection(asol.map(reduce_del)).filter(ee.Filter.notNull(['date','stat']))
# Code for exporting Delhi AOD data to GDrive as a batch process.
task = ee.batch.Export.table.toDrive(
collection = reduced_del,
description = 'aod_stat_del',
folder = 'satellite',
fileFormat = 'CSV',
selectors = ['date','stat']
)
task.start()
4. Importing the data as time series and plotting
We download the data from Google Drive to the local machine and import it using Pandas. The AOD values are reported on the scale of 0.001. We rescale the values to a scale of 1 by dividing all the values by 1000. As there are multiple observations for the same day, we take their average to arrive at single value for each day.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import matplotlib.pyplot as plt
aerosoldf = pd.read_csv('C:UsersADMINDesktopsatellitegeemapdataaod_stat_del.csv')
# Define a function to transform the dataframe into time series
# Rescale AOD values from 0.001 to 1 by dividing the numbers by 1000.
# Group the observations by day and take the average AOD value
def tseries(df):
df['date_formatted'] = pd.to_datetime(df['date']).dt.normalize()
df.stat = df.stat/1000
df = df.groupby('date_formatted').agg(mean_stat = pd.NamedAgg(column = 'stat', aggfunc = 'mean'))
return df
aerosoldf = tseries(aerosoldf)
# we create a figure with pyplot and set the dimensions to 15 x 7
fig, ax = plt.subplots(figsize=(15,7))
# we'll create the plot by setting our dataframe to the data argument
ax.plot(aerosoldf.index, aerosoldf.mean_stat)
# we'll set the labels and title
ax.set_ylabel('aerosol optical depth',fontsize=20)
ax.set_xlabel('date',fontsize=20)
ax.set_title('Daily Aerosol Optical Depth Delhi (2000 - 2020)',fontsize=20);
plt.show()

5. Data processing
The daily time series chart looks very messy. We have plotted a large span of data (daily data for 20 years) with a purpose of validating the claims of Climate Change. It may be noticed that the number of occasions when AOD exceeds the value of 1 have increased in the recent periods implying an increased frequency of severe pollution incidents. Let us count the number of days in each year when AOD was greater than 1.
# Retrieve year information from date index
aerosoldf['year'] = aerosoldf.index.year
#Define AOD threshold as 1
aerosoldf['threshold'] = 1
#Create a binary column which takes value 1 when AOD exceeds threshold and 0 otherwise
aerosoldf['severe'] = np.where(aerosoldf['mean_stat'] > aerosoldf['threshold'], 1,0)
#Group the observations and count the number of days having severe pollution in each year
aerosol = aerosoldf.groupby('year').agg(obs_days = pd.NamedAgg(column = 'severe', aggfunc = 'count'),
severe_days = pd.NamedAgg(column = 'severe', aggfunc = 'sum'))
# As the number of observed days are not uniform across years, we calculate the percentage of severe days and multiply by 365 to make it comparable across years
aerosol['severe_days_annual'] = aerosol['severe_days'] / aerosol['obs_days'] * 365
aerosol

# Function to find slope and intercept of trendline
def trendline(df):
# get coeffs of linear fit
slope, intercept, r_value, p_value, std_err = stats.linregress(df.index,df['severe_days_annual'])
pct = intercept + slope * df.index
return(pct)
#Plot
fig = plt.figure(figsize=(15,6))
ax = fig.add_subplot(111)
ax.plot(aerosol.index.astype(int), aerosol.severe_days_annual, linestyle="dashed", marker = 'o')
ax.plot(aerosol.index.astype(int), trendline(aerosol), color = 'red', marker='o')
ax.set_xticks(range(2000,2021))
# we'll set the labels and title
ax.set_title('Number of days with severe aerosol concentration in Delhi',fontsize=20);
plt.grid()
plt.show()

6. Conclusion
The trendline graph shows that Delhi air is severely polluted with high levels of aerosols for around 4 months in a year. Number of severe days in a year has increased from 70 days in year 2000 to 110 days in year 2020. Incidents of severe pollution have become almost 60% more frequent than they used to be two decades ago. That’s alarming!! This underscores the need for meaningful dialogues and plan of action at national and international levels as effects of pollution are transnational.
Disclaimer: Views expressed in this blog are personal.
References
- https://earthdata.nasa.gov/earth-observation-data/near-real-time/hazards-and-disasters/air-quality
- MCD19A2.006: Terra & Aqua MAIAC Land Aerosol Optical Depth Daily 1km distributed by NASA’s Land Processes Distributed Active Archive Center (LP DAAC)
- Wu, Q., (2020). geemap: A Python package for interactive mapping with Google Earth Engine. The Journal of Open Source Software, 5(51), 2305. https://doi.org/10.21105/joss.02305
- Gorelick, N., Hancher, M., Dixon, M., Ilyushchenko, S., Thau, D., & Moore, R. (2017). Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sensing of Environment.