Why Combine Satellite and Ground Data?

AirQo’s low-cost air quality monitors offer high-frequency, on-the-ground measurements of PM2.5 and PM10. However, they provide limited spatial coverage. On the flip side, satellite AOD data offers rich spatial resolution, but it's not always available daily or at every location due to cloud cover and atmospheric conditions.

By bringing both worlds together:

We fill spatial gaps in monitoring.
Gain contextual insights into local pollution trends.
Improve modeling and forecasting for policymakers, scientists, and communities.

The Data Processing Pipeline

Here's a look under the hood at how this integration works.

1. AirQo Data Preparation

Using daily data from AirQo sensors (including calibrated PM2.5/PM10 values, site coordinates, and timestamps), we:

Convert timestamps into daily aggregates.
Organize data by site, date, and location.

2. MODIS AOD Retrieval via Google Earth Engine

We query MODIS’s MCD19A2 GRANULES for Optical_Depth_047 around each site’s coordinates and timestamp.

To handle MODIS data sparsity:

We implement a dynamic radius search starting at 1km and expanding up to 10km until a valid AOD value is found.

3. Parallel Processing for Speed

To scale efficiently across many records, we use Python’s ThreadPoolExecutor to parallelize AOD retrieval, dramatically speeding up the pipeline.

4. Saving and Filtering Results

Only records with valid AOD data are kept, and the final enriched dataset is saved as a CSV, complete with:

Mean PM2.5/PM10
AOD values
Search radius used
Site metadata

import ee
import geemap
import pandas as pd
from tqdm.notebook import tqdm
import numpy as np
import warnings
import os
import logging
from concurrent.futures import ThreadPoolExecutor

# Suppress warnings
warnings.filterwarnings('ignore')

# Configure logging
logging.basicConfig(level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s')

# Function to authenticate and initialize Earth Engine
def initialize_earth_engine():
    try:
        ee.Authenticate()
        ee.Initialize(project='satelliteAOD')
        logging.info("Earth Engine initialized successfully.")
    except Exception as e:
        logging.error(f"Failed to initialize Earth Engine: {e}")
        raise

# Function to load air quality data
def load_data(file_path):
    try:
        data = pd.read_csv(file_path, parse_dates=['timestamp'])
        logging.info(f"Data loaded successfully from {file_path}.")
        return data
    except Exception as e:
        logging.error(f"Failed to load data from {file_path}: {e}")
        raise

# Function to calculate mean AOD with dynamic radius adjustment
def calculate_mean_aod(row, initial_radius=1000, max_radius=10000):
    lat = row['site_latitude']
    lon = row['site_longitude']
    date = row['date']
    point = ee.Geometry.Point(lon, lat)
    date_str = date.strftime('%Y-%m-%d')
    start_date = ee.Date(date_str)
    end_date = start_date.advance(1, 'day')

    # Select MODIS AOD image collection
    modis_aod = ee.ImageCollection('MODIS/061/MCD19A2_GRANULES') \
                .filterDate(start_date, end_date) \
                .filterBounds(point) \
                .select('Optical_Depth_047')

    if modis_aod.size().getInfo() > 0:
        radius = initial_radius
        while radius <= max_radius:
            try:
                aod = modis_aod.mean().rename('AOD')
                mean_aod = aod.reduceRegion(reducer=ee.Reducer.mean(), geometry=point.buffer(radius), scale=1000).get('AOD')
                aod_value = mean_aod.getInfo()
                if aod_value is not None:
                    logging.info(f"AOD found for lat: {lat}, lon: {lon}, date: {date_str} with radius {radius}m.")
                    return aod_value, radius  # Return both AOD value and search radius
                else:
                    radius *= 2  # Double the radius and try again
            except Exception as e:
                logging.error(f"Error calculating AOD for lat: {lat}, lon: {lon}, date: {date_str} with radius {radius}m: {e}")
                radius *= 2
#        logging.warning(f"No valid AOD found within {max_radius}m radius for lat: {lat}, lon: {lon}, date: {date_str}.")
        return None, max_radius  # Return None for AOD and the maximum radius searched
    else:
#        logging.warning(f"No MODIS image found for lat: {lat}, lon: {lon}, date: {date_str}.")
        return None, 0  # Return None for AOD and 0 for radius (no search performed)

# Main function to process the data
def process_data(aq_data):
    # Ensure the 'timestamp' column is in datetime format
    aq_data['timestamp'] = pd.to_datetime(aq_data['timestamp'])

    # Group by date and site to get daily mean values
    aq_data_daily = aq_data.groupby(['site_id', aq_data['timestamp'].dt.date]).agg({
        'site_name': 'first',
        'site_latitude': 'first',
        'site_longitude': 'first',
        'pm2_5_calibrated_value': 'mean',
        'pm2_5_raw_value': 'mean',
        'pm10_calibrated_value': 'mean',
        'city': 'first',
        'country': 'first'
    }).reset_index().rename(columns={'timestamp': 'date'})

    # Ensure the 'date' column is in datetime format
    aq_data_daily['date'] = pd.to_datetime(aq_data_daily['date'])

    # Use ThreadPoolExecutor to parallelize the AOD calculations
    with ThreadPoolExecutor(max_workers=8) as executor:
        results = list(tqdm(executor.map(calculate_mean_aod, [row for _, row in aq_data_daily.iterrows()]),
                            total=aq_data_daily.shape[0]))

    # Unpack results into AOD values and search radii
    aod_values, search_radii = zip(*results)

    # Add AOD values and search radii to the DataFrame
    aq_data_daily['Mean_AOD'] = aod_values
    aq_data_daily['AOD_Search_Radius'] = search_radii

    # Extract day, month, and year from date
    aq_data_daily['Day'] = aq_data_daily['date'].dt.day
    aq_data_daily['Month'] = aq_data_daily['date'].dt.month

    return aq_data_daily

# Function to create the output directory if it doesn't exist
def ensure_directory_exists(file_path):
    directory = os.path.dirname(file_path)
    if not os.path.exists(directory):
        os.makedirs(directory)
        logging.info(f"Created directory: {directory}")

# Function to save the processed data
def save_data(data, output_path):
    # Ensure the output directory exists
    ensure_directory_exists(output_path)

    # Select only the relevant columns for the output
    data = data[[ 'site_id','date', 'site_name','site_latitude', 'site_longitude',  'pm2_5_calibrated_value', 'pm2_5_raw_value',
                 'pm10_calibrated_value','city', 'country','Month', 'Day','Mean_AOD', 'AOD_Search_Radius'  ]]

    data = data.dropna(subset=['Mean_AOD'])
    data.to_csv(output_path, index=False)
    logging.info(f"Data saved to {output_path}")

# Main script
if __name__ == "__main__":
    try:
        # Initialize Earth Engine
        initialize_earth_engine()

        # Mount Google Drive
        from google.colab import drive
        drive.mount('/content/drive', force_remount=True)

        # Define cities and years
        cities = ['Fort_Portal']
        years = [2020]

        for city in cities:
            for year in years:
                try:
                    # Define file paths
                    data_path = f'/content/drive/MyDrive/LUD/data/{year}/{city}_{year}_sites_dataset.csv'

                    # Load air quality data
                    aq_data = load_data(data_path)

                    # Process data to calculate mean AOD and extract day, month, and year
                    processed_data = process_data(aq_data)

                    # Define output path
                    output_path = f'/content/drive/MyDrive/LUD/MEAN_NEAREST_AOD/{city}/{city}_air_quality_with_mean_aod_{year}.csv'

                    # Save the processed data
                    save_data(processed_data, output_path)

                    logging.info(f"Mean AOD calculation complete for {city} in {year}.")

                except Exception as e:
#                    logging.error(f"An error occurred while processing data for {city} in {year}: {e}")
                    print(f"Error processing data for {city} in {year}: {e}")

    except Exception as e:
        logging.error(f"An error occurred in the main script: {e}")
        print(f"An error occurred in the main script:: {e}")

Bridging the Skies and Streets: Integrating Stationary AirQo Sensors with Satellite AOD Data for Cleaner Air Insights

Overview

Details

Bridging the Skies and Streets: Integrating Stationary AirQo Sensors with Satellite AOD Data for Cleaner Air Insights

Overview

Details