Display AQI for a given LAT/LNG using data from https://fire.airnow.gov (AirNow, AirSis, Purple Air)

airnow.py

#!/usr/bin/python3
# Prints current PM2.5 concentrations based on fire.airnow.gov map
# pip3 install requests pytz geobuf
# usage: LAT=37.9053745 LNG=-122.3048239 MAX_FEATURES=6 MAX_DIST=11 python3 airnow.py
# optional, for AQI instead of raw concentration (µg/m³): pip3 install python-aqi  and pass AQI=1
import collections
import datetime
import math
import os
import requests
import pytz

SOURCES = [
    'https://s3-us-west-2.amazonaws.com/airfire-data-exports/maps/geobuf/purple_air_epa_qc.pbf',
    'https://s3-us-west-2.amazonaws.com/airfire-data-exports/maps/geobuf/airnow_PM2.5_latest10.pbf',
    'https://s3-us-west-2.amazonaws.com/airfire-data-exports/maps/geobuf/airsis_PM2.5_latest10.pbf',
    #'https://s3-us-west-2.amazonaws.com/airfire-data-exports/dev/test_pa/purple_air_epa_qc.geojson',
    #'https://s3-us-west-2.amazonaws.com/airfire-data-exports/monitoring/v1/geojson/airnow_PM2.5_latest10.geojson',
]
LAT = float(os.getenv('LAT'))
LNG = float(os.getenv('LNG'))
MAX_FEATURES = int(os.getenv('MAX_FEATURES', 5))
MAX_DIST = float(os.getenv('MAX_DIST') or 100)
TZ = os.getenv('TZ') or 'America/Los_Angeles'

def distance(origin, destination):
    """
    Calculate the Haversine distance.

    Parameters
    ----------
    origin : tuple of float
        (lat, long)
    destination : tuple of float
        (lat, long)

    Returns
    -------
    distance_in_km : float

    Examples
    --------
    >>> origin = (48.1372, 11.5756)  # Munich
    >>> destination = (52.5186, 13.4083)  # Berlin
    >>> round(distance(origin, destination), 1)
    504.2
    """
    lat1, lon1 = origin
    lat2, lon2 = destination
    radius = 6371  # km

    dlat = math.radians(lat2 - lat1)
    dlon = math.radians(lon2 - lon1)
    a = (math.sin(dlat / 2) * math.sin(dlat / 2) +
         math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) *
         math.sin(dlon / 2) * math.sin(dlon / 2))
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
    d = radius * c

    return d

def colprint(lines):
    widths = {}
    for l in lines:
        for i, c in enumerate(l):
            w = len(c)
            ow = widths.get(i, 0)
            w = max(ow, w)
            widths[i] = w

    for l in lines:
        for i, c in enumerate(l):
            ow = widths[i]
            w = len(c)
            print(c, ' '*(ow-w), end=" ")
        print()

def conc_or_aqi(conc):
    if os.getenv('AQI'):
        try:
            import aqi
            aqi_val = float(aqi.to_iaqi(aqi.POLLUTANT_PM25, conc, algo=aqi.ALGO_EPA))
            return '{:.0f} AQI'.format(aqi_val)
        except ImportError:
            pass
    return "{0:5.1f} µg/m³".format(conc)

def main():
    lines = []
    features = []
    for url in SOURCES:
        resp = requests.get(url)
        if url.endswith('.pbf'):
            import geobuf
            js = geobuf.decode(resp.content)
        else:
            js = resp.json()

        for f in js['features']:
            lng, lat = f['geometry']['coordinates']
            dist = distance((LAT, LNG), (lat, lng))
            if dist <= MAX_DIST:
                features.append((dist, f))

    # closest first
    features.sort(key=lambda f: f[0])

    # which sources did we receive data for?
    features_by_source = collections.defaultdict(list)
    for dist, f in features:
        features_by_source[f['properties']['dataSource']].append((dist, f))
    unique_sources = len(features_by_source)

    # truncate number of features to show
    show_features = features[:MAX_FEATURES]

    # but show at least one feature from every available source
    show_sources = set(f['properties']['dataSource'] for _, f in show_features)
    for source in features_by_source:
        if source not in show_sources:
            show_features = show_features[:-1] + [features_by_source[source][0]]
            show_sources.add(source)

    for d, f in show_features:
        p = f['properties']
        t = datetime.datetime.strptime(p['lastValidUTCTime'], "%Y-%m-%d %H:%M:%S")
        t = pytz.utc.localize(t).astimezone(pytz.timezone(TZ))
        t = t.strftime("%Y-%m-%d %H:%M:%S")
        lines.append([
            "({d:5.2f} km)".format(d=d),
            t,
            p['dataSource'],
            p['siteName'],
            conc_or_aqi(p['PM2.5_1hr']),
        ])

    colprint(lines)

if __name__ == "__main__":
    main()