读取数据
=============================================
从文件中读取DAS数据
------------------------------
``read函数`` 用于读取多类型的DAS数据文件,支持的所有数据格式见下表:
+--------+------+------+----------------------------------------------------------------------------------------------------------------------------------+
| 格式 | 读取 | 写入 | 数据来源 |
+========+======+======+==================================================================================================================================+
| HDF5 | √ | √ | OptaSense、Silixa、Febus、AP Sensing、ASN、Sintela、Aragón Photonics、T8 Sensors、智地感知、AI4EPS、INGV、JAMSTEC和FORESEE数据集 |
+--------+------+------+----------------------------------------------------------------------------------------------------------------------------------+
| TDMS | √ | √ | Silixa和中国科学院半导体研究所 |
+--------+------+------+----------------------------------------------------------------------------------------------------------------------------------+
| SEG-Y | √ | √ | OptaSens和Silixa |
+--------+------+------+----------------------------------------------------------------------------------------------------------------------------------+
| PICKLE | √ | √ | 以二进制保存的daspy.Section或numpy.ndarray类实例 |
+--------+------+------+----------------------------------------------------------------------------------------------------------------------------------+
| NPY | √ | √ | 以二进制保存的numpy.ndarray类实例 |
+--------+------+------+----------------------------------------------------------------------------------------------------------------------------------+
在默认情况下,函数将输出包含数据和元数据的 ``Section`` 类实例:
>>> from daspy import read
>>> sec = read() # 读取示例波形
设置 ``ch1`` 或/和 ``ch2`` 可限制读取的信道范围。可设置 ``output_type='array'`` 以输出数据为 ``numpy.array`` 格式,此时将以 ``dict`` 格式输出元数据:
>>> data, metadata = read(output_type='array', ch1=2700, ch2=2800) # 设置信道范围
在仅需要元数据时,可以设置 ``headonly=True`` 以节省读取时间,此时返回的数据将是一个与原始数据大小相同的全零数组。当文件名的后缀不能准确指示数据格式时,请设置参数 ``ftype`` 为'h5'、'tdms'、'segy'等
自定义读取函数
------------------------------
如果DASPy不支持读取你的数据,可以自定义一个读取函数。尽管函数仅输入文件名 ``fname`` 参数即可完成读取,但我们推荐用户允许函数输入 ``headonly`` 、 ``ch1`` 和 ``ch2`` 参数以保证完全的兼容性:
>>> import h5py
>>> from daspy import DASDateTime
>>>
>>> def read_new_h5(fname, headonly=False, **kwargs):
>>> with h5py.File(fname, 'r') as h5_file:
>>> start_channel = h5_file.attrs['channel_start']
>>> end_channel = h5_file.attrs['channel_end']
>>> ch1 = kwargs.pop('ch1', start_channel)
>>> ch2 = kwargs.pop('ch2', end_channel)
>>> dch = kwargs.pop('dch', 1)
>>> if headonly:
>>> data = np.zeros_like(h5_file['raw_data'])
>>> else:
>>> data = h5_file['raw_data'][ch1-start_channel:ch2-start_channel:dch]
>>> dx = h5_file.attrs['channel spacing m']
>>> metadata = {'dx': dx, 'fs': h5_file.attrs['sampling rate Hz'],
>>> 'gauge_length': h5_file.attrs['GL m'],
>>> 'start_channel': start_channel + ch1,
>>> 'start_distance': (start_channel + ch1) * dx,
>>> 'start_time': DASDateTime.strptime(h5_file.attrs['starttime'], '%Y-%m-%dT%H:%M:%S.%f'),
>>> 'scale': h5_file.attrs['scale factor to strain']}
>>> return data, metadata
``metadata`` 中所必须的关键字是 ``dx`` 和 ``fs`` (若文件元数据中不包含可暂时设置为 ``None`` ),其余参数可按需读取。之后该读取函数便可作为 ``ftype`` 参数输入 ``read`` 函数或用于构建 ``Collection`` 类(详见 :doc:`处理连续数据` ):
>>> sec = read(fname, ftype=read_new_h5, headonly=True)
>>> coll = Collection('../data/*.h5', ftype=read_new_h5)
从其他包的类实例转换
------------------------------
DASPy支持数据从 `ObsPy `_ 、 `DASCore `_ 和 `lightguide `_ 向DASPy流转:
>>> from daspy import Section
>>> sec_obspy = Section.from_obspy_stream(st) # 将obspy.core.stream.Stream实例转换为daspy.section实例
>>> sec_dascore = Section.from_dascore_patch(patch) # 将dascore.core.patch.Patch实例转换为daspy.section实例
>>> sec_lightguide = Section.from_lightguide_blast(blast) # 将lightguide.blast.Blast实例转换为daspy.section实例
Section类
------------------------------
DASPy的功能可以通过调用函数(面向过程)实现,也可以通过类的方法(面向对象)实现。我们推荐面向对象的方式,对DAS数据体的类 ``Section`` 类进行操作。DASPy中的所有面向用户的函数和类的方法都同名,且为名词(如 ``normalization`` 和 ``Section.normalization`` )。
``Section`` 类的属性保存了该DAS数据以及元数据,其中波形数据 ``data`` ,道间距 ``dx`` 、采样率 ``fs`` 是必须设置的(但可以暂时设置为 ``None`` );起始道号 ``start_channel`` 、起始距离 ``start_distance`` 、起始时间 ``start_time`` 在未设置时默认为0;事件起始时间 ``origin_time`` 、标距长度 ``gauge_length`` 、数据类型 ``data_type`` 、数据尺度 ``scale`` 、阵列几何 ``geometry`` 、拐点道号 ``turning_channels`` 、其他头段信息 ``headers`` 可不设置。
访问元数据:
>>> print(sec)
data: shape(500, 5000)
dx: 1 m
fs: 100.0 Hz
start_channel: 2500
start_distance: 2520 m
start_time: 2016-03-21 07:37:30.532309+00:00
origin_time: 2016-03-21 07:37:10.535000+00:00
data_type: strain rate
此外,数据尺寸 ``shape`` ,信道数量 ``nch`` 采样点数 ``nt`` ,结束道号 ``end_channel`` ,结束距离 ``end_distance`` 和结束时间 ``end_time`` 将自动计算,可作为属性被调用。
DASDateTime类
------------------------------
DASPy创建了 ``DASDateTime`` 类,以表示数据的时间信息,包括起始时间 ``start_time`` 、结束时间 ``end_time`` 和事件起始时间 ``origin_time`` 。
``DASDateTime`` 是 ``datetime.DateTime`` 类的子类,继承了 ``datetime.DateTime`` 的所有方法:
>>> from daspy import DASDateTime
>>> DASDateTime.strptime('2021-03-19T1:52:23', '%Y-%m-%dT%H:%M:%S')
DASDateTime(2021, 03, 19, 1, 52, 23)
DASPy内置有当地时区 ``local_tz`` 和utc时区 ``utc`` 用于指定时区:
>>> from daspy import utc, local_tz
>>> DASDateTime.fromtimestamp(1616089943, tz=utc)
DASDateTime(2021, 3, 18, 17, 52, 23, tzinfo=datetime.timezone.utc)
如需使用其他时区可以用 ``datetime.timezone(datetime.timedelta(hours=h))`` 等方法建立。
可以使用 ``local`` 、 ``utc`` 和 ``remove_tz`` 方法转换时区或去除时区信息:
>>> time = DASDateTime.strptime('2021-03-19T1:52:23Z', '%Y-%m-%dT%H:%M:%S%z')
>>> time.local()
DASDateTime(2021, 3, 19, 9, 52, 23, tzinfo=datetime.timezone(datetime.timedelta(seconds=28800)))
>>> time.remove_tz()
DASDateTime(2021, 3, 19, 1, 52, 23)
除了父类本身支持的 ``datetime.datetime`` 和 ``datetime.timedelta`` 之间的加减运算以外, ``DASDateTime`` 还支持输入数字和可迭代对象 ``Iterable`` 来计算加和减,所有时间差都以秒(s)为单位表示,并且会自动处理时区未指定的问题:
>>> DASDateTime(2021, 3, 24, 14, 28, 0, 0) + 100
DASDateTime(2021, 3, 24, 14, 29, 40)
>>> DASDateTime(2021, 3, 24, 14, 28, 0, 0) + [10, 20, 30]
[DASDateTime(2021, 3, 24, 14, 28, 10), DASDateTime(2021, 3, 24, 14, 28, 20), DASDateTime(2021, 3, 24, 14, 28, 30)]
>>> DASDateTime(2021, 3, 24, 14, 28, 0, 0) - 100
DASDateTime(2021, 3, 24, 14, 26, 20)
>>> DASDateTime(2021, 3, 24, 14, 28, 0, 0) - DASDateTime(2021, 3, 19, 1, 52, 23)
477337.0
有需要时可以将 ``DASDateTime`` 实例转换为父类 ``datetime.datetime`` 实例:
>>> DASDateTime(2021, 3, 19, 1, 52, 23).convert_to_datetime()
datetime.datetime(2021, 3, 19, 1, 52, 23)