wps的domain转为shp矢量

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:# 针对d01if i == 0:parent_dx = dxparent_dy = dyparent_ref_lon = x_centerparent_ref_lat = y_centergrid_ratio = parent_grid_ratio[i]# 计算d01的左下角坐标d01_start_x = x_center - ((e_we[i] - 1) / 2) * parent_dxd01_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))# 针对d02else:parent_domain_idx = parent_id[i] - 1grid_ratio = parent_grid_ratio[i]parent_dx = dxparent_dy = dydx = dx / grid_ratiody = dy / grid_ratioparent_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_dxd02_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:# 针对d03parent_domain_idx = parent_id[i] - 1grid_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_dxd03_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_xgrid_start_lat = d_start_y# 计算右上角坐标grid_end_lon = grid_start_lon + (e_we - 1) * dxgrid_end_lat = grid_start_lat + (e_sn - 1) * dy# 各方向的坐标west = grid_start_lonsouth = grid_start_lateast = grid_end_lonnorth = grid_end_latgrid_center_lat = (south + north) / 2grid_center_lon = (west + east) / 2return 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)# 绘制每个嵌套网格的范围并添加到 GeoDataFramefor 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)# 创建一个临时的 GeoDataFrametemp_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)# 保存为 shapefilegdf_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:# 针对d01if i == 0:parent_dx = dxparent_dy = dyparent_ref_lon = x_centerparent_ref_lat = y_centergrid_ratio = parent_grid_ratio[i]# 计算d01的左下角坐标d01_start_x = x_center - ((e_we[i] - 1) / 2) * parent_dxd01_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))# 针对d02else:parent_domain_idx = parent_id[i] - 1grid_ratio = parent_grid_ratio[i]parent_dx = dxparent_dy = dydx = dx / grid_ratiody = dy / grid_ratioparent_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_dxd02_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:# 针对d03parent_domain_idx = parent_id[i] - 1grid_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_dxd03_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_xgrid_start_lat = d_start_y# 计算右上角坐标grid_end_lon = grid_start_lon + (e_we - 1) * dxgrid_end_lat = grid_start_lat + (e_sn - 1) * dy# 各方向的坐标west = grid_start_lonsouth = grid_start_lateast = grid_end_lonnorth = grid_end_latgrid_center_lat = (south + north) / 2grid_center_lon = (west + east) / 2return 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)# 绘制每个嵌套网格的范围并添加到 GeoDataFramefor 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)# 创建一个临时的 GeoDataFrametemp_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)# 保存为 shapefilegdf_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)# 读取shapefileshapefile_path_level1 = r'E:\ruiduobao\数据和代码\行政区划\jiangsu.shp'shapefile_path_level2 = r'E:\ruiduobao\数据和代码\行政区划\xuzhou.shp'shapefile_path_level3 = r'E:\ruiduobao\数据和代码\行政区划\xuzhouxian.shp'# 加载shapefile到 GeoDataFramegdf_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范围为shpoutput_shapefile_path = r'E:\ruiduobao\行政区划\wrf_domains_平面.shp'output_domains_to_shapefile(domains, proj_params, output_shapefile_path)if __name__ == '__main__':main()