wps的domain转为shp矢量

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

简介

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

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