Commit c3f6c217 authored by Mijian Xu's avatar Mijian Xu 😷
Browse files
parent 45662c3a
*.pyc
bin/updateCatalog
bin/bqmail
bin/bqmail_continue
bin/searchDMC
doc/.texpadtmp/bqmail.aux
doc/.texpadtmp/bqmail.log
doc/.texpadtmp/bqmail.out
doc/.texpadtmp/bqmail.synctex.gz
doc/.texpadtmp/bqmail.toc
bin/bqmail_conti
*/__pycache__/
.idea/
*.egg-info
build/
dist/
.vscode/
.DS_Store
*/.DS_Store
This diff is collapsed.
bqmail
===========
Python scripts to request waveform data of events from IRIS
# BQMail
bqmail: Send request mail of events data or custom data to IRIS DMC.
bqmail\_conti: Send request mail of continious data to IRIS DMC.
**BQMail** is a Python module for sending mails to apply for seismic data from the IRIS DMC. It is a front-end API of the BREQ_fast.
searchDMC: Search available stations in IRIS DMC.
## Installation
**BQMail** can currently run on Linux and MAC OSX. **BQMail** is running and testing on Python 3.7.
### Installation via PyPI
```
pip install bqmail
```
uodateCatalog: Update CMT catalog automaticly.
### Installation from source code
The latest version of the **BQMail** is available on Gitlab:
```
git clone https://git.nju.edu.cn/xumi1993/bqmail2.0.git bqmail
```
Then you can install this version:
```
cd bqmail
pip install .
```
## A quick example:
```python
from bqmail.mail import BQMail
from obspy import UTCDateTime
bq = BQMail('xxx@163.com', server='smtp.163.com', password='xxx', username='bqmail')
bq.query_events(starttime=UTCDateTime(2017, 1, 1), endtime=UTCDateTime(2018, 1, 1),
minmagnitude=5.5, catalog='GCMT')
bq.query_stations(network='CB', station='LZH')
bq.send_mail(time_before=0, time_after=1000)
```
\ No newline at end of file
/opt/bqmail/download_seed.py
\ No newline at end of file
#!/usr/bin/env python
#
#Author: Mijian Xu at NJU
#
#Revision History:
# 2014/11/06
# 2015/01/05
# 2015/02/11
# 2015/04/29
# 2015/05/01
# 2015/09/26
# 2015/11/06
# 2017/09/11
#
import datetime
import os, re
import sys, getopt
import time
from obspy import taup, UTCDateTime
import distaz
from util import sendmail, generatemsg, wsfetch
import configparser
config = configparser.ConfigParser()
def Usage():
print('Usage:')
print('python bqmail.py -Nnetwork -Sstation -b -Bsec_begin/sec_end [-Cchannel] [-Plat/lon/phase] [-Llocation] [-cdatetimefile] [-Fformat] [-Mmagmin/magmax] head.cfg')
print('-N -- Network.')
print('-S -- Station.')
print('-b -- Limit to events occurring on or after the specified start time.\n'
' Date and time format: YYYY-MM-DDThh:mm:ss (e.g., 1997-01-31T12:04:32)\n'
' YYYY-MM-DD (e.g., 1997-01-31)')
print('-e -- Limit to events occurring on or before the specified end time\n'
' with the same date and time format as \"-b\".')
print('-B -- Time before/after origal time of events in seconds.')
print('-C -- Channel (e.g., ?H?, HHZ, BH?). Default: BH?')
print('-P -- specify the lat/lon of station and require data by phase. e.g., 20/100/SKS')
print('-L -- Location identifier.')
print('-c -- Directory of date time file. format: "2015,01,04,1,0,0 2015,01,04,10,0,0"')
print('-F -- File format (SEED or miniseed). Default: SEED')
print('-M -- Magnitude range.')
print('head.cfg -- Config file.')
print('Example: bqmail -NCB -SNJ2 -b2015-2-3 -e2015-4-3 -P32.05/118.85/P -B-200/1000 head.cfg')
print(' bqmail -NIC -SBJT -b2015-2-3T00:12:23 -e2015-4-3 -B-100/600 -L10 -Fminiseed head.cfg')
head = ''
argv = sys.argv[1:]
for o in argv:
if os.path.isfile(o):
head = o
argv.remove(o)
break
try:
opts,args = getopt.getopt(argv, "hN:S:C:b:e:B:L:c:F:P:M:")
except:
print('Arguments are not found!')
Usage()
sys.exit(1)
if opts == []:
Usage()
sys.exit(1)
isph = 0
isyrange = 1
chan = "BH?"
fformat = "seed"
loca = ''
magmin = 0
magmax = 10
for op, value in opts:
if op == "-N":
network = value
elif op == "-S":
station = value
elif op == "-b":
starttime = value
elif op == "-e":
endtime = value
elif op == "-c":
datetimefile = value
isyrange = 0
elif op == "-B":
timerange = value
elif op == "-C":
chan = value
elif op == "-L":
loca = value
elif op == "-F":
fformat = value
elif op == "-P":
stla = float(value.split('/')[0])
stlo = float(value.split('/')[1])
phase = value.split('/')[2]
isph = 1
mod=taup.TauPyModel(model='iasp91')
elif op == "-M":
magmin = float(value.split('/')[0])
magmax = float(value.split('/')[1])
elif op == "-h":
Usage()
sys.exit(1)
else:
Usage()
sys.exit(1)
if head == '':
print("Head file are not exist!")
Usage()
sys.exit(1)
ops = [op for op, value in opts]
if '-b' in ops and '-e' in ops:
isyrange = 1
elif '-c' in ops:
isyrange = 0
else:
print('-b and -e must be set at same time otherwise please set -c as custom')
sys.exit(1)
if isyrange:
datemin = UTCDateTime(starttime)
datemax = UTCDateTime(endtime)
event=[]
config.read(head)
eventlst = config.get("lst", "eventlst")
NAME = config.get("info", "NAME")
INST = config.get("info", "INST")
EMAIL = config.get("info", "EMAIL")
MEDIA = config.get("info", "MEDIA")
server = config.get('smtp', 'server')
passwd = config.get('smtp', 'password')
port = config.getint('smtp', 'port')
ALTERNATEMEDIA = MEDIA
if fformat.lower() == 'seed':
recipient = 'breq_fast@iris.washington.edu'
elif fformat.lower() == 'miniseed':
recipient = 'miniseed@iris.washington.edu'
else:
print('Invalid file format!')
sys.exit(1)
if isyrange:
LABEL = 'IRIS_'+datemin.strftime('%Y')+"_"+datemax.strftime('%Y')+"_"+network+"_"+station
else:
LABEL = 'IRIS_'+network+"_"+station
if not isyrange:
EVENT = open(datetimefile, 'r')
for evenum in EVENT.readlines():
evenum = evenum.strip('\n')
evenum_sp = re.split('\W|\s',evenum)
date_beg = datetime.datetime(int(evenum_sp[0]),int(evenum_sp[1]),int(evenum_sp[2]),int(evenum_sp[3]),int(evenum_sp[4]),int(evenum_sp[5]))
date_end = datetime.datetime(int(evenum_sp[6]),int(evenum_sp[7]),int(evenum_sp[8]),int(evenum_sp[9]),int(evenum_sp[10]),int(evenum_sp[11]))
event.append([date_beg.strftime('%Y %m %d %H %M %S'), date_end.strftime('%Y %m %d %H %M %S')])
else:
trange_sp = timerange.split('/')
btime = float(trange_sp[0])
etime = float(trange_sp[1])
cat = wsfetch('IRIS', starttime=datemin, endtime=datemax, minmagnitude=magmin, maxmagnitude=magmax)
for evt in cat:
if isph:
dis = distaz.distaz(stla, stlo, evt[1], evt[2]).delta
arr = mod.get_travel_times(source_depth_in_km=evt[3], distance_in_degree=dis, phase_list=[phase])
if len(arr) != 0:
arr_time = evt[0] + arr[0].time
date = arr_time + btime
dateend = arr_time + etime
else:
date = evt[0] + btime
dateend = evt[0] + etime
event.append([date.strftime('%Y %m %d %H %M %S'), dateend.strftime('%Y %m %d %H %M %S')])
if not event:
print('No events found in the range')
msg = generatemsg(NAME, INST, EMAIL, MEDIA, ALTERNATEMEDIA, LABEL)
for row in event:
msg += station+' '+network+' '+row[0]+' '+row[1]+' 1 '+chan+' '+loca+'\n'
try:
if server == '' or passwd == '':
sendmail(EMAIL, recipient, msg)
else:
sendmail(EMAIL, recipient, msg, server=server, passwd=passwd, port=port)
time.sleep(4)
print("Successful sending the mail of " + network + "." + station + " to IRIS DMC!!!")
except Exception as e:
print('ERROR in sending mail\n{}'.format(e))
sys.exit(1)
from bqmail import query
\ No newline at end of file
import math
import numpy as np
def sind(deg):
rad = math.radians(deg)
return math.sin(rad)
def cosd(deg):
rad = math.radians(deg)
return math.cos(rad)
def tand(deg):
rad = math.radians(deg)
return math.tan(rad)
def cotd(deg):
rad = math.radians(deg)
return math.cos(rad) / math.sin(rad)
def asind(x):
rad = math.asin(x)
return math.degrees(rad)
def acosd(x):
rad = math.acos(x)
return math.degrees(rad)
def atand(x):
rad = math.atan(x)
return math.degrees(rad)
def latlon_from(lat1,lon1,azimuth,gcarc_dist):
lat2 = asind ((sind (lat1) * cosd (gcarc_dist)) + (cosd (lat1) * sind (gcarc_dist) * cosd (azimuth)))
if ( cosd (gcarc_dist) >= (cosd (90 - lat1) * cosd (90 - lat2))):
lon2 = lon1 + asind (sind (gcarc_dist) * sind (azimuth) / cosd (lat2))
else:
lon2 = lon1 + asind (sind (gcarc_dist) * sind (azimuth) / cosd (lat2)) + 180
return lat2, lon2
def km2deg(kilometers):
return kilometers / 111.19
def deg2km(degree):
return degree * 111.19
class distaz:
"""
c Subroutine to calculate the Great Circle Arc distance
c between two sets of geographic coordinates
......@@ -71,67 +71,69 @@ class distaz:
and I vaguely remember a perl port. Long live distaz!
"""
def __init__(self, lat1, lon1, lat2, lon2):
def __init__(self, lat1, lon1, lat2, lon2):
self.stalat = lat1
self.stalon = lon1
self.evtlat = lat2
self.evtlon = lon2
'''
if (lat1 == lat2) and (lon1 == lon2):
self.delta = 0.0
self.az = 0.0
self.baz = 0.0
return
rad=2.*math.pi/360.0
'''
rad = 2. * math.pi / 360.0
"""
c
c scolat and ecolat are the geocentric colatitudes
c as defined by Richter (pg. 318)
c
c Earth Flattening of 1/298.257 take from Bott (pg. 3)
c
c
c scolat and ecolat are the geocentric colatitudes
c as defined by Richter (pg. 318)
c
c Earth Flattening of 1/298.257 take from Bott (pg. 3)
c
"""
sph=1.0/298.257
sph = 1.0 / 298.257
scolat=math.pi/2.0 - math.atan((1.-sph)*(1.-sph)*math.tan(lat1*rad))
ecolat=math.pi/2.0 - math.atan((1.-sph)*(1.-sph)*math.tan(lat2*rad))
slon=lon1*rad
elon=lon2*rad
scolat = math.pi / 2.0 - np.arctan((1. - sph) * (1. - sph) * np.tan(lat1 * rad))
ecolat = math.pi / 2.0 - np.arctan((1. - sph) * (1. - sph) * np.tan(lat2 * rad))
slon = lon1 * rad
elon = lon2 * rad
"""
c
c a - e are as defined by Bullen (pg. 154, Sec 10.2)
c These are defined for the pt. 1
c
"""
a=math.sin(scolat)*math.cos(slon)
b=math.sin(scolat)*math.sin(slon)
c=math.cos(scolat)
d=math.sin(slon)
e=-math.cos(slon)
g=-c*e
h=c*d
k=-math.sin(scolat)
a = np.sin(scolat) * np.cos(slon)
b = np.sin(scolat) * np.sin(slon)
c = np.cos(scolat)
d = np.sin(slon)
e = -np.cos(slon)
g = -c * e
h = c * d
k = -np.sin(scolat)
"""
c
c aa - ee are the same as a - e, except for pt. 2
c
"""
aa=math.sin(ecolat)*math.cos(elon)
bb=math.sin(ecolat)*math.sin(elon)
cc=math.cos(ecolat)
dd=math.sin(elon)
ee=-math.cos(elon)
gg=-cc*ee
hh=cc*dd
kk=-math.sin(ecolat)
aa = np.sin(ecolat) * np.cos(elon)
bb = np.sin(ecolat) * np.sin(elon)
cc = np.cos(ecolat)
dd = np.sin(elon)
ee = -np.cos(elon)
gg = -cc * ee
hh = cc * dd
kk = -np.sin(ecolat)
"""
c
c Bullen, Sec 10.2, eqn. 4
c
"""
delrad=math.acos(a*aa + b*bb + c*cc)
self.delta=delrad/rad
delrad = np.arccos(a * aa + b * bb + c * cc)
self.delta = delrad / rad
"""
c
c Bullen, Sec 10.2, eqn 7 / eqn 8
......@@ -141,13 +143,18 @@ class distaz:
c Calculate baz this way to avoid quadrant problems
c
"""
rhs1=(aa-d)*(aa-d)+(bb-e)*(bb-e)+cc*cc - 2.
rhs2=(aa-g)*(aa-g)+(bb-h)*(bb-h)+(cc-k)*(cc-k) - 2.
dbaz=math.atan2(rhs1,rhs2)
if (dbaz<0.0):
dbaz=dbaz+2*math.pi
self.baz=dbaz/rad
rhs1 = (aa - d) * (aa - d) + (bb - e) * (bb - e) + cc * cc - 2.
rhs2 = (aa - g) * (aa - g) + (bb - h) * (bb - h) + (cc - k) * (cc - k) - 2.
dbaz = np.arctan2(rhs1, rhs2)
dbaz_idx = np.where(dbaz < 0.0)[0]
if len(dbaz_idx) != 0:
if isinstance(dbaz, (int, float)):
dbaz += 2 * math.pi
else:
dbaz[dbaz_idx] += 2 * math.pi
self.baz = dbaz / rad
"""
c
c Bullen, Sec 10.2, eqn 7 / eqn 8
......@@ -155,22 +162,65 @@ class distaz:
c pt. 2 is unprimed, so this is technically the az
c
"""
rhs1=(a-dd)*(a-dd)+(b-ee)*(b-ee)+c*c - 2.
rhs2=(a-gg)*(a-gg)+(b-hh)*(b-hh)+(c-kk)*(c-kk) - 2.
daz=math.atan2(rhs1,rhs2)
if daz<0.0:
daz=daz+2*math.pi
self.az=daz/rad
rhs1 = (a - dd) * (a - dd) + (b - ee) * (b - ee) + c * c - 2.
rhs2 = (a - gg) * (a - gg) + (b - hh) * (b - hh) + (c - kk) * (c - kk) - 2.
daz = np.arctan2(rhs1, rhs2)
daz_idx = np.where(daz < 0.0)[0]
if len(daz_idx) != 0:
if isinstance(daz, (int, float)):
daz += 2 * math.pi
else:
daz[daz_idx] += 2 * math.pi
self.az = daz / rad
"""
c
c Make sure 0.0 is always 0.0, not 360.
c
"""
if (abs(self.baz-360.) < .00001):
self.baz=0.0
if (abs(self.az-360.) < .00001):
self.az=0.0
idx = np.where(np.abs(self.baz - 360.) < .00001)[0]
if len(idx) != 0:
if isinstance(self.baz, float):
self.baz = 0.0
else:
self.baz[idx] = 0.0
idx = np.where(np.abs(self.baz) < .00001)[0]
if len(idx) != 0:
if isinstance(self.baz, float):
self.baz = 0.0
else:
self.baz[idx] = 0.0
idx = np.where(np.abs(self.az - 360.) < .00001)[0]
if len(idx) != 0:
if isinstance(self.az, float):
self.az = 0.0
else:
self.az[idx] = 0.0
idx = np.where(np.abs(self.az) < .00001)[0]
if len(idx) != 0:
if isinstance(self.az, float):
self.az = 0.0
else:
self.az[idx] = 0.0
la_idx = np.where(lat1 == lat2)[0]
lo_idx = np.where(lon1 == lon2)[0]
idx = la_idx[np.where(la_idx == lo_idx)[0]]
if len(idx) != 0:
if isinstance(self.delta, float):
self.delta = 0.
else:
self.delta[idx] = 0.
if isinstance(self.az, float):
self.az = 0.
else:
self.az[idx] = 0.
if isinstance(self.baz, float):
self.baz = 0.
else:
self.baz[idx] = 0.
def getDelta(self):
return self.delta
......@@ -182,7 +232,14 @@ class distaz:
return self.baz
def degreesToKilometers(self):
return self.delta * 111.19
#distaz = DistAz(0, 0, 1,1)
#print "%f %f %f" % (distaz.getDelta(), distaz.getAz(), distaz.getBaz())
return self.delta * 111.19
# distaz = DistAz(0, 0, 1,1)
# print "%f %f %f" % (distaz.getDelta(), distaz.getAz(), distaz.getBaz())
if __name__ == '__main__':
ela = 1
elo = 1
sla = 2
slo = 1
da = distaz(ela, elo, sla, slo)
print(da.baz)
from bqmail.query import Query
from obspy import UTCDateTime
from obspy.taup import TauPyModel
from obspy.clients.iris import Client
import smtplib
from email.mime.text import MIMEText
from email.header import Header
from .distaz import distaz
model = TauPyModel()
cld = Client()
def connectsmtp(server, port, sender, password):
smtpObj = smtplib.SMTP_SSL(server, port)
smtpObj.login(sender, password)
return smtpObj
def loginmail(sender, server='localhost', password='', port=465, test_num=5):
test_it = 1
while test_it <= test_num:
if test_it > 1:
print('Try to link to {} in {} times'.format(server, test_it))
try:
smtpObj = connectsmtp(server, port, sender, password)
except:
test_it += 1
continue
else:
break
if test_it > test_num:
raise ConnectionError("Error in linking {}".format(server))
else:
return smtpObj
def sendmail(smtpObj, sender, contents,
recipient='breq_fast@iris.washington.edu'):
msg = MIMEText(contents, 'text')
msg['Subject'] = Header('batch request seismic data', 'utf-8')
msg['From'] = Header(sender)
msg['To'] = Header(recipient)
try:
smtpObj.sendmail(sender, recipient, msg.as_string())