wrf模擬的domain圖繪製

xiaofeifeixd發表於2021-03-17

wrf模擬的區域繪製，domain圖，利用python的cartopy庫繪製模擬區域

參考Liang Chen的draw_wrf_domian.py這個程式碼, 出處python畫wrf模式的模擬區域

創新點

區別於Liange程式碼的地方在於使用cartopy庫，替換了basemap庫, 方便在最新的python版本下使用。
初學cartopy,使用cartopy根據距離繪製影像是比較難想到的一點, 是在建立投影物件的那裡設定的你敢信？
picture-bed

具體程式碼（使用cartopy)

"""
Author: Forxd
Time: 2021-3-17
Purpose: read in namelist.wps , draw wrf domain and plot some station
"""
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

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="blue",fc="none", lw=0.6, alpha=1)
    plt.gca().add_patch(poly)

def create_map(info):
    """建立一個包含青藏高原區域的Lambert投影的底圖

    Returns:
        ax: 座標圖物件
    """
    ## --建立畫圖空間
    proj = ccrs.PlateCarree()  # 建立座標系

    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']


    # print(true_lat1)

    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,
                )
    ## 建立座標系
    fig = plt.figure(figsize=(4, 4), dpi=500)  # 建立頁面
    ax = fig.add_axes([0.1,0.1,0.8,0.8], projection=proj_lambert)  

    ## 讀取青藏高原地形檔案
    Tibet = cfeat.ShapelyFeature(
        Reader('/home/fengxiang/Data/shp_tp/Tibet.shp').geometries(),
        proj, edgecolor='k',
        facecolor='none', alpha=0.9
        )

    ## 將青藏高原地形檔案加到地圖中區
    ax.add_feature(Tibet, linewidth=0.5, zorder=2)

    ## --設定網格屬性, 不畫預設的標籤
    gl=ax.gridlines(draw_labels=True,linestyle=":",linewidth=0.3 ,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.xformatter = LONGITUDE_FORMATTER  #使橫座標轉化為經緯度格式            
    gl.yformatter = LATITUDE_FORMATTER                                        
    gl.xlocator=mticker.FixedLocator(np.arange(60,120,10))                              
    gl.ylocator=mticker.FixedLocator(np.arange(10,60,10)) 
    gl.xlabel_style={'size':4}#修改經緯度字型大小                             
    gl.ylabel_style={'size':4}
    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)
    print(false_northing)
    print(false_easting)
    return ax
    


def get_information(flnm):
    """根據namelist.wps檔案，獲取地圖的基本資訊

    Args:
        flnm ([type]): [description]

    Returns:
        [type]: [description]
    """
    ## getting namelist.wps domain information
    name_dict={}

    with open(flnm) 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])

    

    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
    

    dict_return = {
                    "dx":dx,
                    "dy":dy,
                    "max_dom":max_dom,
                    "parent_grid_ratio":parent_grid_ratio,
                    "j_parent_start":j_parent_start,
                    "i_parent_start":i_parent_start,
                    "e_sn":e_sn,
                    "e_we":e_we,
                    'ref_lat':ref_lat,
                    'ref_lon':ref_lon,
                    'true_lat1':truelat1,
                    'true_lat2':truelat2,
                    'parent_grid_ratio':parent_grid_ratio,
                }
    return dict_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[3]*1, lat[3]*1., "d02")


def draw_station():
    station = {'TingRi':{'lat':28.6,'lon':87.0},
                'NaQu':{'lat':31.4, 'lon':92.0},
                'LaSa':{'lat':29.6, 'lon':91.1},
                'TuoTuohe':{'lat':34.2, 'lon':92.4},
                'GaiZe':{'lat':32.3, 'lon':84.0},
                'ShenZha':{'lat':30.9, 'lon':88.7},
                'ShiQuanhe':{'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='red',
                 transform=ccrs.PlateCarree(),
                 linewidth=0.1,s=10) 
    ## 給站點加註釋
    for i in range(len(x)):
        print(x[i])
        plt.text(x[i]-2, y[i]+0.5, station_name[i], 
                transform=ccrs.PlateCarree(),
                fontdict={'size':5,}
                )

if __name__ == '__main__':
    file_folder="./"
    file_name="namelist.wps"
    flnm=file_folder+file_name

    info = get_information(flnm) # 獲取namelist.wps檔案資訊
    ax = create_map(info)  # 在domain1區域內，新增地理資訊，建立底圖
    draw_d02(info)  # 繪製domain2區域
    draw_station()  # 將站點位置繪製到圖上
    plt.title('d01', loc='left')
    plt.savefig("domain.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")

Python繪製直方圖 Pygal模擬擲骰子
2018-04-09
Python直方圖
unity3d表格模擬繪製
2012-10-31
Unity3D
Android模擬器繪製實現概述
2019-02-17
Android
使用joinjs繪製流程圖（五）-流程圖繪製
2024-04-29
JS流程圖
java：繪製圖形
2022-11-23
Java
PyQtGraph之多圖繪製
2024-06-13
QT
繪製雷達圖
2024-07-17
審批流程圖怎樣繪製？繪製流程圖方法有哪些
2020-03-25
流程圖
Python例項：僅繪製圖例而不繪製實際的圖形
2024-06-05
Python
流程圖製作: BPMN流程圖線上繪製
2019-06-04
流程圖
流程圖繪製的操作過程
2017-08-18
流程圖
ios股票K線圖的繪製
2017-01-06
iOS
Canvas 繪製雷達圖
2019-05-28
Canvas
使用css繪製圖形
2019-03-13
CSS
如何繪製Wardley地圖？
2021-05-31
地圖
SVG繪製餅狀圖
2018-10-15
SVG
css繪製特殊圖形
2018-12-11
CSS
EazyDraw 向量圖繪製工具
2022-03-31
ECharts 餅圖繪製教程
2021-10-11
Echarts
PyQtGraph繪製折線圖
2024-06-13
QT
繪製雙Y軸圖
2024-07-17
AnyChart繪製折線圖
2024-05-30
canvas 繪製圖案是重疊繪製而不是重置
2018-09-29
Canvas
canvas繪製圖案是重疊繪製而不是重置
2018-09-29
Canvas
使用java繪圖類Graphics繪製圓圈
2016-03-28
Java繪圖
繪製流程圖的軟體下載
2019-09-17
流程圖
繪製流程圖的簡單軟體
2019-10-28
流程圖
等值線圖的Python繪製方法
2018-10-27
Python
業務流程圖的繪製流程分享(一)
2012-06-03
流程圖
資料流圖的繪製要點
2012-05-23
Excalidraw：繪製圖形的新利器
2024-03-19
Canvas(3)---繪製餅狀圖
2020-04-24
Canvas
繪製圖表 go-echarts
2021-07-09
GoEcharts
Python pyecharts繪製餅圖
2021-12-16
PythonEcharts
MATLAB 繪製折線圖
2020-12-23
Matlab
Qt之自繪製餅圖
2016-11-05
QT
使用CreateJS繪製圖形
2015-09-28
JS
如何繪製業務流程圖(二)
2012-06-26
流程圖

wrf模擬的domain圖繪製

wrf模擬的區域繪製，domain圖，利用python的cartopy庫繪製模擬區域

創新點

具體程式碼（使用cartopy)

具體程式碼(使用basemap)

相關文章