基于气象数据ERA5的Python绘制风场

最新推荐文章于 2024-04-24 21:49:12 发布

keila7

最新推荐文章于 2024-04-24 21:49:12 发布

阅读量3.5k

点赞数 1

文章标签： python numpy 深度学习

本文链接：https://blog.csdn.net/weixin_43536702/article/details/127426357

版权

import cartopy.mpl.ticker as cticker
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
from cartopy.mpl.ticker import LongitudeFormatter ,LatitudeFormatter
import cartopy.crs as ccrs
import cartopy.feature as cfeature
p=r'202289.nc'
data=xr.open_dataset(p).sel(time=slice("2022","2022"))
# 地面风场属性字段是u/v10 高空气压风场改成u/v
u=data.u10
v=data.v10

w=np.sqrt(u*u+v*v)
lon=data.longitude.data
lat=data.latitude.data
def make_map(ax, title,box,xstep,ystep):
    # set_extent  set crs
    ax.set_extent(box, crs=ccrs.PlateCarree())
    ax.coastlines(scale)  # set coastline resolution
    # set coordinate axis
    ax.set_xticks(np.arange(box[0], box[1], 2),crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(box[2], box[3], 2),crs=ccrs.PlateCarree())
    ax.xaxis.set_major_formatter(cticker.LongitudeFormatter())
    #经度0不加标识
    ax.yaxis.set_major_formatter(cticker.LatitudeFormatter())
    ax.set_title(title, fontsize=15, loc='center') 
    return ax
fig=plt.figure(figsize=(25,35))
x,y=np.meshgrid(lon,lat)
box1 = [105, 125, 30, 45]  
scale = '50m'            
xstep, ystep = 10, 10
cmap=plt.get_cmap('Spectral_r')#'RdYlBu_r' 配色网上有很多参考 可以搜cmap配色
# 可以改成任意需要的月份
titl=['Aug','Sept']

# =============================================================================
# print(x[::5,::5].shape)
# print(u.data)
# print(v.data)
# print(w.data[0].shape)
# =============================================================================
#几个月份循环内的数字就改成几
for i in range(2):
    print(f"i{i}")
    proj=ccrs.PlateCarree(central_longitude=115)
    ax=fig.add_subplot(4,2,i+1,projection=proj)
    make_map(ax,str(titl[i]),box1,45,45)
#这个代码是读数据绘制箭头的 具体粗细啥的可以找文档调整属性
    cb=ax.quiver(x[::5,::5],y[::5,::5],u.data[i,:,:][::5,::5],v.data[i,:,:][::5,::5],pivot='mid',\
    width=0.0018,scale=150,transform=ccrs.PlateCarree(),color='k',angles='xy',zorder=1)
    # cp=ax.contour(lon,lat,w.data[i],zorder=0,transform=ccrs.PlateCarree(),cmap=cmap,levels=np.arange(0,31,1),extend='both')
    cp=ax.contourf(lon,lat,w.data[i],zorder=0,transform=ccrs.PlateCarree(),cmap=cmap,levels=np.arange(0,31,1),extend='both')
    cbar = fig.colorbar(cp,pad=0.01,label='The wind speed (m/s)')
plt.show()

！特别需要注意的是，我用ArcGIS读了nc数据之后，发现2022年6-9月，存储的格式不一样了，所以可能只能一个月一个月下载数据来绘制（我菜我不会写）

nc数据下载网址：Copernicus Climate Data Store |

我这链接是高空气压的，可以按需自己搜其他数据集来改。

！！有报错就是你catopy库没装好，建议anaconda官方提供的方式下载。