Added NZ catalog reader

csep/__init__.py

@@ -289,6 +289,54 @@ def query_bsi(start_time, end_time, min_magnitude=2.50,
 
     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,
+              verbose=True,
+              apply_filters=False, **kwargs):
+    """
+    Access GNS Science catalog through web service
+
+    Args:
+        start_time: datetime object of start of catalog
+        end_time: datetime object for end of catalog
+        min_magnitude: minimum magnitude to query
+        min_latitude:  maximum magnitude to query
+        max_latitude: max latitude of bounding box
+        min_longitude: min latitude of bounding box
+        max_longitude: max longitude of bounding box
+        max_depth: maximum depth of the bounding box
+        verbose (bool): print catalog summary statistics
+
+    Returns:
+        :class:`csep.core.catalogs.CSEPCatalog
+    """
+
+    # Timezone should be in UTC
+    t0 = time.time()
+    eventlist = readers._query_gns(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)
+    t1 = time.time()
+    gns = catalogs.CSEPCatalog(data=eventlist, date_accessed=utc_now_datetime())
+    if apply_filters:
+        try:
+            gns = gns.filter().filter_spatial()
+        except CSEPCatalogException:
+            gns = gns.filter()
+
+    if verbose:
+        print("Downloaded catalog from GNS Science NZ (GNS) with following parameters")
+        print("Start Date: {}\nEnd Date: {}".format(str(gns.start_time), str(gns.end_time)))
+        print("Min Latitude: {} and Max Latitude: {}".format(gns.min_latitude, gns.max_latitude))
+        print("Min Longitude: {} and Max Longitude: {}".format(gns.min_longitude, gns.max_longitude))
+        print("Min Magnitude: {}".format(gns.min_magnitude))
+        print(f"Found {gns.event_count} events in the gns catalog.")
+
+    return gns
 
 def load_evaluation_result(fname):
     """ Load evaluation result stored as json file

csep/core/catalogs.py

@@ -15,7 +15,8 @@
 from csep.utils.stats import min_or_none, max_or_none
 from csep.utils.calc import discretize
 from csep.core.regions import CartesianGrid2D
-from csep.utils.comcat import SummaryEvent
+import csep.utils.comcat as comcat
+import csep.utils.geonet as geonet
 from csep.core.exceptions import CSEPSchedulerException, CSEPCatalogException, CSEPIOException
 from csep.utils.calc import bin1d_vec
 from csep.utils.constants import CSEP_MW_BINS
@@ -286,7 +287,7 @@ def _get_catalog_as_ndarray(self):
         if isinstance(self.catalog[0], (list, tuple)):
             for i, event in enumerate(self.catalog):
                 catalog[i] = tuple(event)
-        elif isinstance(self.catalog[0], SummaryEvent):
+        elif isinstance(self.catalog[0], (comcat.SummaryEvent, geonet.SummaryEvent)):
             for i, event in enumerate(self.catalog):
                 catalog[i] = (event.id, datetime_to_utc_epoch(event.time),
                               event.latitude, event.longitude, event.depth, event.magnitude)

csep/utils/geonet.py

@@ -0,0 +1,180 @@
+
+# python imports
+from datetime import datetime, timedelta
+from urllib import request
+from urllib.parse import  urlencode
+import json
+
+
+
+class SummaryEvent(object):
+    """Wrapper around summary feature as returned by GeoNet GeoJSON search results.
+    """
+
+    def __init__(self, feature):
+        """Instantiate a SummaryEvent object with a feature.
+        See summary documentation here:
+        https://api.geonet.org.nz/#quakes
+        Args:
+            feature (dict): GeoJSON feature as described at above URL.
+        """
+        self._jdict = feature.copy()
+    @property
+    def url(self):
+        """GeoNet URL.
+        Returns:
+            str: GeoNet URL
+        """
+        url_template= "https://www.geonet.org.nz/earthquake/"
+        return url_template + self._jdict['properties']['publicid']
+
+    @property
+    def latitude(self):
+        """Authoritative origin latitude.
+        Returns:
+            float: Authoritative origin latitude.
+        """
+        return self._jdict['geometry']['coordinates'][1]
+
+    @property
+    def longitude(self):
+        """Authoritative origin longitude.
+        Returns:
+            float: Authoritative origin longitude.
+        """
+        return self._jdict['geometry']['coordinates'][0]
+
+    @property
+    def depth(self):
+        """Authoritative origin depth.
+        Returns:
+            float: Authoritative origin depth.
+        """
+        return self._jdict['properties']['depth']
+
+    @property
+    def id(self):
+        """Authoritative origin ID.
+        Returns:
+            str: Authoritative origin ID.
+        """
+        ## Geonet has eventId or publicid within the properties dict
+        try:
+            return self._jdict['properties']['publicid']
+        except:
+            return self._jdict['properties']['eventId']
+
+    @property
+    def time(self):
+        """Authoritative origin time.
+        Returns:
+            datetime: Authoritative origin time.
+        """
+        from obspy import UTCDateTime
+        time_in_msec = self._jdict['properties']['origintime']
+        # Convert the times
+        if isinstance(time_in_msec, str):
+            event_dtime = UTCDateTime(time_in_msec)
+            time_in_msec = event_dtime.timestamp * 1000
+        time_in_sec = time_in_msec // 1000
+        msec = time_in_msec - (time_in_sec * 1000)
+        dtime = datetime.utcfromtimestamp(time_in_sec)
+        dt = timedelta(milliseconds=msec)
+        dtime = dtime + dt
+        return dtime
+
+    @property
+    def magnitude(self):
+        """Authoritative origin magnitude.
+        Returns:
+            float: Authoritative origin magnitude.
+        """
+        return self._jdict['properties']['magnitude']
+
+    def __repr__(self):
+        tpl = (self.id, str(self.time), self.latitude,
+               self.longitude, self.depth, self.magnitude)
+        return '%s %s (%.3f,%.3f) %.1f km M%.1f' % tpl
+
+
+def gns_search(
+    starttime=None, 
+    endtime=None,
+    minlatitude=-47, 
+    maxlatitude=-34,
+    minlongitude=164, 
+    maxlongitude=180,
+    minmagnitude=2.95,
+    maxmagnitude=None,
+    maxdepth=45.5,
+    mindepth=None):
+
+    """Search the Geonet database for events matching input criteria.
+    This search function is a wrapper around the Geonet Web API described here:
+    https://quakesearch.geonet.org.nz/
+
+    Note:   
+        Geonet has limited search parameters compered to ComCat search parameters, 
+        hence the need for a new function 
+    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.
+        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.
+        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.
+    Returns:
+        list: List of dictionary with event info.
+    """
+    # getting the inputargs must be the first line of the method!
+
+    TIMEFMT = '%Y-%m-%dT%H:%M:%S'
+    TIMEOUT = 120  # how long do we wait for a url to return?
+    try:
+        newargs = {}
+        newargs["bbox"] = f'{minlongitude},{minlatitude},{maxlongitude},{maxlatitude}'
+        newargs["minmag"] = f'{minmagnitude}'
+        newargs["maxdepth"] = f'{maxdepth}'
+        newargs["startdate"] = starttime.strftime(TIMEFMT) 
+        newargs["enddate"] = endtime.strftime(TIMEFMT)
+        if maxmagnitude is not None:
+            newargs["maxmag"] = f'{maxmagnitude}'
+        if mindepth is not None:
+            newargs["mindepth"] = f'{mindepth}'
+
+
+        paramstr = urlencode(newargs)
+        template = "https://quakesearch.geonet.org.nz/geojson?"
+        url = template + '&' + paramstr
+        # print(url)
+        try:
+            fh = request.urlopen(url, timeout=TIMEOUT)
+            data = fh.read().decode('utf8')
+            fh.close()
+            jdict = json.loads(data)
+            events = []
+            for feature in jdict['features']:
+                events.append(SummaryEvent(feature))
+        except Exception as msg:
+            raise Exception(
+                'Error downloading data from url %s.  "%s".' % (url, msg))
+
+        return events
+    except ValueError as e:
+        if len(e.args) > 0 and  'Invalid isoformat string' in e.args[0]:
+            print("Check the input date format. It should follow YYYY-MM-DD \
+                    and is should not be empty")

csep/utils/readers.py

@@ -13,6 +13,7 @@
 # PyCSEP imports
 from csep.utils.time_utils import strptime_to_utc_datetime, strptime_to_utc_epoch, datetime_to_utc_epoch
 from csep.utils.comcat import search
+from csep.utils.geonet import gns_search
 from csep.core.regions import QuadtreeGrid2D
 from csep.core.exceptions import CSEPIOException
 
@@ -705,6 +706,30 @@ def _query_bsi(start_time, end_time, min_magnitude=2.50,
                        starttime=start_time, endtime=end_time, **extra_bsi_params)
 
     return eventlist
+
+# Adding GNS catalog reader
+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, extra_gns_params=None):
+    """
+    Queries GNS catalog.
+    :return: csep.core.Catalog object
+    """
+    extra_gns_params = extra_gns_params or {}
+    geonet_host = 'service.geonet.org.nz'
+    extra_gns_params.update({'host': geonet_host, 'limit': 15000, 'offset': 0})
+    # get eventlist from Comcat
+    eventlist = gns_search(minmagnitude=min_magnitude,
+                       minlatitude=min_latitude, 
+                       maxlatitude=max_latitude,
+                       minlongitude=min_longitude, 
+                       maxlongitude=max_longitude,
+                       maxdepth=max_depth,
+                       starttime=start_time, 
+                       endtime=end_time)
+    return eventlist   
+
 def _parse_datetime_to_zmap(date, time):
         """ Helping function to return datetime in zmap format.
 

tests/artifacts/geonet/vcr_search.yaml

@@ -0,0 +1,56 @@
+interactions:
+- request:
+    body: null
+    headers:
+      Connection:
+      - close
+      Host:
+      - quakesearch.geonet.org.nz
+      User-Agent:
+      - Python-urllib/3.7
+    method: GET
+    uri: https://quakesearch.geonet.org.nz/geojson?&bbox=164%2C-47%2C180%2C-34&minmag=4&maxdepth=45.5&startdate=2023-02-01T00%3A00%3A00&enddate=2023-02-10T00%3A00%3A00
+  response:
+    body:
+      string: '{"type":"FeatureCollection","features":[{"type":"Feature","geometry":{"type":"Point","coordinates":[175.7054138,-37.58729935]},"properties":{"publicid":"2023p087955","eventtype":"earthquake","origintime":"2023-02-02T13:02:43.493Z","modificationtime":"2023-02-17T00:21:34.988Z","depth":5.955107212,"magnitude":4.795154627,"magnitudetype":"MLv","evaluationmethod":"LOCSAT","evaluationstatus":"confirmed","evaluationmode":"manual","earthmodel":"iasp91","usedphasecount":35,"usedstationcount":25,"minimumdistance":0.04634241387,"azimuthalgap":62.78924561,"magnitudeuncertainty":0.2653529139,"originerror":0.3295216371,"magnitudestationcount":110}}]}'
+    headers:
+      Accept-Ranges:
+      - bytes
+      Access-Control-Allow-Methods:
+      - GET
+      Access-Control-Allow-Origin:
+      - '*'
+      Age:
+      - '0'
+      Connection:
+      - close
+      Content-Length:
+      - '645'
+      Content-Security-Policy:
+      - 'base-uri ''none''; frame-ancestors ''self''; object-src ''none''; default-src
+        ''none''; script-src ''self''; connect-src ''self'' https://*.geonet.org.nz
+        https://*.google-analytics.com https://*.analytics.google.com https://*.googletagmanager.com;
+        frame-src ''self'' https://www.youtube.com https://www.google.com; form-action
+        ''self'' https://*.geonet.org.nz; img-src ''self'' *.geonet.org.nz data: https://*.google-analytics.com
+        https://*.googletagmanager.com; font-src ''self'' https://fonts.gstatic.com;
+        style-src ''self'';'
+      Content-Type:
+      - application/vnd.geo+json
+      Date:
+      - Tue, 28 Mar 2023 16:20:00 GMT
+      Referrer-Policy:
+      - no-referrer
+      Strict-Transport-Security:
+      - max-age=63072000
+      Vary:
+      - Accept-Encoding
+      X-Content-Type-Options:
+      - nosniff
+      X-Frame-Options:
+      - DENY
+      X-Xss-Protection:
+      - 1; mode=block
+    status:
+      code: 200
+      message: OK
+version: 1