EPA air quality: NO2 in Phoenix New York San Diego 2016 2017-07-29 Jul

0-EPA-air-quality.md

EPA air quality: NO2 and ozone in 3 US cities in 2016

The Environmental Protection Agency collects data on NO2, ozone, wind etc. in hundreds of US cities, and puts it in files like daily_TEMP_2015.csv on the web site http://aqsdr1.epa.gov/aqsweb/aqstmp/airdata .

Plots

A leading question: how can one set reasonable limits for pollutants that vary so much ?

The files below:

epa_download  -- download daily_TEMP_2016.zip ...
epa_air_smallcsv.py  -- big EPA files -> small csv
2016-epa-ndays-yearavs.tsv  -- city -> nr. days and yearly averages
pr.py  -- for debugging
zutil.py

First look at epa_download -help. This downloads .zip files from the url to the current directory with curl; it also links daily_44602* to daily_no2* and daily_44201* to daily_ozone* .

epa_air_smallcsv.py reads the big .zip files and merges them, e.g.
in: daily_no2_2016.zip daily_ozone_2016.zip daily_TEMP_2016.zip daily_WIND_2016.zip
out: 2016-epa-no2-ozone-temp-wind.csv, ~ 80k rows and 7 columns:

city        state       date        no2     ozone   temp    wind
Albuquerque New-Mexico  2016-01-01  10.9    0.0224  30.3    2.59
Albuquerque New-Mexico  2016-01-02  15.4    0.0208  31.3    2.62
...
Yukon       Oklahoma    2016-12-31  5       NaN     43.4    7.3

Here "NaN", not-a-number, marks missing data for a given city and day. I drop cities with < 350 days, and rows with more than one NaN, i.e. 2 or more of NO2 ozone TEMP WIND missing.

(There's more than one way to skin a cat, split up data; this is adequate for my purpose, which is to just plot the data to see its day-to-day variability. Why so variable ? No idea; comments welcome.)

Missing days

It looks as though many days are missing in the EPA data for no2, ozone, TEMP and WIND, in many cities. It's entirely possible that, as a novice with pandas, I've made mistakes. But simple checks like

egrep -w Chicago daily_TEMP_2016.csv | less

show data for "Chicago-Naperville-Elgin", but not "Chicago".

The file 2016-epa-ndays-yearavs.tsv below gives the numbers of days that have data, and the yearly averages, for the cities in the output .csv:

n_no2 n_ozone n_temp n_wind av_no2 av_ozone av_temp av_wind city state
366 366 366 366  19.6 0.0306 75.5 2.8  Phoenix Arizona 
...
362 362   0 362  15.2 0.0278 NaN 4.77  Chicago Illinois
...

Notes

Interpolate ? One could fill in missing data from nearby days and nearby cities, e.g. with scipy KDTree and inverse-distance weighting . But be warned: interpolation can affect later analysis. If you do interpolate, measure how good it is: throw out some known cities or days, interpolate there, and compare the interpolated to the real values.

Fwiw, Weibull plots of NO2 and ozone don't look that linear; nor, surprisingly, does wind. Would a mixture model, normal + heavy tail, make sense physically ?

Comments are welcome, test cases most welcome.

cheers
-- denis

Last change: 29 July 2017

0-epa_air_smallcsv.py

#!/usr/bin/env python2
""" python this.py [year=2015 to=2016]
    read_epa_download no2 ozone TEMP WIND .zip or .csv
        in download_dirs from epa_download
    write e.g. tmp/2016-epa-no2-ozone-temp-wind.csv  79k rows, 7 columns
        city  state  date  no2  ozone  temp  wind
"""

from __future__ import division
import numpy as np
import pandas as pd

from pr import pr
import zutil as nu

__version__ = "2017-07-29 july  denis-bz-py t-online de"

    # read these columns only, and rename them --
epa_columns_dict = {
"City Name":    "city",
"State Name":   "state",
"Date Local":   "date",
"Arithmetic Mean":  "av",  # hour / 8 hour
# "1st Max Value":    "max1",
"Parameter Name":   "param",  # drop
}

epa_url = "http://aqsdr1.epa.gov/aqsweb/aqstmp/airdata"
    # (download with requests ? mttiw)
download_dirs = "epa/download $epa"  # <-- where you run epa_download

#...............................................................................
def read_epa_download( csv_or_zip, columns=epa_columns_dict, nrows=1e8, dropneg=False ):
    """ read e.g. daily_TEMP_2015.csv or .zip
        the given columns only, and rename them
        WIND: speed only, drop direction
        city and state names with "-", New-York
    """
        # units: NO2 ppb, ozone ppm, temp F, wind knots = .51444 m/s
        # NO2 1 ppb = 1.91 ug/m3, Ozone 1 ppb = 2.0 ug/m3
        # read_csv e.g. daily_no2_2015.zip with a single csv works

    df = pd.read_csv( csv_or_zip, usecols=columns.keys(),
            parse_dates=["Date Local"], infer_datetime_format=True,
            index_col=False, nrows=nrows )
    df.rename( columns=columns, inplace=True )
    df = df[ df.city != "Not in a city" ]
    if "param" in df.columns:
        df = df[ df.param != "Wind Direction - Resultant" ]  # wind speed only
        df = df.drop( "param", axis=1 )
    if dropneg and "av" in df.columns:
        # drop no2 < 0 in Grapevine Midlothian Tyler TX ?
        df = df.drop( df.av < 0 )
    df.city = df.city.str.replace( r" \(.*", "" )  # Dentsville (Dents) -> Dentsville
    df.city = df.city.str.replace( ".", "" )  # St. Louis
    df.city = df.city.str.replace( " ", "-" )  # New-York, for awk
    df.state = df.state.str.replace( " ", "-" )
    df = df.sort_values( by=["city", "date"] )
    df.index = np.arange( len(df) )
    nr, nc = df.shape
    print "\nread_epa_csv %s: %d rows, %d columns: %s" % (
            csv_or_zip, nr, nc, " ".join( df.columns ))
    if "av" in df.columns:
        q = np.r_[0, 10, 25, 50, 75, 90, 100]
        print "Daily averages, percentiles %s: %s" % (q, np.percentile( df.av, q=q ))

    return df


def epa_zipfilename( what, year, suffix=".zip", dirs=download_dirs ):
    basename = "daily_%s_%s%s" % (what, year, suffix)
    return nu.findfile( basename, dirs=dirs )  # else IOError


def write_csv( out, df, csvheader="", **kw ):
    """ pandas df.to_csv with csvheader """
    nr, nc = df.shape
    print "\nsaving to %s  %d rows, %d columns" % (out, nr, nc)
    nu.mkdirpart( out )
    with open( out, "w" ) as f:
        if csvheader:
            f.write( csvheader )  # grr not gist.github
        df.to_csv( f, sep="\t", float_format="%.3g", na_rep="NaN", index=None, **kw )

#...............................................................................
def merge_epa_csvs( year=2015, what="no2 ozone TEMP WIND", mindays=350, maxnan=1,
                nrows=1e8, verbose=1 ):
    """ for no2 ozone ...:
            read_epa_download( findfile )
        -> merged df with columns: city  state  date  no2  ozone  temp  wind
    """
    print "\n-- merge_epa_csvs  mindays %d  maxnan %d --" % (mindays, maxnan)
    whats = what.split()
    avdf_list = []

    for what in whats:
        print "\n", 80 * "-"
        csv = epa_zipfilename( what, year )
        what = what.lower()
        df = read_epa_download( csv, nrows=nrows, dropneg=(what != "temp") )

            # city, state, day -> daily average
        avdf = df.groupby( ["city", "state", "date"] ).mean()

            # drop cities with < 350 days
        if what != "temp":  # except Chicago, no temp
            gr = avdf.groupby( ["city", "state"] )
            nbefore = len(gr)
            avdf = gr.filter( lambda x: x.size >= mindays )
            nafter = len( avdf.groupby( ["city", "state"] ))
            print "dropped %d of %d cities with < %d days %s" % (
                    nbefore - nafter, nbefore, mindays, what )

        avdf = avdf.rename( columns={"av" : what} )
        avdf_list.append( avdf )
        if verbose:
            pr( "-- %s avdf --" % what, avdf )

    merge = pd.concat( avdf_list, axis=1 )

        # drop rows with > maxnan NaNs --
    nbefore = len(merge)
    merge = merge.dropna( thresh = len(whats) - maxnan )
    nafter = len(merge)
    print "dropped %d of %d rows with > %d NaNs " % (
            nbefore - nafter, nbefore, maxnan )
    merge = merge.reset_index()

    return merge, avdf_list


def counts_avs( df ):
    """ -> df: ndays, year avs no2 ozone temp wind """
    gr = df.drop( "date", axis=1 ) .groupby([ "city", "state" ])
    counts = gr.count()
    avs = gr.mean()
    counts.columns = "n_no2 n_ozone n_temp n_wind ".split()
    avs.columns = "av_no2 av_ozone av_temp av_wind ".split()
    navs = pd.concat( [counts, avs], axis=1 ) .reset_index() \
        .iloc[ :, [2,3,4,5,6,7,8,9,0,1] ]  # ndays avs city state
    return navs


#...............................................................................
if __name__ == "__main__":
    import sys

    np.set_printoptions( threshold=20, edgeitems=10, linewidth=140,
            formatter = dict( float=nu.numstr ))
    pd.set_option( "display.width", 140, "display.max_rows", 6,
            "display.precision", 3 )

    #...........................................................................
    year = 2016
    to = 0
    mindays = 350  # drop cities in no2 ozone temp wind with < __ days
    maxnan = 1  # drop rows with > __ NaN s
    nrows = 1e8
    save = "tmp"
    verbose = 1

    # to change these params, run this.py a=1 b=None 'c = "str"'  in sh or ipython
    for arg in sys.argv[1:]:
        exec( arg )

    print "\n", 80 * "="
    print " ".join(sys.argv)
    if to == 0:
        to = year
    print "year: %d to %d" % (year, to)

    from_ = "# %s\n#\n" % nu.From( __file__ )  # pwd, time
    csvheader = \
"""# NO2 ppb, ozone ppm, temp F, wind knots
# (NO2 1 ppb = 1.91 ug/m3,  Ozone 1 ppb = 2.0 ug/m3,  1 knot = .51444 m/s)
""" + from_

    #...........................................................................
        # read, merge, write_csv --
    for year in xrange( year, to+1 ):
        merge, avdfs = merge_epa_csvs( year=year,
                mindays=mindays, maxnan=maxnan, nrows=nrows, verbose=verbose )
        ndays_yearavs = counts_avs( merge )
        if verbose:
            pr( "-- merge_epa_csvs --", merge )
            pr( "-- ndays_yearavs --", ndays_yearavs )
        if save:
            out = "%s/%d-epa-no2-ozone-temp-wind.csv" % (save, year)
            write_csv( out, merge, csvheader=csvheader )

            out = "%s/%d-epa-ndays-yearavs.tsv" % (save, year)
            write_csv( out, ndays_yearavs, csvheader=from_ )  # sort-ndays

2016-epa-ndays-yearavs.tsv

epa_download

#!/bin/bash

url=http://aqsdr1.epa.gov/aqsweb/aqstmp/airdata

case $1 in "" | -h* | --h* )
    cat <<!
epa-download 2015 2016:
    curl daily_no2_2015.zip ... daily_WIND_2016.zip
		from $url
    link daily_44602_2015.zip to daily_no2_2015.zip  and 44201 to ozone

Notes:
    Many days are missing, see EPA-air-quality.md ff.

    pandas.read_csv( daily_no2_2015.zip with a single csv ) works.
!
    exit 1
esac

year=$1
yearto=${2-$year}

for year in `seq $year $yearto`
do
	for zip in daily_{42602,44201,TEMP,WIND}_$year.zip
	do
    if [[ ! -f $zip ]]; then
        curl --remote-time --show-error -O $url/$zip
    fi
	done

	ln daily_{44201,ozone}_$year.zip  2> /dev/null  # not mv, test -f
	ln daily_{42602,no2}_$year.zip  2> /dev/null
done

pr.py

#!/usr/bin/python
""" pr( ... x= y= ... ), probj()  for testing """
# must be thousands such
# use with https://docs.python.org/2.7/library/logging.html ?

from __future__ import division
import numpy as np

__version__ = "2016-03-06 mar  denis"

#...............................................................................
def pr( *args, **kwargs ):
    """ pr( "stage 42:", ntop=ntop, df=df ... )
        probj each arg / kwarg, summary line / ndarray
        apply liberally
    """
    print ""
    for arg in args:
        probj( arg )
    for key, val in sorted( kwargs.items() ):
            # sorted, else random not caller's order
        probj( val, key )
    # print ""

#...............................................................................
def probj( x, nm="" ):
    """ print any ? obj:
        string scalar list tuple dict: 1 summary line "nm: list len 42"
        pandas Series DataFrame with pd.set_option
        numpy ndarray with user's np.set_printoptions
    """
        # sure to be buggy
    if nm:
        print "-- %s:" %  nm ,

    if x is None  or _isstr( x ):
        print x
        return
    if np.isscalar( x ) \
    or (hasattr( x, "ndim" ) and x.ndim == 0):  # np.array( 3 )
        print "%.6g" % x
        return

    if hasattr( x, "name" )  and x.name is not None:  # pandas
        print x.name ,

        # type / classname, shape / len --
    t = getattr( x, "__class__", type(x) )
    print t.__name__ ,
    if hasattr( x, "shape" ):
        print x.shape
        n = x.shape[0]
    else:
        try:
            n = len(x)
            print "len %d" % n
        except (TypeError, AttributeError):  # len() of unsized object ?
            n = np.NaN
            print ""

    if isinstance( x, (tuple, list) ):
        return
    if isinstance( x, dict ):
        print "  keys:", sorted( x.keys() [:10] )
        return

        # pandas DataFrame, Series etc.
        # (to print as ndarray, pr( df.values ))
    if hasattr( x, "head" ):
        print x  # with user's pd.set_option( max_rows max_cols ... )
        return

    if hasattr( x, "values" ):
        x = x.values
    if hasattr( x, "dtype" ):  # np array  kind O ?
        print x  # with user's np.set_printoptions
    print ""


def _isstr( x ):
    """ basestring or np.array( "str" ) """
    return isinstance( x, basestring ) \
    or np.issubdtype( getattr( x, "dtype", "i4"), np.string_ )  # ?


#...............................................................................
if __name__ == "__main__":
    from collections import namedtuple

    np.set_printoptions( threshold=100, edgeitems=10, linewidth=140,
        formatter = dict( float = lambda x: "%.2g" % x ))  # float arrays %.2g
    # pd.set_option( "display.width", 140, "display.precision", 2 )

    class C:
        pass
    c = C()
    Namedtuple = namedtuple( "Namedtuple", "x y" )

    pr(
    adict   = { 1:2, 3:4 },
    alist   = [1, 2],
    array0d = np.array( 3 ),
    none    = None,
    arraynone = [None],
    arraystr = np.array("string"),
    array2  = np.array([ "string", 3 ]),
    eye     = np.eye( 3 ) * np.pi,
    pi      = np.pi,
    arraypi = np.array( np.pi ),
    s       = "string",

    C=C,
    c=c,

    Namedtuple = Namedtuple,
    anamedtuple = Namedtuple( 1, 2 ),
    )

    # also ~/py/etc/pandas/etc/z-groupby.py

zutil.py

#!/usr/bin/env python2
""" zutil.py: Bag asum findfile ints numstr vmap """

from __future__ import division
import numpy as np

class Bag( dict ):
    """ a dict with d.key short for d["key"], aka Dotdict
        d = Bag( k=v ... / **dict / dict.items() / [(k,v) ...] )  just like dict
    """
        # np.savez( **bag )
        # bag = Bag( np.load( .npz ).items() )

    def __init__(self, *args, **kwargs):
        dict.__init__( self, *args, **kwargs )
        self.__dict__ = self

    def __getnewargs__(self):  # for cPickle.dump( d, file, protocol=-1)
        return tuple(self)

Dotdict = Bag

#...............................................................................
def asum( X ):
    """ array summary: "shape dtype min av max [density]" """
    if not hasattr( X, "dtype" ) or X.dtype.kind not in "fiu":  # float int uint
        return str(X)
    if hasattr( X, "todense" ):  # issparse
        sparsetype = type(X).__name__ + " "  # csr_matrix etc.
        density = " of the %.3g %% non-0" % (
                100. * X.nnz / np.prod( X.shape ))
    else:
        sparsetype = density = ""
    dtype = str(X.dtype) .replace( "float64", "" )
    return "%s %s%s  min av max %.3g %.3g %.3g %s" % (
            X.shape, sparsetype, dtype, X.min(), X.mean(), X.max(), density )

def bincount( y ):
    return np.bincount( (y - y.min()) .astype(int) )

def copy( x ):
    return x if x is None  or np.isscalar(x) \
        else np.copy( x )

def findfile( filename, dirs="data ../data $SCIKIT_LEARN_DATA $webdata " ):
    """ -> first dir/file found in a list of dirs
        IOError if none
    """
    from os.path import expanduser, expandvars, isfile, join
    if isinstance( dirs, basestring ):
        dirs = dirs.split()
    for dir in [""] + dirs:  # "": /abspath ~/name name
        dirfile = expanduser( expandvars( join( dir, filename )))
        if isfile( dirfile ):
            return dirfile
    raise IOError( "file \"%s\" not found in folders %s" % (filename, dirs) )

def ints( X ):
    return np.round(X).astype(int)  # NaN Inf -> - maxint

def maxnorm( X, axis=None ):
    return np.linalg.norm( X, ord=np.inf, axis=axis )

def minavmax( X ):
    return "%.3g %.3g %.3g" % (
            X.min(), X.mean(), X.max() )

def mkdirpart( filename ):
    """ tmp/file: mkdir tmp """
    import os
    dir = os.path.dirname( filename )
    if dir and not os.path.isdir( dir ):
        os.makedirs( dir )

def numstr( x, fmt="%.2g" ):
    """ %.2g but 123 not 1.2e+02 """
        # np.set_printoptions( ... formatter = dict( float=ut.numstr ))
    if isinstance( x, basestring ) or x is None:
        return x
    if np.isscalar(x):
        if abs(x) < 1e-10:      return "0"
        if abs(x) < 1e-6:       return "%.1g" % x
        if 100 <= abs(x) < 1e5: return "%.0f" % x
        return                         fmt % x
    x = np.asanyarray(x)
    assert x.ndim == 1, x.ndim  # NotImplemented
    return " ".join([ numstr( _, fmt ) for _ in x ])

def numstr3( x ):
    return numstr( x, fmt="%.3g" )

_pvec = 3
_pvecscale = 0  # > 0: ints( scale * x )

def pvec( x ):
    """ -> str( short vecs ), smax( long ) """
    if x is None or isinstance( x, basestring ):
        return x
    if np.isscalar(x):
        return "%.3g" % x
    if len(x) <= _pvec:
        if _pvecscale > 0:
            x = ints( _pvecscale * np.asarray(x) )
        return str( x )  # np.printoptions
    else:
        return "%.3g" % smax( x )

def smax( x ):
    """ signed max, to see oscillation """
    return x if np.isscalar(x) or x.ndim == 0 \
        else x[ np.fabs(x) .argmax() ]

def vmap( f, X ):
    return np.vstack( map( f, X )) .squeeze()

#...............................................................................
def From( pyfile ):
    """ -> "callersfile pwd time" for tagging output files """
    import os
    pyfile = pyfile.split("/") [-1]
    pwd = os.getenv("PWD") .replace( "/Users/", "~", 1 )  # mac -L not -P
    return "from %s  run in %s  %s" % (pyfile, pwd, isoday())

def isoday():
    import time
    return time.strftime( "%Y-%m-%d %h %H:%M" )  # 2011-11-15 Nov 12:06