GeophyAI

Python与地球物理数据处理

0%

Madagascrar之RSF文件输入与输出

发表于 分类于 Valine：
简要介绍了如何从RSF文件中读取数据、参数以及如何将数据和参数格式化输出到RSF文件中。

Madagascrar中RSF文件读取与写入

本文主要介绍如何在Python脚本中读取由执行scons命令所生成的Madagascrar的.rsf文件，以及如何将numpy数据以及相关参数写入.rsf文件。

在Python中读取RSF文件

这里我们将RSF文件读取封装在一个简单的函数中，该函数使用了m8r中的Input模块来完成对数据和参数的读取，我们只需要将文件名输入到该函数中，即可获得与文件名同名的方法，并可分别通过name.data和name.pars来访问name.rsf中的数据和参数。
函数原型：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
class m8r_read():
    def __init__(self, file_path):
        self.__file_path__ = file_path
    def __read__(self, ):
        self.Input = m8r.Input(self.__file_path__)
        self.pars = self.Input.file.pars
        self.data = self.Input.read()
    def __plot__(self, cmap = plt.cm.gray):
        assert self.data.ndim>1
        plt.figure()
        vmin, vmax = np.percentile(self.data, [2, 98])
        plt.imshow(self.data.T,
                   vmin=vmin, vmax=vmax, 
                   cmap = cmap,
                   aspect = 'auto')
        plt.show()
class warpper():
    m8r.no_swig()
    def __init__(self, root_path):
        """
        root_path: The root path of data.
        """
        self.__root_path__ = root_path
    def __set_attr_by_name__(self, name, postfix = 'rsf'):
        """
        name: The name of rsf file that need to be loaded.
        """
        path = os.path.join(self.__root_path__, name+'.'+postfix)
        self.__setattr__(name, m8r_read(path))
        self.__getattribute__(name).__read__()

我们以rsf/su/rsflab16/seismic.rsf为例，展示如何使用warpper读取RSF文件。

首先，我们通过warp = warpper(file_path)来实例化该类，其中file_path为seismic.rsf文件的路径，然后调用该类的__set_attr_by_name__方法，即warp.__set_attr_by_name__('seismic')，此时已经创建了与seismic.rsf文件同名的warp.seismic方法，该文件中的数据即存储在numpy数组warp.seismic.data中，参数则存储在warp.seismic.pars中。

1
2
3
4
5
6
warp = warpper(file_path)
warp.__set_attr_by_name__('seismic')
# 数据
print(warp.seismic.data.shape)
# 参数
print(warp.seismic.pars)

实质上，数据和参数读取分别调用了m8r.Input.read()和m8r.Input.pars，这里只是将其封装在了新的类中，用户仍旧可以通过warp.seismic.Input.xxx来使用原m8r.Input中的其它方法，例如调用warp.seismic.grey()来显示并存储等。

Numpy数据写入RSF文件

Madagascrar中的部分功能需要使用.rsf文件作为输入参数，由于暂时未找到直接将numpy.ndarray数据作为参数输入m8r.func()中的方法，因此我们先将numpy.ndarray数据格式化写入.rsf文件，然后再调用函数从.rsf文件中读取数据。

数据准备

假设一炮集记录共55道，每道有1500个采样点，将二维炮集记录写入.rsf文件时，需要保证第一个维度为道，第二维度为每道的采样点，此数据形成的numpy数组shape应为(55, 1500)。

文件写入

首先，通过m8r.Output(file_path)创建一个与C中文件指针相当的函数，其中file_path为.rsf文件的路径（包含文件名），然后分别通过此类的put和write方法向文件指针中写入参数与数据，最后调用close方法关闭文件指针。代码示例如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
data.shape = (55, 1500)
type(data) = numpy.ndarray
data.dtype = np.float32

file = m8r.Output('data.rsf')
file.put('n1', data.shape[1])
file.put('n2', data.shape[0])
# 时间采样间隔
file.put('d1', 0.004)
# 空间采样间隔
file.put('d2', 1.)
file.write(data)
file.close()

将numpy格式炮集数据写入.rsf文件时，值得注意的有以下几点：

  1. numpy数据的第一和第二个维度分别为道和时间
  2. n1和n2分别表示时间采样点数和道数；
  3. d1和d2分别代表时间采样间隔和道间距。

结论

由于m8r中的函数不支持numpy数据作为参数输入（可能是由于源码中没有嵌入numpy格式），因此我们需要先将数据写入到.rsf文件中，然后再将文件路径作为参数传入到函数中，完成相关功能的实现。