Merge pull request #139 from mhearne-usgs/countryfix

Countryfix

Author: emthompson-usgs

bin/getcsv

@@ -3,7 +3,9 @@ import argparse
 import sys
 
 from libcomcat.search import search, count
-from libcomcat.utils import maketime
+from libcomcat.utils import (maketime, check_ccode,
+                             get_country_bounds, filter_by_country,
+                             BUFFER_DISTANCE_KM)
 from libcomcat.dataframes import (get_detail_data_frame,
                                   get_summary_data_frame)
 
@@ -97,6 +99,18 @@ def get_parser():
                         metavar=('lonmin', 'lonmax', 'latmin', 'latmax'),
                         dest='bounds', type=float, nargs=4,
                         help=helpstr)
+
+    country_str = '''Specify three character country code and earthquakes 
+    from inside country polygon (50m resolution) will be returned. Earthquakes 
+    in the ocean likely will NOT be returned.'''
+    parser.add_argument('--country', help=country_str)
+
+    buffer_str = '''Use in conjunction with --country. Specify a buffer in km
+    around country border where events will be selected.
+    '''
+    parser.add_argument('--buffer', help=buffer_str,
+                        type=int, default=BUFFER_DISTANCE_KM)
+
     helpstr = 'Search radius in KM (use instead of bounding box)'
     parser.add_argument('-r', '--radius', dest='radius',
                         metavar=('lat', 'lon', 'rmax'),
@@ -180,10 +194,17 @@ def main():
     parser = get_parser()
     args = parser.parse_args()
 
+    tsum = (args.bounds is not None) + \
+        (args.radius is not None) + (args.country is not None)
+    if tsum != 1:
+        print('Please specify a bounding box, radius, or country code.')
+        sys.exit(1)
+
     latitude = None
     longitude = None
     radiuskm = None
     lonmin = latmin = lonmax = latmax = None
+    bounds = None
     if args.radius:
         latitude = args.radius[0]
         longitude = args.radius[1]
@@ -196,6 +217,17 @@ def main():
             lonmin -= 360
     else:
         lonmin, lonmax, latmin, latmax = None, None, None, None
+        bounds = (lonmin, lonmax, latmin, latmax)
+
+    if args.country:
+        ccode = args.country
+        if not check_ccode(ccode):
+            curl = 'https://en.wikipedia.org/wiki/ISO_3166-1_alpha-2'
+            fmt = '%s is not a valid ISO 3166 country code. See %s for the list.'
+            tpl = (ccode, curl)
+            print(fmt % tpl)
+            sys.exit(1)
+        bounds = get_country_bounds(ccode, args.buffer)  # this returns a list
 
     minmag = 0.0
     maxmag = 9.9
@@ -204,7 +236,46 @@ def main():
         maxmag = args.magRange[1]
 
     if args.getCount:
-        nevents = count(starttime=args.startTime,
+        if isinstance(bounds, tuple) or bounds is None:
+            nevents = count(starttime=args.startTime,
+                            endtime=args.endTime,
+                            updatedafter=args.after,
+                            minlatitude=latmin,
+                            maxlatitude=latmax,
+                            minlongitude=lonmin,
+                            maxlongitude=lonmax,
+                            latitude=latitude,
+                            longitude=longitude,
+                            maxradiuskm=radiuskm,
+                            catalog=args.catalog,
+                            contributor=args.contributor,
+                            maxmagnitude=maxmag,
+                            minmagnitude=minmag,
+                            producttype=args.limitByProductType,
+                            verbose=args.verbose)
+        else:
+            for lonmin, lonmax, latmin, latmax in bounds:
+                nevents = 0
+                nevents += count(starttime=args.startTime,
+                                 endtime=args.endTime,
+                                 updatedafter=args.after,
+                                 minlatitude=latmin,
+                                 maxlatitude=latmax,
+                                 minlongitude=lonmin,
+                                 maxlongitude=lonmax,
+                                 latitude=latitude,
+                                 longitude=longitude,
+                                 maxradiuskm=radiuskm,
+                                 catalog=args.catalog,
+                                 contributor=args.contributor,
+                                 maxmagnitude=maxmag,
+                                 minmagnitude=minmag,
+                                 producttype=args.limitByProductType,
+                                 verbose=args.verbose)
+        print('There are %i events matching input criteria.' % nevents)
+        sys.exit(0)
+    if isinstance(bounds, tuple) or bounds is None:
+        events = search(starttime=args.startTime,
                         endtime=args.endTime,
                         updatedafter=args.after,
                         minlatitude=latmin,
@@ -219,31 +290,34 @@ def main():
                         maxmagnitude=maxmag,
                         minmagnitude=minmag,
                         producttype=args.limitByProductType,
+                        host=args.host,
                         verbose=args.verbose)
-        print('There are %i events matching input criteria.' % nevents)
-        sys.exit(0)
-
-    if args.bounds and args.radius:
-        print('Please specify either a bounding box OR radius search.')
-        sys.exit(1)
-
-    events = search(starttime=args.startTime,
-                    endtime=args.endTime,
-                    updatedafter=args.after,
-                    minlatitude=latmin,
-                    maxlatitude=latmax,
-                    minlongitude=lonmin,
-                    maxlongitude=lonmax,
-                    latitude=latitude,
-                    longitude=longitude,
-                    maxradiuskm=radiuskm,
-                    catalog=args.catalog,
-                    contributor=args.contributor,
-                    maxmagnitude=maxmag,
-                    minmagnitude=minmag,
-                    producttype=args.limitByProductType,
-                    host=args.host,
-                    verbose=args.verbose)
+    else:
+        events = []
+        for i, tbounds in enumerate(bounds):
+            lonmin, lonmax, latmin, latmax = tbounds
+            if args.verbose:
+                fmt = 'Checking bounds %i of %i for %s...\n'
+                tpl = (i + 1, len(bounds), ccode)
+                sys.stderr.write(fmt % tpl)
+            tevents = search(starttime=args.startTime,
+                             endtime=args.endTime,
+                             updatedafter=args.after,
+                             minlatitude=latmin,
+                             maxlatitude=latmax,
+                             minlongitude=lonmin,
+                             maxlongitude=lonmax,
+                             latitude=latitude,
+                             longitude=longitude,
+                             maxradiuskm=radiuskm,
+                             catalog=args.catalog,
+                             contributor=args.contributor,
+                             maxmagnitude=maxmag,
+                             minmagnitude=minmag,
+                             producttype=args.limitByProductType,
+                             host=args.host,
+                             verbose=args.verbose)
+            events += tevents
 
     if not len(events):
         print('No events found matching your search criteria. Exiting.')
@@ -272,6 +346,9 @@ def main():
         col_list.remove(column)
     df = df[first_columns + col_list]
 
+    if args.country:
+        df = filter_by_country(df, ccode, buffer_km=args.buffer)
+
     if args.verbose:
         sys.stderr.write('Created table...saving %i records to %s.\n' %
                          (len(df), args.filename))

bin/getproduct

@@ -8,12 +8,20 @@ from datetime import datetime, timedelta
 
 # local imports
 from libcomcat.search import search, get_event_by_id
-from libcomcat.utils import maketime, makedict
+from libcomcat.utils import (maketime, makedict, check_ccode,
+                             get_country_bounds, filter_by_country,
+                             BUFFER_DISTANCE_KM)
 from libcomcat.classes import VersionOption
 
 # third party imports
 import numpy as np
 import pandas as pd
+from shapely.geometry import Point
+
+
+class MyFormatter(argparse.RawTextHelpFormatter,
+                  argparse.ArgumentDefaultsHelpFormatter):
+    pass
 
 
 def _get_product_from_detail(detail, product, contents, folder,
@@ -91,7 +99,7 @@ def get_parser():
     %(prog)s finite-fault .mr -o ~/tmp/chile -b -76.509 -49.804  -67.72 -17.427 -s 2007-01-01 -e 2016-05-01 -m 6.5 9.9
     '''
     parser = argparse.ArgumentParser(
-        description=desc, formatter_class=argparse.RawDescriptionHelpFormatter)
+        description=desc, formatter_class=MyFormatter)
     # positional arguments
     parser.add_argument('product', metavar='PRODUCT',
                         help='The name of the desired product (shakemap, dyfi, etc.)')
@@ -103,6 +111,18 @@ def get_parser():
     parser.add_argument('-b', '--bounds', metavar=('lonmin', 'lonmax', 'latmin', 'latmax'),
                         dest='bounds', type=float, nargs=4,
                         help='Bounds to constrain event search [lonmin lonmax latmin latmax]')
+
+    country_str = '''Specify three character country code and earthquakes 
+    from inside country polygon (50m resolution) will be returned. Earthquakes 
+    in the ocean likely will NOT be returned.'''
+    parser.add_argument('--country', help=country_str)
+
+    buffer_str = '''Use in conjunction with --country. Specify a buffer in km
+    around country border where events will be selected.
+    '''
+    parser.add_argument('--buffer', help=buffer_str,
+                        type=int, default=BUFFER_DISTANCE_KM)
+
     parser.add_argument('-r', '--radius', dest='radius', metavar=('lat', 'lon', 'rmax'), type=float,
                         nargs=3, help='Search radius in KM (use instead of bounding box)')
     parser.add_argument('-s', '--start-time', dest='startTime', type=maketime,
@@ -149,8 +169,10 @@ def main():
                                  args.source, list_only=args.list_only)
         sys.exit(0)
 
-    if args.bounds and args.radius:
-        print('Please specify either a bounding box OR radius search.')
+    tsum = (args.bounds is not None) + \
+        (args.radius is not None) + (args.country is not None)
+    if tsum != 1:
+        print('Please specify a bounding box, radius, or country code.')
         sys.exit(1)
 
     latitude = None
@@ -172,6 +194,7 @@ def main():
     else:
         endtime = args.endTime
 
+    bounds = None
     if args.radius:
         latitude = args.radius[0]
         longitude = args.radius[1]
@@ -182,38 +205,86 @@ def main():
         # fix longitude bounds when crossing dateline
         if lonmin > lonmax and lonmax >= -180:
             lonmin -= 360
+    else:
+        lonmin, lonmax, latmin, latmax = None, None, None, None
+        bounds = (lonmin, lonmax, latmin, latmax)
+
+    if args.country:
+        ccode = args.country
+        if not check_ccode(ccode):
+            curl = 'https://en.wikipedia.org/wiki/ISO_3166-1_alpha-2'
+            fmt = '%s is not a valid ISO 3166 country code. See %s for the list.'
+            tpl = (ccode, curl)
+            print(fmt % tpl)
+            sys.exit(1)
+        bounds = get_country_bounds(ccode, args.buffer)  # this returns a list
 
     minmag = 0.0
     maxmag = 9.9
     if args.magRange:
         minmag = args.magRange[0]
         maxmag = args.magRange[1]
 
-    events = search(starttime=starttime,
-                    endtime=endtime,
-                    updatedafter=args.after,
-                    minlatitude=latmin,
-                    maxlatitude=latmax,
-                    minlongitude=lonmin,
-                    maxlongitude=lonmax,
-                    latitude=latitude,
-                    longitude=longitude,
-                    maxradiuskm=radiuskm,
-                    catalog=args.catalog,
-                    contributor=args.contributor,
-                    producttype=args.product,
-                    eventtype=args.eventType,
-                    maxmagnitude=maxmag,
-                    minmagnitude=minmag,
-                    host=args.host)
+    if isinstance(bounds, tuple) or bounds is None:
+        events = search(starttime=starttime,
+                        endtime=endtime,
+                        updatedafter=args.after,
+                        minlatitude=latmin,
+                        maxlatitude=latmax,
+                        minlongitude=lonmin,
+                        maxlongitude=lonmax,
+                        latitude=latitude,
+                        longitude=longitude,
+                        maxradiuskm=radiuskm,
+                        catalog=args.catalog,
+                        contributor=args.contributor,
+                        producttype=args.product,
+                        eventtype=args.eventType,
+                        maxmagnitude=maxmag,
+                        minmagnitude=minmag,
+                        host=args.host)
+    else:
+        events = []
+        for i, tbounds in enumerate(bounds):
+            lonmin, lonmax, latmin, latmax = tbounds
+            tevents = search(starttime=starttime,
+                             endtime=endtime,
+                             updatedafter=args.after,
+                             minlatitude=latmin,
+                             maxlatitude=latmax,
+                             minlongitude=lonmin,
+                             maxlongitude=lonmax,
+                             latitude=latitude,
+                             longitude=longitude,
+                             maxradiuskm=radiuskm,
+                             catalog=args.catalog,
+                             contributor=args.contributor,
+                             producttype=args.product,
+                             eventtype=args.eventType,
+                             maxmagnitude=maxmag,
+                             minmagnitude=minmag,
+                             host=args.host)
+            events += tevents
 
     if not len(events):
         print('No events found matching your search criteria. Exiting.')
         sys.exit(0)
 
+    if args.country:
+        ids = [event.id for event in events]
+        lats = [event.latitude for event in events]
+        lons = [event.longitude for event in events]
+        df = pd.DataFrame({'id': ids, 'latitude': lats, 'longitude': lons})
+        df2 = filter_by_country(df, ccode, buffer_km=args.buffer)
+        events = [event for event in events if event.id in df2['id'].unique()]
+
     for event in events:
         if not event.hasProduct(args.product):
             continue
+        if args.country:
+            elon = event.longitude
+            elat = event.latitude
+            point = Point(elon, elat)
         try:
             detail = event.getDetailEvent(includesuperseded=True)
         except Exception as e:

libcomcat/data/ne_50m_admin_0_countries.cpg

@@ -0,0 +1 @@
+UTF-8

libcomcat/data/ne_50m_admin_0_countries.dbf

@@ -0,0 +1 @@
+GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]