爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 3092|回复: 3

[经验总结] python作气象站点资料的风速矢量图

[复制链接]

新浪微博达人勋

发表于 2023-5-9 18:27:40 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 menseye 于 2023-5-9 18:27 编辑

最近,我想作气象站点资料的风速矢量图,从而与模式结果对比,但是网上只找到根据气象站点资料作散点图的教程,还要利用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()



风矢量图

风矢量图

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


密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2023-5-9 20:40:20 | 显示全部楼层
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   #这两行需要手动设置
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2023-5-10 16:31:55 | 显示全部楼层
edwardli 发表于 2023-5-9 20:40
1、一维lon、lat生成二维网格lons、lats并无不妥。
但,实际有lon、lat足矣。

感谢老师,可能是我之前考虑的太复杂了,我是按照处理再分析资料的思路来处理站点CSV数据,现在看来确实没有必要作二维数组。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-6-30 15:30:55 | 显示全部楼层
感谢 风向转换这里没想那么仔细就卡了 刚好看到你帖子{:5_213:}{:5_213:}{:5_213:}{:5_213:}
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表