python--循环绘制ERA5风场的空间分布图
使用python封装绘图函数循环绘制ERA5风场资料的空间分布图
通常,在处理气象海洋资料时,经常会绘制风场的空间分布图进行简单分析,而常常需要连续绘制多天,并将多张子图绘制到同一个图片中,因此这就需要用到循环绘图。
- 同时考虑到下载的ERA5风场资料的经度排列顺序是-180~180°,这里也简单进行了经度转换,将其转换为0 ~ 360的排列顺序。
- 根据每个子图的数据,将选取的时间也在循环中加上
- 考虑到绘制全球的处理时间较久,这里自由选取任意经纬度进行绘制
- 由于每张子图的填色范围是固定的,所以统一绘制colorbar,至于图片最下端
- 箭头为风场风量,填色为风速大小
下面附上代码:
# -*- coding: utf-8 -*-
"""
Created on %(date)s
@author: %(jixianpu)s
Email : 211311040008@hhu.edu.cn
introduction : keep learning althongh walk slowly
"""
from matplotlib.ticker import AutoMinorLocator, MultipleLocator, FuncFormatter
import numpy as np
import pandas as pd
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
import matplotlib.pyplot as plt
import xarray as xr
import glob
from scipy import signal
path=r'G:/ERA5_200406/'
str1='uwnd.200406.nc'
str2='vwnd.200406.nc'
# b, a = signal.butter(3, [(2/10)/(1/0.25),(2/3)/(1/0.25)], 'bandpass')
######===================
def make_map(ax, title):
# set_extent set crs
ax.set_extent(box, crs=ccrs.PlateCarree())
land = cfeature.NaturalEarthFeature('physical',
'land',
scale,
edgecolor='grey',
facecolor='grey'
,zorder=2
)
ax.add_feature(land) # set land color
ax.coastlines(scale) # set coastline resolution
# set coordinate axis
ax.set_xticks(np.arange(box[0],box[1] xstep, xstep),crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(box[2], box[3] ystep, ystep),crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label =False))#经度0不加标识
ax.yaxis.set_major_formatter(LatitudeFormatter())
# plt.tick_params(labelsize=25)
ax.set_title(title, fontsize=20, loc='left',pad=12)
ax.yaxis.set_minor_locator(AutoMinorLocator(5))
ax.xaxis.set_minor_locator(AutoMinorLocator(10))
ax.tick_params(which='minor',
direction='out', length=4,width=0.59,
right=True, top=True)
ax.tick_params(which='major',
direction='out',
length=8,
width=0.99,
pad=3,
labelsize=10,
bottom=True, left=True, right=True, top=True)
return ax
# # prepare
box = [100, 180, -30, 0]
scale = '50m'
xstep, ystep = 20, 10
# cmap=plt.get_cmap('RdYlBu_r')#'RdYlBu_r'
# =======================data============================================
day=(np.arange(19,31))
da = xr.open_dataset(path str1).sortby("latitude", ascending= True)
da2=xr.open_dataset(path str2).sortby("latitude", ascending= True)
#####################reverse################################################
lon_name = 'longitude' # whatever name is in the data
da['longitude_adjusted'] = xr.where(da[lon_name] < 0, da[lon_name]60,\
da[lon_name])
da = (
da
.swap_dims({lon_name: 'longitude_adjusted'})
.sel(**{'longitude_adjusted': sorted(da.longitude_adjusted)})
.drop(lon_name))
da = da.rename({'longitude_adjusted': lon_name})
#####################reverse################################################
lon_name = 'longitude' # whatever name is in the data
da2['longitude_adjusted'] = xr.where(da2[lon_name] < 0, da2[lon_name]60,\
da2[lon_name])
da2 = (
da2
.swap_dims({lon_name: 'longitude_adjusted'})
.sel(**{'longitude_adjusted': sorted(da2.longitude_adjusted)})
.drop(lon_name))
da2 = da2.rename({'longitude_adjusted': lon_name})
#####################reverse################################################
fig=plt.figure(figsize=(14,14))
levels=np.arange(0,21,1)
# wspace 调整水平的
plt.subplots_adjust(hspace=0.1)
plt.tight_layout()
count=0
for i in day:
count =1
print(count)
tim_s='2004-06-' str(i) '-00'
# longitude=slice(0,181)
u=da.u.sel(level=slice(850,850),latitude=slice(-30,0),longitude=slice(100,180),time=slice(tim_s,tim_s))[0][0]
v=da2.v.sel(level=slice(850,850),latitude=slice(-30,0),longitude=slice(100,180),time=slice(tim_s,tim_s))[0][0]
w=np.sqrt(u*u v*v)
lon=u.longitude.data
lat=u.latitude.data
x,y=np.meshgrid(lon,lat)
t=u.time
step=10
# =======================plot============================================
ax=fig.add_subplot(4,3,count,projection=ccrs.PlateCarree(central_longitude=180))
make_map(ax,pd.to_datetime(t.data).strftime('%Y_%m_%d_%H:00'))
ax.quiver(x[::step,::step],y[::step,::step],u.data[::step,::step],v.data[::step,::step],pivot='mid',\
width=0.003,scale=200,headlength=4,headwidth=4,
transform=ccrs.PlateCarree(),color='k',angles='xy'
,zorder=2)
# ax.set_aspect('auto')
plot=ax.contourf(x,y,w,transform=ccrs.PlateCarree(),extend='both',levels=levels
,zorder=1)
# cb=fig.colorbar(plot,ax=ax,shrink=0.8,pad=0.05,aspect=15)
# cb.ax.tick_params(labelsize=20)
# cb.ax.set_title('$°C$',fontsize=20,loc='right')
# if count==6:
# break
ax3=fig.add_axes([0.25,0.1,0.5,0.015]) # 0.25控制距离左边的距离,0.01控制距离下面的距离,最后两位控制color的长度和厚度
cb=fig.colorbar(plot,cax=ax3,shrink=0.9,pad=0.04,ticks=[0,5,10,15,20],aspect=15,orientation='horizontal')
cb.ax.tick_params(labelsize=10)
cb.ax.set_title('$m/s$',fontsize=15)
plt.show()
# fig.savefig('wind_850hpa-2004-06-19_30.png',format='png',dpi=100)
绘制结果如下图所示:
这篇好文章是转载于:学新通技术网
- 版权申明: 本站部分内容来自互联网,仅供学习及演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系,请提供相关证据及您的身份证明,我们将在收到邮件后48小时内删除。
- 本站站名: 学新通技术网
- 本文地址: /boutique/detail/tanhgagbaa
系列文章
更多
同类精品
更多
-
photoshop保存的图片太大微信发不了怎么办
PHP中文网 06-15 -
《学习通》视频自动暂停处理方法
HelloWorld317 07-05 -
word里面弄一个表格后上面的标题会跑到下面怎么办
PHP中文网 06-20 -
Android 11 保存文件到外部存储,并分享文件
Luke 10-12 -
photoshop扩展功能面板显示灰色怎么办
PHP中文网 06-14 -
微信公众号没有声音提示怎么办
PHP中文网 03-31 -
excel下划线不显示怎么办
PHP中文网 06-23 -
excel打印预览压线压字怎么办
PHP中文网 06-22 -
TikTok加速器哪个好免费的TK加速器推荐
TK小达人 10-01 -
怎样阻止微信小程序自动打开
PHP中文网 06-13