- Clone an Existing Dataless SEED File （复制现有的无数据SEED文件）
以下代码展示了如何复制一个现有的无数据SEED文件（dataless.seed.BW_RNON）然后使用他作为一个模板为新的台站建立DatalessSEED文件。
首先，我们要做好必要的输入接口然后读取现有的DatalessSEED文件：
>>>from obspy import UTCDateTime
>>>from obspy.io.xseed import Parser
>>>p = Parser("https://examples.obspy.org/dataless.seed.BW_RNON")
>>>blk = p.blockettes
现在我们可以调整仅在DatalessSEED文件开头中出现的信息：
>>>blk[50][0].network_code = 'BW'
>>>blk[50][0].station_call_letters = 'RMOA'
>>>blk[50][0].site_name = "Moar Alm, Bavaria, BW-Net"
>>>blk[50][0].latitude = 47.761658
>>>blk[50][0].longitude = 12.864466
>>>blk[50][0].elevation = 815.0
>>>blk[50][0].start_effective_date = UTCDateTime("2006-07-18T00:00:00.000000Z")
>>>blk[50][0].end_effective_date = ""
>>>blk[33][1].abbreviation_description = "Lennartz LE-3D/1 seismometer"
之后我们还要调整所有三个通道分量的信息：
>>>mult = len(blk[58])/3 >>>for i, cha in enumerate(['Z', 'N', 'E']): blk[52][i].channel_identifier = 'EH%s' % cha blk[52][i].location_identifier = '' blk[52][i].latitude = blk[50][0].latitude blk[52][i].longitude = blk[50][0].longitude blk[52][i].elevation = blk[50][0].elevation blk[52][i].start_date = blk[50][0].start_effective_date blk[52][i].end_date = blk[50][0].end_effective_date blk[53][i].number_of_complex_poles = 3 blk[53][i].real_pole = [-4.444, -4.444, -1.083] blk[53][i].imaginary_pole = [+4.444, -4.444, +0.0] blk[53][i].real_pole_error = [0, 0, 0] blk[53][i].imaginary_pole_error = [0, 0, 0] blk[53][i].number_of_complex_zeros = 3 blk[53][i].real_zero = [0.0, 0.0, 0.0] blk[53][i].imaginary_zero = [0.0, 0.0, 0.0] blk[53][i].real_zero_error = [0, 0, 0] blk[53][i].imaginary_zero_error = [0, 0, 0] blk[53][i].A0_normalization_factor = 1.0 blk[53][i].normalization_frequency = 3.0 # stage sequence number 1, seismometer gain blk[58][i*mult].sensitivity_gain = 400.0 # stage sequence number 2, digitizer gain blk[58][i*mult+1].sensitivity_gain = 1677850.0 # stage sequence number 0, overall sensitivity blk[58][(i+1)*mult-1].sensitivity_gain = 671140000.0
注意：这个例子中没有设置FIR系数。
最后，我们可以将调整后的DatalessSEED卷写入一个新文件：
>>>p.write_seed("dataless.seed.BW_RMOA")
- Export Seismograms to MATLAB（导出数据到MATLAB）
下面的例子展示了如何使用python读取波形文件并保存每条Trace的结果Stream对象到一个MATLAB的MAT格式文件中。这些数据可以通过MATLAB使用load函数读取：
from scipy.io import savemat
import obspy
st = obspy.read("https://examples.obspy.org/BW.BGLD..EH.D.2010.037")
for i, tr in enumerate(st):
mdict = {k: str(v) for k, v in tr.stats.iteritems()}
mdict['data'] = tr.data
savemat("data-%d.mat" % i, mdict)
-
Export Seismograms to ASCII（导出数据为ASCII文件）
- Built-in 内置格式
你可以使用Obspy的write（）函数直接输出Stream对象的波形数据为任何ASCII格式。
>>>from obspy.core import read
>>>stream = read('https://examples.obspy.org/RJOB20090824.ehz')
>>>stream.write('outfile.ascii', format='SLIST')
下面是现在支持的ASCII格式：
- SLIST——使用下列简单列表代表的ASCII时间序列格式：
TIMESERIES BW_RJOB__EHZ_D, 6001 samples, 200 sps, 2009-08-24T00:20:03.000000, SLIST, INTEGER, 288 300 292 285 265 287 279 250 278 278 268 258
- TSPAIR——时间—样本对ASCII格式
TIMESERIES BW_RJOB__EHZ_D, 6001 samples, 200 sps, 2009-08-24T00:20:03.000000, TSPAIR, INTEGER, 2009-08-24T00:20:03.000000 288 2009-08-24T00:20:03.005000 300 2009-08-24T00:20:03.010000 292 2009-08-24T00:20:03.015000 285 2009-08-24T00:20:03.020000 265 2009-08-24T00:20:03.025000 287
- SH_ASC——Seismic Handler支持的ASCII格式
DELTA: 5.000000e-03 LENGTH: 6001 START: 24-AUG-2009_00:20:03.000 COMP: Z CHAN1: E CHAN2: H STATION: RJOB CALIB: 1.000000e+00 2.880000e+02 3.000000e+02 2.920000e+02 2.850000e+02 2.650000e+02 2.870000e+02 2.790000e+02 2.500000e+02
下面的脚本示例如何使用自定义头转换每条Trace或者地震记录文件到ASCII文件。波形数据会乘上给定的calibration校准因子，然后使用Numpy的savetxt（）函数保存。
"""
USAGE: export_seismograms_to_ascii.py in_file out_file calibration
"""
from __future__ import print_function
import sys
import numpy as np
import obspy
try:
in_file = sys.argv[1]
out_file = sys.argv[2]
calibration = float(sys.argv[3])
except:
print(__doc__)
raise
st = obspy.read(in_file)
for i, tr in enumerate(st):
f = open("%s_%d" % (out_file, i), "w")
f.write("# STATION %s\n" % (tr.stats.station))
f.write("# CHANNEL %s\n" % (tr.stats.channel))
f.write("# START_TIME %s\n" % (str(tr.stats.starttime)))
f.write("# SAMP_FREQ %f\n" % (tr.stats.sampling_rate))
f.write("# NDAT %d\n" % (tr.stats.npts))
np.savetxt(f, tr.data * calibration, fmt="%f")
f.close()
- Anything to MiniSEED（转换任意文件格式为MiniSEED）
下列程序展示如何转换任何数据到MiniSEED格式。例子中将几行气象站的输出写入MiniSEED文件。准确的元信息（starttime、sampling_rate、station名称等）也被编码。（注意：必须是MiniSEED标准允许的编码）。如果你想发送日志信息或者通过SeedLink协议输出气象站或其他任何东西，那么转换任意ASCII为MiniSEED极有帮助。
from __future__ import print_function
import numpy as np
from obspy import UTCDateTime, read, Trace, Stream
weather = """
00.0000 0.0 ??? 4.7 97.7 1015.0 0.0 010308 000000
00.0002 0.0 ??? 4.7 97.7 1015.0 0.0 010308 000001
00.0005 0.0 ??? 4.7 97.7 1015.0 0.0 010308 000002
00.0008 0.0 ??? 4.7 97.7 1015.4 0.0 010308 000003
00.0011 0.0 ??? 4.7 97.7 1015.0 0.0 010308 000004
00.0013 0.0 ??? 4.7 97.7 1015.0 0.0 010308 000005
00.0016 0.0 ??? 4.7 97.7 1015.0 0.0 010308 000006
00.0019 0.0 ??? 4.7 97.7 1015.0 0.0 010308 000007
"""
# Convert to NumPy character array
data = np.fromstring(weather, dtype='|S1')
# Fill header attributes
stats = {'network': 'BW', 'station': 'RJOB', 'location': '',
'channel': 'WLZ', 'npts': len(data), 'sampling_rate': 0.1,
'mseed': {'dataquality': 'D'}}
# set current time
stats['starttime'] = UTCDateTime()
st = Stream([Trace(data=data, header=stats)])
# write as ASCII file (encoding=0)
st.write("weather.mseed", format='MSEED', encoding=0, reclen=256)
# Show that it worked, convert NumPy character array back to string
st1 = read("weather.mseed")
print(st1[0].data.tostring())
- Beachball Plot（绘制沙滩球图）
下面展示如何创建震源机制的图形表示：
from obspy.imaging.beachball import beachball mt = [0.91, -0.89, -0.02, 1.78, -1.55, 0.47] beachball(mt, size=200, linewidth=2, facecolor='b') mt2 = [150, 87, 1] beachball(mt2, size=200, linewidth=2, facecolor='r') mt3 = [-2.39, 1.04, 1.35, 0.57, -2.94, -0.94] beachball(mt3, size=200, linewidth=2, facecolor='g')
简单的Basemap绘制可以通过Inventory或Catalog对象的内置方法执行：Inventory.plot（），Catalog.plot()。
from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
from obspy import read_inventory, read_events
# Set up a custom basemap, example is taken from basemap users' manual
fig, ax = plt.subplots()
# setup albers equal area conic basemap
# lat_1 is first standard parallel.
# lat_2 is second standard parallel.
# lon_0, lat_0 is central point.
m = Basemap(width=8000000, height=7000000,
resolution='c', projection='aea',
lat_1=40., lat_2=60, lon_0=35, lat_0=50, ax=ax)
m.drawcoastlines()
m.drawcountries()
m.fillcontinents(color='wheat', lake_color='skyblue')
# draw parallels and meridians.
m.drawparallels(np.arange(-80., 81., 20.))
m.drawmeridians(np.arange(-180., 181., 20.))
m.drawmapboundary(fill_color='skyblue')
ax.set_title("Albers Equal Area Projection")
# we need to attach the basemap object to the figure, so that obspy knows about
# it and reuses it
fig.bmap = m
# now let's plot some data on the custom basemap:
inv = read_inventory()
inv.plot(fig=fig, show=False)
cat = read_events()
cat.plot(fig=fig, show=False, title="", colorbar=False)
plt.show()
-
- 确定地点的带有Beachball的Basemap plot
下面的例子展示如何在basemap图中绘制一些台站的beachball图。该例需要安装basemap包。演示的SRTM文件可以从here下载。我们的SRTM前几行是这样的：
ncols 400 nrows 200 xllcorner 12°40'E yllcorner 47°40'N xurcorner 13°00'E yurcorner 47°50'N cellsize 0.00083333333333333 NODATA_value -9999 682 681 685 690 691 689 678 670 675 680 681 679 675 671 674 680 679 679 675 671 668 664 659 660 656 655 662 666 660 659 659 658 ....
import gzip
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from obspy.imaging.beachball import beach
# read in topo data (on a regular lat/lon grid)
# (SRTM data from: http://srtm.csi.cgiar.org/)
with gzip.open("srtm_1240-1300E_4740-4750N.asc.gz") as fp:
srtm = np.loadtxt(fp, skiprows=8)
# origin of data grid as stated in SRTM data file header
# create arrays with all lon/lat values from min to max and
lats = np.linspace(47.8333, 47.6666, srtm.shape[0])
lons = np.linspace(12.6666, 13.0000, srtm.shape[1])
# create Basemap instance with Mercator projection
# we want a slightly smaller region than covered by our SRTM data
m = Basemap(projection='merc', lon_0=13, lat_0=48, resolution="h",
llcrnrlon=12.75, llcrnrlat=47.69, urcrnrlon=12.95, urcrnrlat=47.81)
# create grids and compute map projection coordinates for lon/lat grid
x, y = m(*np.meshgrid(lons, lats))
# Make contour plot
cs = m.contour(x, y, srtm, 40, colors="k", lw=0.5, alpha=0.3)
m.drawcountries(color="red", linewidth=1)
# Draw a lon/lat grid (20 lines for an interval of one degree)
m.drawparallels(np.linspace(47, 48, 21), labels=[1, 1, 0, 0], fmt="%.2f",
dashes=[2, 2])
m.drawmeridians(np.linspace(12, 13, 21), labels=[0, 0, 1, 1], fmt="%.2f",
dashes=[2, 2])
# Plot station positions and names into the map
# again we have to compute the projection of our lon/lat values
lats = [47.761659, 47.7405, 47.755100, 47.737167]
lons = [12.864466, 12.8671, 12.849660, 12.795714]
names = [" RMOA", " RNON", " RTSH", " RJOB"]
x, y = m(lons, lats)
m.scatter(x, y, 200, color="r", marker="v", edgecolor="k", zorder=3)
for i in range(len(names)):
plt.text(x[i], y[i], names[i], va="top", family="monospace", weight="bold")
# Add beachballs for two events
lats = [47.751602, 47.75577]
lons = [12.866492, 12.893850]
x, y = m(lons, lats)
# Two focal mechanisms for beachball routine, specified as [strike, dip, rake]
focmecs = [[80, 50, 80], [85, 30, 90]]
ax = plt.gca()
for i in range(len(focmecs)):
b = beach(focmecs[i], xy=(x[i], y[i]), width=1000, linewidth=1)
b.set_zorder(10)
ax.add_collection(b)
plt.show()
几个Note：
- GDAL包允许直接读取GeoTiff到numpy ndarray。
>>> geo = gdal.Open("file.geotiff")
>>> x = geo.ReadAsArray()
- 可从ASTER得到Geo Tiff高度数据。
- 可添加阴影/照明。查看plotmap_shaded.py例程获取更多信息。
- 带有beachball的全球Basemap
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.basemap import Basemap from obspy.imaging.beachball import beach m = Basemap(projection='cyl', lon_0=142.36929, lat_0=38.3215, resolution='c') m.drawcoastlines() m.fillcontinents() m.drawparallels(np.arange(-90., 120., 30.)) m.drawmeridians(np.arange(0., 420., 60.)) m.drawmapboundary() x, y = m(142.36929, 38.3215) focmecs = [0.136, -0.591, 0.455, -0.396, 0.046, -0.615] ax = plt.gca() b = beach(focmecs, xy=(x, y), width=10, linewidth=1, alpha=0.85) b.set_zorder(10) ax.add_collection(b) plt.show()
