<jupyter_start><jupyter_text>Executive Summary<jupyter_code># This analysis looks at raw wind data from HPN from 2006-18 as recorded by the ASOS station on the field
# Data is recorded every hour from 2006-15 and every 5 minutes from 2016-18

# Data shows that during the last few years winds have favored 16 consistently more than they have during earlier periods
# In aggrigate during periods when the airport is most active 16 was recently favored ~25% more than 10+ years ago

# For 2013, 2015 and 2018 (years of the HMMH study) 16 winds increased from 28% to 34, and increase of 21%
# The above is in line with the increase in jet traffic observed by the study on 16 (+24% for 'day' and +18% for 'night')

# Based on the trend of winds favoring 16 along with the observed increases in actual traffic on 16 it is reasonable to conclude that changes in wind patterns is a primary, if not the primary, cause of increased traffic on the 16 approach path in recent years
# The code below will pull all the raw data from source, and conduct analysis and produce charts as reviewed by the HPN AAB<jupyter_output><empty_output><jupyter_text>Import Required Packages<jupyter_code>from __future__ import print_function
import json
import time
import datetime
from urllib.request import urlopen

import os
import math
import matplotlib
import pandas as pd
import openpyxl<jupyter_output><empty_output><jupyter_text>Define Functions  Data Fetching Service<jupyter_code># This service goes to the site below and automatically downloads data flat files into the /data folder

SERVICE = "http://mesonet.agron.iastate.edu/cgi-bin/request/asos.py?"

# Inspired by https://github.com/akrherz/iem/blob/master/scripts/asos/iem_scraper_example.py

# Number of attempts to download data
MAX_ATTEMPTS = 6

def download_data(uri):
    """Fetch the data from the IEM
    The IEM download service has some protections in place to keep the number
    of inbound requests in check.  This function implements an exponential
    backoff to keep individual downloads from erroring.
    Args:
      uri (string): URL to fetch
    Returns:
      string data
    """
    attempt = 0
    while attempt < MAX_ATTEMPTS:
        try:
            data = urlopen(uri, timeout=300).read().decode("utf-8")
            if data is not None and not data.startswith("ERROR"):
                return data
        except Exception as exp:
            print("download_data(%s) failed with %s" % (uri, exp))
            time.sleep(5)
        attempt += 1

    print("Exhausted attempts to download, returning empty data")
    return ""

def fetchASOS(stations, startts, endts):
    service = SERVICE + "data=all&tz=Etc/UTC&format=comma&latlon=yes&"

    service += startts.strftime("year1=%Y&month1=%m&day1=%d&")
    service += endts.strftime("year2=%Y&month2=%m&day2=%d&")
    
    for station in stations:
        uri = "%s&station=%s" % (service, station)
        print("Downloading: %s" % (station,))
        data = download_data(uri)
        outfn = "%s_%s_%s_%s.txt" % (
            'Data/',
            station,
            startts.strftime("%Y%m%d%H%M"),
            endts.strftime("%Y%m%d%H%M"),
        )
        out = open(outfn, "w")
        out.write(data)
        out.close()
        print("Wrote output: %s" % (outfn,))
        return str(outfn)<jupyter_output><empty_output><jupyter_text>Fetch for Defined Years and Airports<jupyter_code># Uses the above service to fetch data. Takes a list of airport IDs and years as inputs.

def importASOS(years, airports):
    df = pd.DataFrame()
    
    for airport in airports:
        print("Processing Airport: %s" % (airport,))
        for year in years:
            outfn = fetchASOS([airport], datetime.datetime(year, 1, 1), datetime.datetime(year, 12, 31))
            df = pd.concat([df, pd.read_csv(outfn, skiprows=5)])
            
            # Drop unused columns
            df = df[['station', 'valid', 'drct', 'sknt', 'vsby', 'skyc1', 'skyl1']]
            
            # Set the valid column as a timestamp
            df['valid'] = pd.to_datetime(df['valid'])
            
            # Add a column for the year
            df['year'] = df['valid'].dt.year
    
    return df<jupyter_output><empty_output><jupyter_text>Fetch Data for Analysis<jupyter_code>airports = ['HPN']
years = [2006, 2007, 2008, 2009, 2010, 2011, 2012 ,2013, 2014, 2015, 2016, 2017, 2018, 2019]

#airports = ['HPN', 'LGA']
#years = [2013,2015,2018]

df = importASOS(years,airports)

# Pre-calculate some parameters to aid in later filtering and grouping
df['valid'] = pd.to_datetime(df['valid'])
df['hour'] = df['valid'].dt.hour
df['day'] = df['valid'].dt.dayofweek
df['month'] = df['valid'].dt.month
df.reset_index(drop=True, inplace=True)
# We quickly check the number of datapoints recorded for each year
# Note that in 2016 the data goes from being recorded every hour to every 5 minutes

HPNg = df.groupby('year').count().station
HPNg = pd.DataFrame(HPNg)
HPNg.unstack().plot(kind='bar', stacked=True).legend(loc='center left', bbox_to_anchor=(1, 0.5))<jupyter_output><empty_output><jupyter_text>Define Wind Roses<jupyter_code># The wind rose below is defined around the main 16/34 runway at HPN. 
# ASOS data is 'true' and runway head headings are magnetic. HPN deviation is 13W.
# Wind readings are reported in 10 degree bins so 10 degress are subtracted from the ASOS readings to allign with runway orientation.

def HPN_Wind(airport, winddir):
    rwy = ''
    
    if airport == 'HPN':
        # A small number of readings are missing (likely due to station maintenance). These are omitted from the analysis
        if winddir == 'M':
            rwy = 'M'

        # The ASOS data records calm winds as having a heading of 0.00; North winds have a heading of 360.00
        if winddir in ['0.00']:
            rwy = 'Calm'
            
        # These wind directions would be direct crosswinds on the 16-34 runway    
        if winddir in ['240.00', '60.00']:
            rwy = 'Crosswind'

        # These wind directions would provide a headwind component favoring use of runway 16
        if winddir in ['70.00', '80.00', '90.00', '100.00', '110.00', '120.00', '130.00', '140.00', '150.00', '160.00', '170.00', '180.00', '190.00', '200.00', '210.00', '220.00', '230.00']:
            rwy = '16'

        # These wind directions would provide a headwind component favoring use of runway 34
        if winddir in ['250.00', '260.00', '270.00', '280.00', '290.00', '300.00', '310.00', '320.00', '330.00', '340.00', '350.00', '360.00', '10.00', '20.00', '30.00', '40.00', '50.00']:
            rwy = '34'
    else:
        rwy = 'NO ROSE'
            
    return(rwy)<jupyter_output><empty_output><jupyter_text>Analysis  Determine Rwy Preference For Each Period<jupyter_code># Remove rows with missing wind direction data
df = df[df['drct'] != 'M']
for i in df.index:
    if df.at[i, 'station'] == 'HPN':
        df.at[i, 'Raw_Rwy_Preference'] = HPN_Wind(df.at[i, 'station'], df.at[i, 'drct'])
    else:
        df.at[i, 'Raw_Rwy_Preference'] = 'OTHER AIRPORT'
# View the raw data with runway preference added
df<jupyter_output><empty_output><jupyter_text>Native Runway Preference for HPN  All Times<jupyter_code># Looking at all the data, we can see that the last several years have seen 16 winds much higher than historical levels going back to 2006

HPNg = df[df['station']=='HPN'].groupby(['year', 'Raw_Rwy_Preference']).count().station
HPNg = HPNg.groupby(level=0).apply(lambda x:100 * x / float(x.sum()))
HPNg = pd.DataFrame(HPNg)
HPNg.unstack().plot(kind='bar', stacked=True, title='HPN Wind Direction by Favored Runway (All Times)', figsize=(7,7)).legend(loc='center left', bbox_to_anchor=(1, 0.5), title='Legend')
HPNg.to_excel("GraphData/HPN_Wind_All_Times.xlsx")<jupyter_output><empty_output><jupyter_text>All Times - By Hour of Day (UTC)  All Times - By Month<jupyter_code># Looking at the above by hour of day we see that 16 winds are most common during daylight hours, peaking around 20Z (2 or 3 PM local)
# Those afternoon periods also have very few periods with calm winds

HPNg = df[df['station']=='HPN'].groupby(['hour', 'Raw_Rwy_Preference']).count().station
HPNg = HPNg.groupby(level=0).apply(lambda x:100 * x / float(x.sum()))
HPNg = pd.DataFrame(HPNg)
HPNg.unstack().plot(kind='bar', stacked=True, title='HPN Wind Direction by Favored Runway (All Times; By Hour of Day UTC)', figsize=(7,7)).legend(loc='center left', bbox_to_anchor=(1, 0.5), title='Legend')
HPNg.to_excel("GraphData/HPN_Wind_All_Times-By_Hour.xlsx") 
# Looking at the above by month we see some seasonality with spring into summer having more 16 winds than cooler months of the year

HPNg = df[df['station']=='HPN'].groupby(['month', 'Raw_Rwy_Preference']).count().station
HPNg = HPNg.groupby(level=0).apply(lambda x:100 * x / float(x.sum()))
HPNg = pd.DataFrame(HPNg)
HPNg.unstack().plot(kind='bar', stacked=True, title='HPN Wind Direction by Favored Runway (All Times; By Month)', figsize=(7,7)).legend(loc='center left', bbox_to_anchor=(1, 0.5), title='Legend')
HPNg.to_excel("GraphData/HPN_Wind_All_Times-By_Month.xlsx")<jupyter_output><empty_output><jupyter_text>All Times - By Day of Week (0 = Monday, 6 = Sunday)<jupyter_code># Looking across the 14 years we see little consistent seasonality by day of week (as expected)
# As we will see later, certain years can develop localized patterns that do favor shorter term preferences (e.g. weekends)

HPNg = df[df['station']=='HPN'].groupby(['day', 'Raw_Rwy_Preference']).count().station
HPNg = HPNg.groupby(level=0).apply(lambda x:100 * x / float(x.sum()))
HPNg = pd.DataFrame(HPNg)
HPNg.unstack().plot(kind='bar', stacked=True, title='HPN Wind Direction by Favored Runway (All Times; By Day of Week)', figsize=(7,7)).legend(loc='center left', bbox_to_anchor=(1, 0.5), title='Legend')
HPNg.to_excel("GraphData/HPN_Wind_All_Times-By_DayofWeek.xlsx")<jupyter_output><empty_output><jupyter_text>Excluding Overnight Periods<jupyter_code># All the above data looks at all hours of the day. But few flights arrive during the overnight hours (5Z - 12Z)
# The rest of the analysis below will exclude the overnight periods and focus on the times when the airport is generally active

dfControlled = df[(df['hour'] > 11) | (df['hour'] < 4)]
# Looking at the overal trend when overnight periods are excluded we again see recent years with high 16 winds

HPNg = dfControlled[dfControlled['station']=='HPN'].groupby(['year', 'Raw_Rwy_Preference']).count().station
HPNg = HPNg.groupby(level=0).apply(lambda x:100 * x / float(x.sum()))
HPNg = pd.DataFrame(HPNg)
HPNg.unstack().plot(kind='bar', stacked=True, title='HPN Wind Direction by Favored Runway (Excl Overnight; By Year)', figsize=(7,7)).legend(loc='center left', bbox_to_anchor=(1, 0.5), title='Legend')
HPNg.to_excel("GraphData/HPN_Wind_Excl_Overnight-By_Year.xlsx")<jupyter_output><empty_output><jupyter_text>Excluding Overnight - 'Outdoor Months' (Apr, May, Jun, Jul, Aug, Sep, Oct)<jupyter_code>dfControlledOutdoors = dfControlled[(dfControlled['month'] <= 10) & (dfControlled['month'] >= 4)]
# Here we look just at 'outdoor months' where people are most likly to be outside during the daytime

HPNg = dfControlledOutdoors[dfControlledOutdoors['station']=='HPN'].groupby(['year', 'Raw_Rwy_Preference']).count().station
HPNg = HPNg.groupby(level=0).apply(lambda x:100 * x / float(x.sum()))
HPNg = pd.DataFrame(HPNg)
HPNg.unstack().plot(kind='bar', stacked=True, title='HPN Wind Direction by Favored Runway (Outdoor Months, Excl Overnight)', figsize=(7,7)).legend(loc='center left', bbox_to_anchor=(1, 0.5), title='Legend')
HPNg.to_excel("GraphData/HPN_Wind_Excl_Overnight_OutdoorMonths-By_Year.xlsx")<jupyter_output><empty_output><jupyter_text>'Outdoor Months' by Hour of Day (UTC)<jupyter_code>dfOutdoors = df[(df['month'] <= 10) & (df['month'] >= 4)]
# Looking by hour of day we see that overnight periods have much less 16 winds and >20% calm winds... during afternoons there are few calm winds and 16 is the most common wind

HPNg = dfOutdoors[dfOutdoors['station']=='HPN'].groupby(['hour', 'Raw_Rwy_Preference']).count().station
HPNg = HPNg.groupby(level=0).apply(lambda x:100 * x / float(x.sum()))
HPNg = pd.DataFrame(HPNg)
HPNg.unstack().plot(kind='bar', stacked=True, title='HPN Wind Direction by Favored Runway (Outdoor Months; By Hour of Day)', figsize=(7,7)).legend(loc='center left', bbox_to_anchor=(1, 0.5), title='Legend')
HPNg.to_excel("GraphData/HPN_Wind_OutdoorMonths-By_Hour.xlsx")<jupyter_output><empty_output><jupyter_text>Excluding Overnight - 'Outdoor Months' (Apr, May, Jun, Jul, Aug, Sep, Oct) Sat & Sun Only<jupyter_code>dfControlledOutdoorsWeekend = dfControlled[(dfControlled['month'] <= 10) & (dfControlled['month'] >= 4) & ((dfControlled['day'] == 5) | (dfControlled['day'] == 6))]
# Although across the whole dataset weekday has little impact, just looking at weekends during outdoor months (likely the most sensitive time) does show some significant yearly changes

HPNg = dfControlledOutdoorsWeekend[dfControlledOutdoorsWeekend['station']=='HPN'].groupby(['year', 'Raw_Rwy_Preference']).count().station
HPNg = HPNg.groupby(level=0).apply(lambda x:100 * x / float(x.sum()))
HPNg = pd.DataFrame(HPNg)
HPNg.unstack().plot(kind='bar', stacked=True, title='HPN Wind Direction by Favored Runway (Excl Overnight; Outdoor Months; Weekends Only)', figsize=(7,7)).legend(loc='center left', bbox_to_anchor=(1, 0.5), title='Legend')
HPNg.to_excel("GraphData/HPN_Wind_Excl_Overnight_WarmMonths_Weekends-By_Year.xlsx")<jupyter_output><empty_output><jupyter_text>Excluding Overnight - 'Indoor Months' (Nov, Dec, Jan, Feb, Mar)<jupyter_code>dfControlledIndoors = dfControlled[(dfControlled['month'] >= 11) | (dfControlled['month'] <= 3)]
HPNg = dfControlledIndoors[dfControlledIndoors['station']=='HPN'].groupby(['year', 'Raw_Rwy_Preference']).count().station
HPNg = HPNg.groupby(level=0).apply(lambda x:100 * x / float(x.sum()))
HPNg = pd.DataFrame(HPNg)
HPNg.unstack().plot(kind='bar', stacked=True, title='HPN Wind Direction by Favored Runway (Excl Overnight; Indoor Months)', figsize=(7,7)).legend(loc='center left', bbox_to_anchor=(1, 0.5), title='Legend')
HPNg.to_excel("GraphData/HPN_Wind_Excl_Overnight_IndoorMonths-By_Year.xlsx")<jupyter_output><empty_output><jupyter_text>Excluding Overnight - 'Indoor Months' (Nov, Dec, Jan, Feb, Mar) Sat & Sun Only<jupyter_code>dfControlledIndoorsWeekend = dfControlled[((dfControlled['month'] >= 11) | (dfControlled['month'] <= 3)) & ((dfControlled['day'] == 5) | (dfControlled['day'] == 6))]
# The same analysis can be done looking at weekends during 'indoor' months; Of interest is that while 2016 had the most 16 winds in the summer weekends it had the least 16 winds on weekends during indoor months

HPNg = dfControlledIndoorsWeekend[dfControlledIndoorsWeekend['station']=='HPN'].groupby(['year', 'Raw_Rwy_Preference']).count().station
HPNg = HPNg.groupby(level=0).apply(lambda x:100 * x / float(x.sum()))
HPNg = pd.DataFrame(HPNg)
HPNg.unstack().plot(kind='bar', stacked=True, title='HPN Wind Direction by Favored Runway (Excl Overnight; Indoor Months; Weekends Only)', figsize=(7,7), color=['#ff7f0e', '#2ca02c', '#d62728', '#9467bd']).legend(loc='center left', bbox_to_anchor=(1, 0.5), title='Legend')
HPNg.to_excel("GraphData/HPN_Wind_Excl_Overnight_IndoorMonths_Weekends-By_Year.xlsx")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13663829510753037, 'fraction_numerical': 0.034110843999740105, 'mean_word_length': 4.724060989215322, 'pattern_counts': {'":': 0, '<': 72, '<?xml version=': 0, '>': 72, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '223939', 'license': 'MIT', 'n_tokens_mistral': 5628, 'n_tokens_neox': 5150, 'n_words': 1487, 'path': 'HPN Wind Rose Analysis.ipynb', 'repo_name': 'nthartman/HPN_Wind_Analysis', 'size': '15391'}