读取数据 ============================================= 从文件中读取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)