import matplotlib.pyplot as plt
%matplotlib inline

"""
Explanation: x lines of Python
Download global seismic data after a nuclear test
North Korea made another nuclear test this morning just after 0330 UTC. Inspired by Steven Gibbons's tweets this morning, I wanted to download the data for myself, in as few lines of Python as possible.
Here we go! First some preliminaries:
End of explanation
"""


from obspy.clients.fdsn import Client

client = Client("IRIS")

"""
Explanation: We'll use the awesome ObsPy library. Following the instructions in the ObsPy documentation, and using IRIS as the data source (per Steve's tweets)...
End of explanation
"""


from obspy import UTCDateTime

t = UTCDateTime("2017-09-03_03:30:00")
st = client.get_waveforms("IC", "MDJ", "00", "BHZ", t, t + 5*60)
st.plot()  

"""
Explanation: We want station MDJ in eastern China in the IC network.
End of explanation
"""


st

"""
Explanation: This st thing is the ObsPy stream.
End of explanation
"""


inventory = client.get_stations(network="IC", station="MDJ")
inventory.plot()
plt.show()

inventory

"""
Explanation: We requested the broadband channel, but some stations have lots of others:

LHZ: long period ~1 Hz, with Z for vertical, or N or E for horizontal.
BHZ: broadband 10–80 Hz, with Z for vertical, or N or E for horizontal.
HHZ: broadband 80–250 Hz, with Z for vertical, or N or E for horizontal.
Some stations have loads of others, see Appendix A in the SEED Manual.

Where is the station?
To do this you need basemap and cartopy in your environment:
conda install basemap
conda install cartopy
End of explanation
"""


t1 = UTCDateTime("2016-09-09_00:30:00")
st1 = client.get_waveforms("IC", "MDJ", "00", "BHZ", t1, t1 + 5*60)
st1.plot()  

"""
Explanation: How does the test compare to previous tests?
Let's compare to the test in September 2016:
End of explanation
"""


label = "{network}.{station}.{location}.{channel}".format(**st.traces[0].meta)

plt.figure(figsize=(15,5))
plt.plot(st.traces[0], label="3 Sept 2017")
plt.plot(st1.traces[0], label="9 Sept 2016")
plt.text(-200, 2.5e6, label, size=20)
plt.legend()
plt.grid()
plt.show()

"""
Explanation: We can plot the tests together:
End of explanation
"""


st = client.get_waveforms("IC", "MDJ", "*", "HHZ", t, t + 15*60)
st.plot()  

z = st.traces[0]
NFFT = 2048
Fs = st.traces[0].meta.sampling_rate  # Sample rate in Hz

plt.figure(figsize=(15, 5))
Pxx, freqs, bins, im = plt.specgram(z, NFFT=NFFT, Fs=Fs, noverlap=int(0.9*NFFT), cmap='viridis', vmin=-50, vmax=50)
plt.show()

"""
Explanation: Clearly a substantially more powerful explosion.
Frequency analysis
Let's get a little more data, 15 minutes' worth, and we'll query the high-resolution channel, HHZ.
End of explanation
"""


