wps的domain转为shp矢量

news2024/10/9 16:22:32

wps的namelist制作、python出图和转矢量

简介

wps（WRF Preprocessing System）是中尺度数值天气预报系统WRF(Weather Research and Forecasting)的预处理系统。

wrf模型示意图

wps的安装地址在GitHub上：https://github.com/wrf-model/WPS

下载完成后，就可以进行WPS安装，教程请查看官网：https://www2.mmm.ucar.edu/wrf/OnLineTutorial/compilation_tutorial.php

在wps安装成功后，会得到三个编译软件分别为geogrid、ungrib和metgrid。分别的功能是：

  • geogrid处理下垫面数据；

  • ungrib是将不同格式的气象数据转为统一格式

  • metgrid会将气象数据统一"裁剪"下垫面数据的网格上。

这三个exe文件的控制参数都是由namelist.wps进行控制的。

我今天想介绍的主要内容为根据namelist.input反推domain的空间位置并导出shp矢量，其他东西只会简单介绍。

namelist参数含义

不要去看汉化后的博客，官网有详细解释的：https://www2.mmm.ucar.edu/wrf/users/namelist_best_prac_wps.html#io_form_geogrid

由于我们要根据namelist.wps反推domain的空间分布，由几个特别重要的的参数，官网没有详细解释，我这里做一下补充：

（1）每一个domain都是由格网组成的，每层的格网大小由dx、dy和parent_grid_ratio决定。

请注意，e_we和e_sn代表的是格的端点数，它们如果为n，则对应的格子数目为n-1，如图所示（e_we不会只有3，只是距离）：

（2）e_we和e_sn都是指在该层范围内分辨率下的格子数目，比如304就是d02每行拥有303个格子。

（3）i_parent_start是经度的指标，一般用X表示，数字代表的子域d02左下角和d01的左下角的x方向的水平方向格子数。特别需要注意，

i_parent_start所代表的格子数量是d01的，不是d02。

j_parent_start是纬度的指标，一般用Y标识，是垂直方向的格子数。

domain坐标计算的原理

（1）d01层的平面坐标计算

ref_lon,ref_lat是d01的中心位置，此外，它还有个特殊的作用：在wps中所有网格都是在一个投影的平面坐标系中的，(ref_lon,ref_lat)代表了在平面坐标系下的原点（0，0）。如果我们以d01举例，求它的四个顶点的坐标，我们需要明确，在兰伯特投影下的平面坐标系单位是米，d01的格子横向长度为dx，d01的格子纵向长度为dy。现在我们就可以求出d01的四个顶点在投影的平面坐标系下的坐标。

（2）d02、d03层的平面坐标计算

上面一步，我们已经求得了d01的四个顶点的坐标值。

d02的坐标，我们需要根据d01的左下角的坐标求得。

由此，我们获得了d02的四个顶点在投影下的平面坐标。

在获得d02的左下角坐标之后，按照相同方法，我们可以获得d03的平面坐标。

自此，我们明白了wps的namelist.wps的平面坐标计算原理，开始写代码。计算坐标点的代码如下：

# 初始化网格域
def initialize_domains(e_we, e_sn, i_parent_start, j_parent_start, parent_id, parent_grid_ratio, dx, dy, ref_lat, ref_lon, truelat1, truelat2):
    domains = []

    # 坐标系详细信息
    # 兰伯特 坐标系信息
    proj_params = {
        'proj': 'lcc',
        'lat_1': truelat1,
        'lat_2': truelat2,
        'lat_0': ref_lat,
        'lon_0': ref_lon,
        'x_0': 0,
        'y_0': 0,
        'datum': 'WGS84'
    }

    # 创建平面坐标系
    p_lcc = Proj(proj_params)

    # 计算d01的中心点平面坐标系的坐标
    x_center, y_center = p_lcc(ref_lon, ref_lat)

    print("x_center, y_center",x_center, y_center)

    for i in range(len(e_we)):
        if parent_id[i] == 1:
            # 针对d01
            if i == 0:
                parent_dx = dx
                parent_dy = dy

                parent_ref_lon = x_center
                parent_ref_lat = y_center

                grid_ratio = parent_grid_ratio[i]
                # 计算d01的左下角坐标
                d01_start_x = x_center - ((e_we[i] - 1) / 2) * parent_dx
                d01_start_y = y_center - ((e_sn[i] - 1) / 2) * parent_dy

                # 计算d01的平面坐标
                domains.append(compute_boundaries(parent_ref_lat, parent_ref_lon, e_we[i], e_sn[i], dx, dy, i_parent_start[i], j_parent_start[i], parent_dx, parent_dy, d01_start_y, d01_start_x, grid_ratio))
            # 针对d02
            else:
                parent_domain_idx = parent_id[i] - 1
                grid_ratio = parent_grid_ratio[i]
                parent_dx = dx
                parent_dy = dy
                dx = dx / grid_ratio
                dy = dy / grid_ratio
                parent_ref_lat = domains[parent_domain_idx][2]
                parent_ref_lon = domains[parent_domain_idx][0]

                # 计算d02的左下角坐标
                d02_start_x = domains[parent_domain_idx][0] + (i_parent_start[i] - 1) * parent_dx
                d02_start_y = domains[parent_domain_idx][2] + (j_parent_start[i] - 1) * parent_dy

                # 计算d02的经纬度
                domains.append(compute_boundaries(parent_ref_lat, parent_ref_lon, e_we[i], e_sn[i], dx, dy, i_parent_start[i], j_parent_start[i], parent_dx, parent_dy, d02_start_y, d02_start_x, grid_ratio))
        else:
            # 针对d03
            parent_domain_idx = parent_id[i] - 1
            grid_ratio = parent_grid_ratio[i]

            parent_dx = domains[parent_domain_idx][6]
            parent_dy = domains[parent_domain_idx][7]

            parent_ref_lat = domains[parent_domain_idx][2]
            parent_ref_lon = domains[parent_domain_idx][0]
            # 计算d03的左下角坐标
            d03_start_x = domains[parent_domain_idx][0] + (i_parent_start[i] - 1) * parent_dx
            d03_start_y = domains[parent_domain_idx][2] + (j_parent_start[i] - 1) * parent_dy

            # 计算d03的经纬度
            domains.append(compute_boundaries(parent_ref_lat, parent_ref_lon, e_we[i], e_sn[i], dx / grid_ratio, dy / grid_ratio, i_parent_start[i], j_parent_start[i], parent_dx, parent_dy, d03_start_y, d03_start_x, grid_ratio))

    return domains, proj_params

# 计算网格的边界
def compute_boundaries(ref_lat, ref_lon, e_we, e_sn, dx, dy, i_start, j_start, parent_dx, parent_dy, d_start_y, d_start_x, grid_ratio):
    # 计算子网格的左下角坐标
    grid_start_lon = d_start_x
    grid_start_lat = d_start_y

    # 计算右上角坐标
    grid_end_lon = grid_start_lon + (e_we - 1) * dx
    grid_end_lat = grid_start_lat + (e_sn - 1) * dy

    # 各方向的坐标
    west = grid_start_lon
    south = grid_start_lat
    east = grid_end_lon
    north = grid_end_lat

    grid_center_lat = (south + north) / 2
    grid_center_lon = (west + east) / 2

    return west, east, south, north, grid_center_lat, grid_center_lon, dx, dy

python出图

# 绘制地图
def plot_map(domains, gdf_level1, gdf_level2, gdf_level3, truelat1, truelat2, stand_lon, proj_params):
    proj = ccrs.LambertConformal(central_longitude=stand_lon, standard_parallels=(truelat1, truelat2))

    fig, ax = plt.subplots(figsize=(10, 10), subplot_kw={'projection': proj})

    # 绘制shapefile背景
    gdf_level1.to_crs(proj).plot(ax=ax, edgecolor='blue', facecolor='none', linewidth=1)  # 蓝色边框,空心
    gdf_level2.to_crs(proj).plot(ax=ax, edgecolor='red', facecolor='none', linewidth=1)   # 红色边框,空心
    gdf_level3.to_crs(proj).plot(ax=ax, edgecolor='black', facecolor='#43A7EE', linewidth=1)  # 黑色边框,RGB(67,167,238)填充

    # 创建坐标系对象
    crs_wgs84 = CRS.from_epsg(4326)  # 使用 EPSG 代码 4326 表示 WGS84 地理坐标系
    crs_lcc = CRS(proj_params)
    transformer = Transformer.from_crs(crs_lcc, crs_wgs84, always_xy=True)

    # 绘制每个嵌套网格的范围
    for i, (west, east, south, north, _, _, _, _) in enumerate(domains):
        # 将网格边界转换为经纬度
        # 左下
        west_lon_ZUO, south_lat_ZUO = transformer.transform(west, south)
        # 左上
        west_lon_ZUO2, north_lat_ZUO2 = transformer.transform(west, north)
        # 右上
        east_lon_YOU2, north_lat_YOU2 = transformer.transform(east, north)
        # 右下
        east_lon_YOU, south_lat_YOU = transformer.transform(east, south)

        # 打印经纬度范围
        print(f"Domain {i+1} bounds (west, east, south, north):")
        print(f"Longitude: {west_lon_ZUO} to {east_lon_YOU2}")
        print(f"Latitude: {south_lat_ZUO} to {north_lat_ZUO2}")

        # 使用 Polygon 创建每个网格的多边形,按照逆时针顺序连接点
        vertices = [(west_lon_ZUO, south_lat_ZUO), (east_lon_YOU, south_lat_YOU), (east_lon_YOU2, north_lat_YOU2), (west_lon_ZUO2, north_lat_ZUO2)]
        polygon = Polygon(vertices)
        ax.add_geometries([polygon], crs=ccrs.PlateCarree(), edgecolor='red' if i == 0 else 'blue' if i == 1 else 'green', facecolor='none', linewidth=2, label=f'Domain {i+1}')
        
        # 计算标注的位置(使用多边形的右上角点)
        lon_label, lat_label = east_lon_YOU2, north_lat_YOU2

        # 添加标注,并调整标注的位置
        ax.text(lon_label, lat_label, f'd0{i+1}', color='red', fontsize=12, ha='right', va='top', transform=ccrs.PlateCarree())

    # 转换 d01 的边界坐标到地理坐标
    west_d01, south_d01 = transformer.transform(domains[0][0], domains[0][2])
    east_d01, north_d01 = transformer.transform(domains[0][1], domains[0][3])

    print("west_d01, south_d01, east_d01, north_d01", west_d01, south_d01, east_d01, north_d01)
    # 设置显示范围,在经度和纬度方向上各自添加边距
    lon_margin = (east_d01 - west_d01) * 0.1  # 经度方向上的边距为d01经度范围的10%
    lat_margin = (north_d01 - south_d01) * 0.1  # 纬度方向上的边距为d01纬度范围的10%
    ax.set_extent([west_d01 - lon_margin, east_d01 + lon_margin, south_d01 - lat_margin, north_d01 + lat_margin], crs=ccrs.PlateCarree())

    # 添加海岸线和网格线
    ax.gridlines(draw_labels=True)

    plt.title('WRF Domains')
    plt.show()

我们加入研究区的矢量如下，出图效果如下：

namelist转矢量shp

既然我们已经知道d01到d03的坐标点，按照逆时针把矢量点串联起来，获得shp矢量。

# 输出d01到d03范围为shp
def output_domains_to_shapefile(domains, proj_params, output_shapefile_path):
    # 创建一个空的 GeoDataFrame,用于存储域的范围
    gdf_domains = gpd.GeoDataFrame(columns=['geometry', 'domain_id'], crs="EPSG:4326")
    
    # 创建坐标系对象
    crs_wgs84 = CRS.from_epsg(4326)  # 使用 EPSG 代码 4326 表示 WGS84 地理坐标系
    crs_lcc = CRS(proj_params)
    transformer = Transformer.from_crs(crs_lcc, crs_wgs84, always_xy=True)

    # 绘制每个嵌套网格的范围并添加到 GeoDataFrame
    for i, (west, east, south, north, _, _, _, _) in enumerate(domains):
        # 将网格边界转换为经纬度
        # 左下
        west_lon_ZUO, south_lat_ZUO = transformer.transform(west, south)
        # 左上
        west_lon_ZUO2, north_lat_ZUO2 = transformer.transform(west, north)
        # 右上
        east_lon_YOU2, north_lat_YOU2 = transformer.transform(east, north)
        # 右下
        east_lon_YOU, south_lat_YOU = transformer.transform(east, south)

        # 使用 Polygon 创建每个网格的多边形,按照逆时针顺序连接点
        vertices = [(west_lon_ZUO, south_lat_ZUO), (east_lon_YOU, south_lat_YOU), (east_lon_YOU2, north_lat_YOU2), (west_lon_ZUO2, north_lat_ZUO2)]
        polygon = Polygon(vertices)
        
        # 创建一个临时的 GeoDataFrame
        temp_gdf = gpd.GeoDataFrame([{'geometry': polygon, 'domain_id': f'Domain {i+1}'}], crs="EPSG:4326")
        
        # 使用 pd.concat 将临时的 GeoDataFrame 添加到主要的 GeoDataFrame 中
        gdf_domains = pd.concat([gdf_domains, temp_gdf], ignore_index=True)

    # 保存为 shapefile
    gdf_domains.to_file(output_shapefile_path, driver='ESRI Shapefile')

特别注意，我们需要最后生成的shp是wgs84坐标系，所以需要把平面坐标转回为wgs84坐标系。

完整代码

import re
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import geopandas as gpd
import pandas as pd
from shapely.geometry import Polygon
from pyproj import Proj,transform
from pyproj import CRS, Transformer


# 提取单个参数的函数
def get_param(pattern, content, index=0):
    match = re.search(pattern, content)
    if match:
        return float(match.group(1))
    else:
        raise ValueError(f'Parameter {pattern} not found in namelist.wps')

# 提取多个参数的函数
def get_params(pattern, content):
    match = re.search(pattern, content)
    if match:
        return [int(x) for x in match.group(1).split(',')]
    else:
        raise ValueError(f'Parameter {pattern} not found in namelist.wps')

# 读取并解析namelist.wps文件
def parse_namelist(namelist_path):
    with open(namelist_path, 'r') as file:
        namelist_content = file.read()

    dx = get_param(r'dx\s*=\s*(\d+)', namelist_content)
    dy = get_param(r'dy\s*=\s*(\d+)', namelist_content)
    ref_lat = get_param(r'ref_lat\s*=\s*([-+]?\d*\.\d+|\d+)', namelist_content)
    ref_lon = get_param(r'ref_lon\s*=\s*([-+]?\d*\.\d+|\d+)', namelist_content)
    e_we = get_params(r'e_we\s*=\s*([\d,\s]+)', namelist_content)
    e_sn = get_params(r'e_sn\s*=\s*([\d,\s]+)', namelist_content)
    i_parent_start = get_params(r'i_parent_start\s*=\s*([\d,\s]+)', namelist_content)
    j_parent_start = get_params(r'j_parent_start\s*=\s*([\d,\s]+)', namelist_content)
    parent_id = get_params(r'parent_id\s*=\s*([\d,\s]+)', namelist_content)
    parent_grid_ratio = get_params(r'parent_grid_ratio\s*=\s*([\d,\s]+)', namelist_content)
    truelat1 = get_param(r'truelat1\s*=\s*([-+]?\d*\.\d+|\d+)', namelist_content)
    truelat2 = get_param(r'truelat2\s*=\s*([-+]?\d*\.\d+|\d+)', namelist_content)
    stand_lon = get_param(r'stand_lon\s*=\s*([-+]?\d*\.\d+|\d+)', namelist_content)

    return dx, dy, ref_lat, ref_lon, e_we, e_sn, i_parent_start, j_parent_start, parent_id, parent_grid_ratio, truelat1, truelat2, stand_lon

# 初始化网格域
def initialize_domains(e_we, e_sn, i_parent_start, j_parent_start, parent_id, parent_grid_ratio, dx, dy, ref_lat, ref_lon, truelat1, truelat2):
    domains = []

    # 坐标系详细信息
    # 兰伯特 坐标系信息
    proj_params = {
        'proj': 'lcc',
        'lat_1': truelat1,
        'lat_2': truelat2,
        'lat_0': ref_lat,
        'lon_0': ref_lon,
        'x_0': 0,
        'y_0': 0,
        'datum': 'WGS84'
    }

    # 创建平面坐标系
    p_lcc = Proj(proj_params)

    # 计算d01的中心点平面坐标系的坐标
    x_center, y_center = p_lcc(ref_lon, ref_lat)

    print("x_center, y_center",x_center, y_center)

    for i in range(len(e_we)):
        if parent_id[i] == 1:
            # 针对d01
            if i == 0:
                parent_dx = dx
                parent_dy = dy

                parent_ref_lon = x_center
                parent_ref_lat = y_center

                grid_ratio = parent_grid_ratio[i]
                # 计算d01的左下角坐标
                d01_start_x = x_center - ((e_we[i] - 1) / 2) * parent_dx
                d01_start_y = y_center - ((e_sn[i] - 1) / 2) * parent_dy

                # 计算d01的平面坐标
                domains.append(compute_boundaries(parent_ref_lat, parent_ref_lon, e_we[i], e_sn[i], dx, dy, i_parent_start[i], j_parent_start[i], parent_dx, parent_dy, d01_start_y, d01_start_x, grid_ratio))
            # 针对d02
            else:
                parent_domain_idx = parent_id[i] - 1
                grid_ratio = parent_grid_ratio[i]
                parent_dx = dx
                parent_dy = dy
                dx = dx / grid_ratio
                dy = dy / grid_ratio
                parent_ref_lat = domains[parent_domain_idx][2]
                parent_ref_lon = domains[parent_domain_idx][0]

                # 计算d02的左下角坐标
                d02_start_x = domains[parent_domain_idx][0] + (i_parent_start[i] - 1) * parent_dx
                d02_start_y = domains[parent_domain_idx][2] + (j_parent_start[i] - 1) * parent_dy

                # 计算d02的经纬度
                domains.append(compute_boundaries(parent_ref_lat, parent_ref_lon, e_we[i], e_sn[i], dx, dy, i_parent_start[i], j_parent_start[i], parent_dx, parent_dy, d02_start_y, d02_start_x, grid_ratio))
        else:
            # 针对d03
            parent_domain_idx = parent_id[i] - 1
            grid_ratio = parent_grid_ratio[i]

            parent_dx = domains[parent_domain_idx][6]
            parent_dy = domains[parent_domain_idx][7]

            parent_ref_lat = domains[parent_domain_idx][2]
            parent_ref_lon = domains[parent_domain_idx][0]
            # 计算d03的左下角坐标
            d03_start_x = domains[parent_domain_idx][0] + (i_parent_start[i] - 1) * parent_dx
            d03_start_y = domains[parent_domain_idx][2] + (j_parent_start[i] - 1) * parent_dy

            # 计算d03的经纬度
            domains.append(compute_boundaries(parent_ref_lat, parent_ref_lon, e_we[i], e_sn[i], dx / grid_ratio, dy / grid_ratio, i_parent_start[i], j_parent_start[i], parent_dx, parent_dy, d03_start_y, d03_start_x, grid_ratio))

    return domains, proj_params

# 计算网格的边界
def compute_boundaries(ref_lat, ref_lon, e_we, e_sn, dx, dy, i_start, j_start, parent_dx, parent_dy, d_start_y, d_start_x, grid_ratio):
    # 计算子网格的左下角坐标
    grid_start_lon = d_start_x
    grid_start_lat = d_start_y

    # 计算右上角坐标
    grid_end_lon = grid_start_lon + (e_we - 1) * dx
    grid_end_lat = grid_start_lat + (e_sn - 1) * dy

    # 各方向的坐标
    west = grid_start_lon
    south = grid_start_lat
    east = grid_end_lon
    north = grid_end_lat

    grid_center_lat = (south + north) / 2
    grid_center_lon = (west + east) / 2

    return west, east, south, north, grid_center_lat, grid_center_lon, dx, dy

# 绘制地图
def plot_map(domains, gdf_level1, gdf_level2, gdf_level3, truelat1, truelat2, stand_lon, proj_params):
    proj = ccrs.LambertConformal(central_longitude=stand_lon, standard_parallels=(truelat1, truelat2))

    fig, ax = plt.subplots(figsize=(10, 10), subplot_kw={'projection': proj})

    # 绘制shapefile背景
    gdf_level1.to_crs(proj).plot(ax=ax, edgecolor='blue', facecolor='none', linewidth=1)  # 蓝色边框,空心
    gdf_level2.to_crs(proj).plot(ax=ax, edgecolor='red', facecolor='none', linewidth=1)   # 红色边框,空心
    gdf_level3.to_crs(proj).plot(ax=ax, edgecolor='black', facecolor='#43A7EE', linewidth=1)  # 黑色边框,RGB(67,167,238)填充

    # 创建坐标系对象
    crs_wgs84 = CRS.from_epsg(4326)  # 使用 EPSG 代码 4326 表示 WGS84 地理坐标系
    crs_lcc = CRS(proj_params)
    transformer = Transformer.from_crs(crs_lcc, crs_wgs84, always_xy=True)

    # 绘制每个嵌套网格的范围
    for i, (west, east, south, north, _, _, _, _) in enumerate(domains):
        # 将网格边界转换为经纬度
        # 左下
        west_lon_ZUO, south_lat_ZUO = transformer.transform(west, south)
        # 左上
        west_lon_ZUO2, north_lat_ZUO2 = transformer.transform(west, north)
        # 右上
        east_lon_YOU2, north_lat_YOU2 = transformer.transform(east, north)
        # 右下
        east_lon_YOU, south_lat_YOU = transformer.transform(east, south)

        # 打印经纬度范围
        print(f"Domain {i+1} bounds (west, east, south, north):")
        print(f"Longitude: {west_lon_ZUO} to {east_lon_YOU2}")
        print(f"Latitude: {south_lat_ZUO} to {north_lat_ZUO2}")

        # 使用 Polygon 创建每个网格的多边形,按照逆时针顺序连接点
        vertices = [(west_lon_ZUO, south_lat_ZUO), (east_lon_YOU, south_lat_YOU), (east_lon_YOU2, north_lat_YOU2), (west_lon_ZUO2, north_lat_ZUO2)]
        polygon = Polygon(vertices)
        ax.add_geometries([polygon], crs=ccrs.PlateCarree(), edgecolor='red' if i == 0 else 'blue' if i == 1 else 'green', facecolor='none', linewidth=2, label=f'Domain {i+1}')
        
        # 计算标注的位置(使用多边形的右上角点)
        lon_label, lat_label = east_lon_YOU2, north_lat_YOU2

        # 添加标注,并调整标注的位置
        ax.text(lon_label, lat_label, f'd0{i+1}', color='red', fontsize=12, ha='right', va='top', transform=ccrs.PlateCarree())

    # 转换 d01 的边界坐标到地理坐标
    west_d01, south_d01 = transformer.transform(domains[0][0], domains[0][2])
    east_d01, north_d01 = transformer.transform(domains[0][1], domains[0][3])

    print("west_d01, south_d01, east_d01, north_d01", west_d01, south_d01, east_d01, north_d01)
    # 设置显示范围,在经度和纬度方向上各自添加边距
    lon_margin = (east_d01 - west_d01) * 0.1  # 经度方向上的边距为d01经度范围的10%
    lat_margin = (north_d01 - south_d01) * 0.1  # 纬度方向上的边距为d01纬度范围的10%
    ax.set_extent([west_d01 - lon_margin, east_d01 + lon_margin, south_d01 - lat_margin, north_d01 + lat_margin], crs=ccrs.PlateCarree())

    # 添加海岸线和网格线
    ax.gridlines(draw_labels=True)

    plt.title('WRF Domains')
    plt.show()

# 输出d01到d03范围为shp
def output_domains_to_shapefile(domains, proj_params, output_shapefile_path):
    # 创建一个空的 GeoDataFrame,用于存储域的范围
    gdf_domains = gpd.GeoDataFrame(columns=['geometry', 'domain_id'], crs="EPSG:4326")
    
    # 创建坐标系对象
    crs_wgs84 = CRS.from_epsg(4326)  # 使用 EPSG 代码 4326 表示 WGS84 地理坐标系
    crs_lcc = CRS(proj_params)
    transformer = Transformer.from_crs(crs_lcc, crs_wgs84, always_xy=True)

    # 绘制每个嵌套网格的范围并添加到 GeoDataFrame
    for i, (west, east, south, north, _, _, _, _) in enumerate(domains):
        # 将网格边界转换为经纬度
        # 左下
        west_lon_ZUO, south_lat_ZUO = transformer.transform(west, south)
        # 左上
        west_lon_ZUO2, north_lat_ZUO2 = transformer.transform(west, north)
        # 右上
        east_lon_YOU2, north_lat_YOU2 = transformer.transform(east, north)
        # 右下
        east_lon_YOU, south_lat_YOU = transformer.transform(east, south)

        # 使用 Polygon 创建每个网格的多边形,按照逆时针顺序连接点
        vertices = [(west_lon_ZUO, south_lat_ZUO), (east_lon_YOU, south_lat_YOU), (east_lon_YOU2, north_lat_YOU2), (west_lon_ZUO2, north_lat_ZUO2)]
        polygon = Polygon(vertices)
        
        # 创建一个临时的 GeoDataFrame
        temp_gdf = gpd.GeoDataFrame([{'geometry': polygon, 'domain_id': f'Domain {i+1}'}], crs="EPSG:4326")
        
        # 使用 pd.concat 将临时的 GeoDataFrame 添加到主要的 GeoDataFrame 中
        gdf_domains = pd.concat([gdf_domains, temp_gdf], ignore_index=True)

    # 保存为 shapefile
    gdf_domains.to_file(output_shapefile_path, driver='ESRI Shapefile')

def main():
    # 读取namelist.wps文件
    read_path = r"E:\ruiduobao\namelis设置\namelist.wps"
    dx, dy, ref_lat, ref_lon, e_we, e_sn, i_parent_start, j_parent_start, parent_id, parent_grid_ratio, truelat1, truelat2, stand_lon = parse_namelist(read_path)

    # 初始化网格域
    domains, proj_params = initialize_domains(e_we, e_sn, i_parent_start, j_parent_start, parent_id, parent_grid_ratio, dx, dy, ref_lat, ref_lon, truelat1, truelat2)

    # 读取shapefile
    shapefile_path_level1 = r'E:\ruiduobao\数据和代码\行政区划\jiangsu.shp'
    shapefile_path_level2 = r'E:\ruiduobao\数据和代码\行政区划\xuzhou.shp'
    shapefile_path_level3 = r'E:\ruiduobao\数据和代码\行政区划\xuzhouxian.shp'

    # 加载shapefile到 GeoDataFrame
    gdf_level1 = gpd.read_file(shapefile_path_level1)
    gdf_level2 = gpd.read_file(shapefile_path_level2)
    gdf_level3 = gpd.read_file(shapefile_path_level3)

    # 绘制地图
    plot_map(domains, gdf_level1, gdf_level2, gdf_level3, truelat1, truelat2, stand_lon, proj_params)

    # 输出d01到d03范围为shp
    output_shapefile_path = r'E:\ruiduobao\行政区划\wrf_domains_平面.shp'
    output_domains_to_shapefile(domains, proj_params, output_shapefile_path)

if __name__ == '__main__':
    main()

最后，我们把生成的namelist.wps的矢量放到GIS软件中，就可以任意编辑了：

总结

这个代码看起来很简单，但我实际上搞了快两天才弄懂里面的原理，尴尬，故写一篇技术博客方便以后自己查阅。我主要遇到以下问题:

(1)一开始我是用wgs84的经纬度去计算的各个domain的空间位置的，但实际上会有偏移，因为每度随着纬度的不同是会变化的，需要放到平面坐标系中才会有正确的结果。

(2)我刚开始是计算每一个domain的中心点，但实际上这是比较傻的方法，因为i_parent_start等是从左下角开始计数的。

参考

https://github.com/wrf-model/WPS

https://www2.mmm.ucar.edu/wrf/OnLineTutorial/compilation_tutorial.php

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.coloradmin.cn/o/1863223.html

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

循环神经网络——RNN

循环神经网络——RNN

循环神经网络 在之前NLP基础章节-语言模型中我们介绍了 n n n 元语法，其中单词 x t x_t xt​ 在时间步 t t t 的条件概率仅取决于前面 n n n 个单词，若是想要将之前单词的影响也加入那么模型参数数量会指数级增长。但是可能之前的单词存在重要的信息…
阅读更多...
Linux-笔记 高级I/O操作

Linux-笔记 高级I/O操作

前言 I/O（Input/Output，输入/输出）是计算机系统中的一个重要组成部分，它是指计算机与 外部世界之间的信息交流过程。I/O 操作是计算机系统中的一种基本操作，用于向外部设备（如 硬盘、键盘、鼠标、网络等&am…
阅读更多...
服务器数据恢复—异常断电导致RAID6阵列中磁盘出现坏扇区的数据恢复案例

服务器数据恢复—异常断电导致RAID6阵列中磁盘出现坏扇区的数据恢复案例

服务器存储数据恢复环境： 一台存储中有一组由12块SAS硬盘组建的RAID6磁盘阵列，划分为一个卷，分配给几台Vmware ESXI主机做共享存储。该卷中存放了大量Windows虚拟机，这些虚拟机系统盘是统一大小，数据盘大小不确定&…
阅读更多...
服务器硬件及RAID配置

服务器硬件及RAID配置

目录 一、RAID磁盘阵列 1.概念 2.RAID 0 3.RAID 1 4.RAID 5 5.RAID 6 6.RAID 10 二、阵列卡 1.简介 2.缓存 三、创建 1.创建RAID 0 2.创建RAID 1 3.创建RAID 5 4.创建RAID 10 四、模拟故障 一、RAID磁盘阵列 1.概念 （1）是Redundant Array …
阅读更多...
求任意方阵每行，每列，两对角线上元素之和

求任意方阵每行，每列，两对角线上元素之和

注：其中对角线，我们可以分为正副两个，正：左上角指向右下角，副：右上角指向左下角 //这里我们以阶层为5为例子进行代码的实现 #define N 5 void arr_diagonal(int arr[N][N]) {int sum1 0, sum2 0, sum 0…
阅读更多...
Js逆向爬虫基础篇

Js逆向爬虫基础篇

这里写自定义目录标题 逆向技巧断点一 、请求入口定位1. 关键字搜索2. 请求堆栈3. hook4. JSON.stringify 二、响应入口定位：1. 关键字搜索2. hook3. JSON.parse 逆向技巧 断点 普通断点 条件断点 日志断点 XHR断点 一 、请求入口定位 1. 关键字搜索 key关…
阅读更多...
C++ | Leetcode C++题解之第198题打家劫舍

C++ | Leetcode C++题解之第198题打家劫舍

题目&#xff1a; 题解&#xff1a; class Solution { public:int rob(vector<int>& nums) {if (nums.empty()) {return 0;}int size nums.size();if (size 1) {return nums[0];}int first nums[0], second max(nums[0], nums[1]);for (int i 2; i < size; …
阅读更多...
10.XSS绕过之htmlspecialchars()函数

10.XSS绕过之htmlspecialchars()函数

XSS绕过之htmlspecialchars()函数 首先可以测试一下是否将字符被转移成html实体&#xff0c;输入字符测试 1111"<>$点击提交 查看页面元素代码&#xff0c;发现单引号不变&#xff0c;可以利用 重新输入攻击代码&#xff0c;用单引号闭合前面的&#xff0c;进…
阅读更多...
深圳大学 软件测试作业 #2

深圳大学 软件测试作业 #2

声明&#xff1a;本人上课摆烂选手&#xff0c;稍微听了下&#xff0c;答案仅供参考。 ———————— 1. 考虑下面这个代码&#xff0c;并回答以下的问题。 (a) 请画出上面代码的控制流程图。(20分) (b) 请画出上面代码的数据流程图。(10分) (c) 找出每个变量的定义使…
阅读更多...
SpringBoot整合Mybatis并实现数据库增删改查

SpringBoot整合Mybatis并实现数据库增删改查

写在前面 Mybatis一个基于Java的持久层框架&#xff0c;它通过XML或注解的方式&#xff0c;将SQL语句和Java方法进行映射&#xff0c;使得开发者可以轻松地进行数据库操作。下面我会演示mybatis的配置与使用并实现数据库的增删改查。 1.准备测试数据 使用mybatis实现对数据库…
阅读更多...
Java银系统/超市收银系统/智慧新零售/ERP进销存管理/线上商城/h5/小程序

Java银系统/超市收银系统/智慧新零售/ERP进销存管理/线上商城/h5/小程序

>>>系统简述&#xff1a; 神点收银系统支持B2B2C多商户模式&#xff0c;系统基于前后端分离的架构&#xff0c;后端采用Java SpringBoot Mysql Mybatis Plus&#xff0c;前端基于当前流行的Uniapp、Element UI&#xff0c;支持小程序、h5。架构包含&#xff1a;会员端…
阅读更多...
AI智能写作工具,AI写作助手大全

AI智能写作工具,AI写作助手大全

随着人工智能技术的快速发展&#xff0c;AI智能写作工具助手已成为学术研究、内容创作和商业文案等领域的重要辅助工具。它们不仅能够提高写作效率&#xff0c;还能激发创意灵感&#xff0c;为各行各业的专业人士提供了强大的支持。下面小编将为大家全面介绍目前市场上备受瞩目…
阅读更多...
Mac(M1芯片)安装多个jdk，Mac卸载jdk

Mac(M1芯片)安装多个jdk，Mac卸载jdk

1.jdk下载 oracle官方链接&#xff1a;oracle官方下载链接 2.安装 直接下一步&#xff0c;下一步就行 3.查看是否安装成功 出现下图内容表示安装成功。 4.配置环境变量 open -e .bash_profile 路径建议复制过去 #刷新环境变量 source ~/.bash_profile 5.切换方法 6.jdk…
阅读更多...
HTML+CSS 彩色浮雕按钮

HTML+CSS 彩色浮雕按钮

效果演示 实现了一个彩色按钮特效&#xff0c;包括一个按钮&#xff08;button&#xff09;和一个前景色&#xff08;::before&#xff09;。按钮具有四种不同的颜色&#xff0c;当鼠标悬停在按钮上时&#xff0c;前景色会出现渐变效果&#xff0c;并且按钮的颜色、文本阴影和边…
阅读更多...
【研究】AI大模型需要什么样的硬件？

【研究】AI大模型需要什么样的硬件？

关注AI大模型 x 硬件的两条思路 从22年11月OpenAI推出ChatGPT至今&#xff0c;我们看到Chatbot应用的能力不断增强&#xff0c;从最初的文字问答&#xff0c;迅速向具有自主记忆、推理、规划和执行的全自动能力的AI Agent发展。我们认为端侧智能是大模型发展的重要分支。建议投…
阅读更多...
昇思25天学习打卡营第二天|张量

昇思25天学习打卡营第二天|张量

张量 Tensor 张量&#xff08;Tensor&#xff09;是一个可用来表示在一些矢量、标量和其他张量之间的线性关系的多线性函数&#xff0c;这些线性关系的基本例子有内积、外积、线性映射以及笛卡儿积。其坐标在 &#x1d45b;&#x1d45b; 维空间内&#xff0c;有  &#x1…
阅读更多...
北尔Beijer软件iXDeveloper2触摸屏和使用说明手侧

北尔Beijer软件iXDeveloper2触摸屏和使用说明手侧

北尔Beijer软件iXDeveloper2触摸屏和使用说明手侧
阅读更多...
Python笔记 文件的写，追加，备份操作

Python笔记 文件的写，追加，备份操作

一、文件的写操作 案例演示&#xff1a; # 1.打开文件 f open(python.txt,w)# 2.文件写入 f.write(hello world)# 3.内容刷新 f.flush() 注意&#xff1a; 直接调用write&#xff0c;内容并为真正的写入文件&#xff0c;二十会积攒在程序的内存中&#xff0c;称之为缓冲区…
阅读更多...
SpringBoot控制反转和依赖注入

SpringBoot控制反转和依赖注入

目录 一、内聚和耦合 二、分层解耦 三、具体实现 四、bean的组件扫描 五、bean注入 一、内聚和耦合 在了解分层解耦的概念之前我们我们要去先了解一下内聚和耦合。内聚&#xff1a;通常将的是软件中各个模块之间的功能联系。耦合衡量软件各个模块之间的依赖、关联的程度。一…
阅读更多...
Lua网站开发之文件表单上传

Lua网站开发之文件表单上传

这个代码示例演示如何上传文件或图片&#xff0c;获取上传信息及保存文件到本地。 local fw require("fastweb") local request require("fastweb.request") local response require("fastweb.response") local cjson require("cjson&q…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023