python作气象站点资料的风速矢量图

menseye

最近，我想作气象站点资料的风速矢量图，从而与模式结果对比，但是网上只找到根据气象站点资料作散点图的教程，还要利用arcgis的方法，但是我不会用arcgis，遂自己琢磨出方法，我的思路是：先将站点资料里的一维的经度和纬度转为二维数组，然后将风速风向转为u,v分量，再将一维的u,v转变为对角阵，最后就可以画出风速矢量图了，代码如下：

####气象站观测数据可视化##########
from metpy.plots import SkewT
import math
import os
import numpy as np
import pandas as pd
#from datetime import datetime
import xarray as xr
#from scipy.ndimage import gaussian_filter
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import metpy.calc as mpcalc
from metpy.units import units
#from metpy.units import units
import  netCDF4 as nc
import metpy.constants as constants
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
#from mpl_toolkits.basemap import Basemap

plt.rcParams['font.sans-serif']=['SimHei'] #显示中文标签
plt.rcParams['axes.unicode_minus']=False   #这两行需要手动设置

#################打开地图文件，大家可以用shp文件，效果一样######################
with open(r'D:\code\china-geospatial-data-GB2312\china-geospatial-data-GB2312\CN-border-La.gmt') as src:
    context = src.read()
    blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]
    borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]
####################################################################


###风速风向转为UV分量#########
def wsd2uv(ws, wd):

    wd = 270.0- wd
    wd = wd /180 *np.pi
    x = ws * np.cos(wd)
    y = ws * np.sin(wd)
    return(x, y)

####打开站点数据csv文件##########
data = pd.read_csv('D:\\DATA\\observation\\test_python.csv',sep=',',header=1,
                   names=['WIN_S_INST',        'STATION_ID_C',        'LON','WIN_D_INST','WIN_S_AVG_2MI',
                          'HOUR','WIN_D_AVG_10MI','WIN_S_AVG_10MI',        'PRE_1H','RHU',        'MON','DAY',        'TEM','LAT'])

ws=data['WIN_S_INST'].values
wd=data['WIN_D_INST'].values
wd=wd.astype(float)            ###将风向转变为float类型
wd[(wd==999017)]=None    ###设置缺测值为空值
ws[(ws==999017)]=None
lon=data['LON'].values
lat=data['LAT'].values
lons,lats=np.meshgrid(lon,lat)       ###将经纬度数组化   

(u,v)=wsd2uv(ws, wd)     ###调用函数
u=np.diag(u)                  ###
u[(u==0.000)]=None
# u[(u==None)]=0.
v=np.diag(v)
v[(v==0.000)]=None   

#####画图部分###########
datacrs=ccrs.PlateCarree()
fig=plt.figure(figsize=(9,9))
plt.rc('font',family='Times New Roman')
plt.rcParams['font.sans-serif']=['SimHei'] #显示中文标签
plt.rcParams['axes.unicode_minus']=False   #这两行需要手动设置
ax=fig.add_axes([0.06,0.05,0.8,0.8],projection=datacrs)
# rejion=[113,121,35,43]
rejion=[113.7,118.5,37.3,41.2]

ax.set_extent(rejion,datacrs)
clevers=np.arange(0,3000,100)


cf=ax.scatter(lon,lat,s =20,c=data["RHU"],cmap='jet',transform=ccrs.PlateCarree())  ####画出站点资料的相对湿度散点图，这儿是用的一维经纬度
cax = fig.add_axes([0.90,0.05,0.02,0.8]) # bar的举行描述（归一化后的尺度
cb = plt.colorbar(cf,cax=cax)
########画省界，可用shp文件代替##########################
for line in borders:
    ax.plot(line[0::2], line[1::2], '-', color='gold',linewidth=0.8,transform=ccrs.Geodetic())
#############################################

uv=ax.quiver(lons,lats, u, v, color='r',units='dots',scale=0.04, transform=datacrs) ########使用二维数组的经纬度和对角阵后的u,v分量
plt.quiverkey(uv,0.93,1.01,2,label='2m/s',transform=datacrs)
  


ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.set_title('obs_wind',fontsize=20)
gl=ax.gridlines(xlocs=np.arange(113,122,1),ylocs=np.arange(35,44,1),draw_labels=True,
linewidth=1, color='k', alpha=0.5, linestyle='--',crs=datacrs)
gl.xformatter = LONGITUDE_FORMATTER  #x轴设为经度格式
gl.yformatter = LATITUDE_FORMATTER  #y轴设为纬度格式
gl.top_labels = False  #关闭顶端标签
gl.right_labels = False  #关闭右侧标签
gl.xlabel_style = {'fontsize':14}
gl.ylabel_style = {'fontsize':14}
plt.show()

风矢量图

图中风矢量与相对湿度的散点重合，基本说明该方法的正确性。
若有问题，欢迎批评指正！

edwardli

1、一维lon、lat生成二维网格lons、lats并无不妥。
但，实际有lon、lat足矣。

2、风向风速→u、v分解
自己写没有问题，但会出问题。建议使用
from metpy.units import units
from metpy.calc import wind_components
df_w['u'], df_w['v'] = wind_components((df_w['ws'].values * units('m/s')),df_w['wd'].values * units.degree)
不是代码行数问题，是“不再重复造轮子”的问题。
3、所谓的“对角矩阵”完全是多此一举，只要是quiver（x,y,u,v），每个量都是一维矩阵array-like。

4、pd.read_csv，默认的分隔符就是逗号，无需赋sep关键字参数。



##### 绘制极大风速【因为只有这个需要“方向”！！！】
# from metpy.plots import StationPlot
# stationplot = StationPlot(ax, df_w['lon'], df_w['lat'], transform=proj,fontsize=4.5)
# stationplot.plot_barb(df_w['u'],df_w['v'],transform=proj,barb_increments={'half':2, 'full':4, 'flag':20},
#                       sizes=dict(spacing=0.15,height=0.35,width=0.3),linewidth=0.5,barbcolor=list(dfw['color'].values))
ax.barbs(df_w['lon'], df_w['lat'], df_w['u'], df_w['v'], transform = proj,
         barb_increments = {'half':2, 'full':4, 'flag':20},
         sizes = dict(spacing=0.1,height=0.3,width=0.3),
         linewidth = 0.3, barbcolor=  list(df_w['color']))

5、另from metpy.plots import SkewT
import math
import os等没有使用的包要舍弃。

6、另另，没有中文输出，此两行可简（shan）化（chu）

plt.rcParams['font.sans-serif']=['SimHei'] #显示中文标签
plt.rcParams['axes.unicode_minus']=False   #这两行需要手动设置

menseye

edwardli 发表于 2023-5-9 20:40
1、一维lon、lat生成二维网格lons、lats并无不妥。
但，实际有lon、lat足矣。

感谢老师，可能是我之前考虑的太复杂了，我是按照处理再分析资料的思路来处理站点CSV数据，现在看来确实没有必要作二维数组。

guanqiby

感谢 风向转换这里没想那么仔细就卡了 刚好看到你帖子{:5_213:}{:5_213:}{:5_213:}{:5_213:}