#!/usr/bin/env python
from obspy.clients.fdsn import Client
import numpy as np
from obspy.core import UTCDateTime, read, Stream

# specify webservice
client = Client("RESIF")

# request station metadata at channel level
station_chan = client.get_stations(network="FR,RA",station="*",channel="HN?",level="channel")
# plot stations
station_chan.plot()

# number of stations for FR network
nsta_FR = station_chan[0].selected_number_of_stations
# number of stations for RA network
nsta_RA = station_chan[1].selected_number_of_stations

# get coordinates of every station
FR_coords = np.empty([nsta_FR, 2])
RA_coords = np.empty([nsta_RA, 2])

for k in range(nsta_FR):
    FR_coords[k,0] = station_chan[0][k].latitude
    FR_coords[k,1] = station_chan[0][k].longitude

for m in range(nsta_RA):
    RA_coords[m,0] = station_chan[1][m].latitude
    RA_coords[m,1] = station_chan[1][m].longitude

# get station names
FR_sta = ["" for x in range(nsta_FR)]
for p in range(nsta_FR):
    FR_sta[p] = station_chan[0][p].code

# get event metadata of events within radius=0.5 degrees of stations
starttime = UTCDateTime("2002-01-01")
endtime = UTCDateTime("2002-03-01")
event_radius = 0.5 # degrees
for p in range(nsta_FR):
    # create empty Stream object
    st = Stream()
    # check if station was active during requested time frame
    act = station_chan[0][p].is_active(starttime=starttime, endtime=endtime)
    if act is True:
        lat = FR_coords[p,0] # station latitude
        lon = FR_coords[p,1] # station longitude
        client = Client("IRIS")

        try:
            FR_events = client.get_events(starttime=starttime, endtime=endtime, latitude=lat, longitude=lon, maxradius=event_radius) # produces a Catalog object
            FR_events.write("%s_events.xml" %FR_sta[p], format="QUAKEML")
            n_events = FR_events.count() # number of events found
            print '%i event(s) found for station %s' %(n_events,FR_sta[p])

            # request waveform data from origin time minus 1 minute to origin time plus 10 minutes:
            client = Client("RESIF")                        
            for j in range(n_events):
                # event origin time
                otime = UTCDateTime(FR_events[j].origins[0].time)
                try: 
                    st += client.get_waveforms("FR",FR_sta[p],"*","HN?",otime-60, otime+10*60)
                    print 'Waveform data found for station %s event %i!' %(FR_sta[p],j+1)
                except Exception:
                    print 'No waveform data found for station %s event %i...' %(FR_sta[p],j+1)
                    continue

            FR_events.plot(projection="local",title="Event locations %s" %FR_sta[p]) 
            st.plot()
            # write data to file
            st.write("%s_data.MSEED" %FR_sta[p], format="MSEED")
        except Exception:
            print 'No events found for station %s...' %FR_sta[p]
            continue
    else:
        print 'Station %s was not active in requested time frame' %FR_sta[p]

# request station metadata at response level
client = Client("RESIF")
station_resp = client.get_stations(network="FR,RA", station="*", channel="HN?", level="response")

# plot responses for network FR
station_resp.plot_response(min_freq=0.001,network="FR",station="GRN")