Last active
September 22, 2023 13:52
-
-
Save bradlipovsky/994dbc715bc6be7cc048829df3ebbc28 to your computer and use it in GitHub Desktop.
extremely simplified eq location
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import numpy as np | |
import matplotlib.pyplot as plt | |
import datetime | |
import geopy.distance | |
T0 = datetime.datetime(2023,9,16,12,31,0).timestamp() # test origin time | |
v = 1 # surface wave speed | |
def G(lon,lat,m,vv): | |
#m = (x0,y0,t0) is the source parameter vector | |
# return np.sqrt( (xx-m[0])**2 + (yy-m[1])**2 ) / vv | |
out = [] | |
for lo,la in zip(lon,lat): | |
this = geopy.distance.geodesic( (la,lo) , (m[1],m[0]) ).km | |
out.append(this) | |
return out | |
# predict travel times at sites: IU.KBS,IU.KONO,IU.SFJD,NS.JNW,DK.SCO | |
# (the following vectors are arranged in that order) | |
lon_obs = np.array([11.9385,9.5946,-50.62076,-8.42367,-21.949699]) | |
lat_obs = np.array([78.9154,59.6521,66.9961,71.02734,70.485603]) | |
tobs = np.array([datetime.datetime(2023,9,16,12,40,0).timestamp(), | |
datetime.datetime(2023,9,16,12,45,0).timestamp(), | |
datetime.datetime(2023,9,16,12,40,0).timestamp(), | |
datetime.datetime(2023,9,16,12,38,0).timestamp(), | |
datetime.datetime(2023,9,16,12,36,0).timestamp()]) | |
# now loop over test event locations | |
nlon=20 | |
nlat=20 | |
lons = np.linspace(-30,-20,nlon) | |
lats = np.linspace(70,82,nlat) | |
residuals = np.zeros((nlon,nlat)) | |
for i,lon_test in enumerate(lons): | |
for j,lat_test in enumerate(lats): | |
tt = G(lon_obs,lat_obs,[lon_test,lat_test],v) | |
residuals[i,j] = np.mean(np.abs(tt - (tobs - T0))) | |
fig,ax=plt.subplots() | |
c=plt.contourf(lons,lats,residuals) | |
plt.colorbar(c,label='Misfit (seconds)') | |
ax.axis('equal') | |
plt.xlabel('Longitude') | |
plt.ylabel('Latitude') | |
plt.show() |
Author
bradlipovsky
commented
Sep 22, 2023
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment