Merge pull request #234 from Lachlan-Adams/develop

Addition of Source Specific Station Term (SSST) earthquake relocation algorithm.

seismic/catalogue/download_catalogues/download_GA/eatws/format_data.py

@@ -6,7 +6,7 @@
 from obspy.geodetics import FlinnEngdahl
 from pprint import pprint
 
-def azimuth(lon1, lat1, lon2, lat2, units='degrees'):
+def azimuth(lon1, colat1, lon2, colat2, units='degrees'):
     """
     Function to calculate the azimuth from (lon1, lat1) to (lon2, lat2).
     
@@ -16,14 +16,14 @@ def azimuth(lon1, lat1, lon2, lat2, units='degrees'):
     lon1 : float
         Longitude of point 1.
         
-    lat1 : float
-        Latitude of point 1.
+    colat1 : float
+        Colatitude of point 1.
         
     lon2 : float
         Longitude of point 2.
     
-    lat2 : float
-        Latitude of point 2.
+    colat2 : float
+        Colatitude of point 2.
         
     units : string (optional)
         'degrees' or 'radians' describing the units in which lon1, lat1, lon2,
@@ -33,30 +33,31 @@ def azimuth(lon1, lat1, lon2, lat2, units='degrees'):
     Returns
     -------
     azim : float
-        Azimuth from (lon1, lat1) to (lon2, lat2).
+        Azimuth from (lon1, colat1) to (lon2, colat2).
         
         
     """
     if units == 'degrees':
-        colat1 = 90 - lat1
-        colat2 = 90 - lat2
         degrad = np.pi/180.0
         a = lon1*degrad
         b = colat1*degrad
         x = lon2*degrad
         y = colat2*degrad
     else:
         a = lon1
-        b = np.pi/2 - lat1
+        b = colat1
         x = lon2
-        y = np.pi/2 - lat2
+        y = colat2
     #end if
-    azim = np.arctan(np.sin(x - a)/(np.sin(b)*np.cos(y)/np.sin(y) - \
+    azim = np.arctan(np.sin(x - a)/(np.sin(b)/np.tan(y) - \
                                     np.cos(b)*np.cos(x - a)))
-    if lon2 > lon1 and colat2 < colat1: pass
-    elif colat2 > colat1: azim = azim + np.pi
-    elif lon2 < lon1 and colat2 < colat1: azim = azim + 2*np.pi
-    if units == 'degrees': azim = azim/degrad
+
+    offset1 = np.pi*(colat1 < colat2)
+    offset2 = np.pi*(np.logical_and(colat1 == colat2, b > np.pi/2.0))
+
+    azim = (azim + offset1 + offset2 + 2*np.pi) % (2*np.pi)
+    if units == 'degrees': 
+        azim = azim/degrad
     return azim
 #end func
 
@@ -85,9 +86,9 @@ def get_arrivals_and_picks(event_data):
             # horizontal_slowness=None,
             # horizontal_slowness_errors=None,
             backazimuth=azimuth(station_information["properties"]["longitude"],
-                                station_information["properties"]["latitude"],
+                                90.0-station_information["properties"]["latitude"],
                                 event_data["event_details"]["properties"]["longitude"],
-                                event_data["event_details"]["properties"]["latitude"]),
+                                90.0-event_data["event_details"]["properties"]["latitude"]),
             backazimuth_errors=QuantityError(),
             # slowness_method_id=None,
             # onset=None,
@@ -111,9 +112,9 @@ def get_arrivals_and_picks(event_data):
             phase=station_information["properties"]["phase"],
             # time_correction=None,
             azimuth=azimuth(event_data["event_details"]["properties"]["longitude"],
-                            event_data["event_details"]["properties"]["latitude"],
+                            90.0-event_data["event_details"]["properties"]["latitude"],
                             station_information["properties"]["longitude"],
-                            station_information["properties"]["latitude"]),
+                            90.0-station_information["properties"]["latitude"]),
             distance=station_information["properties"]["distance"],
             # takeoff_angle=None,
             takeoff_angle_errors=QuantityError(),

seismic/catalogue/download_catalogues/download_ISC/README.md

@@ -1,20 +1,32 @@
 # Description
 This script will download data from the ISC web services portal.
 
+
 # download_isc.py
 This script will download catalogues of events from the ISC web services portal, and produce .xml files in the QuakeML format.
 
-Arguments:
+
+Arguments
+---------
 - bbox: Bounding box, specified as "min_lon min_lat max_lon max_lat".
+
 - start_date: Start date in format YYYY-MM-DD.
+
 - start_time: Start time in format hh:mm:ss.
+
 - end_date: End date in format YYYY-MM-DD.
+
 - end_time: End time in format hh:mm:ss.
+
 - mag_range: Magnitude range, specified as "min_mag max_mag".
+
 - output_path: Directory in which to save output xml files.
+
 - split_by_day: If true, one xml file will be produced for each day in the time span, or if else, each month in the time span.
 
-Usage:
+
+Usage
+-----
 >
 > python download_isc.py --output_path ~ --start_date YYYY-MM-DD --start_time hh:mm:ss --end_date YYYY-MM-DD --end_time hh:mm:ss --bbox min_lon min_lat max_lon max_lat --mag_range min_mag max_mag --split_by_day True
 >

seismic/catalogue/download_catalogues/download_ISC/download_isc.py

@@ -1,9 +1,10 @@
 """
-Description:
-    This script will download a catalogue of events from the ISC database given
-    various input criteria by the user. An output file will be generated for 
-    each month or each day in the requested time span (depending on argument 
-    "--split_by_day".
+Description
+-----------
+This script will download a catalogue of events from the ISC database given
+various input criteria by the user. An output file will be generated for each 
+month or each day in the requested time span (depending on argument 
+"--split_by_day").
     
 Developer: Lachlan Adams 
 Contact: lachlan.adams@ga.gov.au or lachlan.adams.1996@outlook.com

seismic/catalogue/download_catalogues/download_USGS/README.md

@@ -1,33 +1,52 @@
 # Description
 These scripts will download catalogues of events from the USGS web services portal.
 
+
 # download_usgs_hypocentres.py
 This script will download hypocentres for events on the USGS database.
 
-Arguments:
+
+Arguments
+---------
 - bbox: Bounding box, specified as "min_lon min_lat max_lon max_lat".
+
 - start_date: Start date in format YYYY-MM-DD.
+
 - start_time: Start time in format hh:mm:ss.
+
 - end_date: End date in format YYYY-MM-DD.
+
 - end_time: End time in format hh:mm:ss.
+
 - mag_range: Magnitude range, specified as "min_mag max_mag".
-- output_path: Directory in which to save output xml files.
+
+- output_path: Directory in which to save output files.
+
 - output_format: Desired save file format. Possible values are 'csv', 'xml', 'text'.
 
-Usage:
+
+Usage
+-----
 >
 > python download_usgs_hypocentres.py --output_path ~ --start_date YYYY-MM-DD --start_time hh:mm:ss --end_date YYYY-MM-DD --end_time hh:mm:ss --bbox min_lon min_lat max_lon max_lat --mag_range min_mag max_mag --output_format format
 >
 
+
 # download_usgs_quakeml.py
 This script will download quakeml files for events on the USGS database. Events on this database are stored in .zip files spanning a week each. There is no way to filter the returned events except by origin time. The script connects to ftp://hazards.cr.usgs.gov/ and extracts events from directory /NEICPDE/quakeml/.
 
-Arguments:
+
+Arguments
+---------
 - start_date: Start date in format YYYY-MM-DD.
+
 - end_date: End date in format YYYY-MM-DD.
+
 - output_path: Directory in which to save output xml files.
 
-Usage:
+
+Usage
+-----
 >
 > python download_usgs_catalogue.py --output_path ~ --start_date YYYY-MM-DD --end_date YYYY-MM-DD
 >

seismic/catalogue/download_catalogues/download_USGS/download_usgs_hypocentres.py

@@ -1,9 +1,8 @@
 """
-Description:
-    This script will download a catalogue of events from the ISC database given
-    various input criteria by the user. An output file will be generated for 
-    each month or each day in the requested time span (depending on argument 
-    "--split_by_day".
+Description
+-----------
+This script will download a catalogue of events from the USGS database given 
+various input criteria by the user. 
     
 Developer: Lachlan Adams 
 Contact: lachlan.adams@ga.gov.au or lachlan.adams.1996@outlook.com

seismic/catalogue/download_catalogues/download_USGS/download_usgs_quakeml.py

@@ -1,9 +1,10 @@
 """
-Description:
-    This script will download a catalogue of events from the ISC database given
-    various input criteria by the user. An output file will be generated for 
-    each month or each day in the requested time span (depending on argument 
-    "--split_by_day".
+Description
+-----------
+This script will download a catalogue of events from the USGS database given 
+various input criteria by the user. The USGS ftp server doesn't allow filtering 
+of the requested data so all events within the specified time range are 
+downloaded. A single event is stored in each .xml file returned.
     
 Developer: Lachlan Adams 
 Contact: lachlan.adams@ga.gov.au or lachlan.adams.1996@outlook.com

seismic/catalogue/merge_catalogues_python_new/Convert_Catalogues_To_CSV.py

@@ -146,61 +146,7 @@ def read_station_kml(path):
     #end if
     return dct
 #end func
-
-def azimuth(lon1, lat1, lon2, lat2, units='degrees'):
-    """
-    Function to calculate the azimuth from (lon1, lat1) to (lon2, lat2).
-    
-    
-    Parameters
-    ----------
-    lon1 : float
-        Longitude of point 1.
-        
-    lat1 : float
-        Latitude of point 1.
-        
-    lon2 : float
-        Longitude of point 2.
-    
-    lat2 : float
-        Latitude of point 2.
-        
-    units : string (optional)
-        'degrees' or 'radians' describing the units in which lon1, lat1, lon2,
-        lat2 are input. Default is 'degrees'.
-        
-        
-    Returns
-    -------
-    azim : float
-        Azimuth from (lon1, lat1) to (lon2, lat2).
-        
-        
-    """
-    if units == 'degrees':
-        colat1 = 90 - lat1
-        colat2 = 90 - lat2
-        degrad = np.pi/180.0
-        a = lon1*degrad
-        b = colat1*degrad
-        x = lon2*degrad
-        y = colat2*degrad
-    else:
-        a = lon1
-        b = np.pi/2 - lat1
-        x = lon2
-        y = np.pi/2 - lat2
-    #end if
-    azim = np.arctan(np.sin(x - a)/(np.sin(b)*np.cos(y)/np.sin(y) - \
-                                    np.cos(b)*np.cos(x - a)))
-    if lon2 > lon1 and colat2 < colat1: pass
-    elif colat2 > colat1: azim = azim + np.pi
-    elif lon2 < lon1 and colat2 < colat1: azim = azim + 2*np.pi
-    if units == 'degrees': azim = azim/degrad
-    return azim
-#end func
-
+
 def ang_dist(lon1, lat1, lon2, lat2, units='degrees'):
     """
     Function to calculate the angular distance from (lon1, lat1) to 

seismic/catalogue/merge_catalogues_python_new/Filter_Catalogues.py

@@ -64,30 +64,31 @@ def azimuth(lon1, lat1, lon2, lat2, units='degrees'):
     Returns
     -------
     azim : float
-        Azimuth from (lon1, lat1) to (lon2, lat2).
+        Azimuth from (lon1, colat1) to (lon2, colat2).
         
         
     """
     if units == 'degrees':
-        colat1 = 90 - lat1
-        colat2 = 90 - lat2
         degrad = np.pi/180.0
         a = lon1*degrad
-        b = colat1*degrad
+        b = (90.0-lat1)*degrad
         x = lon2*degrad
-        y = colat2*degrad
+        y = (90.0-lat2)*degrad
     else:
         a = lon1
-        b = np.pi/2 - lat1
+        b = np.pi/2.0 - lat1
         x = lon2
-        y = np.pi/2 - lat2
+        y = mp.pi/2.0 - lat2
     #end if
-    azim = np.arctan(np.sin(x - a)/(np.sin(b)*np.cos(y)/np.sin(y) - \
+    azim = np.arctan(np.sin(x - a)/(np.sin(b)/np.tan(y) - \
                                     np.cos(b)*np.cos(x - a)))
-    if lon2 > lon1 and colat2 < colat1: pass
-    elif colat2 > colat1: azim = azim + np.pi
-    elif lon2 < lon1 and colat2 < colat1: azim = azim + 2*np.pi
-    if units == 'degrees': azim = azim/degrad
+
+    offset1 = np.pi*(lat1 > lat2)
+    offset2 = np.pi*(np.logical_and(lat1 == lat2, b > np.pi/2.0))
+
+    azim = (azim + offset1 + offset2 + 2*np.pi) % (2*np.pi)
+    if units == 'degrees': 
+        azim = azim/degrad
     return azim
 #end func
 

seismic/catalogue/merge_catalogues_python_new/Merge_Catalogues.py

@@ -16,7 +16,8 @@
 def reorder_list_by_date(lst):
     """
     Takes a list of events and picks, extracts the event rows, sorts them by 
-    origin time, and returns the sorted list with attached picks.
+    origin time, and returns the sorted list with attached picks. Also converts
+    depth to kilometres so it can be used by the pick harvester routine.
     
     
     Parameters
@@ -74,7 +75,8 @@ def reorder_list_by_date(lst):
         if row[0][0] == '#':
             i = i + 1
             row[-1] = i
-            row[-3] = str('smi:local/' + str(i))
+            row[-3] = str(i) #str('smi:local/' + str(i))
+            row[9] = str(float(row[9])/1e3)
         #end if
     #end for
 

seismic/catalogue/merge_catalogues_python_new/README.md

@@ -136,4 +136,4 @@ Usage
 -----
 >
 > python3 Order_csv_output_by_date.py --input_file ~ --output_file ~
->
+>