<jupyter_start><jupyter_text>Download all the data needed for the projectDownloaded files will be stored in downloads/ for faster processing next timeDownloading the first time may take much longer due to the source files being rate-limitedData will be minimally converted to dataframes and saved for the next workbook Note: if you got this from my github, you don't have to run thisThe data has all been saved in the dataframes/ directory in both the (slow and safe) csv.gz format and the (fast and dangerous) .pickle.gz format Run time: 1 minute (if data already downloaded)  *  * * *  * *  Get Hourly Electric Demand dataOperational Data is provided publicly by Southwest Power Pool, Inc.<jupyter_code>import urllib.request
import os

downloads_dir = 'downloads/'
years_to_load = ['2015', '2016', '2017', '2018', '2019']

sourceurl = 'ftp://pubftp.spp.org/Operational_Data/HourlyLoad/%s/%s.zip'

for year in years_to_load:
    if os.path.exists(downloads_dir + year + '.zip') == False: 
        print ("Trying to retrieve hourly demand for " + year, end=' ... ', flush=True)
        urllib.request.urlretrieve(sourceurl%(year, year), downloads_dir + year + ".zip")
        print ("done")
    else:
        print (downloads_dir + year + ".zip already downloaded")

from zipfile import ZipFile
import pandas as pd
import numpy as np
import io

try: 
    del demand_df
except: 
    pass

for year in years_to_load:
    with ZipFile(downloads_dir + year + '.zip', 'r') as zipObj:
        # Get list of files names in zip
        listOfiles = zipObj.namelist()
        # Iterate over the list of file names in given list & print them
        for elem in listOfiles:
#            print(elem)
            if (elem.split('.')[-1:][0] == 'csv'):
                df = pd.read_csv(io.BytesIO( zipObj.read(elem)) ) 
                try: 
                    demand_df = demand_df.append(df)
                except: 
                    demand_df = df

demand_df = demand_df.reset_index(drop=True)
demand_df.head()
# Yikes, all my column names start with spaces. Fix that real quick:
renamedict={}
for i in demand_df.columns: 
    if i.startswith(' '):
        renamedict[i] = i[1:]
demand_df.rename(columns=renamedict, inplace=True)
demand_df.columns
demand_df.count()
# drop rows where MarketHour is all blanks (https://stackoverflow.com/questions/29314033/drop-rows-containing-empty-cells-from-a-pandas-dataframe)

demand_df['MarketHour'].replace('', np.nan, inplace=True)
demand_df.dropna(subset=['MarketHour'], inplace=True)
demand_df.count()

# also there are a handful of bad rows where CSWS and others are blank
demand_df['CSWS'].replace('', np.nan, inplace=True)
demand_df.dropna(subset=['CSWS'], inplace=True)
demand_df.count()<jupyter_output><empty_output><jupyter_text>What is the last MarketHour really?  And what is that silly time format? I am amazed; I have spent weeks of my life dealing with bad timestamp formats and this one line of code fixes them<jupyter_code>pd.to_datetime(['2021-09-01 17:00', '16-Jul-1970 3:30PM', '9/11/2001 8:46AM', '3/4/5 6:7' ])
# wow.
demand_df['MarketHour'] = pd.to_datetime(demand_df.MarketHour)
demand_df.sort_values(by='MarketHour', inplace=True, ignore_index=True)
demand_df.tail()<jupyter_output><empty_output><jupyter_text>MarketHour is in an industry format known as Hour Ending GMT.So the demand reported for 2020-01-01 06:00:00 is summed over the last hour of 2019, which ends at the start of the first hour of 2020.  This is an inconvenient format because grouping to a higher time level like day or month is difficult. So, I want to subtract an hour from each MarketHour and call it HourBeginUTC so we are all clear on what it means, and all the hours for a given day match that day.<jupyter_code>demand_df.MarketHour -= pd.Timedelta("1 hour")
demand_df.rename(columns={"MarketHour":"HourBeginUTC"}, inplace=True)

# Convert UTC to US/Central
# OMG this is an ugly hack please refactor correctly
demand_df['HourBeginCPT']=demand_df.HourBeginUTC.dt.tz_localize('utc').dt.tz_convert('US/Central')
demand_df.tail()<jupyter_output><empty_output><jupyter_text>I see some duplicate rows, due to the way the files are stored on the source site.De-duplicate:<jupyter_code>demand_df.drop_duplicates(subset='HourBeginUTC', keep='first', ignore_index=True, inplace=True)
demand_df.tail()
print ("Expecting 5 years of hours, or about ", 365.25*5*24, "observations")
demand_df.count()<jupyter_output>Expecting 5 years of hours, or about  43830.0 observations<jupyter_text>The original analysis explained footprint-wide demand based on temperature.  Add a "SPP" column that is the fooprint wide demand:<jupyter_code>demand_df['SPP'] = demand_df.sum(axis=1, skipna=True, level=None, numeric_only=True, min_count=0)
demand_df.head()<jupyter_output><empty_output><jupyter_text>One company WAUE has fewer rows; it probably joined at a later date. See if we can tell when:<jupyter_code>demand_df[['HourBeginUTC', 'WAUE']].dropna()
# That's still a lot of data; just go with the uniform footprint by dropping all N/A columns for now.
demand_df = demand_df.dropna().reset_index(drop=True)
demand_df.head()<jupyter_output><empty_output><jupyter_text>It is always useful to visualize the data to see if it makes sense<jupyter_code># from https://pandas.pydata.org/pandas-docs/dev/getting_started/intro_tutorials/09_timeseries.html
# which has other good tips

import matplotlib.pyplot as plt
fig, axs = plt.subplots(figsize=(12, 4))
demand_df.groupby(demand_df["HourBeginCPT"].dt.hour)["SPP"].mean().plot(
    kind='bar', rot=0, ax=axs
)
plt.xlabel("CPT Hour of the day");  # custom x label using matplotlib
plt.ylabel("Mean Demand by Hour, MWh");<jupyter_output><empty_output><jupyter_text>Save out the hourly demand dataframe<jupyter_code># that's pretty clean. Save it off for further verification. 

# CSV has problems with datatypes and indexes:
demand_df.to_csv('dataframes/demand_df.csv.gz', compression='infer', index=False)

# Pickle works perfectly if you don't have to change Python versions (such as PC to Google Colab)
# demand_df.to_pickle('dataframes/demand_df.pickle.gz', compression='infer')

# Trying Feather
#demand_df.to_feather('dataframes/demand_df.feather.zstd', compression='zstd')

# Trying JSON
# demand_df.to_json('dataframes/demand_df.json.gz', compression='infer')<jupyter_output><empty_output><jupyter_text>The original analysis predicted total daily demand, not total hourly demand. To compare with the orignal, I'll sum to the daily level too and save that<jupyter_code># this didn't work; I need to truncate the hour to the day level before the group by 
# demand_df.groupby(['HourBeginUTC'])['SPP'].sum().reset_index()

# maybe use The resample() method is similar to a groupby operation:

#    it provides a time-based grouping, by using a string (e.g. M, 5H,…) that defines the target frequency

#    it requires an aggregation function such as mean, max,…

# this still doesn't work: 
# daily_sum = demand_df.resample("M").sum()

# but I sure can cheat with SQLite: https://pythonspeed.com/articles/indexing-pandas-sqlite/
import sqlite3
con=sqlite3.connect("")
demand_df.to_sql('demand_df', con, if_exists='replace', index=False)

daily_demand_df=pd.read_sql("""
select datetime( HourBeginCPT, \'start of day\') as opday, 
sum(SPP) as sum_spp_load
from demand_df
group by datetime( HourBeginCPT, \'start of day\')
""", con)
daily_demand_df['opday'] = pd.to_datetime(daily_demand_df.opday)

daily_demand_df.head()

# todo:  replace this with one line of Pandas
# from https://pandas.pydata.org/pandas-docs/dev/getting_started/intro_tutorials/09_timeseries.html
# which has other good tips

import matplotlib.pyplot as plt
fig, axs = plt.subplots(figsize=(12, 4))
daily_demand_df.groupby(daily_demand_df["opday"].dt.month)["sum_spp_load"].mean().plot(
    kind='bar', rot=0, ax=axs
)
plt.xlabel("Month");  # custom x label using matplotlib
plt.ylabel("Mean Daily Demand");<jupyter_output><empty_output><jupyter_text>WAUE started 10/1/2015 06:00. So we need to either drop rows before then, or feed another weather station into the model.  We probably have enough data wihout those values before 10/1/2015.   Save out the daily demand dataframe<jupyter_code># save daily_demand_df 

daily_demand_df.to_csv('dataframes/daily_demand_df.csv.gz', compression='infer', index=False)

# daily_demand_df.to_pickle('dataframes/daily_demand_df.pickle.gz', compression='infer')
# Trying Feather
# daily_demand_df.to_feather('dataframes/daily_demand_df.feather.zstd', compression='zstd')
# Trying JSON
# daily_demand_df.to_json('dataframes/daily_demand_df.json.gz', compression='infer')<jupyter_output><empty_output><jupyter_text>*  * * *  * *  Now look for weather station dataThese daily weather datasets are documeted at https://www.ncei.noaa.gov/data/global-historical-climatology-network-daily/doc/  Get a list of all weather stations with available data<jupyter_code>sourceurl = 'https://www.ncei.noaa.gov/data/global-historical-climatology-network-daily/doc/ghcnd-stations.txt'
save_file = downloads_dir + 'ghcnd-stations.txt'
if os.path.exists(save_file) == False: 
    print ("Trying to retrieve " + sourceurl + " to " + save_file, end=' ... ', flush=True)
    urllib.request.urlretrieve(sourceurl, save_file)
    print ("done")
else:
    print (save_file + " already downloaded")
    
# load station file from NOAA site; it is fixed-width formatted.

st_df = pd.read_fwf(save_file, 
                 widths=[11,1,8,1,9,1,6,1,2,1,30,1,3,1,3,1,5])
    
# rename columns and drop spaces
st_df.rename(columns={'ACW00011604' : 'station',
                   '17.1167' : 'latitude',
                   '-61.7833' : 'longitude',
                   '10.1' : 'altitude?',
                   'ST JOHNS COOLIDGE FLD' : 'station name',
                   'Unnamed: 16' : 'zip code',
                   'Unnamed: 8' : 'state code',
                   }, inplace=True)

st_df = st_df[['station', 'station name', 'latitude', 'longitude', 'state code', 'zip code']]<jupyter_output><empty_output><jupyter_text>We are interested only in weather stations that could impact demand in the utility footprint*. Very roughly, this is a box from latitudes 31 to 100, and longitudes -106 to -92.  If the station is outside that box, we can not consider it a predictor of demand in SPP.  This excludes most of the stations in the list.(*) I thought this was true but it may not be.  For instance, Little Rock is consistently one of the most predictive stations, as well as others outside the untility.  The electric grid is interconnected and demand may be reflecting those outside connections.<jupyter_code>st_df = st_df.loc[(st_df['longitude'] >=  -106) &
        (st_df['longitude'] <=  -92)  & 
        (st_df['latitude'] >=  31)  & 
        (st_df['latitude'] <=  100) &
        (st_df['station'].str.contains("USW"))
       ]
st_df<jupyter_output><empty_output><jupyter_text>Make the station names label-friendlySince I will be using station name as a label for predictor variables, I need to define a function to transform the station name as listed into something that can be referenced as a label.<jupyter_code># define a lookup from the station ID to the station name

def station_name(station): 
    df = st_df[st_df.station == station][['state code','station name']]
#    print (str(df.values))
    return '_'.join(df.values[0]).replace(' ', '_')
    

# The initial analysis used only these stations: 

for station in [ 
    'USW00013967',  # Oklahoma City OK
    'USW00013996',  # Topeka KS 
    'USW00023047',  # Amarillo TX
    'USW00013968',  # Tulsa OK
    'USW00014939',  # Lincoln NE 
    'USW00024011',  # Bismark ND
    ]: 
    print (f'Station {station} is {station_name(station)}')<jupyter_output>Station USW00013967 is OK_OKLAHOMA_CITY_WILL_ROGERS_AP
Station USW00013996 is KS_TOPEKA_MUNI_AP
Station USW00023047 is TX_AMARILLO
Station USW00013968 is OK_TULSA_INTL_AP
Station USW00014939 is NE_LINCOLN_MUNI_AP
Station USW00024011 is ND_BISMARCK<jupyter_text>Download selected weather stations.<jupyter_code># Temporarily load them ALL for a possible genetic-algorithm search for the best ones 

stations_to_load = st_df.station

sourceurl = "https://www.ncei.noaa.gov/data/global-historical-climatology-network-daily/access/%s.csv"

for station in stations_to_load: 
    absolute_url = sourceurl%station
    save_file = downloads_dir + station + ".csv"
    
    if os.path.exists(save_file) == False: 
        print ("Trying to retrieve " + absolute_url + " to " + save_file, end=' ... ', flush=True)
        urllib.request.urlretrieve(absolute_url, save_file)
        print ("done")
    else:
        print (save_file + " already downloaded")<jupyter_output>downloads/USW00003013.csv already downloaded
downloads/USW00003017.csv already downloaded
downloads/USW00003024.csv already downloaded
downloads/USW00003026.csv already downloaded
downloads/USW00003027.csv already downloaded
downloads/USW00003028.csv already downloaded
downloads/USW00003030.csv already downloaded
downloads/USW00003031.csv already downloaded
downloads/USW00003047.csv already downloaded
downloads/USW00003054.csv already downloaded
downloads/USW00003055.csv already downloaded
downloads/USW00003063.csv already downloaded
downloads/USW00003067.csv already downloaded
downloads/USW00003068.csv already downloaded
downloads/USW00003072.csv already downloaded
downloads/USW00003078.csv already downloaded
downloads/USW00003082.csv already downloaded
downloads/USW00003083.csv already downloaded
downloads/USW00003086.csv already downloaded
downloads/USW00003088.csv already downloaded
downloads/USW00003089.csv already downloaded
downloads/USW00003090.csv already downloaded
downloads/[...]<jupyter_text>If the quality of the data for a weather station is good, append it to a single weather dataframe<jupyter_code>try: 
    del raw_weather_df
except: 
    pass

for station in stations_to_load:
    save_file = downloads_dir + station + ".csv"
    
    # print ("reading " + save_file)
    df = pd.read_csv(save_file, low_memory=False,  parse_dates=["DATE"]) 
    
    # trim here such that DATE is greater or equal to our first load year
    df = df.loc[(df['DATE'] >= pd.to_datetime(years_to_load[0]))]

    # and not past our last year to load 
    df = df.loc[(df['DATE'] <= (pd.to_datetime(years_to_load[-1]) + pd.Timedelta(days=366)))]
    
    
    if 'TMIN' not in df.columns:  # we can't really use this one; we at least need TMIN and TMAX
        print (f'file {save_file} has no TMIN attribute; skipping')
        continue
        
    if 'TAVG' not in df.columns: 
        df['TAVG'] = (df.TMIN+df.TMAX)/2.0
        
    # if TAVG is a column but not populated, populate it: 
    df['TAVG'] = np.where(df['TAVG'].isnull(), (df['TMIN'] + df['TMAX']) / 2.0, df['TAVG'] )
    
    good_station_name = station_name(station)
    df['STATION'] = good_station_name
    
    # cut it down to just columns we wish to keep
    df = df[['STATION', 'DATE', 'TMAX', 'TMIN', 'TAVG' ]]
    
    if df.dropna().count().min() >= 1827:  # the specific number of observations I am expecting... FRAGILE
        print ("Adding " + good_station_name)
    else: 
        print ("Station " + good_station_name + " has too many NA values; skipping")
        continue
    
    if 'raw_weather_df' not in vars():  # first time through  
        raw_weather_df = df
    else:  # append to existing df
        raw_weather_df = raw_weather_df.append(df)  
  
raw_weather_df.groupby(['STATION']).size()
raw_weather_df[['TMAX', 'TMIN', 'TAVG']].describe()<jupyter_output><empty_output><jupyter_text>Extract DATE as a timestamp, and TMIN/TMAX/TAVG as degrees CAlso choose here which attributes to keep. The old analysis only used "heating degree days" and "cooling degree days", which was calculated as difference of (tmax+tmin)/2 from 65F.   What exactly are "Degree Days"? Here is a very good paper on using Degree Days as a metric for energy use estimation, along with some of the problems of this technique:  https://www.energylens.com/articles/degree-days<jupyter_code># Drop observations that don't have a TMIN weather component
weather_df = raw_weather_df[['STATION', 'DATE', 'TMIN', 'TMAX', 'TAVG']].dropna(subset=['TMIN'])

# covert the date to a datetime
weather_df['DATE'] = pd.to_datetime(weather_df.DATE)

weather_df['TMIN'] /= 10.0 # convert from 1/10th of degrees C to degrees C
weather_df['TMAX'] /= 10.0
weather_df['TAVG'] /= 10.0

# convert degrees C to degrees F to calculate HDD around 65.0 degrees F, which is commen in the US per the paper above
def c_to_f(x): 
    return (x * 9.0 / 5.0 + 32.0)

# add values of HDD and CDD consistent with the old analysis (65 F)
weather_df['CDD'] = weather_df['TAVG'].apply(lambda x: c_to_f(x) - 65.0 if c_to_f(x) > 65.0 else 0)  
weather_df['HDD'] = weather_df['TAVG'].apply(lambda x: 65.0 - c_to_f(x) if c_to_f(x) < 65.0 else 0)  

weather_df.sort_values(by='DATE', inplace=True, ignore_index=True)<jupyter_output><empty_output><jupyter_text>Show a sample of the combined weather DF<jupyter_code>weather_df.sample(n=20)
weather_df.describe()
# This is what the old analysis used the average of TMIN,TMAX as TAVG.
# Use it because all stations have it
# weather_df['TAVG'] = (weather_df['TMIN'] + weather_df['TMAX']) / 2.0<jupyter_output><empty_output><jupyter_text>Plot weather TAVG over time to make sure it looks OK - are degrees right? Are seasons reasonable?<jupyter_code># todo: have one plot per station, add trendlines
weather_df[['DATE', 'TAVG', 'CDD', 'HDD']].set_index('DATE').sort_index().plot(figsize=(16,8))
# from https://pandas.pydata.org/pandas-docs/dev/getting_started/intro_tutorials/09_timeseries.html
# which has other good tips

import matplotlib.pyplot as plt
fig, axs = plt.subplots(figsize=(12, 4))
weather_df.groupby(weather_df["DATE"].dt.month)["TAVG"].mean().plot(
    kind='bar', rot=0, ax=axs
)
plt.xlabel("Month");  # custom x label using matplotlib
plt.ylabel("Mean Temperature");

weather_df.head()
weather_df.dtypes<jupyter_output><empty_output><jupyter_text>Save out the weather dataframe<jupyter_code># save for later
weather_df.to_csv('dataframes/weather_df.csv.gz', compression="infer", index=False)
#weather_df.to_pickle('dataframes/weather_df.pickle.gz', compression="infer")

# Trying Feather
# weather_df.to_feather('dataframes/weather_df.feather.zstd', compression='zstd')

# Trying JSON
# weather_df.to_json('dataframes/weather_df.json.gz', compression='infer')<jupyter_output><empty_output><jupyter_text>*  * * *  * *  Get list of holidays on which people might demand a different quantity of electric powerusing https://github.com/dr-prodigy/python-holidays installed via "conda install holidays"<jupyter_code>from datetime import date

import holidays
from datetime import datetime

# ignore the days people still go to work:  MLK, Wash, Columb, Veterans; load is probably unchanged these days
# alternatively: one-hot encode these days
us_holidays = holidays.UnitedStates()

holiday_list = []

for date, name in sorted(holidays.US(state='OK', years=[int(i) for i in years_to_load]).items()):
    if name not in ("<NAME>", 
                    "Washington's Birthday", 
                    "Veterans Day (Observed)",
                    "Veterans Day",
                    "Columbus Day",
                   ):
          holiday_list.append([date, name])

# many people don't work Christmas Eve; try adding that?
holiday_list.append([datetime.strptime('2015-12-24', '%Y-%m-%d').date(), 'Christmas Eve'])
holiday_list.append([datetime.strptime('2016-12-24', '%Y-%m-%d').date(), 'Christmas Eve'])
holiday_list.append([datetime.strptime('2017-12-24', '%Y-%m-%d').date(), 'Christmas Eve'])
holiday_list.append([datetime.strptime('2018-12-24', '%Y-%m-%d').date(), 'Christmas Eve'])
holiday_list.append([datetime.strptime('2019-12-24', '%Y-%m-%d').date(), 'Christmas Eve'])
            
holiday_df = pd.DataFrame(holiday_list, columns = ['date', 'holiday']).sort_values('date')
holiday_df['date'] = pd.to_datetime(holiday_df.date)
holiday_df.set_index('date', inplace=True)

holiday_df.head()
import re
# clip the word "Observed"
holiday_df["holiday"].replace(' \(Observed\)','',regex=True, inplace = True)

holiday_df.head()<jupyter_output><empty_output><jupyter_text>Save out the holiday dataframe<jupyter_code># save the holiday df
holiday_df.to_csv('dataframes/holiday_df.csv.gz', compression="infer")
#holiday_df.to_pickle('dataframes/holiday_df.pickle.gz', compression="infer")

# Trying Feather
#holiday_df.to_feather('dataframes/holiday_df.feather.zstd', compression='zstd')
# Trying JSON
# holiday_df.to_json('dataframes/holiday_df.json.gz', compression='infer')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11704490584258812, 'fraction_numerical': 0.028343795267986478, 'mean_word_length': 4.562986838571045, 'pattern_counts': {'":': 1, '<': 87, '<?xml version=': 0, '>': 88, 'https://': 10, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 63, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '403479', 'license': 'MIT', 'n_tokens_mistral': 7476, 'n_tokens_neox': 6791, 'n_words': 1909, 'path': '01_get_data.ipynb', 'repo_name': 'k5dru/CPSC_5383', 'size': '20730'}