obspy中文教程（四）

目录

Clone an Existing Dataless SEED File（复制现有的无数据SEED文件）

Export Seismograms to MATLAB（导出数据到MATLAB）

Anything to MiniSEED（转换任意文件格式为MiniSEED）

Beachball Plot（绘制沙滩球图）

---

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 Plots

设置自定义投影的Basemap plot

简单的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：

1. GDAL包允许直接读取GeoTiff到numpy ndarray。

>>> geo = gdal.Open("file.geotiff")  

>>> x = geo.ReadAsArray() 

1. 可从ASTER得到Geo Tiff高度数据。
2. 可添加阴影/照明。查看plotmap_shaded.py例程获取更多信息。
   1. 带有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()