In [1]:

    

%matplotlib ipympl
import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd
import scipy
import seismic_pandas
import seispy
home = os.environ["HOME"]


    

In [2]:

    

with pd.HDFStore(os.path.join(home,
                              "Projects",
                              "San-Jacinto",
                              "proc",
                              "v2.10",
                              "NLLoc",
                              "loc",
                              "FANG16+CVMH-15.1",
                              "San-Jacinto.h5")) as store:
    df_NLL = pd.concat([store[key] for key in store if not key == "/END_origin"], ignore_index=True)
df = df_NLL


    

In [3]:

    

origin = seispy.coords.as_geographic([33.4991, -116.522, 0])
length, width = 75, 15
strike = -45
plt.close("all")
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
bm = seispy.mapping.Basemap()
bm.scatter(df["lon"], df["lat"], 
           s=0.1, 
           c=df["depth"], 
           cmap=plt.get_cmap("gnuplot_r"), 
           zorder=2,
           alpha=0.1)
df["N"] = df["E"] = df["D"] = np.nan
df.loc[:, ["N", "E", "D"]] = seispy.coords.as_geographic(df[["lat", "lon", "depth"]]
                                                        ).to_cartesian(
                                                        ).to_ned(origin=origin
                                                        ).rotate(np.radians(strike))
df = df[(df["N"].abs() < length)
       &(df["E"].abs() < width)
       &(df["errz"] < 5)]
bm.scatter(df["lon"], df["lat"], 
           s=0.1, 
           c=df["depth"], 
           cmap=plt.get_cmap("gnuplot_r"), 
           zorder=2,
           alpha=0.1)

rect_ned = seispy.coords.as_ned([[-length, -width, 0],
                                 [-length,  width, 0],
                                 [ length,  width, 0],
                                 [ length, -width, 0],
                                 [-length, -width, 0]],
                                origin=origin
                               )
rect_geo = rect_ned.rotate(-np.radians(strike)).to_cartesian().to_geographic()
bm.plot(rect_geo[:,1], rect_geo[:,0], "w")

# line_ned = seispy.coords.as_ned([[-length, 0, 0],
#                                  [length,  0, 0]],
#                                 origin=seispy.coords.as_geographic([lat0, lon0, depth0])
#                                 )
# line_geo = line_ned.rotate(-np.radians(strike)).to_cartesian().to_geographic()
# bm.plot(line_geo[:,1], line_geo[:,0], "k")

# bm.scatter([lon0], [lat0], c="w", zorder=3)


    

    




/Users/malcolcw/Local/anaconda3/envs/py36/lib/python3.6/site-packages/matplotlib/__init__.py:942: MatplotlibDeprecationWarning: nbagg.transparent is deprecated and ignored. Use figure.facecolor instead.
  mplDeprecation)






    






 
 










    
Out[3]:





[<matplotlib.lines.Line2D at 0x1a1d3f0898>]

In [ ]:

    

plt.close("all")
fig = plt.figure(figsize=(9, 2))
ax = fig.add_subplot(1, 1, 1, aspect=1, facecolor="0.5")
ax.scatter(df["N"], 
           df["D"],
           c=df["errz"],
           cmap=plt.get_cmap("gnuplot_r"),
           vmin=0,
           vmax=4,
           s=0.1,
           alpha=0.2)
ax.invert_yaxis()
ax.set_xlim(-length, length)
ax.set_ylim(25, 0)


    

In [ ]:

    

plt.close("all")

fig = plt.figure()
ax = fig.add_subplot(2, 1, 1, aspect=1)
stat, x_edge, y_edge, binnum = scipy.stats.binned_statistic_2d(df["N"], df["D"], df["errz"])
qmesh = ax.pcolormesh(x_edge, y_edge, stat,
                      vmin=0,
                      vmax=5)
fig.colorbar(qmesh, 
             ax=ax,
             shrink=0.5)
ax.invert_yaxis()

ax = fig.add_subplot(2, 1, 2, aspect=1)
stat, x_edge, y_edge, binnum = scipy.stats.binned_statistic_2d(df["N"], df["D"], df["errz"], statistic=np.median)
qmesh = ax.pcolormesh(x_edge,y_edge, stat,
                      vmin=0,
                      vmax=5)
fig.colorbar(qmesh, 
             ax=ax,
             shrink=0.5)
ax.invert_yaxis()


    

In [ ]: