<jupyter_start><jupyter_text>Introduction on Using Python to access and analyse GeoNet's Tide Gauge Data In this tutorial we will look at geting waveform data form Tide Gauges. By the end of this notebook you should be able to make a graph like the one below.   Table of Contents  1. Introduction 2. Set up 3. Get the Location of all  Stations with Tide Gauge 4. Look more at one Station 5. Get Waveform Data  &nbsp; 5. 1 Getting Near real-time waveform data 6. Save Waveform Data  1. Introduction In this tutorial we will be learning how to use Python to be able to access the data recorded by GeoNets network of water level sensors’. Then we will look at some ways we can interact with this data in Python.  You will need to install the packages obspy (https://github.com/obspy/obspy/wiki) and basemap (https://matplotlib.org/basemap/) useing `pip install (package name)` in the terminal / command prompt This tutorial assumes that you have basic knowledge of Python.We will be using a package called obspy (https://github.com/obspy/obspy/wiki)  About Geonet tide gauge data Geonet has 3 types of data from tide gauges.Then raw data of the pressure sensorThe data corrected to heightThe data that has been de-tided  Geonet also samples at 2 rates. 10 Hz and 1 Hz. The 3 letter channel code tells you all of this.All tide gauges are in the form (B/L)T(Z/H/T)The first letter is the sample rate. B = 10 Hz and L = 1 Hz, The second letter say that is a water level sensor The third letter say what kind of data it is, Z = raw, H = corrected to height and T = de-tided data For a full list see https://www.geonet.org.nz/data/supplementary/channels   2. Set Up  Import obspy<jupyter_code>from obspy.clients.fdsn import Client as FDSNClient # Import a module to retrieve data from FDSN Web Services 
from obspy import UTCDateTime # Import a module to manage time series format and rename as UTCDateTime
from obspy import read_inventory # Imports module to read stattion data<jupyter_output><empty_output><jupyter_text>For more information on obspy run `help(obspy)`<jupyter_code>#import matplotlib
#matplotlib.use('TkAgg')<jupyter_output><empty_output><jupyter_text>`matplotlib` is the package used to make plots in python. The line `matplotlib.use('TkAgg')` will make all plots open in an interactive window, known as a Tkinter window, so you can zoom in on plots and save them as a png. (The `Tk` is to make it in a Tkinter window and the `Agg` stands for the Anti-Grain Geometry rendering engine that is used to make high-quality images). This is a needed line of code in Python 2 to stop plotting errors but is not need in Python 3.  Set up URL's<jupyter_code>client = FDSNClient("GEONET")
clinet_nrt = FDSNClient('https://service-nrt.geonet.org.nz')<jupyter_output><empty_output><jupyter_text>GeoNet has two FDSN services For waveform dataan archive service holds verified data starting 7 days after collection, and a near real-time service holds data that could be incomplete from the last 8 days  3. Get the Location of all  Stations with Tide Gauges  Set the start and end time of the data you want to select<jupyter_code>starttime = UTCDateTime("2011-03-11T17:00:00.000") # initial time
endtime = UTCDateTime("2011-03-12T07:00:00.000") # end time<jupyter_output><empty_output><jupyter_text>Note that time must be in UTC (Universal Time Coordinated) which is 13 or 12 hours behind nzt (New Zealand Time) (depending on daylight savings time).  Also Note that GMT or Greenwich mean time is the same time as UTC.    In this case we are looking at the 2011 Tōhoku earthquake and tsunami which happened on the 11 March 2011 at 14:46 JST and if it efected new zeeland it would have done so about 12 hours after it happened.     There are online converters if needed.    You can add to the time in number of seconds eg. `stattime + 60` would be the time 2011-03-11 17:01:00.000  Get all stations that have a tide gauge<jupyter_code>Inventory = client.get_stations(starttime = starttime,endtime = endtime, channel="?T?" )<jupyter_output>/home/mkortink/anaconda3/envs/GeoNet/lib/python3.7/site-packages/obspy/io/stationxml/core.py:98: UserWarning: The StationXML file has version 1, ObsPy can read versions (1.0, 1.1). Proceed with caution.
  version, ", ".join(READABLE_VERSIONS)))<jupyter_text>`get_stations` gets all stations that where running at some point between `starttime` and `endtime` and hold to any other conditions put in the `()`such as channel, station, latitude and longitude with a maxradiusAll channels with a code in the form (B/L)T(Z/H/T) are tide gauges The ?s say it will look for data with a channel code in the format (any letter)T(any letter) Therefor the code above finds all stations which have a channel code in the format (any letter)T(any letter)   Print all stations that have a tide gauge<jupyter_code>Inventory<jupyter_output><empty_output><jupyter_text>Plot all stations that have a tide gauge<jupyter_code>Inventory.plot(projection="local")<jupyter_output><empty_output><jupyter_text>`projection="local"` makes a plot that is a map and has all the data points in it. More information can be found at https://matplotlib.org/basemap/  4. Look more at one Station In this section we will learn how to get more information on a station including the instrument response information.  Get more data on the Raoul Island station<jupyter_code>Raoul_Island = client.get_stations(starttime = starttime,endtime = endtime,station="RBCT", level="response" )<jupyter_output>/home/mkortink/anaconda3/envs/GeoNet/lib/python3.7/site-packages/obspy/io/stationxml/core.py:98: UserWarning: The StationXML file has version 1, ObsPy can read versions (1.0, 1.1). Proceed with caution.
  version, ", ".join(READABLE_VERSIONS)))<jupyter_text>Having seen the Raoul Island station on the map and its station code, I can get more data about this station. We can use the parameter `level` to pick the "level" of detail in data we want. The “level” of the data is about how much on each station you get. `level="response"` is the  most complete "level" of data we can get from `get_stations` and includes the instrument response data.`station="RBCT"` say that the code get data for the station with the station code RBCT  Print Raoul Island station metadata<jupyter_code>channels = Raoul_Island [0][0]
print(channels)<jupyter_output>Station RBCT (Raoul Island Boat Cove)
	Station Code: RBCT
	Channel Count: 18/54 (Selected/Total)
	2009-05-18T00:00:00.000000Z - 
	Access: open 
	Latitude: -29.28, Longitude: -177.89, Elevation: 0.0 m
	Available Channels:
		RBCT.40.BTZ, RBCT.40.BTH, RBCT.40.BTT, RBCT.40.LTZ, RBCT.40.LTH, 
		RBCT.40.LTT, RBCT.40.VTZ, RBCT.40.VTH, RBCT.40.VTT, RBCT.41.BTZ, 
		RBCT.41.BTH, RBCT.41.BTT, RBCT.41.LTZ, RBCT.41.LTH, RBCT.41.LTT, 
		RBCT.41.VTZ, RBCT.41.VTH, RBCT.41.VTT<jupyter_text>Remember the 3 letter channel code.All tide gauges are in the form (B/L)T(Z/H/T)The first letter is the sample rate. B = 10 Hz and L = 1 Hz, The second letter say that is a water level sensor The third letter say what kind of data it is, Z = raw, H = corrected to height and T = de-tided data Note that `Raoul_Island [0][0]` says get the first network's first station.  Print Raoul Island Channel metadata including instrument response data for channels at 10Hz<jupyter_code>for i in range (len(channels)): # make a lop over the channels data
        if i==1 or i==4 or i==7 : # only print the 3 channel we want to look more into
            print(channels[i])
            print(channels[i].response)
            print("")<jupyter_output>Channel 'BTZ', Location '41' 
	Time range: 2010-10-27T03:00:00.000000Z - 2017-03-26T22:00:00.000000Z
	Latitude: -29.28, Longitude: -177.89, Elevation: 0.0 m, Local Depth: 0.0 m
	Azimuth: 0.00 degrees from north, clockwise
	Dip: 0.00 degrees down from horizontal
	Channel types: CONTINUOUS, GEOPHYSICAL
	Sampling Rate: 10.00 Hz
	Sensor (Description): Pressure Sensor (Druck PTX)
	Response information available
Channel Response
	From m () to count ()
	Overall Sensitivity: 86402.7 defined at 1.000 Hz
	3 stages:
		Stage 1: PolynomialResponseStage from m to V, gain: 0.206
		Stage 2: CoefficientsTypeResponseStage from V to count, gain: 419430
		Stage 3: FIRResponseStage from count to count, gain: 1

Channel 'VTH', Location '41' 
	Time range: 2010-10-27T03:00:00.000000Z - 2017-03-26T22:00:00.000000Z
	Latitude: -29.28, Longitude: -177.89, Elevation: 0.0 m, Local Depth: 0.0 m
	Azimuth: 0.00 degrees from north, clockwise
	Dip: 0.00 degrees down from horizontal
	Channel types: CONTINUOUS, GEOPHYSICA[...]<jupyter_text>Note that the 2nd, 5th and 8th channels are printed as the have the channel code in the form BT(Z/H/T) In other words we know these 3 channels are tide gauges that have a sample rate of 10Hz.  5. Get Waveform Data  Get Waveform Data from 10Hz Channels<jupyter_code>RBCT= client.get_waveforms('NZ',"RBCT", "40", "BT?", starttime,endtime)<jupyter_output><empty_output><jupyter_text>The parameters are in the order network, station, location, channel, starttime, and endtime.The ? is a wildcard character for only a single character, so we will get data with channel codes BTH, BTT and BTZ  Check basic metadata from miniSEED header<jupyter_code>print(RBCT[0].stats)<jupyter_output>network: NZ
               station: RBCT
              location: 40
               channel: BTH
             starttime: 2011-03-11T16:59:59.969541Z
               endtime: 2011-03-12T06:59:59.969541Z
         sampling_rate: 10.0
                 delta: 0.1
                  npts: 504001
                 calib: 1.0
_fdsnws_dataselect_url: http://service.geonet.org.nz/fdsnws/dataselect/1/query
               _format: MSEED
                 mseed: AttribDict({'dataquality': 'D', 'number_of_records': 2328, 'encoding': 'STEIM2', 'byteorder': '>', 'record_length': 512, 'filesize': 3575808})
            processing: ['ObsPy 1.2.2: trim(endtime=UTCDateTime(2011, 3, 12, 6, 59, 59, 969541)::fill_value=None::nearest_sample=True::pad=False::starttime=UTCDateTime(2011, 3, 11, 16, 59, 59, 969541))']<jupyter_text>Plot this data<jupyter_code>RBCT.plot()<jupyter_output><empty_output><jupyter_text>Note the differences in the 3 plots and remember the 3rd letter say what kind of data it is, Z = raw, H = corrected to height and T = de-tided data.This can be seen clearly when plotted  Print what channels we have saved under RBCT<jupyter_code>print(RBCT)<jupyter_output>3 Trace(s) in Stream:
NZ.RBCT.40.BTH | 2011-03-11T16:59:59.969541Z - 2011-03-12T06:59:59.969541Z | 10.0 Hz, 504001 samples
NZ.RBCT.40.BTT | 2011-03-11T16:59:59.969541Z - 2011-03-12T06:59:59.969541Z | 10.0 Hz, 504001 samples
NZ.RBCT.40.BTZ | 2011-03-11T16:59:59.969541Z - 2011-03-12T06:59:59.969541Z | 10.0 Hz, 504001 samples<jupyter_text>Data is in order BTH,BTT,BTZ so is in order data corrected to height, de-tided, then raw data.  Split the RBCT data up into its channels<jupyter_code>corrected = RBCT[0]
detided = RBCT[1]
raw = RBCT[2]<jupyter_output><empty_output><jupyter_text>5. 1 Getting Near real-time waveform data  Set new start and end time of the data you want to look at<jupyter_code>t1 = UTCDateTime.now()-777600 #9 days ago
t2 = UTCDateTime.now()-518400 #6 days ago<jupyter_output><empty_output><jupyter_text>Note that time must be in UTC (Universal Time Coordinated) which is 13 or 12 hours behind nzt (New Zealand Time) (depending on daylight savings time).There are online converters if needed.`UTCDateTime.now()-777600` say take the time the computer says it is now, convert it to UTC and then takes away 777,600 seconds from it, which is $60$(minutes)$*60$(seconds)$*24$(hours)$*9$(days).   Try both near real-time and  archive service to get waveform data<jupyter_code>try:
    stnrt = clinet_nrt.get_waveforms('NZ',"RBCT", "40", "BTH", t1 , t2)#try near real-time
except:
    print('nrt client not successful')
    
try:
    starc = client.get_waveforms('NZ',"RBCT", "40", "BTH", t1 , t2)#try achive
except:
    print('arc client not successful')<jupyter_output><empty_output><jupyter_text>Remember that GeoNet has two FDSN services An archive service holds verified data starting 7 days after collectionA near real-time service holds data that could be incomplete from the last 8 daysAs The time between 9 and 6 days overlaps this time the code tries to get the data from both services  Merge this data<jupyter_code>st = stnrt #rename near real-time data 
st += starc #add archive data to this 
st.merge(fill_value = 'interpolate') #merge this data<jupyter_output><empty_output><jupyter_text>If data is not merged then it is harder to plot and manipulate. `st += starc` is the same as `st = st + starc`. `fill_value = 'interpolate'` says that if there is a missing value then it will be filled with a "best guess" data point by looking at the two data points next to it.   Plot this data<jupyter_code>st.plot()<jupyter_output><empty_output><jupyter_text>Trim this data<jupyter_code>t = UTCDateTime.now()-604800 #7 days ago
st.trim(t,t+86400)<jupyter_output><empty_output><jupyter_text>`st.trim` is used to "zoom in" on part of the data that you already have. In this example the code is looking at data from 7 days after now till 8 days after now ($60$(minutes)$*60$(seconds)$*24$(hours)$*7$(days)$= 604800$ and $60$(minutes)$*60$(seconds)$*24$(hours)$= 86400$ or 1 day)  Plot the trimmed data<jupyter_code>st.plot()<jupyter_output><empty_output><jupyter_text>6. Save Waveform Data  Save this data<jupyter_code>corrected.write('corrected.csv', format='TSPAIR')
detided.write('detided.csv', format='TSPAIR')
raw.write('raw.csv', format='TSPAIR')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0898537947248549, 'fraction_numerical': 0.05355962089486445, 'mean_word_length': 4.560457516339869, 'pattern_counts': {'":': 0, '<': 83, '<?xml version=': 0, '>': 84, 'https://': 6, 'lorem ipsum': 0, 'www.': 1, 'xml': 2}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '521791', 'license': 'Apache-2.0', 'n_tokens_mistral': 4978, 'n_tokens_neox': 4377, 'n_words': 1704, 'path': 'Tide_Gauge_Data/GeoNet_Tide_Gauge_Data.ipynb', 'repo_name': 'megankortink/data-tutorials', 'size': '13611'}