mkdocs/doc/meshkernelpy/examples/06_mesh2d_refinement_gridded_samples_gebco.md

基于不同水深精度的 GEBCO 网格采样数据进行二维网格加密

本教程简要介绍如何使用 GEBCO 全球数据集中的网格采样数据进行网格加密。

同时演示以下操作：

1. 在指定范围内生成网格。
2. 将生成的网格保存为 UGrid 文件。
3. 从大型 NetCDF 文件中读取水深数据。
4. 根据均匀间距水深采样数据（较快）或非均匀间距水深采样数据（较慢）进行网格加密。

返回示例目录

以下保留原笔记本的代码和已保存输出；转换过程中未重新执行代码。

导入所需的库，并关闭所有图窗。

import meshkernel
import numpy as np
import matplotlib.pyplot as plt

plt.close("all")

1. 在指定范围内生成网格

# 通用设置
lon_min, lon_max = -1, -0.2
lat_min, lat_max = 49.1, 49.6
lon_res, lat_res = 0.1, 0.1
figsize = (10, 4)
crs = "EPSG:4326"

"""
Make a regular (potentially rotated) rectilinear grid. First generate a curvilinear grid than convert the curvilinear grid into unstructured grid. The steps are the following:
- curvilinear_compute_uniform_on_extension, see the following notebook: https://github.com/Deltares/MeshKernelPy/blob/AddCurvilinearGridSupport/docs/examples/04_curvilineargrid_basics.ipynb
- curvilinear_convert_to_mesh2d: https://github.com/Deltares/MeshKernelPy/blob/118cb4953c4e95d5b18ed283bb37f391134b2bb2/meshkernel/meshkernel.py#L1399 
"""

# 创建 MakeGridParameters 实例并设置参数值
make_grid_parameters = meshkernel.MakeGridParameters()
make_grid_parameters.origin_x = lon_min
make_grid_parameters.origin_y = lat_min
make_grid_parameters.upper_right_x = lon_max
make_grid_parameters.upper_right_y = lat_max
make_grid_parameters.block_size_x = lon_res
make_grid_parameters.block_size_y = lat_res


mk2 = meshkernel.MeshKernel(projection=meshkernel.ProjectionType.SPHERICAL)
mk2.curvilinear_compute_rectangular_grid_on_extension(make_grid_parameters)
mk2.curvilinear_convert_to_mesh2d()  # 转换为 UGrid/二维网格

mesh2d = mk2.mesh2d_get()
fig, ax = plt.subplots()
mesh2d.plot_edges(ax)

2. 将生成的网格保存为 UGrid 文件

可以使用以下代码将生成的网格保存为 UGrid 文件。执行这些代码前，需要安装 UGrid Python 包。

ugrid_installed = False
if ugrid_installed:
    from ugrid import UGrid, UGridMesh2D

    mesh2d_ugrid = UGrid.from_meshkernel_mesh2d_to_ugrid_mesh2d(
        mesh2d=mesh2d, name="mesh2d", is_spherical=True
    )

    attribute_dict = {
        "name": "Unknown projected",
        "epsg": np.array([4326], dtype=int),
        "grid_mapping_name": "Unknown projected",
        "longitude_of_prime_meridian": np.array([0.0], dtype=float),
        "semi_major_axis": np.array([6378137.0], dtype=float),
        "semi_minor_axis": np.array([6356752.314245], dtype=float),
        "inverse_flattening": np.array([6356752.314245], dtype=float),
        "EPSG_code": "EPSG:4326",
        "value": "value is equal to EPSG code",
    }

    with UGrid("./gebco_mesh2d_net.nc", "w+") as ug:
        # 1. 定义新的二维网格
        topology_id = ug.mesh2d_define(mesh2d_ugrid)
        # 3. 写入新的二维网格
        ug.mesh2d_put(topology_id, mesh2d_ugrid)
        # 3. 将坐标参考系统添加到文件
        ug.variable_int_with_attributes_define("wgs84", attribute_dict)
        # 4. 添加数据约定（全局属性）
        conventions = {
            "institution": "Deltares",
            "references": "Unknown",
            "source": "Unknown Unknown. Model: Unknown",
            "history": "Created on 2017-11-27T18:05:09+0100, Unknown",
            "Conventions": "CF-1.6 UGRID-1.0/Deltares-0.8",
        }
        ug.attribute_global_define(conventions)

3. 从大型 NetCDF 文件中读取水深数据

可以使用 xarray 从大型数据集中读取水深数据，它支持仅加载选定区域的数据。当按非均匀间距处理经纬度数组时，双线性插值的速度较慢。

# 选择并绘制水深数据
uniform_grid_spacing = False
if not uniform_grid_spacing:
    import xarray as xr

    file_nc_bathy = r"p:\metocean-data\open\GEBCO\2021\GEBCO_2021.nc"
    data_bathy = xr.open_dataset(file_nc_bathy)
    data_bathy_sel = data_bathy.sel(
        lon=slice(lon_min - 1 / 4, lon_max + 1 / 4),
        lat=slice(lat_min - 1 / 4, lat_max + 1 / 4),
    )

    lon_np = data_bathy_sel.lon.to_numpy().flatten().astype("float")
    lat_np = data_bathy_sel.lat.to_numpy().flatten().astype("float")
    values_np = data_bathy_sel.elevation.to_numpy().flatten().astype("float32")

另一种方式是从 ASCII 文件读取数据，并假设采样间距均匀。

def read_asc_file(file_path, dtype=np.float32):
    """读取 ASC 文件，返回文件头和 NumPy 数组形式的数据。
    参数：
            file_path (str)：文件路径。
    返回值：
            header：ASCII 文件头。
            data：以双精度 NumPy 数组表示的 ASCII 数据。
    """

    header = {}
    data = []

    with open(file_path, "r") as file:
        # 读取文件头信息
        for _ in range(6):
            line = file.readline().strip().split()
            header[line[0]] = float(line[1])

        # 读取数据值
        for line in file:
            data_row = [float(value) for value in line.strip().split()]
            data.insert(0, data_row)  # 将该行插入到开头

        # 将数据展平
        data = np.array(data).flatten().astype(dtype)

    return header, data

4. 基于网格采样数据进行加密：水深类型为 np.float32 或 np.int16

将文件中的水深读取为 np.float32。

header, values_np = read_asc_file("./data_examples/gebco.asc", dtype=np.float32)

将文件中的水深读取为 np.int16。

header, values_np = read_asc_file("./data_examples/gebco.asc", dtype=np.int16)

采样网格属性。

num_x = int(header["ncols"])
num_y = int(header["nrows"])
x_origin = header["xllcenter"]
y_origin = header["yllcenter"]

网格采样数据。

gridded_samples = meshkernel.GriddedSamples(
    num_x=num_x,
    num_y=num_y,
    x_origin=x_origin,
    y_origin=y_origin,
    cell_size=0.0041666666666,
    values=values_np,
)

另一种方式是通过非均匀间距的 x、y 坐标数组生成网格采样数据。

if not uniform_grid_spacing:
    gridded_samples = meshkernel.GriddedSamples(
        x_coordinates=lon_np,
        y_coordinates=lat_np,
        num_x=len(lon_np),
        num_y=len(lat_np),
        values=values_np,
    )

定义网格加密参数。

mesh_refinement_parameters = meshkernel.MeshRefinementParameters(
    refine_intersected=False,
    use_mass_center_when_refining=False,
    min_edge_size=500,
    refinement_type=meshkernel.RefinementType.WAVE_COURANT,
    connect_hanging_nodes=True,
    account_for_samples_outside_face=False,
    max_refinement_iterations=3,
    smoothing_iterations=5,
    max_courant_time=120.0,
    directional_refinement=0,
)

执行加密，此处使用双线性插值。

mk2.mesh2d_refine_based_on_gridded_samples(
    gridded_samples=gridded_samples,
    mesh_refinement_params=mesh_refinement_parameters,
    use_nodal_refinement=True,
)

绘制加密后的网格。

# 放大绘图，以观察零散的海岸线
fig1 = plt.figure(figsize=(16, 12))
ax1 = fig1.add_subplot(111)

mesh2d_grid2 = mk2.mesh2d_get()
mesh2d_grid2.plot_edges(ax1, linewidth=1)