【代码】水深测量数据可视化及转格式

ytkz2024-03-292024-09-09

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

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

这次，展示水深测量数据的可视化及ASCII文件转换为TIF文件。

数据处理过程

我们的任务是

1.将ASCIIGrid格式的地形数据进行可视化。

2.将ASCIIGrid格式的地形数据转换为TIFF格式。

ASCIIGrid是一种网格数据格式，它以ASCII文本形式存储规则网格的单元格值。

以下是对这个任务的数据处理过程的详细描述：

读取ASCIIGrid格式的地形数据：首先，我们需要读取存储在ASCIIGrid格式文件中的地形数据。ASCIIGrid是一种常见的地理空间数据格式，它将地形数据存储为一个二维数组，每个元素代表一个地理位置的海拔高度。我们可以使用numpy的loadtxt函数来读取这个文件。由于ASCIIGrid文件的前几行包含了关于数据的元信息，我们需要跳过这些行（前6行），只读取实际的数据部分。
可视化地形数据：读取数据后，我们可以使用matplotlib库来可视化地形数据。这可以帮助我们理解数据的分布和特征。我们使用imshow函数来创建一个图像，其中每个像素的颜色代表一个地理位置的海拔高度。我们还添加了一个颜色条，以显示不同颜色对应的海拔高度。（视觉效果与一般的DEM数据相似）
转换数据格式：最后，我们将地形数据从ASCIIGrid格式转换为TIFF格式。TIFF是一种常用的图像格式，它可以存储大量的元数据，并且支持多种数据类型和颜色空间。我们使用gdal库来创建一个新的TIFF文件，并将数据写入这个文件。我们需要指定文件的高度、宽度、数据类型、坐标参考系统（CRS）和转换矩阵。

以上就是我们处理地形数据的过程。这个过程涵盖了数据读取、可视化和格式转换三个主要步骤，每个步骤都是为了实现我们的最终目标：理解地形数据并将其转换为我们需要的格式。

可视化部分的代码

import matplotlib.pyplot as plt
from osgeo import gdal, osr
import numpy as np
def visualize_asciigrid(asc_file):
    # 读取ASCIIGrid文件
    # 需要跳过前6行，因为这些行包含了文件头
    data = np.loadtxt(asc_file, skiprows=6)

    # 将NoData值（在这个例子中是-9999）替换为NaN，以便在可视化中忽略它们
    data[data == -9999] = np.nan

    # 用黑体显示中文, 解决中文乱码问题
    import matplotlib
    matplotlib.rcParams['font.sans-serif'] = ['SimHei']

    # 使用imshow显示数据
    plt.imshow(data, cmap='jet')
    plt.colorbar(label='Value')
    plt.title('水深网格')
    plt.show()
visualize_asciigrid('bathymetry.asc')

可视化结果如下：

格式转换部分的代码

def read_and_georeference_asc(asc_file, output_file):
    # 读取头部信息
    with open(asc_file) as f:
        ncols = int(f.readline().split()[1])
        nrows = int(f.readline().split()[1])
        xllcorner = float(f.readline().split()[1])
        yllcorner = float(f.readline().split()[1])
        cellsize = float(f.readline().split()[1])
        no_data = float(f.readline().split()[1])

    # 读取数据
    data = np.loadtxt(asc_file, skiprows=6)

    # 创建输出文件
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(output_file, ncols, nrows, 1, gdal.GDT_Float32)

    # 设置地理参考
    dataset.SetGeoTransform((xllcorner, cellsize, 0, yllcorner + nrows * cellsize, 0, -cellsize))

    # 写入数据
    dataset.GetRasterBand(1).WriteArray(data)
    dataset.GetRasterBand(1).SetNoDataValue(no_data)

    dataset = None  # 关闭文件