【代码】XYZ格式地形(水深)数据转换为ASCIIGrid格式

水深测量数据，算是一种遥感数据，只是与平常的航天航空影像数据数据的格式、形式不一样。

这是最常用的水深测量方法，特别是在大型水体（如海洋和大湖）中。声纳（SONAR）是一种利用声波在水中传播和反射的技术，通过发送声波并测量其返回时间来计算水深。多普勒效应可以用来测量海底的粗糙度和其它特性。这种方法可以生成大面积的高分辨率水深图。

这次，展示水深测量数据的一般处理过程。

数据的格式和特点

本项目处理的数据是XYZ格式的地形数据。XYZ格式是一种常见的地理数据格式，其中每一行包含一个点的X、Y和Z坐标。在地形数据中，X和Y通常代表地理坐标（例如经度和纬度），Z代表地形高度或深度。

XYZ格式的特点包括：

  • 简单：每行包含一个数据点，每个数据点包含三个字段（X、Y和Z）。
  • 灵活：可以表示任何三维空间的数据点。
  • 通用：被许多地理信息系统（GIS）和数据处理工具所支持。

然而，XYZ格式的一个缺点是它不直接支持网格数据。在许多应用中，我们需要将散点数据转换为规则的网格，以便进行进一步的分析和可视化。这就是本项目的目的。

image-20240329150243476

数据处理过程

我们的任务是将XYZ格式的地形数据转换为ASCIIGrid格式。ASCIIGrid是一种网格数据格式，它以ASCII文本形式存储规则网格的单元格值。

转换过程包括以下步骤：

  1. 读取XYZ文件。
  2. 计算数据的最小和最大X和Y坐标。
  3. 根据给定的分辨率，为X和Y坐标创建bins。
  4. 使用numpy的histogram2d函数，根据X和Y的bins创建一个2D直方图（即网格）。这个函数会计算每个网格单元格中的点的数量和Z值的总和。
  5. 在每个网格单元格中，将Z值的总和除以点的数量，得到平均Z值。这就是我们的网格数据。
  6. 将没有数据的网格单元格（即没有点的单元格和只包含NaN值的单元格）的值设置为NoData值（在本例中为-9999）。
  7. 沿Y轴翻转网格，以纠正方向。
  8. 将网格数据写入ASCIIGrid文件。

代码


import numpy as np

def xyz_to_asciigrid(xyz_file, asc_file, resolution=1.0):
    # 读取XYZ文件（忽略第四个字段）
    xyz_data = np.loadtxt(xyz_file, usecols=(0, 1, 2))

    # 获取最小和最大坐标
    x_min, y_min = np.min(xyz_data[:, 0]), np.min(xyz_data[:, 1])
    x_max, y_max = np.max(xyz_data[:, 0]), np.max(xyz_data[:, 1])

    # 为x和y创建bins
    x_bins = np.arange(x_min, x_max + resolution, resolution)
    y_bins = np.arange(y_min, y_max + resolution, resolution)

    # 创建一个2D直方图来表示网格
    grid, _, _ = np.histogram2d(xyz_data[:, 1], xyz_data[:, 0], bins=(y_bins, x_bins), weights=xyz_data[:, 2])
    count, _, _ = np.histogram2d(xyz_data[:, 1], xyz_data[:, 0], bins=(y_bins, x_bins))

    # 在count不为零的地方，将网格除以count
    grid[count != 0] /= count[count != 0]

    # 用NoData值替换NaN值和count为零的单元格
    grid[np.isnan(grid) | (count == 0)] = -9999

    # 沿Y轴翻转网格以纠正方向
    grid = np.flipud(grid)

    # 写入ASCIIGrid文件
    nrows, ncols = grid.shape
    header = f"ncols {ncols}\nnrows {nrows}\nxllcorner {x_min}\nyllcorner {y_min}\ncellsize {resolution}\nNODATA_value -9999\n"
    with open(asc_file, 'w') as file:
        file.write(header)
        np.savetxt(file, grid, fmt="%f")

if __name__ == '__main__':
    xyz_file = r'bathymetry.xyz'
    asc_file = r'bathymetry.asc'
    xyz_to_asciigrid(xyz_file, asc_file)

bathymetry.asc结果展示：

image-20240329150415643

打开bathymetry.asc文件，文件前六行是数据的属性，

从第7行到最后一行，是一个网格数据，代表水深信息。其中-9999是无效值。

整个工程的文件如下：

image-20240329150323872