<jupyter_start><jupyter_text>Interpolating profiler data to a regular grid---<jupyter_code># Imports
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np

from scipy.interpolate import griddata
sns.set_style('whitegrid')
profiler_data = pd.read_csv('../../coastal_upwelling_output/profiler_data_2018.csv')
profiler_data.time
profiler_data['time'] = pd.to_datetime(profiler_data['time'])
profiler_hourly = profiler_data.resample('H', on='time').mean().dropna(how='all').reset_index()
# Plot hourly-averaged METBK data
fig, ax = plt.subplots(1, 1, sharex=True, figsize=(12,8))
ax.plot(profiler_hourly['time'], profiler_hourly['seawater_pressure'], linewidth=2)
ax.set_xlabel('Time', size=24)
ax.set_ylabel('Seawater Pressure', size=24)
ax.set_title('Oregon Offshore profiler hourly averages', size=24)
ax.invert_yaxis()<jupyter_output><empty_output><jupyter_text>Since there are two significant gaps in the data, the results of the interpolation will include a large ratio of "made up" data points. It's recommended to add the original gaps back to the data after interpolation to avoid feeding models large chunks of fabricated data, according to: https://www.mathworks.com/matlabcentral/answers/363643-use-griddata-ignoring-gaps-in-data. See the [misc_code](https://github.com/dgumustel/coastal-upwelling/blob/main/extras/misc_code.ipynb) notebook for the start and end dates of some of the gaps in the profiler data.   Using `griddata()` To start, interpolate a single day's worth of data. I arbitrarily chose December 1.<jupyter_code>df = profiler_data

df_resampled = df.resample('T', on='time').mean().dropna(how='all').reset_index()
df_resampled['time']
mask = (df_resampled['time'][:] > '2018-12-01') & (df_resampled['time'] <= '2018-12-02')
df_resampled = df_resampled.loc[mask].reset_index(drop=True)<jupyter_output><empty_output><jupyter_text>The only way I could get `griddata` to work was by converting the times to integers.<jupyter_code>df_resampled['time_'] = df_resampled['time'].values.astype(object)
df_resampled['time_']<jupyter_output><empty_output><jupyter_text>However, theme numbers are enormous and leads to computation overflow errors, so I divide them into smaller integers and can reverse the division and type change later to get the human-readable times back for creating the final dataframe.<jupyter_code># data coordinates and values
x = df_resampled.time_.astype(np.int64) / 100000000000 
y = df_resampled.seawater_pressure 
z = df_resampled.seawater_temperature 
# target grid to interpolate to
xi = np.arange(min(x), max(x), 60) # intervals of 60 seconds
yi = np.arange(min(y), max(y), 1) # intervals of 1 meter
xi.shape
yi.shape
xi,yi = np.meshgrid(xi,yi)
# interpolate
zi_linear = griddata((x,y),z,(xi,yi),method='linear', rescale=True)
zi_linear.shape<jupyter_output><empty_output><jupyter_text>Plot the original data as scatter points to demonstrate exactly what interpolation does.<jupyter_code>fig, ax = plt.subplots(figsize=(12,6))

ax.invert_yaxis()
ax.grid()

sc = ax.scatter(x,y,c=z,edgecolors='k')
ax.set_xlabel('Time', size=16)
ax.set_ylabel('Pressure (dbar)', size=16)
ax.set_title('Temperature (degC)', size=16)
cb = fig.colorbar(sc,ax=ax)
cb.set_label(label='Temperature (deg C)', size=16)
cb.ax.tick_params(labelsize=16)

plt.tight_layout()
plt.savefig('../figures/dec-01_raw.png')<jupyter_output><empty_output><jupyter_text>These are the measurements collected by the shallow profiler resampled to 1-minute averages. Next, we can see how using the linear interpolation method fills in the gaps between measurements.<jupyter_code>fig, ax = plt.subplots(figsize=(12,6))

cs2 = ax.contourf(xi,yi,zi_linear)
plt.colorbar(cs2, ax=ax)
ax.scatter(x, y, c=z,edgecolors='k')
ax.invert_yaxis()
# ax.set_title("Using method='linear'",fontsize=16)
ax.set_xlabel('Time', size=16)
ax.set_ylabel('Pressure (dbar)', size=16)
ax.set_title('Temperature (degC)', size=16)
# cb = fig.colorbar(sc,ax=ax)
# cb = plt.colorbar(cs2, ax=ax)
cb.set_label(label='Temperature (deg C)', size=16)
cb.ax.tick_params(labelsize=16)

plt.tight_layout()
plt.savefig('../figures/dec-01_griddata.png')<jupyter_output><empty_output><jupyter_text>See how most of the white spaces are filled in? This produces a grid of data that is aligned along pressure and time, so we can turn this into a pandas dataframe with time as the index and different pressure as the columns. Then we can use each pressure column as individual features in models, and even do PCA to reduce dimensionality.  Interpolating 1 monthNow do the same for a full month of data to see how it performs; I chose the month of December for this one.<jupyter_code>df_resampled = df
df_resampled['time']
mask = (df_resampled['time'][:] > '2018-12-01') & (df_resampled['time'] <= '2018-12-31')
df_resampled = df_resampled.loc[mask].reset_index(drop=True)
df_resampled['time_'] = df_resampled['time'].values.astype(object)
df_resampled['time_']
# data coordinates and values
x = df_resampled.time_.astype(np.int64) / 100000000000
y = df_resampled.seawater_pressure 
z = df_resampled.seawater_temperature 
# target grid to interpolate to
xi = np.arange(min(x), max(x), 60)
yi = np.arange(min(y), max(y), 1)
xi,yi = np.meshgrid(xi,yi)
xi.shape
yi.shape
# interpolate
zi_linear = griddata((x,y),z,(xi,yi),method='linear', rescale=True)
fig, ax = plt.subplots(figsize=(16,8))

cf = ax.contourf(xi,yi,zi_linear)
plt.colorbar(cf, ax=ax)
# ax2.plot(x,y,'k.')
ax.set_xlabel('time')
ax.set_ylabel('depth')
ax.invert_yaxis()
ax.set_title("Using method='linear'",fontsize=16)

plt.savefig('../figures/december_gridded.png')<jupyter_output><empty_output><jupyter_text>Interpolating the full datasetNow do it for all of the 2018 data.<jupyter_code>df_resampled = df
df_resampled['time_'] = df_resampled['time'].values.astype(object)
df_resampled['time_']
# data coordinates and values
x = df_resampled.time_.astype(np.int64) / 100000000000 
y = df_resampled.seawater_pressure 
z = df_resampled.seawater_temperature 
# target grid to interpolate to
xi = np.arange(min(x), max(x), 60)
yi = np.arange(min(y), max(y), 1)
xi,yi = np.meshgrid(xi,yi)
xi.shape
yi.shape
# interpolate
zi_linear = griddata((x,y),z,(xi,yi),method='linear', rescale=True)
fig, ax = plt.subplots(figsize=(24,8))

cf = ax.contourf(xi,yi,zi_linear, levels=10)
cb = plt.colorbar(cf, ax=ax)
# ax.plot(x,y,'k.') # uncomment to plot locations of real measurements over interpolated values

ax.set_xlabel('Time',fontsize=16)
ax.set_ylabel('Pressure (dbar)',fontsize=16)
ax.invert_yaxis()
ax.set_title("2018 Interpolated Profiler Data",fontsize=24)
plt.tick_params(axis='both', which='major', labelsize=16)

cb.set_label(label='Temperature (deg C)', size=16)
cb.ax.tick_params(labelsize=16)

plt.savefig('../figures/2018_gridded.png')<jupyter_output><empty_output><jupyter_text>You can change the number of levels shown in the plot above if you'd like to see the water column structure in more detail. Also note two sections of data (around 1.536 and 1.538 on the x-axis) - data was missing here prior to interpolation, so `griddata()` filled these areas in according to the measurements on either side of the gap. These are data points we might consider removing prior to modeling.  --- Create a dataframe from the griddata() results<jupyter_code>z_df = pd.DataFrame(zi_linear, columns=np.arange(min(x), max(x), 60), index=np.arange(min(y), max(y), 1))
z_df = z_df.T
z_df<jupyter_output><empty_output><jupyter_text>A few things to fix here before saving:* Round the pressure columns to fewer decimals* Set time as the index* Convert the index back to human-readable times<jupyter_code># Round the pressure columns to fewer decimals
col_dict = {x:round(x, 2) for x in z_df.columns}
z_df.rename(col_dict, inplace=True, axis=1)<jupyter_output><empty_output><jupyter_text>There are also plenty of nans to clean up. These are mostly located around the edges of the dataset:<jupyter_code>import missingno as msno
msno.matrix(z_df);<jupyter_output><empty_output><jupyter_text>This would look more like a checkerboard if there were whole rows of data missing. Since all the nans are around the edges, we can drop these edge rows and columns to get a clean grid.<jupyter_code>print(len(z_df))
print(*z_df.isna().sum())
f'Drop all columns with more than {len(z_df) * 0.05} nans, then consider dropping the rest as rows.'<jupyter_output><empty_output><jupyter_text>The cell to drop columns below is commented out to prevent running it multiple times. Uncomment it, run it once, and then recomment it to prevent accidental loss of data. However, make sure you don't skip it, or the next cell will completely wipe your dataframe by removing every row, because each row contains nans in the first and last columns.<jupyter_code># z_df.drop(columns=z_df.columns[[0, -1, -2, -3, -4]], inplace=True)
print(*z_df.isna().sum())<jupyter_output>59 21 5 3 3 3 3 3 3 3 3 3 3 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 10 79<jupyter_text>Now the rest of the nans can be cleaned up by dropping just three rows<jupyter_code>z_df.dropna(axis=0, inplace=True)
z_df
1 - len(z_df) / 2403<jupyter_output><empty_output><jupyter_text>Only lost 3.28% of the observations and 5 meters of depth. Now get a datetime column back.<jupyter_code># Get a column of the timestamps to convert them back to human-readable times
z_df['time'] = z_df.index
# Do conversions, then reset the index to be time
z_df['time'] = 100000000000 * z_df['time']
z_df['time'] = pd.to_datetime(z_df['time'])
z_df<jupyter_output><empty_output><jupyter_text>--- Append CUTI data<jupyter_code>CUTI_data = pd.read_csv('../../coastal_upwelling_output/CUTI_daily.csv',
                  parse_dates=[[0,1,2]],
                  infer_datetime_format=True)
CUTI_data.rename(columns={'year_month_day':'time'},inplace=True)
CUTI_data = CUTI_data[CUTI_data['time'].dt.year == 2018][['time', '44N']].reset_index(drop=True)
CUTI_data<jupyter_output><empty_output><jupyter_text>First, append the CUTI data according to time. I couldn't get this to work using z_df.index, so I added a time column back only to delete it again after.<jupyter_code>for month in range(1, 13):
    for day in range(1, 32):
        try:
            indices = z_df[(z_df['time'].dt.month == month) & (z_df['time'].dt.day == day)].index
            cuti_value = CUTI_data.loc[(CUTI_data['time'].dt.month == month) & (CUTI_data['time'].dt.day == day)]['44N'].values[0]
#             print(month, day, cuti_value)
            z_df.loc[indices,'CUTI'] = cuti_value
        except:
            pass<jupyter_output><empty_output><jupyter_text>We can finally turn time back into the index.<jupyter_code>z_df.set_index('time', inplace=True)<jupyter_output><empty_output><jupyter_text>Then create a binary target variable for classification modeling.<jupyter_code>z_df['upwelling'] = z_df['CUTI'].apply(lambda x: 1 if x > 0 else 0)
z_df
z_df['upwelling'].value_counts(normalize=True)<jupyter_output><empty_output><jupyter_text>Save it for modeling!<jupyter_code>z_df.to_csv('../../coastal_upwelling_output/interpolated.csv')<jupyter_output><empty_output><jupyter_text>--- Other griddata optionsIf you're curious to see how `griddata`'s other methods interpolate this data, run the cells below. Warning, they take a few minutes!<jupyter_code>df_resampled = df
df_resampled['time_'] = df_resampled['time'].values.astype(object)
# data coordinates and values
x = df_resampled.time_.astype(np.int64) / 100000000000 # np.random.random(100)
y = df_resampled.seawater_pressure # np.random.random(100)
z = df_resampled.seawater_temperature # np.random.random(100)
# target grid to interpolate to
xi = np.arange(min(x), max(x), 60)
yi = np.arange(min(y), max(y), 1)
xi,yi = np.meshgrid(xi,yi)
# interpolate
zi_nearest = griddata((x,y),z,(xi,yi),method='nearest', rescale=True)
zi_linear = griddata((x,y),z,(xi,yi),method='linear', rescale=True)
zi_cubic = griddata((x,y),z,(xi,yi),method='cubic', rescale=True)
# plot
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(16,16))

cs1 = ax1.contourf(xi,yi,zi_nearest)
plt.colorbar(cs1, ax=ax1)
ax1.set_xlabel('time')
ax1.set_ylabel('depth')
ax1.invert_yaxis()
ax1.set_title("Using method='nearest'",fontsize=16)

cs2 = ax2.contourf(xi,yi,zi_linear)
plt.colorbar(cs2, ax=ax2)
ax2.set_xlabel('time')
ax2.set_ylabel('depth')
ax2.invert_yaxis()
ax2.set_title("Using method='linear'",fontsize=16)

cs3 = ax3.contourf(xi,yi,zi_cubic)
plt.colorbar(cs3, ax=ax3)
ax3.set_xlabel('time')
ax3.set_ylabel('depth')
ax3.invert_yaxis()
ax3.set_title("Using method='cubic'",fontsize=16)

plt.tight_layout();
plt.savefig('../figures/gridded_data.png')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12972389325929354, 'fraction_numerical': 0.039025142680857626, 'mean_word_length': 5.792561550550026, 'pattern_counts': {'":': 0, '<': 86, '<?xml version=': 0, '>': 87, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 5, 'word_list_counts': {'cursed_substrings.json': 31, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '420379', 'license': 'MIT', 'n_tokens_mistral': 5201, 'n_tokens_neox': 4737, 'n_words': 1303, 'path': 'notebooks/.ipynb_checkpoints/06_0_interpolation_griddata-checkpoint.ipynb', 'repo_name': 'dgumustel/coastal-upwelling', 'size': '12966'}