<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Cleaning resampling ISH temperature datasets
Step1:  Constants / Parameters
Step2:  Interpolation and cleanup
Step3:  Deduplication
More precisely, we can only have one value for each time,
otherwise interpolation doesn't make much sense (or work)
Step4:  Outlier removal
Using a deviation from running median/sigam threshold method
Step5:  "by-hand" fixes for particular datasets, hopefully minimal
Step6:  Identify bigger gaps which will get filled day-over-day interpolation
Interpolate based on same hour-of-day across days.
Step7:  Interpolate to produce on-the-hour values
Simple interpolation hour-to-hour
Step8:  Fill the bigger gaps
Step9:  Check that it looks OK...
Plot the temperature data
Step10:  Save final cleaned temperatures
Step11:  Compute the normals
Need the normal (repated so it covers 2 years) for running medfoes on the normals
Not needed for this particular study
<ASSISTANT_TASK:>
Python Code:
# boilerplate includes
import sys
import os
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
#from mpl_toolkits.mplot3d import Axes3D
import pandas as pd
import seaborn as sns
import datetime
import scipy.interpolate
# import re
from IPython.display import display, HTML
%matplotlib notebook
plt.style.use('seaborn-notebook')
pd.set_option('display.max_columns', None)
Explanation: Cleaning resampling ISH temperature datasets
End of explanation
# PARAMETERS (might be overridden by a calling script)
# if not calling from another script (batch), SUBNOTEBOOK_FLAG might not be defined
try:
    SUBNOTEBOOK_FLAG
except NameError:
    SUBNOTEBOOK_FLAG = False
    
# Not calling as a sub-script? define params here
if not SUBNOTEBOOK_FLAG:
    
    # SET PARAMETER VARIABLES HERE UNLESS CALLING USING %run FROM ANOTHER NOTEBOOK
    
    DATADIR = '../data/temperatures/ISD'
    OUTDIR = '../data/temperatures'
    FTPHOST = 'ftp.ncdc.noaa.gov'
    FETCH_STATIONS_LIST_FILE = True
    TEMP_COL = 'AT' # The label of the hourly temperature column we make/output
    # Resampling and interpolation parameters
    # spline order used for converting to on-the-hour and filling small gaps
    BASE_INTERPOLATION_K = 1 # 1 for linear interpolation
    # give special treatment to data gaps longer than...
    POTENTIALLY_PROBLEMATIC_GAP_SIZE = pd.Timedelta('03:00:00')
    # Time range to use for computing normals (30 year, just like NOAA uses)
    NORM_IN_START_DATE = '1986-07-01'
    NORM_IN_END_DATE = '2016-07-01'
    # Time range or normals to output to use when running 'medfoes on normal temperature' (2 years, avoiding leapyears)
    NORM_OUT_START_DATE = '2014-01-01'
    NORM_OUT_END_DATE = '2015-12-31 23:59:59'
    
print("Cleaning temperature data for ",STATION_CALLSIGN)
# Potentially turn interactive figure display off
if SUPPRESS_FIGURE_DISPLAY:
    plt.ioff()
Explanation: Constants / Parameters
End of explanation
# Load the data
fn = "{}_AT.h5".format(STATION_CALLSIGN)
ot = pd.read_hdf(os.path.join(DATADIR,fn), 'table')
Explanation: Interpolation and cleanup
End of explanation
t = ot.copy(deep=True) # not needed, just safety
# just showing the duplicates
tmp = t[t.index.duplicated(keep=False)].sort_index()
print(len(tmp), 'duplicates')
#display(tmp) # decomment to see the list of duplicates
# actually remove duplicates, just keeping the first
# @TCC could somehow try to identify the most reliable or take mean or such
t = t[~t.index.duplicated(keep='first')].sort_index()
Explanation: Deduplication
More precisely, we can only have one value for each time,
otherwise interpolation doesn't make much sense (or work)
End of explanation
# fairly permissive settings
rolling_sigma_window = 24*5 # None or 0 to just use median instead of median/sigma
rolling_median_window = 5
thresh = 1.5 # deviation from media/sigma to trigger removal
multipass = True # cycle until no points removed, or False for not
tin = t
cum_num = 0
while multipass:
    if rolling_sigma_window:
        sigma = t['AT'].rolling(window=rolling_sigma_window, center=True).std()
    else:
        sigma = 1
    diff = (t['AT']-t['AT'].rolling(window=rolling_median_window, center=True).median())/sigma
    outlier_mask = diff.abs() > thresh
    num = np.count_nonzero(outlier_mask)
    cum_num += num
    print("removing {} points".format(num))
    if num == 0:
        break
    # plotting each step
#     ax = t.plot(linestyle='-', marker='*')
#     if np.count_nonzero(outlier_mask) > 0:
#         t[outlier_mask].plot(ax=ax, linestyle='none', marker='o', color='red')
#     diff.abs().plot(ax=ax)
#     if np.count_nonzero(outlier_mask) > 0:
#         diff.abs()[outlier_mask].plot(ax=ax, linestyle='none', marker='o', color='yellow')
    t = t[~outlier_mask]
# plot showing what is being removed
if cum_num > 0:
    fig = plt.figure()
    ax = fig.add_subplot(1,1,1)
    ax = tin[~tin.index.isin(t.index)].plot(ax=ax, linestyle='none', marker='o', color='r', zorder=8)
    ax = tin.plot(ax=ax, linestyle='-', linewidth=1, marker=None, color='red')
    ax = t.plot(ax=ax, linestyle='-', marker='.', color='blue')
    ax.set_ylabel('air temperature [$\degree$ C]')
    ax.legend(['outlier', 'original', 'cleaned'])
    ax.set_title(STATION_CALLSIGN)
# saving figure
# saving
fn = '{}_outlier.png'.format(STATION_CALLSIGN)
fig.savefig(os.path.join(OUTDIR,fn))
#mpld3.save_html(fig, '{}_outler.html'.format(STATION_CALLSIGN))
# Actually apply the outlier removal
ot = t
Explanation: Outlier removal
Using a deviation from running median/sigam threshold method
End of explanation
def remove_spurious_temps(ot, query_op, date1, date2=None, plot=True, inplace=False):
    if date2 is None:
        date2 = date1
    ax = ot.loc[date1:date2].plot(ax=None, linestyle='-', marker='o') # plot
    out_t = ot.drop(ot.loc[date1:date2].query('AT {}'.format(query_op)).index, inplace=inplace)
    if inplace:
        out_t = ot
    out_t.loc[date1:date2].plot(ax=ax, linestyle='-', marker='*') # plot'
    ax.set_title("Remove AT {}, range=[{}:{}]".format(query_op, date1, date2))
    return out_t
STATION_CALLSIGN
if STATION_CALLSIGN == 'KSNA': # KSNA (Orange County)
    # 2016-08-14 to 2016-08-15 overnight has some >0 values when they should be more like 19-20
    remove_spurious_temps(ot, '< 0', '2016-08-14', '2016-08-15', inplace=True)
    
if STATION_CALLSIGN == 'KSFO':
    remove_spurious_temps(ot, '< 0', '1976-07-16', '1976-07-17', inplace=True)
if STATION_CALLSIGN == 'KRIV':
    remove_spurious_temps(ot, '< 0', '1995-11-15', '1995-11-15', inplace=True)
Explanation: "by-hand" fixes for particular datasets, hopefully minimal
End of explanation
# flag the gaps in the original data that are possibly too long for the simple interpolation we did above
gaps_filename = os.path.join(OUTDIR, "{}_AT_gaps.tsv".format(STATION_CALLSIGN))
gaps = ot.index.to_series().diff()[1:]
idx = np.flatnonzero(gaps > POTENTIALLY_PROBLEMATIC_GAP_SIZE)
prob_gaps = gaps[idx]
# save to file for future reference
with open(gaps_filename,'w') as fh:
    # output the gaps, biggest to smallest, to review
    print('#', STATION_CALLSIGN, ot.index[0].isoformat(), ot.index[-1].isoformat(), sep='\t', file=fh)
    print('# Potentially problematic gaps:', len(prob_gaps), file=fh)
    tmp = prob_gaps.sort_values(ascending=False)
    for i in range(len(tmp)):
        rng = [tmp.index[i]-tmp.iloc[i], tmp.index[i]]
        print(rng[0], rng[1], rng[1]-rng[0], sep='\t', file=fh)
if not SUPPRESS_FIGURE_DISPLAY:
    # go ahead and just print it here too
    with open(gaps_filename) as fh:
        for l in fh:
            print(l, end='')
else:
    print('# Potentially problematic gaps:', len(prob_gaps))
Explanation: Identify bigger gaps which will get filled day-over-day interpolation
Interpolate based on same hour-of-day across days.
End of explanation
# Interpolate to get on-the-hour values
newidx = pd.date_range(start=ot.index[0].round('d')+pd.Timedelta('0h'),
                       end=ot.index[-1].round('d')-pd.Timedelta('1s'),
                       freq='1h', tz='UTC')
if True:
    # Simple linear interpolation
    at_interp_func = scipy.interpolate.interp1d(ot.index.astype('int64').values, 
                                           ot['AT'].values, 
                                           kind='linear', 
                                           fill_value=np.nan, #(0,1) 
                                           bounds_error=False)
else:
    # Should be better method, but has some screwy thing using updated data
    at_interp_func = scipy.interpolate.InterpolatedUnivariateSpline(
        ot.index.astype('int64').values, 
        ot['AT'].values, 
        k=BASE_INTERPOLATION_K,
        ext='const')
nt = pd.DataFrame({'AT':at_interp_func(newidx.astype('int64').values)},
                   index=newidx)
Explanation: Interpolate to produce on-the-hour values
Simple interpolation hour-to-hour
End of explanation
# Fill those gaps using day-to-day (at same hour) interpolation
gap_pad = pd.Timedelta('-10m') # contract the gaps a bit so we don't remove good/decent edge values
t = nt.copy(deep=True) # operate on a copy so we can compare with nt
# fill the gap ranges with nan (replacing the default interpolation)
for i in range(len(prob_gaps)):
    rng = [prob_gaps.index[i]-prob_gaps.iloc[i], prob_gaps.index[i]]
    t[rng[0]-gap_pad:rng[1]+gap_pad] = np.nan
# reshape so each row is a whole day's (24) data points
rows = int(t.shape[0]/24)
foo = pd.DataFrame(t.iloc[:rows*24].values.reshape((rows,24)))
# simple linear interpolation
foo.interpolate(metnod='time', limit=24*60, limit_direction='both', inplace=True)
# # Alternative interpolation using running means
# # @TCC not great for very large gaps
# RUNNING_MEAN_WINDOW_SIZE = 3
# while True:
#     # interpolate each column (temp at hour x on each day)
#     # filling nans with values from a windowed running mean
#     foo.fillna(foo.rolling(window=RUNNING_MEAN_WINDOW_SIZE, min_periods=1, center=True).mean(), inplace=True)
#     if not foo.isnull().values.any():
#         break
# reshape back
t = pd.DataFrame({'AT':foo.stack(dropna=False).values}, index=t.index[:rows*24])
Explanation: Fill the bigger gaps
End of explanation
# You can specify a specific range by setting r1 and r2, or None for full range
#r1, r2 = '1952-05-07', '1952-05-23'
r1, r2 = None, None
if r1 is None:
    r1 = t.index[0]
if r2 is None:
    r2 = t.index[-1]
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(ot.loc[r1:r2].index, ot.loc[r1:r2]['AT'], linestyle='none', marker='.', label='raw')
#ax.scatter(ot.loc[r1:r2].index, ot.loc[r1:r2]['AT'], marker='.', label='raw')
ax.plot(nt.loc[r1:r2].index, nt.loc[r1:r2]['AT'], linestyle='-', marker=None, lw=1, label='interpolated')
# ax.plot(t.loc[r1:r2].index, t.loc[r1:r2]['AT'], '-*', lw=1, label='filled')
# @TCC maybe make a single dataframe with the parts I don't want deleted or masked out
for i in range(len(prob_gaps)):
    if i == 0: # only label first segment
        label = 'filled'
    else:
        label = ''
    rng = [tmp.index[i]-tmp.iloc[i], tmp.index[i]]
    ax.plot(t.loc[rng[0]:rng[1]].index, t.loc[rng[0]:rng[1]]['AT'], '.-', lw=1, color='r', label=label)
# # mark the big gaps with vertical lines
# for i in range(len(prob_gaps)):
#     ax.axvline(prob_gaps.index[i]-prob_gaps.iloc[i],
#                c='k', ls=':', lw=0.5)
#     ax.axvline(prob_gaps.index[i],
#                c='k', ls=':', lw=0.5)
    
ax.set_xlim((r1,r2))
ax.set_xlabel('DateTime')
ax.set_ylabel('Temperature [$\degree$C]')
ax.set_title(STATION_CALLSIGN)
ax.legend()
# saving
fig.savefig(os.path.join(OUTDIR, '{}_cleaning.png'.format(STATION_CALLSIGN)))
#mpld3.save_html(fig, '{}_cleaning.html'.format(STATION_CALLSIGN))
Explanation: Check that it looks OK...
Plot the temperature data
End of explanation
outfn = os.path.join(OUTDIR, "{}_AT_cleaned".format(STATION_CALLSIGN))
print("Saving cleaned temp data to:", outfn)
t.to_hdf(outfn+'.h5', 'table', mode='w',
          data_colums=True, complevel=5, complib='bzip2',
          dropna=False)
Explanation: Save final cleaned temperatures
End of explanation
# # Time range to use for computing normals (30 year, just like NOAA uses)
# NORM_IN_START_DATE = '1986-07-01'
# NORM_IN_END_DATE = '2016-07-01'
# # Time range or normals to output to use when running 'medfoes on normal temperature' (2 years, avoiding leapyears)
# NORM_OUT_START_DATE = '2014-01-01'
# NORM_OUT_END_DATE = '2015-12-31 23:59:59'
# %run "Temperature functions.ipynb" # for compute_year_over_year_norm function
# tempnorm = compute_year_over_year_norm(ot,
#                                        NORM_OUT_START_DATE, NORM_OUT_END_DATE,
#                                        NORM_IN_START_DATE, NORM_IN_END_DATE,
#                                        freq='hourly',
#                                        interp_method='linear',
#                                        norm_method='mean')
# # Save as csv for medfoes input
# outfn = os.path.join(OUTDIR, "{}_AT_cleaned_normalsX2.csv".format(STATION_CALLSIGN))
# print("Saving temp normals data to:",outfn)
# tempnorm.to_csv(outfn, index_label='datetime')
# tempnorm.plot()
# Turn iteractive display back on, if we turned it off
if SUPPRESS_FIGURE_DISPLAY:
    plt.ioff()
Explanation: Compute the normals
Need the normal (repated so it covers 2 years) for running medfoes on the normals
Not needed for this particular study
End of explanation
<END_TASK>