wrf模拟的domain图绘制
wrf模拟的区域绘制，domain图，利用python的cartopy库绘制模拟区域
参考Liang Chen的draw_wrf_domian.py这个代码, 出处python画wrf模式的模拟区域
创新点
区别于Liange代码的地方在于使用cartopy库，替换了basemap库, 方便在最新的python版本下使用。
初学cartopy,使用cartopy根据距离绘制图像是比较难想到的一点, 是在创建投影对象的那里设置的你敢信？
具体代码（使用cartopy)
"""
Author: Forxd
Time: 2021-3-17
Purpose: read in namelist.wps , draw wrf domain and plot some station
"""
#!/home/fengxiang/anaconda3/envs/wrfout/bin/python
# -*- encoding: utf-8 -*-
'''
Description:
read in namelist.wps , draw wrf domain and plot some station
-----------------------------------------
Time :2021/03/28 17:28:59
Author :Forxd
Version :1.0
'''
import xarray as xr
import numpy as np
import salem
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.io.shapereader import Reader, natural_earth
import cartopy.feature as cf
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import geopandas
import cmaps
import re
from matplotlib.path import Path
import matplotlib.patches as patches
def draw_screen_poly(lats, lons):
'''
lats: 纬度列表
lons: 经度列表
purpose: 画区域直线
'''
x, y = lons, lats
xy = list(zip(x, y))
# print(xy)
poly = plt.Polygon(xy, edgecolor="black", fc="none", lw=.8, alpha=1)
plt.gca().add_patch(poly)
def create_map(info):
"""创建一个包含青藏高原区域的Lambert投影的底图
Returns:
ax: 坐标图对象
"""
## --创建画图空间
ref_lat = info['ref_lat']
ref_lon = info['ref_lon']
true_lat1 = info['true_lat1']
true_lat2 = info['true_lat2']
false_easting = (info['e_we'][0] - 1) / 2 * info['dx']
false_northing = (info['e_sn'][0] - 1) / 2 * info['dy']
proj_lambert = ccrs.LambertConformal(
central_longitude=ref_lon,
central_latitude=ref_lat,
standard_parallels=(true_lat1, true_lat2),
cutoff=-30,
false_easting=false_easting,
false_northing=false_northing,
)
# proj = ccrs.PlateCarree(central_longitude=ref_lon) # 创建坐标系
proj = ccrs.PlateCarree() # 创建坐标系
## 创建坐标系
fig = plt.figure(figsize=(6, 5), dpi=300) # 创建页面
ax = fig.add_axes([0.08, 0.02, 0.85, 0.95], projection=proj_lambert)
## 读取青藏高原地形文件
Tibet = cfeat.ShapelyFeature(
Reader('/home/fengxiang/Data/shp_tp/Tibet.shp').geometries(),
proj,
edgecolor='black',
lw=1.,
linewidth=1.,
facecolor='none',
alpha=1.)
## 将青藏高原地形文件加到地图中区
ax.add_feature(Tibet, linewidth=0.5, zorder=2)
ax.coastlines()
## --设置网格属性, 不画默认的标签
gl = ax.gridlines(draw_labels=True,
dms=True,
linestyle=":",
linewidth=0.3,
x_inline=False,
y_inline=False,
color='k')
# # gl=ax.gridlines(draw_labels=True,linestyle=":",linewidth=0.3 , auto_inline=True,x_inline=False, y_inline=False,color='k')
## 关闭上面和右边的经纬度显示
gl.top_labels = False #关闭上部经纬标签
# gl.bottom_labels = False
# # gl.left_labels = False
gl.right_labels = False
## 这个东西还挺重要的，对齐坐标用的
gl.rotate_labels = None
gl.xformatter = LONGITUDE_FORMATTER #使横坐标转化为经纬度格式
gl.yformatter = LATITUDE_FORMATTER
gl.xlocator = mticker.FixedLocator(np.arange(55, 130, 10))
gl.ylocator = mticker.FixedLocator(np.arange(10, 55, 5))
gl.xlabel_style = {'size': 10} #修改经纬度字体大小
gl.ylabel_style = {'size': 10}
ax.spines['geo'].set_linewidth(0.6) #调节边框粗细
# ax.set_extent([60, 120, 10, 60], crs=proj)
# ax.set_extent([0, 2237500*2, 0, 1987500*2], crs=proj_lambert)
ax.set_extent([0, false_easting * 2, 0, false_northing * 2],
crs=proj_lambert)
# 标注d01, 这个位置需要根据经纬度手动调整
ax.text(65,
50,
'd01',
transform=ccrs.PlateCarree(),
fontdict={
'size': 14,
})
return ax
def get_information(flnm):
"""根据namelist.wps文件，获取地图的基本信息
Args:
flnm ([type]): [description]
Returns:
[type]: [description]
"""
## 设置正则表达式信息
pattern = {}
pattern['dx'] = 'dx\s*=\s*\d*,'
pattern['dy'] = 'dy\s*=\s*\d*,'
pattern['max_dom'] = 'max_dom\s*=\s*\d\s*,'
pattern[
'parent_grid_ratio'] = 'parent_grid_ratio\s*=\s*\d,\s*\d,\s*\d,\s*\d,'
pattern['j_parent_start'] = 'j_parent_start\s*=\s*\d,\s*\d*,\s*\d*,\s*\d*,'
pattern['i_parent_start'] = 'i_parent_start\s*=\s*\d,\s*\d*,\s*\d*,\s*\d*,'
pattern['e_sn'] = 'e_sn\s*=\s*\d*,\s*\d*,\s*\d*,\s*\d*'
pattern['e_we'] = 'e_we\s*=\s*\d*,\s*\d*,\s*\d*,\s*\d*'
pattern['ref_lat'] = 'ref_lat\s*=\s*\d*.?\d*,'
pattern['ref_lon'] = 'ref_lon\s*=\s*\d*.?\d*,'
pattern['true_lat1'] = 'truelat1\s*=\s*\d*.?\d*,'
pattern['true_lat2'] = 'truelat2\s*=\s*\d*.?\d*,'
f = open(flnm)
fr = f.read()
def get_var(var, pattern=pattern, fr=fr):
"""处理正则表达式得到的数据"""
ff1 = re.search(pattern[var], fr, flags=0)
str_f1 = ff1.group(0)
str1 = str_f1.replace('=', ',')
aa = str1.split(',')
bb = []
for i in aa[1:]:
if i != '':
bb.append(i.strip())
return bb
dic_return = {}
aa = get_var('parent_grid_ratio')
var_list = [
'dx',
'dy',
'max_dom',
'parent_grid_ratio',
'j_parent_start',
'i_parent_start',
'e_sn',
'e_we',
'ref_lat',
'ref_lon',
'true_lat1',
'true_lat2',
]
for i in var_list:
aa = get_var(i)
if i in [
'parent_grid_ratio',
'j_parent_start',
'i_parent_start',
'e_we',
'e_sn',
]:
bb = aa
bb = [float(i) for i in bb]
else:
bb = float(aa[0])
dic_return[i] = bb
return dic_return
def draw_d02(info):
"""绘制domain2
Args:
info ([type]): [description]
"""
max_dom = info['max_dom']
dx = info['dx']
dy = info['dy']
i_parent_start = info['i_parent_start']
j_parent_start = info['j_parent_start']
parent_grid_ratio = info['parent_grid_ratio']
e_we = info['e_we']
e_sn = info['e_sn']
if max_dom >= 2:
### domain 2
# 4 corners 找到四个顶点和距离相关的坐标
ll_lon = dx * (i_parent_start[1] - 1)
ll_lat = dy * (j_parent_start[1] - 1)
ur_lon = ll_lon + dx / parent_grid_ratio[1] * (e_we[1] - 1)
ur_lat = ll_lat + dy / parent_grid_ratio[1] * (e_sn[1] - 1)
lon = np.empty(4)
lat = np.empty(4)
lon[0], lat[0] = ll_lon, ll_lat # lower left (ll)
lon[1], lat[1] = ur_lon, ll_lat # lower right (lr)
lon[2], lat[2] = ur_lon, ur_lat # upper right (ur)
lon[3], lat[3] = ll_lon, ur_lat # upper left (ul)
draw_screen_poly(lat, lon) # 画多边型
## 标注d02
# plt.text(lon[0] * 1+100000, lat[0] * 1. - 225000, "d02", fontdict={'size':14})
plt.text(lon[2] * 1 - 440000,
lat[2] * 1. - 200000,
"d02",
fontdict={'size': 14})
if max_dom >= 3:
### domain 3
## 4 corners
ll_lon += dx / parent_grid_ratio[1] * (i_parent_start[2] - 1)
ll_lat += dy / parent_grid_ratio[1] * (j_parent_start[2] - 1)
ur_lon = ll_lon + dx / parent_grid_ratio[1] / parent_grid_ratio[2] * (
e_we[2] - 1)
ur_lat = ll_lat + dy / parent_grid_ratio[1] / parent_grid_ratio[2] * (
e_sn[2] - 1)
## ll
lon[0], lat[0] = ll_lon, ll_lat
## lr
lon[1], lat[1] = ur_lon, ll_lat
## ur
lon[2], lat[2] = ur_lon, ur_lat
## ul
lon[3], lat[3] = ll_lon, ur_lat
draw_screen_poly(lat, lon)
if max_dom >= 4:
### domain 3
## 4 corners
ll_lon += dx / parent_grid_ratio[1] / parent_grid_ratio[2] * (
i_parent_start[3] - 1)
ll_lat += dy / parent_grid_ratio[1] / parent_grid_ratio[2] * (
j_parent_start[3] - 1)
ur_lon = ll_lon + dx / parent_grid_ratio[1] / parent_grid_ratio[
2] / parent_grid_ratio[3] * (e_we[3] - 1)
ur_lat = ll_lat + dy / parent_grid_ratio[1] / parent_grid_ratio[
2] / parent_grid_ratio[3] * (e_sn[3] - 1)
## ll
lon[0], lat[0] = ll_lon, ll_lat
## lr
lon[1], lat[1] = ur_lon, ll_lat
## ur
lon[2], lat[2] = ur_lon, ur_lat
## ul
lon[3], lat[3] = ll_lon, ur_lat
draw_screen_poly(lat, lon)
def draw_station():
"""画站点
"""
station = {
# 'TR': {
# 'lat': 28.6,
# 'lon': 87.0
# },
'NQ': {
'lat': 31.4,
'lon': 92.0
},
'LS': {
'lat': 29.6,
'lon': 91.1
},
'TTH': {
'lat': 34.2,
'lon': 92.4
},
'GZ': {
'lat': 32.3,
'lon': 84.0
},
'SZ': {
'lat': 30.9,
'lon': 88.7
},
'SQH': {
'lat': 32.4,
'lon': 80.1
},
# 'JinChuan': {
# 'lat': 31.29,
# 'lon': 102.04
# },
# 'JinLong': {
# 'lat': 29.00,
# 'lon': 101.50
# },
}
values = station.values()
station_name = list(station.keys())
# print(type(station_name[0]))
# print(station_name[0])
x = []
y = []
for i in values:
y.append(float(i['lat']))
x.append(float(i['lon']))
## 标记出站点
ax.scatter(x,
y,
color='black',
transform=ccrs.PlateCarree(),
linewidth=0.2,
s=12)
## 给站点加注释
for i in range(len(x)):
# print(x[i])
if station_name[i] == 'LS':
# x[i] = x[i]
y[i] = y[i] - 2.0
if station_name[i] == 'SQH':
x[i] = x[i] - 0.5
plt.text(x[i] - 1,
y[i] + 0.5,
station_name[i],
transform=ccrs.PlateCarree(),
fontdict={
'size': 9,
})
if __name__ == '__main__':
file_folder = "./"
file_name = "namelist.wps_l"
flnm = file_folder + file_name
info = get_information(flnm) # 获取namelist.wps文件信息
# print(info['ref_lat'])
ax = create_map(info) # 在domain1区域内，添加地理信息，创建底图
draw_d02(info) # 绘制domain2区域
plt.savefig("domain_lyh.png")
具体代码(使用basemap)
'''
File name: draw_wrf_domain.py
Author: Liang Chen
E-mail: chenliang@tea.ac.cn
Date created: 2016-12-22
Date last modified: 2021-3-3
##############################################################
Purpos:
this function reads in namelist.wps and plot the wrf domain
'''
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, cm
from matplotlib.colors import LinearSegmentedColormap
import shapefile
from matplotlib.collections import LineCollection
import matplotlib.colors
import sys
import numpy as np
def draw_screen_poly( lats, lons):
'''
lats: 纬度列表
lons: 经度列表
purpose: 画区域直线
'''
x, y = lons, lats
xy = list(zip(x,y))
print(xy)
poly = plt.Polygon( xy, edgecolor="blue",fc="none", lw=2, alpha=1)
plt.gca().add_patch(poly)
sShapeFiles="/home/fengxiang/Data/shp_tp/"
shape_line=['Tibet.shp',]
## setting namelist.wps domain information
file_folder="./"
file_name="namelist.wps"
sfile=file_folder+file_name
name_dict={}
with open(sfile) as fr:
for line in fr:
if "=" in line: # 这里没有考虑注释的那些行吧, 不过wps一般也没人注释就是了
line=line.replace("=","").replace(",","")
name_dict.update({line.split()[0]: line.split()[1:]}) # 这个字典直接可以更新
dx = float(name_dict["dx"][0])
dy = float(name_dict["dy"][0])
max_dom = int(name_dict["max_dom"][0])
print(max_dom)
parent_grid_ratio = list(map(int, name_dict["parent_grid_ratio"]))
i_parent_start = list(map(int, name_dict["i_parent_start"]))
j_parent_start = list(map(int, name_dict["j_parent_start"]))
e_sn = list(map(int, name_dict["e_sn"]))
e_we = list(map(int, name_dict["e_we"]))
ref_lat= float(name_dict["ref_lat"][0]) # 模式区域中心位置
ref_lon= float(name_dict["ref_lon"][0])
truelat1 = float(name_dict["truelat1"][0]) # 和投影相关的经纬度
truelat2 = float(name_dict["truelat2"][0])
# # ## draw map
fig = plt.figure(figsize=(7,6)) # 设置画板大小
#Custom adjust of the subplots
plt.subplots_adjust(left=0.05,right=0.97,top=0.9,bottom=0.1) # 调整画布大小
ax = plt.subplot(111)
m = Basemap(resolution="l", projection="lcc", rsphere=(6370000.0, 6370000.0), lat_1=truelat1, lat_2=truelat2, lat_0=ref_lat, lon_0=ref_lon, width=dx*(e_we[0]-1), height=dy*(e_sn[0]-1))
# m.drawcoastlines()
#m.drawcountries(linewidth=2)
#m.drawcountries()
#m.fillcontinents()
#m.fillcontinents(color=(0.8,1,0.8))
#m.drawmapboundary()
#m.fillcontinents(lake_color="aqua")
#m.drawmapboundary(fill_color="aqua")
### 根据地形文件，画底图
ii=0 # 控制变量
for sr in shape_line:
# print(sr)
r = shapefile.Reader(sShapeFiles+sr) # 读地形文件
shapes = r.shapes()
records = r.records()
for record, shape in zip(records,shapes):
lons,lats = zip(*shape.points)
data = np.array(m(lons, lats)).T
if len(shape.parts) == 1:
segs = [data,]
else:
segs = []
for i in range(1,len(shape.parts)):
index = shape.parts[i-1]
index2 = shape.parts[i]
segs.append(data[index:index2])
segs.append(data[index2:])
lines = LineCollection(segs,antialiaseds=(1,))
# lines.set_facecolors(cm.jet(np.random.rand(1)))
if ii==0:
lines.set_edgecolors('black')
lines.set_linewidth(2)
else:
lines.set_edgecolors('k')
lines.set_linewidth(1)
ax.add_collection(lines)
ii=ii+1
## 画标签
m.drawparallels(np.arange(-90, 90, 10), labels = [1,0,0,0], fontsize=16,dashes=[1,1])
# m.drawmeridians(np.arange(-180, 180, 10), labels = [0,0,0,1], fontsize=16,dashes=[1,1])
print(ref_lat, ref_lon)
## plot center position 画中心点
cenlon= np.arange(max_dom); cenlat=np.arange(max_dom)
cenlon_model=dx*(e_we[0]-1)/2.0
cenlat_model=dy*(e_sn[0]-1)/2.0
cenlon[0], cenlat[0]=m(cenlon_model, cenlat_model, inverse=True)
## 画区域1的中点和标注
plt.plot(cenlon_model,cenlat_model, marker="o", color="gray")
plt.text(cenlon_model*0.8, cenlat_model*1.01, "({cenlat}, {cenlon})".format(cenlat=round(cenlat[0],2), cenlon=round(cenlon[0],2)))
#### draw nested domain rectangle
#### 区域2
#### 画多边形
lon=np.arange(4); lat=np.arange(4)
if max_dom >= 2:
### domain 2
# 4 corners
ll_lon = dx*(i_parent_start[1]-1)
ll_lat = dy*(j_parent_start[1]-1)
ur_lon = ll_lon + dx/parent_grid_ratio[1] * (e_we[1]-1)
ur_lat = ll_lat + dy/parent_grid_ratio[1] * (e_sn[1]-1)
## lower left (ll)
lon[0],lat[0] = ll_lon, ll_lat
## lower right (lr)
lon[1],lat[1] = ur_lon, ll_lat
## upper right (ur)
lon[2],lat[2] = ur_lon, ur_lat
## upper left (ul)
lon[3],lat[3] = ll_lon, ur_lat
print(lat)
print(lon)
draw_screen_poly(lat, lon) # 画多边型
## 标注d02
plt.text(lon[3]*1, lat[3]*1., "d02")
### 区域2画多边形中点
cenlon_model = ll_lon + (ur_lon-ll_lon)/2.0
cenlat_model = ll_lat + (ur_lat-ll_lat)/2.0
cenlon[1], cenlat[1]=m(cenlon_model, cenlat_model, inverse=True)
# plt.plot(cenlon_model, cenlat_model,marker="o") # 这个画的是区域2的中点
# plt.text(cenlon_model*0.8, cenlat_model*1.01, "({cenlat}, {cenlon})".format(cenlat=round(cenlat[1],2), cenlon=round(cenlon[1],2)))
if max_dom >= 3:
### domain 3
## 4 corners
ll_lon += dx/parent_grid_ratio[1]*(i_parent_start[2]-1)
ll_lat += dy/parent_grid_ratio[1]*(j_parent_start[2]-1)
ur_lon = ll_lon +dx/parent_grid_ratio[1]/parent_grid_ratio[2]*(e_we[2]-1)
ur_lat =ll_lat+ dy/parent_grid_ratio[1]/parent_grid_ratio[2]*(e_sn[2]-1)
## ll
lon[0],lat[0] = ll_lon, ll_lat
## lr
lon[1],lat[1] = ur_lon, ll_lat
## ur
lon[2],lat[2] = ur_lon, ur_lat
## ul
lon[3],lat[3] = ll_lon, ur_lat
draw_screen_poly(lat, lon)
plt.text(lon[0]-lon[0]/10,lat[0]-lat[0]/10,"({i}, {j})".format(i=i_parent_start[2], j=j_parent_start[2]))
#plt.plot(lon,lat,linestyle="",marker="o",ms=10)
cenlon_model = ll_lon + (ur_lon-ll_lon)/2.0
cenlat_model = ll_lat + (ur_lat-ll_lat)/2.0
# plt.plot(cenlon,cenlat,marker="o",ms=15)
#print m(cenlon, cenlat)cenlon, cenlat, ll_lon, ll_lat, ur_lon, ur_lat
#print m(cenlon, cenlat,inverse=True)
cenlon[2], cenlat[2]=m(cenlon_model, cenlat_model, inverse=True)
if max_dom >= 4:
### domain 3
## 4 corners
ll_lon += dx/parent_grid_ratio[1]/parent_grid_ratio[2]*(i_parent_start[3]-1)
ll_lat += dy/parent_grid_ratio[1]/parent_grid_ratio[2]*(j_parent_start[3]-1)
ur_lon = ll_lon +dx/parent_grid_ratio[1]/parent_grid_ratio[2]/parent_grid_ratio[3]*(e_we[3]-1)
ur_lat =ll_lat+ dy/parent_grid_ratio[1]/parent_grid_ratio[2]/parent_grid_ratio[3]*(e_sn[3]-1)
## ll
lon[0],lat[0] = ll_lon, ll_lat
## lr
lon[1],lat[1] = ur_lon, ll_lat
## ur
lon[2],lat[2] = ur_lon, ur_lat
## ul
lon[3],lat[3] = ll_lon, ur_lat
draw_screen_poly(lat, lon)
#plt.plot(lon,lat,linestyle="",marker="o",ms=10)
cenlon_model = ll_lon + (ur_lon-ll_lon)/2.0
cenlat_model = ll_lat + (ur_lat-ll_lat)/2.0
# plt.plot(cenlon,cenlat,marker="o",ms=15)
#print m(cenlon, cenlat)cenlon, cenlat, ll_lon, ll_lat, ur_lon, ur_lat
#print m(cenlon, cenlat,inverse=True)
cenlon[3], cenlat[3]=m(cenlon_model, cenlat_model, inverse=True)
## 标注站点
plt.plot(cenlon_model, cenlat_model,marker="o") # 这个画的是区域2的中点
print(cenlon_model/25000, cenlat_model/25000)
# plt.text(cenlon_model*0.8, cenlat_model*1.01, "({cenlat}, {cenlon})".format(cenlat=round(cenlat[1],2), cenlon=round(cenlon[1],2)))
cenlon_model=dx*(e_we[0]-1)/2.0
print(dx)
print(dy)
Tingri={'lat':28.6,'lon':87.0,'name':'Tingri'}
plt.plot(Tingri['lon']*25000, Tingri['lat']*25000,marker="o")
# plt.text(Tingri['lon']*0.8, Tingri['lat']*1.01, "({cenlat}, {cenlon})".format(cenlat=round(cenlat[1],2), cenlon=round(cenlon[1],2)))
plt.savefig("tttt.png")
