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")