Acea Group: Predicting Water Level in The Petrignano Aquifer

path = r'C:\Users\matte\OneDrive\Desktop\GitHub\data\Petrignano\Aquifer_Petrignano.csv'
df = pd.read_csv(path)
df.head()

df.info()

# Parse object type to datetime type
df['date'] = pd.to_datetime(df['date'], format = '%d/%m/%Y')
df.head()

df.isna().sum()

fig, ax = plt.subplots(nrows=5, ncols=1, figsize=(15, 25))

for i, column in enumerate(df.drop('date', axis=1).columns):
    # For now I handle missing value with fillna
    sns.lineplot(x=df['date'], y=df[column].fillna(method='ffill'), ax=ax[i], color='steelblue')
    ax[i].set_title(column, fontsize=17)
    ax[i].set_ylabel(column, fontsize=17)
                      
    ax[i].set_xlim([date(2009, 1, 1), date(2020, 6, 30)])     

fig, ax = plt.subplots(figsize=(12, 5))

sns.lineplot(x=df['date'], y=df['river_hydrometry'], color='steelblue',
             label='original', ax=ax)
sns.lineplot(x=df['date'], y=df['hydro_fill'], color='orange',
             label='interpolate', ax=ax)
ax.set_xlim([date(2012, 5, 1), date(2016, 10, 1)])
plt.show()

# interpolate in this case is quite good, so I'll use it for each of the columns

df.loc[df['depth_to_groundwater']==0, 'depth_to_groundwater'] = np.nan
df.loc[df['river_hydrometry']==0, 'river_hydrometry'] = np.nan

df['river_hydrometry'] = df['river_hydrometry'].interpolate()
df['drainage_volume'] = df['drainage_volume'].interpolate()
df['depth_to_groundwater'] = df['depth_to_groundwater'].interpolate()
df.drop(columns=['hydro_fill', 'river_hydrometry_back_up'], inplace=True)

df.head()

fig, ax = plt.subplots(ncols=1, nrows=3, sharex=True, figsize=(16,12))

sns.lineplot(x=df['date'], y=df['temperature'], color='steelblue', ax=ax[0])
ax[0].set_title('Daily Temperature', fontsize=14)

weekly_temp = df[['date','temperature']].resample('7D', on='date').mean().reset_index(drop=False)
sns.lineplot(x=weekly_temp['date'], y=weekly_temp['temperature'], color='steelblue', ax=ax[1])
ax[1].set_title('Weekly Temperature', fontsize=14)

monthly_temp = df[['date','temperature']].resample('1M', on='date').mean().reset_index(drop=False)
sns.lineplot(x=monthly_temp['date'], y=monthly_temp['temperature'], color='steelblue', ax=ax[2])
ax[2].set_title('Monthly Temperature', fontsize=14)

for i in range(3):
    ax[i].set_xlim([date(2009, 1, 1), date(2020, 6, 30)])
plt.show()

df['year'] = df['date'].dt.year
df['month'] = df['date'].dt.month
df['day'] = df['date'].dt.day
df['day_of_year'] = df['date'].dt.dayofyear
df['quarter'] = df['date'].dt.quarter
df['season'] = df['month'] % 12 // 3 + 1

df[['date', 'year', 'month', 'day', 'day_of_year', 'quarter','season']].head()

df.head()

fig, ax = plt.subplots(ncols=1, nrows=4, sharex=True, figsize=(16,8))
col = 'depth_to_groundwater'

res = seasonal_decompose(df[col], period=52, model='additive',
                         extrapolate_trend = 'freq')

ax[0].set_title(f'{col} Observed')
res.observed.plot(ax= ax[0], legend=False, color='steelblue')

ax[1].set_title(f'{col} Trend')
res.trend.plot(ax=ax[1], legend=False, color='steelblue')

ax[2].set_title(f'{col} Seasonal')
res.seasonal.plot(ax=ax[2], legend=False, color='steelblue')

ax[3].set_title(f'{col} Residual')
res.resid.plot(ax=ax[3], legend=False, color='steelblue')

plt.show()

f, ax = plt.subplots(nrows=1, ncols=2, figsize=(16, 8))

corrmat = df[to_decompose].corr()

sns.heatmap(corrmat, annot=True, vmin=-1, vmax=1, cmap='coolwarm_r', ax=ax[0])
ax[0].set_title('Correlation Matrix of Core Features', fontsize=16)

shifted_cols = [
    'depth_to_groundwater_seasonal',         
    'temperature_seasonal_shift_b_2m',
    'drainage_volume_seasonal_shift_2m', 
    'river_hydrometry_seasonal_shift_3m'
]
corrmat = df[shifted_cols].corr()

sns.heatmap(corrmat, annot=True, vmin=-1, vmax=1, cmap='coolwarm_r', ax=ax[1])
ax[1].set_title('Correlation Matrix of Lagged Features', fontsize=16)


plt.tight_layout()
plt.show()

df[to_decompose].corr().to_csv(r'C:\Users\matte\OneDrive\Desktop\GitHub\data\Petrignano\4 tableau\corr_matrix_normal.csv')
df[to_decompose].corr()

df[shifted_cols].corr().to_csv(r'C:\Users\matte\OneDrive\Desktop\GitHub\data\Petrignano\4 tableau\corr_matrix_shift.csv')
df[shifted_cols].corr()

train_size = int(0.85 * len(df))
test_size = len(df) - train_size

univariate_df = df[['date', 'depth_to_groundwater']].copy()
univariate_df.columns = ['ds', 'y']

train = univariate_df.iloc[:train_size, :]

x_train, y_train = pd.DataFrame(univariate_df.iloc[:train_size, 0]), pd.DataFrame(univariate_df.iloc[:train_size, 1])
x_valid, y_valid = pd.DataFrame(univariate_df.iloc[train_size:, 0]), pd.DataFrame(univariate_df.iloc[train_size:, 1])

print(len(train), len(x_valid))

from prophet import Prophet

# Train the model
model = Prophet()
model.fit(train)

y_hat = model.predict(x_valid)

# Calcuate metrics
score_mae = mean_absolute_error(y_valid, y_hat.tail(test_size)['yhat'])
score_rmse = np.sqrt(mean_squared_error(y_valid, y_hat.tail(test_size)['yhat']))

print(f'RMSE: {score_rmse}')

ml_predict = pd.merge(df[['date', 'depth_to_groundwater']], y_hat[['ds', 'yhat']],
         how='outer', left_on='date', right_on='ds')

ml_predict.drop(columns='ds', inplace=True)
ml_predict.tail()

fig, ax = plt.subplots(figsize=(12,6))

model.plot(y_hat, ax=ax)
sns.lineplot(x=x_valid['ds'], y=y_valid['y'], ax=ax, color='orange', label='Ground truth')

ax.set_title(f'Prediction \n MAE: {score_mae:.2f}, RMSE: {score_rmse:.2f}', fontsize=14)
ax.set_xlabel(xlabel='Date', fontsize=14)
ax.set_ylabel(ylabel='Depth to Groundwater', fontsize=14)

plt.show()