MateoFile.py

# -*- coding: utf-8 -*-
"""Copy of Groningen_FirstInformation_4Hojjat.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/19nt0wujgKGEVGl4BNxoyjHIk5CGbeidr

##
# Import libraries
Python packages to use. Run before running seperate sections
"""
# %%
# Commented out IPython magic to ensure Python compatibility.
# Loading and data transfer functions
# from google_drive_downloader import GoogleDriveDownloader as gdd
from utils import Functions_spatial
# Array and Dataformating
import numpy as np
import h5py
import pandas as pd
# Plotting functions
import matplotlib.pylab as plt
from matplotlib.patches import Wedge, Polygon
# %matplotlib inline
# %%
"""# Load the HDF5 file from Mateo Acosta's google drive"""

# gdd.download_file_from_google_drive(file_id='1revYMCdicTOScv_-OzA8J5DkArxh2mDO',dest_path='./Groningen_Data_Example.hdf5',unzip=True)
f = h5py.File('data/Groningen_Data_Example.hdf5', 'r')
# keys() shows the contents of the file
"""# Extract the data to use variables simpler"""

# extract all the data.
# note that these hdf4 datasets should behave exactly as np arrays

# The reservoir data
x_res = f['x_res']
y_res = f['y_res']
reservoir_outline = f['res_outline']
x_outline = f['x_outline']
y_outline = f['y_outline']
# The original well extraction data
extractions = f['extraction_data']
wells_locations = f['wells_locations_list']
wells_names = f['wells_names']
# The computed pressures from the linearized pressure diffusion model
pressures = f['pressures_list']
x_mesh = f['x_mesh']  # The meshes used for the pressure diffusion model
y_mesh = f['y_mesh']  # The meshes used for the pressure diffusion model
# The computed deformations and coulomb stress change from the BorisBlocks model
deformations = f['deformations_list']
# The meshes used for the BorisBlocks model
x_reservoir_mesh = f['x_reservoir_mesh']
# The meshes used for the BorisBlocks model
y_reservoir_mesh = f['y_reservoir_mesh']
max_coulomb_stress = f['max_coulomb_stress']
# the maximum coulomb stress smoothed with a 6km Gaussian2D kernel calculated 10m below the reservoir.
smoothed_coulomb_stress = f['smoothed_coulomb_stress']

# The simulation given in this example is from 1956 to 2019 for a total of 756 months
y0, y1, dy = 1956, 2019, 12

dates = np.linspace(y0, y1+1, (y1-y0)*dy+1)

"""# Import the seismicity catalog"""

# Downloading the calatog:
#gdd.download_file_from_google_drive(file_id='1kKJbecpDO7km0t8DHL7GZGmcVSrx_XaM', dest_path = ('./catalog2.txt'))
cat = pd.read_csv('data/catalog2.txt', sep='\t',
                  names=['Date', 'Time', 'Mag', 'Depth', 'RDX', 'RDY'])
# creating the dates and times
cat['DateTime'] = pd.to_datetime(
    cat.Date + ' ' + cat.Time, format='%d %m %Y %H %M %S')
temp_index = pd.DatetimeIndex(cat.DateTime)
cat['decDate'] = temp_index.year + temp_index.month/12 + (temp_index.day + (
    temp_index.hour + (temp_index.minute + (temp_index.second/60)/60)/24))/365.25
# filtering the catalog to magnitudes > mc
mc = 1.5
cat = cat[cat.Mag > mc]
cat = cat.reset_index()

"""# Plot a couple things"""

# A couple plots to get you started:

fig, ax = plt.subplots(1, 3, figsize=(18, 5))
# the final map of pressures
quad1 = ax[0].scatter(x_mesh, y_mesh, c=pressures[-1],
                      cmap='viridis', label='Pressure')
ax[0].scatter(wells_locations[:, 0], wells_locations[:, 1],
              c='k', marker='o', s=100, label='Wells')
ax[0].plot(x_outline[:], y_outline[:], color='r')
cb1 = plt.colorbar(quad1, ax=ax[0])
cb1.set_label('Final pressure (MPa)')
ax[1].set_title('Pressure')
# the final map of deformation
quad2 = ax[1].scatter(x_reservoir_mesh, y_reservoir_mesh, 50, deformations[-1] -
                      deformations[0], cmap='inferno', label='Deformations', alpha=1)
ax[1].plot(x_outline[:], y_outline[:], 'k--')
title = ax[1].set_title('Z Displacement (mm)')
cb1 = plt.colorbar(quad2, ax=ax[1])
cb1.set_label('Final vertical displacement (mm)')
# the final map of coulomb stress change
quad3 = ax[2].scatter(np.array(x_res).flatten(), np.array(
    y_res).flatten(), c=smoothed_coulomb_stress[:, :, -1], cmap='coolwarm')
ax[2].plot(x_outline[:], y_outline[:], 'k--')
ax[2].set_title('Maximum coulomb stress change')
cb1 = plt.colorbar(quad3, ax=ax[2])
cb1.set_label('Maximum Coulomb stress change (MPa)')
# the map labels for the three top subplots
for ii in [0, 1, 2]:
    ax[ii].set_xlabel('X (km)')
    ax[ii].set_ylabel('Y (km)')
    ax[ii].set_xlim([230, 270])
    ax[ii].set_ylim([565, 615])
    ax[ii].set_aspect('equal')
# %%
# the time evolution of average stress and pressure
fig, ax = plt.subplots(1, figsize=(18, 5))
ax.plot(dates, np.nanmean(pressures, axis=1), 'b', label='Pressure')
ax2 = ax.twinx()  # second y-axis for the coulomb stress
ax2.plot(dates, np.nanmean(smoothed_coulomb_stress,
         axis=(0, 1)), 'k', label='CSC')
ax.set_xlabel('Time')
ax.set_ylabel('Average reservoir pressure (MPa)')
ax2.set_ylabel('Maximum Coulomb Stress change averaged over reservoir (MPa)')
ax.legend(loc=9)
ax2.legend(loc=0)

# the time evolution of seismicity
fig, ax = plt.subplots(1, figsize=(18, 5))
ax.plot(cat.decDate, cat.index, 'r', label='Cumulative events in catalog')
ax.set_xlabel('Time')
ax.set_ylabel('Cumulative number of events')
# %%
# Functions: Function that discern whether a point is inside the border
# For this project we need an integrator in space
# To do so lets first see how does the data look like in position
# Based on the the Fig 1 (Right hand side)
# Based on this line:
# quad3 = ax[2].scatter(np.array(x_res).flatten(), np.array(y_res).flatten(), c=smoothed_coulomb_stress[:,:,-1], cmap='coolwarm')

# So here np.array(x_res).flatten() and np.array(x_res).flatten() are vector
# But the question is, should we limit ourselves to the outline? or not?
# So np.array(x_res) is an array 97*99 matrix the columns of this matrix are equal and each row is added by the previous plus .5
# so we have 97 number of X
# So np.array(y_res) is an array 97*99 matrix the rows of this matrix are equal and each column is added by the previous plus .5
# Actually, we notice that the data has already been discretized into 0.25km^2, so we do not need to do integration: we can simply sum with weights of (0.25km**2)
# Let's see what is x_out and y_out, you can see it by running the following line:
# We also need a function to check whether the point is inside the Groninger border

# plt.plot(np.array(x_outline),np.array(y_outline))
# Ok let's see what does the stress data looklike


# from matplotlib import cm
# fig, ax = plt.subplots(subplot_kw={"projection": "3d"})

# surf = ax.plot_surface(np.array(x_res), np.array(y_res), np.array(smoothed_coulomb_stress[:,:,700]), cmap=cm.coolwarm,
#                        linewidth=0, antialiased=False)
# Ok, what else do we need?
# We also need a function to calculate the R for the model
#
# In the model, we input the history of a specific point, and predict the R for that specific point.
# So inside that function we need a time integrator that sums integrate an exponential function
# So in the inner most part, we need a function that integrates in time.

# Trying to write the code in a way that has the most readability.
# I should also check the the seismicity rate is sink with stress time or not.
# we should also check the t_b and ds_t to see whether they are match or not.
# So the time for ds is dates
# Now lets check the integrater function
# %%
# Checking the Functions.TimeIntegrator and Functions.FindTimeIndex
# Asig=0.006
# t_a=8700
# Ds_c=0.17
# Ds_t=smoothed_coulomb_stress[45,45,:]
# t_b=1980
# t_end=1985
# i=Functions.FindTimeIndex(dates,t_b)
# j=Functions.FindTimeIndex(dates, t_end)
# #plt.plot(Ds_t)
# plt.plot(dates[i:j+1],np.exp((Ds_t[i:j+1] - Ds_c) / Asig))
# # plt.yscale("log")

# ans=Functions.TimeIntegrator(Ds_t,t_b,Ds_c,dates,Asig,t_end)
# # From the Geometry the integral should be around .35
# %%
# So, now we need a function that uses the integral to calculate the R
# in the formula

# %%

fac = 1
DX, DY, a = np.shape(smoothed_coulomb_stress)
Dx_new = DX//fac
Dy_new = DY//fac

data_counts, years = np.histogram(
    dates, np.linspace(1956, 2019, (2019-1956)+1))
Yearly_Ds = np.zeros((Dx_new, Dy_new, len(years)-1))
stop_i = 0

for ti, t in enumerate(data_counts):
    stop_i += t
    Yearly_Ds[:, :, ti] = np.sum(
        smoothed_coulomb_stress[:, :, stop_i - t:stop_i], axis=2)/(t)


# %%
fac = 1

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

# the final map of pressures
quad1 = ax[0].scatter(np.array(x_res).flatten(), np.array(
    y_res).flatten(), c=np.mean(Yearly_Ds, axis=2), cmap='coolwarm')
ax[0].plot(x_outline[:], y_outline[:], 'k--')
ax[0].set_title('Maximum coulomb stress change')
cb1 = plt.colorbar(quad1, ax=ax[0])
cb1.set_label('Maximum Coulomb stress change (MPa)')
# the final map of deformation
quad2 = ax[1].scatter(np.array(x_res).flatten(), np.array(
    y_res).flatten(), c=smoothed_coulomb_stress[:, :, -1], cmap='coolwarm')
ax[1].plot(x_outline[:], y_outline[:], 'k--')
ax[1].set_title('Maximum coulomb stress change')
cb1 = plt.colorbar(quad2, ax=ax[1])
cb1.set_label('Maximum Coulomb stress change (MPa)')
# the final map of coulomb stress change

# the map labels for the three top subplots
for ii in [0, 1]:
    ax[ii].set_xlabel('X (km)')
    ax[ii].set_ylabel('Y (km)')
    ax[ii].set_xlim([230, 270])
    ax[ii].set_ylim([565, 615])
    ax[ii].set_aspect('equal')

# %%

thresh = 0.2
mask = Functions_spatial.Get_mask(Yearly_Ds, thresh)

avg_ds = np.mean(Yearly_Ds, axis=2)
fig, ax = plt.subplots(1, 2, figsize=(12, 5))

# the final map of pressures
quad1 = ax[0].scatter(np.array(x_res).flatten(), np.array(
    y_res).flatten(), c=avg_ds, cmap='coolwarm')
ax[0].plot(x_outline[:], y_outline[:], 'k--')
ax[0].set_title('Maximum coulomb stress change')
cb1 = plt.colorbar(quad1, ax=ax[0])
cb1.set_label('Maximum Coulomb stress change (MPa)')
# the final map of deformation
quad2 = ax[1].scatter(np.array(x_res).flatten(), np.array(
    y_res).flatten(), c=mask, cmap='coolwarm')
ax[1].plot(x_outline[:], y_outline[:], 'k--')
ax[1].set_title('Maximum coulomb stress change')
cb1 = plt.colorbar(quad2, ax=ax[1])
cb1.set_label('Maximum Coulomb stress change (MPa)')
# the final map of coulomb stress change

# the map labels for the three top subplots
for ii in [0, 1]:
    ax[ii].set_xlabel('X (km)')
    ax[ii].set_ylabel('Y (km)')
    ax[ii].set_xlim([230, 270])
    ax[ii].set_ylim([565, 615])
    ax[ii].set_aspect('equal')

# %%

# spatial averaging of Dst
# Yearly_Ds_combined = np.sum(Yearly_Ds,axis=(0,1))/np.sum(zzz)

Yearly_R0 = np.zeros((2018-1992+1,))


for i in range(1, Yearly_R0.size):

    year = i+1992
    print(year)

    start = np.min(np.where(cat['decDate'] >= year))
    end = np.max(np.where(cat['decDate'] < year+1))

    Yearly_R0[i] = cat['index'][end+1]-cat['index'][start]


plt.show()
plt.plot(range(1992, 2019), Yearly_R0)
plt.show()

# %%

R0 = Yearly_R0
Ds_t = Functions_spatial.Organize_Dst(Yearly_Ds, mask)
Ds_t = Ds_t[-27:, :]*1e6

time = np.array(range(27))
plt.plot(time, R0)
plt.show()
plt.plot(time, Ds_t)

# %% Priors (ranges)
r = 5e-6
t_a = 8700
Asig = 0.006e6
Ds_c = 0.17e6

# Defining ranges of priors
r_min = 4e-7  # 6.2e-7 4.9e-6
r_max = 6e-5  # 2.5e-5 5.1e-6

t_a_min = 6000
t_a_max = 1e4

Asig_min = 0.003e6  # 1e3
Asig_max = 0.01e6  # 1e6

Ds_c_min = 0.1e5  # 0
Ds_c_max = 0.3e6  # .5e6

r_range = r_max-r_min
t_a_range = t_a_max-t_a_min
Asig_range = Asig_max-Asig_min
Ds_c_range = Ds_c_max-Ds_c_min


# u=np.array([r,t_a,Asig,Ds_c])
u_min = np.array([r_min, t_a_min, Asig_min, Ds_c_min])
u_max = np.array([r_max, t_a_max, Asig_max, Ds_c_max])


# %%
N = 1000
M = 4
iters = 300

U0 = Functions_spatial.GetU0_Uniform(N, M, u_min, u_max)

U_history = Functions_spatial.ParticleFilter(
    u_min, u_max, U0, R0, Ds_t, time, iters)

for m in range(4):
    plt.show()
    plt.hist(U_history[:, m], bins=100, range=(u_min[m], u_max[m]))
    plt.show()

# %%
wts = Functions_spatial.FindWeights(U_history, R0, Ds_t, time)
u_pred = np.dot(wts, U_history)

R_pred = Functions_spatial.FindR(Ds_t, u_pred, time)

plt.plot(time, R0, label='Actual')
plt.plot(time, R_pred, label='Predicted')
plt.xlabel("Time (Years)")
plt.ylabel("R")
plt.title("")
plt.legend()
plt.savefig('figs/particle_filer.png', bbox_inches='tight', dpi=400)
plt.show()

# %%