提交 31e7d68b 编辑于 作者: Mijian Xu's avatar Mijian Xu 😷
浏览文件

1.1

上级
此差异已折叠。
bqmail
===========
This is a Python script to request waveform data of events from IRIS
Usage:
python bqmail.py -Nnetwork -Sstation -Yyear1/month1/day1/year2/month2/day2 head.cfg
#!/usr/bin/env python
#
#Author: Mijian Xu at NJU
#
#Revision History:
# 2014/11/06
# 2015/01/05
# 2015/02/11
#
def Usage():
print('Usage:')
print('python bqmail.py -Nnetwork -Sstation -Yyear1/month1/day1/year2/month2/day2 -Cdatetimefile head.cfg')
import datetime
import os, re
import sys, getopt
try:
import configparser
config = configparser.ConfigParser()
except:
import ConfigParser
config = ConfigParser.ConfigParser()
try:
opts,args = getopt.getopt(sys.argv[1:], "hN:S:C:")
except:
print('arguments are not found!')
Usage()
sys.exit(1)
iscustom = 0
isyrange = 0
for op, value in opts:
if op == "-N":
network = value
elif op == "-S":
station = value
elif op == "-Y":
yrange = value
isyrange = 1
elif op == "-C":
datetimefile = value
iscustom = 1
elif op == "-h":
Usage()
sys.exit(1)
else:
Usage()
sys.exit(1)
head = []
for o in sys.argv[1:]:
if os.path.isfile(o):
head = o
break
if head == []:
print("Aruments or head file are not exist!")
Usage()
sys.exit(1)
if isyrange:
y_split = yrange.split('/')
year1 = int(y_split[0])
mon1 = int(y_split[1])
day1 = int(y_split[2])
year2 = int(y_split[3])
mon2 = int(y_split[4])
day2 = int(y_split[5])
datemin=datetime.datetime(year1,mon1,day1)
datemax=datetime.datetime(year2,mon2,day2)
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")
ALTERNATEMEDIA = config.get("info","ALTERNATEMEDIA")
if isyrange:
LABEL = 'IRIS_'+str(year1)+"_"+str(year2)+"_"+network+"_"+station
else:
LABEL = 'IRIS_'+network+"_"+station
if iscustom:
EVENT = open(datetimefile,'r')
for evenum in EVENT:
evenum = evenum.strip('\n')
evenum_sp = re.split(',|\s',evenum)
# print(evenum_sp)
event.append(evenum_sp)
else:
EVENT = open(eventlst,'r+')
for evenum in EVENT:
evenum_split = evenum.split()
year=int(evenum_split[0])
mon=int(evenum_split[1])
day=int(evenum_split[2])
jjj=int(evenum_split[3])
hour=int(evenum_split[4])
min=int(evenum_split[5])
sec=float(evenum_split[6])
lat=float(evenum_split[7])
lon=float(evenum_split[8])
dep=float(evenum_split[9])
mw=float(evenum_split[10])
date = datetime.datetime(year,mon,day,hour,min)
dt = datetime.timedelta(hours=1)
dateend = date + dt
if datemin <= date <= datemax:
event.append([date.strftime('%Y'),date.strftime('%m'),date.strftime('%d'),date.strftime('%H'),date.strftime('%M'),dateend.strftime('%Y'),dateend.strftime('%m'),dateend.strftime('%d'),dateend.strftime('%H'),dateend.strftime('%M')])
OUT = open(network+'_'+station+'.bq','w+')
OUT.write('.NAME '+NAME+'\n')
OUT.write('.INST '+INST+'\n')
OUT.write('.MAIL\n')
OUT.write('.EMAIL '+EMAIL+'\n')
OUT.write('.PHONE\n')
OUT.write('.FAX\n')
OUT.write('.MEDIA '+MEDIA+'\n')
OUT.write('.ALTERNATE MEDIA '+ALTERNATEMEDIA+'\n')
OUT.write('.ALTERNATE MEDIA '+ALTERNATEMEDIA+'\n')
OUT.write('.LABEL '+LABEL+'\n')
OUT.write('.END\n')
if not iscustom:
for row in event:
OUT.write(station+' '+network+' '+row[0]+' '+row[1]+' '+row[2]+' '+row[3]+' '+row[4]+' 00.0 '+row[5]+' '+row[6]+' '+row[7]+' '+row[8]+' '+row[9]+' 00.0 1 BH?\n')
else:
for row in event:
OUT.write(station+' '+network+' '+row[0]+' '+row[1]+' '+row[2]+' '+row[3]+' '+row[4]+' '+row[5]+' '+row[6]+' '+row[7]+' '+row[8]+' '+row[9]+' '+row[10]+' '+row[11]+' 1 BH?\n')
os.system('mail breq_fast@iris.washington.edu <'+network+'_'+station+'.bq')
print("Successful sending the mail of "+network+"."+station+" to IRIS DMC!!!")
os.system('rm '+network+'_'+station+'.bq')
import math
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
c
c Equations take from Bullen, pages 154, 155
c
c T. Owens, September 19, 1991
c Sept. 25 -- fixed az and baz calculations
c
P. Crotwell, Setember 27, 1995
Converted to c to fix annoying problem of fortran giving wrong
answers if the input doesn't contain a decimal point.
H. P. Crotwell, September 18, 1997
Java version for direct use in java programs.
*
* C. Groves, May 4, 2004
* Added enough convenience constructors to choke a horse and made public double
* values use accessors so we can use this class as an immutable
H.P. Crotwell, May 31, 2006
Port to python, thus adding to the great list of languages to which
distaz has been ported from the origin fortran: C, Tcl, Java and now python
and I vaguely remember a perl port. Long live distaz!
"""
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
"""
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
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
"""
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)
"""
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)
"""
c
c Bullen, Sec 10.2, eqn. 4
c
"""
delrad=math.acos(a*aa + b*bb + c*cc)
self.delta=delrad/rad
"""
c
c Bullen, Sec 10.2, eqn 7 / eqn 8
c
c pt. 1 is unprimed, so this is technically the baz
c
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
"""
c
c Bullen, Sec 10.2, eqn 7 / eqn 8
c
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
"""
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
def getDelta(self):
return self.delta
def getAz(self):
return self.az
def getBaz(self):
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())
#!/usr/bin/env python
#
# This is a example of bqmail for request data from network YA.
import os
network = 'YA'
for sta in os.popen('matchsta.py -N'+network+' -R95/110/20/40'):
sta_sp = sta.split()
staname = sta_sp[0]
date_begin = sta_sp[4]
date_end = sta_sp[5]
os.system('bqmail -NYA -S'+staname+' -Y'+date_begin+'/'+date_end+' head.cfg')
[lst]
eventlst = ./EventCMT.dat
[info]
NAME = huihuang
INST = NJU
EMAIL = oldyellow9451@163.com
MEDIA = Electronic (FTP)
ALTERNATEMEDIA = Electronic (FTP)
[lst]
eventlst = /Users/xumj/Documents/GitHub/bqmail/EventCMT.dat
[info]
NAME = MijianXu
INST = NJU
EMAIL = gomijianxu@163.com
MEDIA = Electronic (FTP)
ALTERNATEMEDIA = Electronic (FTP)
#!/usr/bin/env python
#
#
#
import urllib.request as rq
import os
import re
import sys, getopt
def Usage():
print('Usage:')
print('python matchsta.py -Nnetwork -Rlon1/lon2/lat1/lat2 -K > network.lst')
try:
opts,args = getopt.getopt(sys.argv[1:], "hN:R:K")
except:
print('arguments are not found!')
Usage()
sys.exit(1)
islalo = 0
ismap = 0
for op, value in opts:
if op == "-N":
network = value
elif op == "-R":
lat_lon = value
islalo = 1
elif op == "-K":
ismap = 1
elif op == "-h":
Usage()
sys.exit(1)
else:
Usage()
sys.exit(1)
if islalo:
lalo_split = lat_lon.split('/')
lon1 = float(lalo_split[0])
lon2 = float(lalo_split[1])
lat1 = float(lalo_split[2])
lat2 = float(lalo_split[3])
url = 'http://ds.iris.edu/mda/'+network
response = rq.urlopen(url)
html = str(response.read())
find_re = re.compile(r'TITLE.+?</A></TD>.+?</TR>',re.DOTALL)
for info in find_re.findall(html):
find_sta = re.compile('\w+</A></TD>.+?</TR>',re.DOTALL)
sta_info = find_sta.findall(info)
if sta_info == []:
continue
sta_info = sta_info[0]
info_s1 = sta_info.split('<')
staname = info_s1[0]
info_s2 = re.split('</TD>.+?>',sta_info)
stlat = float(info_s2[2])
stlon = float(info_s2[3])
stel = float(info_s2[4])
yrange1 = info_s2[5]
yrange2 = info_s2[6]
if islalo:
if lat1<=stlat<=lat2 and lon1<=stlon<=lon2:
print(staname+' '+str(stlat)+' '+str(stlon)+' '+info_s2[4]+' '+yrange1+' '+yrange2)
else:
print(staname+' '+str(stlat)+' '+str(stlon)+' '+info_s2[4]+' '+yrange1+' '+yrange2)
if ismap:
google = open(network+'.kml','w+')
google.write('<?xml version="1.0" encoding="UTF-8"?><kml xmlns="http://www.google.com/earth/kml/2.0"><NetworkLink><name>Selected stations</name><description>'+network+'</description><Link><href>http://www.iris.edu/cgi-bin/kmlstationinfo/'+network+'</href><refreshMode>onInterval</refreshMode><refreshInterval>86400</refreshInterval></Link></NetworkLink></kml>')
#!/usr/bin/env python
#
#
#
import distaz, math
import urllib.request as rq
import os
import re
import sys, getopt
import glob
def Usage():
print('Usage: python searchlalo.py -Rlon1/lon2/lat1/lat2 -Dlat/lon/dis1/dis2 -K')
print('-R -- Search range.')
print('-D -- Search by distance.')
print('-K -- Output Google Earth kml file.')
try:
opts,args = getopt.getopt(sys.argv[1:], "hR:D:KO:")
except:
print('arguments are not found!')
Usage()
sys.exit(1)
iskml = 0
islalo = 0
for op, value in opts:
if op == "-R":
lat_lon = value
islalo = 1
elif op == "-K":
iskml = 1
elif op == "-D":
lat_lon = value
elif op == "-h":
Usage()
sys.exit(1)
else:
Usage()
sys.exit(1)
lat_lon_split = lat_lon.split('/')
if islalo:
lon1 = lat_lon_split[0]
lon2 = lat_lon_split[1]
lat1 = lat_lon_split[2]
lat2 = lat_lon_split[3]
lalo = lon1+'_'+lon2+'_'+lat1+'_'+lat2
else:
lon = lat_lon_split[0]
lat = lat_lon_split[1]
dis1 = float(lat_lon_split[2])
dis2 = float(lat_lon_split[3])
# print(lon,lat,dis)
# [lat1,lon1] = distaz.latlon_from(float(lat),float(lon),225,dis*math.sqrt(2))
# [lat2,lon2] = distaz.latlon_from(float(lat),float(lon),45,dis*math.sqrt(2))
# print(lon1,lon2,lat1,lat2)
lon1 = str(0)
lat1 = str(-90)
lon2 = str(0)
lat2 = str(90)
lalo = lon+'_'+lat+'_'+lat_lon_split[2]
url = 'http://ds.iris.edu/cgi-bin/xmlstationinfo?minlat='+lat1+'&maxlat='+lat2+'&minlon='+lon1+'&maxlon='+lon2
response = rq.urlopen(url)
html = str(response.read())
find_re = re.compile(r'<station\s.+?"\s/>',re.DOTALL)
for info in find_re.findall(html):
sta_info = re.split('\w+="|"\s+?\w+="|"\s/>',info)
if sta_info == []:
continue
network = sta_info[1]
staname = sta_info[2]
stlat = sta_info[4]
stlon = sta_info[5]
if sta_info[-2] == 'No archive data':
continue
yrange1 = sta_info[-5]
yrange2 = sta_info[-4]
if not islalo:
delta = distaz.distaz(float(lat),float(lon),float(stlat),float(stlon))
if dis1 < delta.delta < dis2:
print(network+' '+staname+' '+stlat+' '+stlon+' '+yrange1+' '+yrange2)
else:
print(network+' '+staname+' '+stlat+' '+stlon+' '+yrange1+' '+yrange2)
if iskml:
google = open('Station_'+lalo+'.kml','w+')
google.write('<?xml version="1.0" encoding="UTF-8"?><kml xmlns="http://www.google.com/earth/kml/2.0"><NetworkLink><name>Selected stations</name><description>Station List</description><Link><href>http://www.iris.edu/cgi-bin/kmlstationinfo?minlat='+lat1+'&amp;maxlat='+lat2+'&amp;minlon='+lon1+'&amp;maxlon='+lon2+'&amp;kmz=1</href><refreshMode>onInterval</refreshMode><refreshInterval>86400</refreshInterval></Link></NetworkLink></kml>')
支持 Markdown
0% .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑!
想要评论请 注册