Skip to content

gCMT Catalog accessor #217

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
May 30, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
126 changes: 90 additions & 36 deletions csep/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,9 @@
'__version__'
]

def load_stochastic_event_sets(filename, type='csv', format='native', **kwargs):

def load_stochastic_event_sets(filename, type='csv', format='native',
**kwargs):
""" General function to load stochastic event sets

This function returns a generator to iterate through a collection of catalogs.
Expand Down Expand Up @@ -106,7 +108,8 @@ def load_stochastic_event_sets(filename, type='csv', format='native', **kwargs):
raise ValueError('format must be either "native" or "csep!')


def load_catalog(filename, type='csep-csv', format='native', loader=None, apply_filters=False, **kwargs):
def load_catalog(filename, type='csep-csv', format='native', loader=None,
apply_filters=False, **kwargs):
""" General function to load single catalog

See corresponding class documentation for additional parameters.
Expand All @@ -122,9 +125,12 @@ def load_catalog(filename, type='csep-csv', format='native', loader=None, apply_
Returns (:class:`~csep.core.catalogs.AbstractBaseCatalog`)
"""


if type not in ('ucerf3', 'csep-csv', 'zmap', 'jma-csv', 'ingv_horus', 'ingv_emrcmt', 'ndk') and loader is None:
raise ValueError("type must be one of the following: ('ucerf3', 'csep-csv', 'zmap', 'jma-csv', 'ndk', 'ingv_horus', 'ingv_emrcmt').")
if type not in (
'ucerf3', 'csep-csv', 'zmap', 'jma-csv', 'ingv_horus',
'ingv_emrcmt',
'ndk') and loader is None:
raise ValueError(
"type must be one of the following: ('ucerf3', 'csep-csv', 'zmap', 'jma-csv', 'ndk', 'ingv_horus', 'ingv_emrcmt').")

# map to correct catalog class, at some point these could be abstracted into configuration file
# this maps a human readable string to the correct catalog class and the correct loader function
Expand Down Expand Up @@ -167,7 +173,8 @@ def load_catalog(filename, type='csep-csv', format='native', loader=None, apply_
if loader is None:
loader = class_loader_mapping[type]['loader']

catalog = catalog_class.load_catalog(filename=filename, loader=loader, **kwargs)
catalog = catalog_class.load_catalog(filename=filename, loader=loader,
**kwargs)

# convert to csep format if needed
if format == 'native':
Expand Down Expand Up @@ -213,12 +220,15 @@ def query_comcat(start_time, end_time, min_magnitude=2.50,
# Timezone should be in UTC
t0 = time.time()
eventlist = readers._query_comcat(start_time=start_time, end_time=end_time,
min_magnitude=min_magnitude,
min_latitude=min_latitude, max_latitude=max_latitude,
min_longitude=min_longitude, max_longitude=max_longitude,
max_depth=max_depth)
min_magnitude=min_magnitude,
min_latitude=min_latitude,
max_latitude=max_latitude,
min_longitude=min_longitude,
max_longitude=max_longitude,
max_depth=max_depth)
t1 = time.time()
comcat = catalogs.CSEPCatalog(data=eventlist, date_accessed=utc_now_datetime(), **kwargs)
comcat = catalogs.CSEPCatalog(data=eventlist,
date_accessed=utc_now_datetime(), **kwargs)
print("Fetched ComCat catalog in {} seconds.\n".format(t1 - t0))

if apply_filters:
Expand All @@ -229,9 +239,13 @@ def query_comcat(start_time, end_time, min_magnitude=2.50,

if verbose:
print("Downloaded catalog from ComCat with following parameters")
print("Start Date: {}\nEnd Date: {}".format(str(comcat.start_time), str(comcat.end_time)))
print("Min Latitude: {} and Max Latitude: {}".format(comcat.min_latitude, comcat.max_latitude))
print("Min Longitude: {} and Max Longitude: {}".format(comcat.min_longitude, comcat.max_longitude))
print("Start Date: {}\nEnd Date: {}".format(str(comcat.start_time),
str(comcat.end_time)))
print(
"Min Latitude: {} and Max Latitude: {}".format(comcat.min_latitude,
comcat.max_latitude))
print("Min Longitude: {} and Max Longitude: {}".format(
comcat.min_longitude, comcat.max_longitude))
print("Min Magnitude: {}".format(comcat.min_magnitude))
print(f"Found {comcat.event_count} events in the ComCat catalog.")

Expand Down Expand Up @@ -266,11 +280,14 @@ def query_bsi(start_time, end_time, min_magnitude=2.50,
t0 = time.time()
eventlist = readers._query_bsi(start_time=start_time, end_time=end_time,
min_magnitude=min_magnitude,
min_latitude=min_latitude, max_latitude=max_latitude,
min_longitude=min_longitude, max_longitude=max_longitude,
min_latitude=min_latitude,
max_latitude=max_latitude,
min_longitude=min_longitude,
max_longitude=max_longitude,
max_depth=max_depth)
t1 = time.time()
bsi = catalogs.CSEPCatalog(data=eventlist, date_accessed=utc_now_datetime(), **kwargs)
bsi = catalogs.CSEPCatalog(data=eventlist,
date_accessed=utc_now_datetime(), **kwargs)
print("Fetched BSI catalog in {} seconds.\n".format(t1 - t0))

if apply_filters:
Expand All @@ -280,19 +297,25 @@ def query_bsi(start_time, end_time, min_magnitude=2.50,
bsi = bsi.filter()

if verbose:
print("Downloaded catalog from Bollettino Sismico Italiano (BSI) with following parameters")
print("Start Date: {}\nEnd Date: {}".format(str(bsi.start_time), str(bsi.end_time)))
print("Min Latitude: {} and Max Latitude: {}".format(bsi.min_latitude, bsi.max_latitude))
print("Min Longitude: {} and Max Longitude: {}".format(bsi.min_longitude, bsi.max_longitude))
print(
"Downloaded catalog from Bollettino Sismico Italiano (BSI) with following parameters")
print("Start Date: {}\nEnd Date: {}".format(str(bsi.start_time),
str(bsi.end_time)))
print("Min Latitude: {} and Max Latitude: {}".format(bsi.min_latitude,
bsi.max_latitude))
print(
"Min Longitude: {} and Max Longitude: {}".format(bsi.min_longitude,
bsi.max_longitude))
print("Min Magnitude: {}".format(bsi.min_magnitude))
print(f"Found {bsi.event_count} events in the BSI catalog.")

return bsi


def query_gns(start_time, end_time, min_magnitude=2.950,
min_latitude=-47, max_latitude=-34,
min_longitude=164, max_longitude=180,
max_depth=45.5,
min_latitude=-47, max_latitude=-34,
min_longitude=164, max_longitude=180,
max_depth=45.5,
verbose=True,
apply_filters=False, **kwargs):
"""
Expand Down Expand Up @@ -338,6 +361,29 @@ def query_gns(start_time, end_time, min_magnitude=2.950,

return gns


def query_gcmt(start_time, end_time, min_magnitude=5.0,
max_depth=None,
catalog_id=None,
min_latitude=None, max_latitude=None,
min_longitude=None, max_longitude=None):

eventlist = readers._query_gcmt(start_time=start_time,
end_time=end_time,
min_magnitude=min_magnitude,
min_latitude=min_latitude,
max_latitude=max_latitude,
min_longitude=min_longitude,
max_longitude=max_longitude,
max_depth=max_depth)

catalog = catalogs.CSEPCatalog(data=eventlist,
name='gCMT',
catalog_id=catalog_id,
date_accessed=utc_now_datetime())
return catalog


def load_evaluation_result(fname):
""" Load evaluation result stored as json file

Expand All @@ -361,7 +407,8 @@ def load_evaluation_result(fname):
evaluation_type = json_dict['type']
except:
evaluation_type = 'default'
eval_result = evaluation_result_factory[evaluation_type].from_dict(json_dict)
eval_result = evaluation_result_factory[evaluation_type].from_dict(
json_dict)
return eval_result


Expand Down Expand Up @@ -404,15 +451,18 @@ class with the region and magnitude members correctly assigned.

# sanity checks
if not os.path.exists(fname):
raise FileNotFoundError(f"Could not locate file {fname}. Unable to load forecast.")
raise FileNotFoundError(
f"Could not locate file {fname}. Unable to load forecast.")
# sanity checks
if loader is not None and not callable(loader):
raise AttributeError("Loader must be callable. Unable to load forecast.")
raise AttributeError(
"Loader must be callable. Unable to load forecast.")
extension = os.path.splitext(fname)[-1][1:]
if extension not in forecast_loader_mapping.keys() and loader is None:
raise AttributeError("File extension should be in ('dat','xml','h5','bin') if loader not provided.")
raise AttributeError(
"File extension should be in ('dat','xml','h5','bin') if loader not provided.")

if extension in ('xml','h5','bin'):
if extension in ('xml', 'h5', 'bin'):
raise NotImplementedError

# assign default loader
Expand All @@ -425,7 +475,8 @@ class with the region and magnitude members correctly assigned.
return forecast


def load_catalog_forecast(fname, catalog_loader=None, format='native', type='ascii', **kwargs):
def load_catalog_forecast(fname, catalog_loader=None, format='native',
type='ascii', **kwargs):
""" General function to handle loading catalog forecasts.

Currently, just a simple wrapper, but can contain more complex logic in the future.
Expand All @@ -444,10 +495,12 @@ def load_catalog_forecast(fname, catalog_loader=None, format='native', type='asc
"""
# sanity checks
if not os.path.exists(fname):
raise FileNotFoundError(f"Could not locate file {fname}. Unable to load forecast.")
raise FileNotFoundError(
f"Could not locate file {fname}. Unable to load forecast.")
# sanity checks
if catalog_loader is not None and not callable(catalog_loader):
raise AttributeError("Loader must be callable. Unable to load forecast.")
raise AttributeError(
"Loader must be callable. Unable to load forecast.")
# factory methods for loading different types of catalogs
catalog_loader_mapping = {
'ascii': catalogs.CSEPCatalog.load_ascii_catalogs,
Expand All @@ -456,17 +509,18 @@ def load_catalog_forecast(fname, catalog_loader=None, format='native', type='asc
if catalog_loader is None:
catalog_loader = catalog_loader_mapping[type]
# try and parse information from filename and send to forecast constructor
if format == 'native' and type=='ascii':
if format == 'native' and type == 'ascii':
try:
basename = str(os.path.basename(fname.rstrip('/')).split('.')[0])
split_fname = basename.split('_')
name = split_fname[0]
start_time = strptime_to_utc_datetime(split_fname[1], format="%Y-%m-%dT%H-%M-%S-%f")
start_time = strptime_to_utc_datetime(split_fname[1],
format="%Y-%m-%dT%H-%M-%S-%f")
# update kwargs
_ = kwargs.setdefault('name', name)
_ = kwargs.setdefault('start_time', start_time)
except:
pass
# create observed_catalog forecast
return CatalogForecast(filename=fname, loader=catalog_loader, catalog_format=format, catalog_type=type, **kwargs)

return CatalogForecast(filename=fname, loader=catalog_loader,
catalog_format=format, catalog_type=type, **kwargs)
132 changes: 132 additions & 0 deletions csep/utils/iris.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,132 @@
# python imports
from datetime import datetime
from urllib import request
from urllib.parse import urlencode

# PyCSEP imports
from csep.utils.time_utils import datetime_to_utc_epoch

HOST_CATALOG = "https://service.iris.edu/fdsnws/event/1/query?"
TIMEOUT = 180


def gcmt_search(format='text',
starttime=None,
endtime=None,
updatedafter=None,
minlatitude=None,
maxlatitude=None,
minlongitude=None,
maxlongitude=None,
latitude=None,
longitude=None,
maxradius=None,
catalog='GCMT',
contributor=None,
maxdepth=1000,
maxmagnitude=10.0,
mindepth=-100,
minmagnitude=0,
offset=1,
orderby='time-asc',
host=None,
verbose=False):
"""Search the IRIS database for events matching input criteria.
This search function is a wrapper around the ComCat Web API described here:
https://service.iris.edu/fdsnws/event/1/

This function returns a list of SummaryEvent objects, described elsewhere in this package.
Args:
starttime (datetime):
Python datetime - Limit to events on or after the specified start time.
endtime (datetime):
Python datetime - Limit to events on or before the specified end time.
updatedafter (datetime):
Python datetime - Limit to events updated after the specified time.
minlatitude (float):
Limit to events with a latitude larger than the specified minimum.
maxlatitude (float):
Limit to events with a latitude smaller than the specified maximum.
minlongitude (float):
Limit to events with a longitude larger than the specified minimum.
maxlongitude (float):
Limit to events with a longitude smaller than the specified maximum.
latitude (float):
Specify the latitude to be used for a radius search.
longitude (float):
Specify the longitude to be used for a radius search.
maxradius (float):
Limit to events within the specified maximum number of degrees
from the geographic point defined by the latitude and longitude parameters.
catalog (str):
Limit to events from a specified catalog.
contributor (str):
Limit to events contributed by a specified contributor.
maxdepth (float):
Limit to events with depth less than the specified maximum.
maxmagnitude (float):
Limit to events with a magnitude smaller than the specified maximum.
mindepth (float):
Limit to events with depth more than the specified minimum.
minmagnitude (float):
Limit to events with a magnitude larger than the specified minimum.
offset (int):
Return results starting at the event count specified, starting at 1.
orderby (str):
Order the results. The allowed values are:
- time order by origin descending time
- time-asc order by origin ascending time
- magnitude order by descending magnitude
- magnitude-asc order by ascending magnitude
host (str):
Replace default ComCat host (earthquake.usgs.gov) with a custom host.
Returns:
list: List of SummaryEvent() objects.
"""

# getting the inputargs must be the first line of the method!
inputargs = locals().copy()
newargs = {}

for key, value in inputargs.items():
if value is True:
newargs[key] = 'true'
continue
if value is False:
newargs[key] = 'false'
continue
if value is None:
continue
newargs[key] = value

del newargs['verbose']

events = _search_gcmt(**newargs)

return events


def _search_gcmt(**_newargs):
"""
Performs de-query at ISC API and returns event list and access date

"""
paramstr = urlencode(_newargs)
url = HOST_CATALOG + paramstr
fh = request.urlopen(url, timeout=TIMEOUT)
data = fh.read().decode('utf8').split('\n')
fh.close()
eventlist = []
for line in data[1:]:
line_ = line.split('|')
if len(line_) != 1:
id_ = line_[0]
time_ = datetime.fromisoformat(line_[1])
dt = datetime_to_utc_epoch(time_)
lat = float(line_[2])
lon = float(line_[3])
depth = float(line_[4])
mag = float(line_[10])
eventlist.append((id_, dt, lat, lon, depth, mag))

return eventlist
Loading