- 积分
- 341
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2019-11-6
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
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()
风矢量图
图中风矢量与相对湿度的散点重合,基本说明该方法的正确性。
若有问题,欢迎批评指正!
|
|