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。分别的功能是：
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) :# 坐标系详细信息 # 兰伯特 坐标系信息 'proj' : 'lcc' ,'lat_1' : truelat1,'lat_2' : truelat2,'lat_0' : ref_lat,'lon_0' : ref_lon,'x_0' : 0 ,'y_0' : 0 ,'datum' : 'WGS84' # 创建平面坐标系 # 计算d01的中心点平面坐标系的坐标 "x_center, y_center" ,x_center, y_center)for i in range(len(e_we)):if
1 :# 针对d01 if i == 0 :# 计算d01的左下角坐标 1 ) / 2 ) * parent_dx1 ) / 2 ) * parent_dy# 计算d01的平面坐标 # 针对d02 else :1 2 ]0 ]# 计算d02的左下角坐标 0 ] + (i_parent_start[i] - 1 ) * parent_dx2 ] + (j_parent_start[i] - 1 ) * parent_dy# 计算d02的经纬度 else :# 针对d03 1 6 ]7 ]2 ]0 ]# 计算d03的左下角坐标 0 ] + (i_parent_start[i] - 1 ) * parent_dx2 ] + (j_parent_start[i] - 1 ) * parent_dy# 计算d03的经纬度 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) :# 计算子网格的左下角坐标 # 计算右上角坐标 1 ) * dx1 ) * dy# 各方向的坐标 2 2 return west, east, south, north, grid_center_lat, grid_center_lon, dx, dypython出图 # 绘制地图 def plot_map (domains, gdf_level1, gdf_level2, gdf_level3, truelat1, truelat2, stand_lon, proj_params) :10 , 10 ), subplot_kw={'projection' : proj})# 绘制shapefile背景 'blue' , facecolor='none' , linewidth=1 ) # 蓝色边框,空心 'red' , facecolor='none' , linewidth=1 ) # 红色边框,空心 'black' , facecolor='#43A7EE' , linewidth=1 ) # 黑色边框,RGB(67,167,238)填充 # 创建坐标系对象 4326 ) # 使用 EPSG 代码 4326 表示 WGS84 地理坐标系 True
# 绘制每个嵌套网格的范围 for i, (west, east, south, north, _, _, _, _) in enumerate(domains):# 将网格边界转换为经纬度 # 左下 # 左上 # 右上 # 右下 # 打印经纬度范围 f"Domain {i+1 } bounds (west, east, south, north):" )f"Longitude: {west_lon_ZUO} to {east_lon_YOU2} " )f"Latitude: {south_lat_ZUO} to {north_lat_ZUO2} " )# 使用 Polygon 创建每个网格的多边形,按照逆时针顺序连接点 'red' if i == 0 else 'blue' if i == 1 else 'green' , facecolor='none' , linewidth=2 , label=f'Domain {i+1 } ' )# 计算标注的位置(使用多边形的右上角点) # 添加标注,并调整标注的位置 f'd0{i+1 } ' , color='red' , fontsize=12 , ha='right' , va='top' , transform=ccrs.PlateCarree())# 转换 d01 的边界坐标到地理坐标 0 ][0 ], domains[0 ][2 ])0 ][1 ], domains[0 ][3 ])"west_d01, south_d01, east_d01, north_d01" , west_d01, south_d01, east_d01, north_d01)# 设置显示范围,在经度和纬度方向上各自添加边距 0.1 # 经度方向上的边距为d01经度范围的10% 0.1 # 纬度方向上的边距为d01纬度范围的10% # 添加海岸线和网格线 True )'WRF Domains' )我们加入研究区的矢量如下，出图效果如下：
namelist转矢量shp 既然我们已经知道d01到d03的坐标点，按照逆时针把矢量点串联起来，获得shp矢量。
# 输出d01到d03范围为shp def output_domains_to_shapefile (domains, proj_params, output_shapefile_path) :# 创建一个空的 GeoDataFrame,用于存储域的范围 'geometry' , 'domain_id' ], crs="EPSG:4326" )# 创建坐标系对象 4326 ) # 使用 EPSG 代码 4326 表示 WGS84 地理坐标系 True )# 绘制每个嵌套网格的范围并添加到 GeoDataFrame for i, (west, east, south, north, _, _, _, _) in enumerate(domains):
# 左下 # 左上 # 右上 # 右下 # 使用 Polygon 创建每个网格的多边形,按照逆时针顺序连接点 # 创建一个临时的 GeoDataFrame 'geometry' : polygon, 'domain_id' : f'Domain {i+1 } ' }], crs="EPSG:4326" )# 使用 pd.concat 将临时的 GeoDataFrame 添加到主要的 GeoDataFrame 中 True )# 保存为 shapefile 'ESRI Shapefile' )特别注意，我们需要最后生成的shp是wgs84坐标系，所以需要把平面坐标转回为wgs84坐标系。
完整代码 import reimport matplotlib.pyplot as pltimport cartopy.crs as ccrsimport geopandas as gpdimport pandas as pdfrom shapely.geometry import Polygonfrom pyproj import Proj,transformfrom pyproj import CRS, Transformer# 提取单个参数的函数 def get_param (pattern, content, index=0 ) :if match:return float(match.group(1 ))else :raise ValueError(f'Parameter {pattern} not found in namelist.wps' )# 提取多个参数的函数 def get_params (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:r'dx\s*=\s*(\d+)' , namelist_content)r'dy\s*=\s*(\d+)' , namelist_content)r'ref_lat\s*=\s*([-+]?\d*\.\d+|\d+)' , namelist_content)r'ref_lon\s*=\s*([-+]?\d*\.\d+|\d+)' , namelist_content)r'e_we\s*=\s*([\d,\s]+)' , namelist_content)r'e_sn\s*=\s*([\d,\s]+)' , namelist_content)r'i_parent_start\s*=\s*([\d,\s]+)' , namelist_content)r'j_parent_start\s*=\s*([\d,\s]+)' , namelist_content)r'parent_id\s*=\s*([\d,\s]+)' , namelist_content)r'parent_grid_ratio\s*=\s*([\d,\s]+)' , namelist_content)r'truelat1\s*=\s*([-+]?\d*\.\d+|\d+)' , namelist_content)r'truelat2\s*=\s*([-+]?\d*\.\d+|\d+)' , namelist_content)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) :# 坐标系详细信息 # 兰伯特 坐标系信息 'proj' : 'lcc' ,'lat_1' : truelat1,'lat_2' : truelat2,'lat_0' : ref_lat,'lon_0' : ref_lon,'x_0' : 0 ,'y_0' : 0 ,'datum' : 'WGS84' # 创建平面坐标系 # 计算d01的中心点平面坐标系的坐标 "x_center, y_center" ,x_center, y_center)for i in range(len(e_we)):if parent_id[i] == 1 :# 针对d01 if i == 0 :# 计算d01的左下角坐标 1 ) / 2 ) * parent_dx1 ) / 2 ) * parent_dy# 计算d01的平面坐标 # 针对d02 else :1 2 ]0 ]# 计算d02的左下角坐标 0 ] + (i_parent_start[i] - 1 ) * parent_dx2 ] + (j_parent_start[i] - 1 ) * parent_dy# 计算d02的经纬度 else :# 针对d03 1 6 ]7 ]2 ]0 ]# 计算d03的左下角坐标 0 ] + (i_parent_start[i] - 1 ) * parent_dx2 ] + (j_parent_start[i] - 1 ) * parent_dy# 计算d03的经纬度 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) :# 计算子网格的左下角坐标 # 计算右上角坐标 1 ) * dx1 ) * dy# 各方向的坐标 2 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) :10 , 10 ), subplot_kw={'projection' : proj})# 绘制shapefile背景 'blue' , facecolor='none' , linewidth=1 ) # 蓝色边框,空心 'red' , facecolor='none' , linewidth=1 ) # 红色边框,空心 'black' , facecolor='#43A7EE' , linewidth=1 ) # 黑色边框,RGB(67,167,238)填充 # 创建坐标系对象 4326 ) # 使用 EPSG 代码 4326 表示 WGS84 地理坐标系 True )# 绘制每个嵌套网格的范围 for i, (west, east, south, north, _, _, _, _) in enumerate(domains):# 将网格边界转换为经纬度 # 左下 # 左上 # 右上 # 右下 # 打印经纬度范围 f"Domain {i+1 } bounds (west, east, south, north):" )f"Longitude: {west_lon_ZUO} to {east_lon_YOU2} " )f"Latitude: {south_lat_ZUO} to {north_lat_ZUO2} " )# 使用 Polygon 创建每个网格的多边形,按照逆时针顺序连接点 'red' if i == 0 else 'blue' if i == 1 else 'green' , facecolor='none' , linewidth=2 , label=f'Domain {i+1 } ' )# 计算标注的位置(使用多边形的右上角点) # 添加标注,并调整标注的位置 f'd0{i+1 } ' , color='red' , fontsize=12 , ha='right' , va='top' , transform=ccrs.PlateCarree())# 转换 d01 的边界坐标到地理坐标 0 ][0 ], domains[0 ][2 ])0 ][1 ], domains[0 ][3 ])"west_d01, south_d01, east_d01, north_d01" , west_d01, south_d01, east_d01, north_d01)# 设置显示范围,在经度和纬度方向上各自添加边距 0.1 # 经度方向上的边距为d01经度范围的10% 0.1 # 纬度方向上的边距为d01纬度范围的10% # 添加海岸线和网格线 True )'WRF Domains' )# 输出d01到d03范围为shp def output_domains_to_shapefile (domains, proj_params, output_shapefile_path) :# 创建一个空的 GeoDataFrame,用于存储域的范围 'geometry' , 'domain_id' ], crs="EPSG:4326" )# 创建坐标系对象 4326 ) # 使用 EPSG 代码 4326 表示 WGS84 地理坐标系 True )# 绘制每个嵌套网格的范围并添加到 GeoDataFrame for i, (west, east, south, north, _, _, _, _) in enumerate(domains):# 将网格边界转换为经纬度 # 左下 # 左上 # 右上 # 右下 # 使用 Polygon 创建每个网格的多边形,按照逆时针顺序连接点 # 创建一个临时的 GeoDataFrame 'geometry' : polygon, 'domain_id' : f'Domain {i+1 } ' }], crs="EPSG:4326" )# 使用 pd.concat 将临时的 GeoDataFrame 添加到主要的 GeoDataFrame 中 True )# 保存为 shapefile 'ESRI Shapefile' )def main () :# 读取namelist.wps文件 r"E:\ruiduobao\namelis设置\namelist.wps" # 初始化网格域 # 读取shapefile r'E:\ruiduobao\数据和代码\行政区划\jiangsu.shp' r'E:\ruiduobao\数据和代码\行政区划\xuzhou.shp' r'E:\ruiduobao\数据和代码\行政区划\xuzhouxian.shp' # 加载shapefile到 GeoDataFrame # 绘制地图 # 输出d01到d03范围为shp r'E:\ruiduobao\行政区划\wrf_domains_平面.shp' if __name__ == '__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
声明： 欢迎转载、转发。气象学家公众号转载信息旨在传播交流，其内容由作者负责，不代表本号观点。文中部分图片来源于网络，如涉及内容、版权和其他问题，请联系小编(微信：qxxjgzh) 处理。
