python气象数据可视化学习记录1——基于ERA5数据画风场和海平面气压填色叠加图

  • 1. 写在前面
  • 2. 图片效果
  • 3. 逐步代码解析
    • 3.1导入库
    • 3.2 读取NC格式数据
    • 3.3 对数据进行加工1
    • 3.4 建立画布和子图，选择投影方式
    • 3.5 （最重要的部分）自定义画图函数
      • 3.5.1 为什么要定义画图函数？
      • 3.5.2 生成二维的风场和海平面气压数据
      • 3.5.3 利用cartopy绘制地图底图
      • 3.5.4 画填色图和风场图
    • 3.6 调用函数，画图
    • 3.7 添加colorbar
    • 3.8 添加图中的text
    • 3.9 图片输出
  • 4. 代码完整版

1. 写在前面

大学时为了下载某些资源申请CSDN账号，11年过去了，一篇博文也没有，惭愧。受小王同学的启发，应该把平时自己的一些东西和学习经验记录下来，也是一种很好的经验。所以有了这第一篇文章。

学习python想了好几年了，但因为自己受NCL影响太多，一直走不了这个坑。近些年来，python在气象和大气环境领域有了很好的发展，自己也有了危机感，也决定有一些改变。记录博文也是为了督促自己，当然如果有人能从中得到收获，那也是更开心的。

因为自己还是新手，很多东西是看书+自己琢磨+官网+各类公众号和气象家园的综合，如果有任何错误和版权问题，欢迎大家指出！

2. 图片效果

最近自己科研工作上有一个需求，要利用ERA5 0.25°的小时资料画风场和海平面气压叠加的图，包含4个子图，分别代表了不同的时间，效果如下：

  • 填色代表海平面气压hPa，箭头为风场

3. 逐步代码解析

3.1导入库

使用到的库主要有

  1. matplotlib用来做图，包含了画图的pyplot,ticker等函数
  2. numpy 简单计算
  3. cartopy 绘地图
  4. netCDF4 读取NC格式数据
import matplotlib.pyplot as plt###引入库包
import numpy as np
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
import netCDF4 as nc
import matplotlib.pyplot as plt

此外，为了控制图内所有字体大小、轴线粗细、字体等信息一致，使图片更加美化，统一进行如下设置：

mpl.rcParams["font.family"] = 'Arial'  #默认字体类型
mpl.rcParams["mathtext.fontset"] = 'cm' #数学文字字体
mpl.rcParams["font.size"] = 12   #字体大小
mpl.rcParams["axes.linewidth"] = 1   #轴线边框粗细（默认的太粗了）

最后，一个python大神告诉我的方法，加入一个魔法函数，可以不使用plt.show()就可以显示图片的方法：

%matplotlib inline

到此，引入库包的步骤就完成啦！

3.2 读取NC格式数据

读取文件，这里用了ERA5的pressure文件来选取850hpa上的uv风场，single文件里则代表了地表的参量。

f1=nc.Dataset('E:/winter-photo-data/sum/data/ERA-data/202010-ERA-single.nc')
f2=nc.Dataset('E:/winter-photo-data/sum/data/ERA-data/202010-ERA-pressure.nc')

读取变量，需要用到msl（海平面气压），lat，lon，time，u，v等信息。
注意：

  1. 如果不知道变量名，可以print(f1)来看下具体有哪些变量。
  2. ERA5中lat是从50->30这样从北往南排列的。
  3. lat,lon 都是一维函数，msl是（时间，lat,lon）三维函数，u850和v850都是（时间，层数，lat，lon）的四位函数，绘图时只需要二维结果，所以后续要对数据进行挑选和加工。
msl=f1.variables["msl"][:]/100   #pa->hpa
lat=f1.variables['latitude'][:]
lon=f1.variables['longitude'][:]
time=f1.variables['time'][:]
u850=f2.variables["u"][:]
v850=f2.variables["v"][:]

3.3 对数据进行加工1

首先，要建立二维的lon2d,和lat2d，后面画填色图的时候经纬度信息要是二维的。

lon2d, lat2d = np.meshgrid(lon, lat)

其次，用布尔索引的方法挑选需要的时间段。

time_all=[(time>=1058760) & (time<=1059096)]   #13-27 Oct
time_20=[(time>=1058928) & (time<=1058951)]   # 20 Oct
time_25=[(time>=1059048) & (time<=1059071)]   # 25 Oct
time_clean=[(time>=1058760) & (time<1058928)]

这样，数据就挑选完成了，接下来将用lon2d,lat2d,msl,u850,v850，还有time等逻辑值来画图啦。

3.4 建立画布和子图，选择投影方式

使用的是最简单的lat,lon的投影，图片大小是10*8，有两行两列4个子图，这样就建立了画布，得到了各个子图的轴ax。

proj = ccrs.PlateCarree(central_longitude=130)  #中国为左
fig = plt.figure(figsize=(10,8),dpi=550)  # 创建画布
ax = fig.subplots(2, 2, subplot_kw={'projection': proj})  # 创建子图

3.5 （最重要的部分）自定义画图函数

3.5.1 为什么要定义画图函数？

要画多张图，可统一定义函数，后面直接传递参数，调用画图即可。不然要重复多次步骤。

3.5.2 生成二维的风场和海平面气压数据

定义的函数名为plot4，传递的参量只有time_range（对应time_all, time_20等布尔逻辑值，用来挑选所需要的时次），ax轴，labels为图左上角的标签。

msl_all等为时间平均后的二维数组，u_all等为时间平均后包含高度层信息的三维数组，后面再直接选择特定的层就行，这里没有选择。

def plot4(time_range,ax,labels):msl_all=np.mean(msl[time_range],axis=0)u_all=np.mean(u850[time_range],axis=0)v_all=np.mean(v850[time_range],axis=0)

3.5.3 利用cartopy绘制地图底图

这里参考了小白学习cartopy画地图的第一天（中国行政区域图，含南海）的博文，添加了经纬度信息，并对中国地区的边界进行了重新的设定。具体不再说明，上面的文章里写的很详细。如果用cartopy画地图的话可以学习下。需要注意的是选择extent=[100,125,30,50]经纬度信息时，一定要跟你画图时二维数组的经纬度信息一致，不然出的图是错的。

#-----------绘制地图-------------------------------------------ax.add_feature(cfeature.LAND.with_scale('50m'))####添加陆地######ax.add_feature(cfeature.COASTLINE.with_scale('50m'),lw=0.3)#####添加海岸线#########ax.add_feature(cfeature.OCEAN.with_scale('50m'))######添加海洋#########-----------添加经纬度---------------------------------------extent=[100,125,30,50]##经纬度范围gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0., color='k', alpha=0.5, linestyle='--')gl.top_labels = False ##关闭上侧坐标显示gl.right_labels = False ##关闭右侧坐标显示gl.xformatter = LONGITUDE_FORMATTER ##坐标刻度转换为经纬度样式gl.yformatter = LATITUDE_FORMATTER gl.xlocator = mticker.FixedLocator(np.arange(extent[0], extent[1]+1, 5))gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3]+1, 5))gl.xlabel_style={'size':10}gl.ylabel_style={'size':10}ax.set_extent(extent)     #显示所选择的区域#----------修改国界，并添加省界-----------------------------with open('C:/Users/hj/.local/share/cartopy/shapefiles/natural_earth/physical/CN-border-La.dat') 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]for line in borders:ax.plot(line[0::2], line[1::2], '-', color='k',lw=0.3, transform=ccrs.Geodetic())

3.5.4 画填色图和风场图

  1. 首先levels定义了填色图色标的范围，1004-1032之间，间隔1个给一个颜色。
  2. ax.contourf绘制了海平面气压的填色图，参数lon2d, lat2d ,msl_all都是二维的数组，colormap选择了’Spectral_r’，这个可以看个人的需求，网上有很多资源。
  3. ax.quiver是绘制风场的，也需要输入lon2d, lat2d，u,v的信息（此时uv为三维数组，需要挑选自己所需要的层数，变成二维的），scale用来调整画出来的风场的大小，值越大，风的线条越小。width是风线条的粗细，headwidth和headlength是箭头粗细和长度。uv和lat,lon都是每隔5个显示，是为了使风场图不那么密集，这个可以根据自己的需求进行调节。
    最重要的是，一定要加上transform=ccrs.PlateCarree()，图片中地图需要转换下，使地图和数据直接做出来的二维图对应上。
  4. 绘制每个子图的标题，pad代表标题离轴的远近。
    返回了填色图和风场图，这样函数就完成了，后面只需要传递参数调用就行。
#-------------------plot---------------------------levels = np.arange(1004, 1032 + 1, 1)cb = ax.contourf(lon2d,lat2d,msl_all, levels=levels, cmap='Spectral_r',transform=ccrs.PlateCarree())cq = ax.quiver(lon2d[::5,::5],lat2d[::5,::5],u_all[16,::5,::5],v_all[16,::5,::5],color='k',scale=75,zorder=10,width=0.005,headwidth=3,headlength=4.5,transform=ccrs.PlateCarree())   #925hpa->16ax.set_title(labels,loc="left",fontsize=12,pad=1)return cb, cq

3.6 调用函数，画图

subplots_adjust调整了图片的大小和子图之间的距离等，使图片更加美观。
后面直接进行4次函数调用，就可以啦。

#------------------------plot-----------------------
plt.subplots_adjust(left=0.15,right=0.85,top=0.8,bottom=0.2,wspace=0.15,hspace=0.2)cb1,cq1=plot4(time_all,ax[0][0],'(a) Mean')
cb2,cq2=plot4(time_clean,ax[0][1],'(b) Clean days')
cb3,cq3=plot4(time_20,ax[1][0],'(c) 10/20')
cb4,cq4=plot4(time_25,ax[1][1],'(d) 10/25')

3.7 添加colorbar

因为是填色图，所以一定要条件colorbar. 因为是4幅图一个色标，所以ax=(ax[0][0],ax[0][1],ax[1][0],ax[1][1])，方向为垂直放置在cb3旁边，tick显示的范围为1004-1032，每隔4个显示。aspect是colorbar的长宽高比例，shrink是代表要显示多长，是百分比。pad是离图片的距离。

cbar=fig.colorbar(cb3,ax=(ax[0][0],ax[0][1],ax[1][0],ax[1][1]),orientation='vertical',ticks=np.arange(1004, 1032 + 4,4),aspect=30,shrink=0.9,pad=0.02)
cbar.ax.tick_params(pad=0.01,length=2,width=0.15)

ax.tick_params里的pad的意思是色标上的数字离色标的距离，length和width是色标tick的长断粗细。

3.8 添加图中的text

为了在图中标出高压和低压中心，需要用ax.text的进行添加。ax为需要添加到哪个图中，0.85和0.87则用百分比的方式决定了点的位置，其他不再赘述，比较简单。一定要添加transform=ax[1][0].transAxes，进行坐标转换。

#--------------------add text----------------------------
ax[1][0].text(0.85,0.87,'L',color='red',fontsize=15,weight='bold',va='bottom',ha='center',transform=ax[1][0].transAxes)
ax[1][1].text(0.92,0.12,'H',color='red',fontsize=15,weight='bold',va='bottom',ha='center',transform=ax[1][1].transAxes)

3.9 图片输出

可以输出eps,pdf等矢量图，也可以有Png这样的栅格图。OK，大功告成！

plt.savefig("Figure2.png")

4. 代码完整版

import matplotlib.pyplot as plt###引入库包
import numpy as np
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
import netCDF4 as nc
import matplotlib.pyplot as pltmpl.rcParams["font.family"] = 'Arial'  #默认字体类型
mpl.rcParams["mathtext.fontset"] = 'cm' #数学文字字体
mpl.rcParams["font.size"] = 12
mpl.rcParams["axes.linewidth"] = 1
%matplotlib inlinef1=nc.Dataset('E:/winter-photo-data/sum/data/ERA-data/202010-ERA-single.nc')
f2=nc.Dataset('E:/winter-photo-data/sum/data/ERA-data/202010-ERA-pressure.nc')
msl=f1.variables["msl"][:]/100   #pa->hpa
lat=f1.variables['latitude'][:]
lon=f1.variables['longitude'][:]
time=f1.variables['time'][:]
u850=f2.variables["u"][:]
v850=f2.variables["v"][:]
lon2d, lat2d = np.meshgrid(lon, lat)#-----------------mean-----------------------------
time_all=[(time>=1058760) & (time<=1059096)]   #13-27 Oct
time_20=[(time>=1058928) & (time<=1058951)]   # 20 Oct
time_25=[(time>=1059048) & (time<=1059071)]   # 25 Oct
time_clean=[(time>=1058760) & (time<1058928)]
proj = ccrs.PlateCarree(central_longitude=130)  #中国为左
fig = plt.figure(figsize=(10,8),dpi=550)  # 创建画布
ax = fig.subplots(2, 2, subplot_kw={'projection': proj})  # 创建子图
def plot4(time_range,ax,labels):msl_all=np.mean(msl[time_range],axis=0)u_all=np.mean(u850[time_range],axis=0)v_all=np.mean(v850[time_range],axis=0)#-----------绘制地图-------------------------------------------ax.add_feature(cfeature.LAND.with_scale('50m'))####添加陆地######ax.add_feature(cfeature.COASTLINE.with_scale('50m'),lw=0.3)#####添加海岸线#########ax.add_feature(cfeature.OCEAN.with_scale('50m'))######添加海洋#########-----------添加经纬度---------------------------------------extent=[100,125,30,50]##经纬度范围gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0., color='k', alpha=0.5, linestyle='--')gl.top_labels = False ##关闭上侧坐标显示gl.right_labels = False ##关闭右侧坐标显示gl.xformatter = LONGITUDE_FORMATTER ##坐标刻度转换为经纬度样式gl.yformatter = LATITUDE_FORMATTER gl.xlocator = mticker.FixedLocator(np.arange(extent[0], extent[1]+1, 5))gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3]+1, 5))gl.xlabel_style={'size':10}gl.ylabel_style={'size':10}ax.set_extent(extent)     #显示所选择的区域#----------修改国界，并添加省界-----------------------------with open('C:/Users/hj/.local/share/cartopy/shapefiles/natural_earth/physical/CN-border-La.dat') 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]for line in borders:ax.plot(line[0::2], line[1::2], '-', color='k',lw=0.3, transform=ccrs.Geodetic())#-------------------plot---------------------------levels = np.arange(1004, 1032 + 1, 1)cb = ax.contourf(lon2d,lat2d,msl_all, levels=levels, cmap='Spectral_r',transform=ccrs.PlateCarree())cq = ax.quiver(lon2d[::5,::5],lat2d[::5,::5],u_all[16,::5,::5],v_all[16,::5,::5],color='k',scale=75,zorder=10,width=0.005,headwidth=3,headlength=4.5,transform=ccrs.PlateCarree())   #925hpa->16ax.set_title(labels,loc="left",fontsize=12,pad=1)return cb, cq
#------------------------plot-----------------------
plt.subplots_adjust(left=0.15,right=0.85,top=0.8,bottom=0.2,wspace=0.15,hspace=0.2)cb1,cq1=plot4(time_all,ax[0][0],'(a) Mean')
cb2,cq2=plot4(time_clean,ax[0][1],'(b) Clean days')
cb3,cq3=plot4(time_20,ax[1][0],'(c) 10/20')
cb4,cq4=plot4(time_25,ax[1][1],'(d) 10/25')cbar=fig.colorbar(cb3,ax=(ax[0][0],ax[0][1],ax[1][0],ax[1][1]),orientation='vertical',ticks=np.arange(1004, 1032 + 4,4),aspect=30,shrink=0.9,pad=0.02)
cbar.ax.tick_params(pad=0.01,length=2,width=0.15)#--------------------add text----------------------------ax[1][0].text(0.85,0.87,'L',color='red',fontsize=15,weight='bold',va='bottom',ha='center',transform=ax[1][0].transAxes)
ax[1][1].text(0.92,0.12,'H',color='red',fontsize=15,weight='bold',va='bottom',ha='center',transform=ax[1][1].transAxes)plt.savefig("Figure2.png")

python气象数据可视化学习记录1——基于ERA5数据画风场和海平面气压填色叠加图相关推荐

  1. Python气象数据可视化学习笔记5——基于cartopy绘制contour并对中国地区进行白化（包含南海）

    基于cartopy绘制contour并对中国地区进行白化(包含南海) 1. 写在前面 利用cartopy画填色图已经掌握,这一篇主要记录了在填色的基础上叠加白化.主要参考了气象家园的两篇帖子,并进行了 ...

  2. python气象数据可视化学习笔记6——利用python地图库cnmaps绘制地图填色图并白化

    文章目录 1. 效果图 2. cnmaps简介及安装 2.1 写在前面 2.2 cnmaps简介和安装 3. 导入库 4. 定义绘图函数 4.1 使用get_adm_maps返回地图边界 4.2 ax ...

  3. Python数据可视化学习（初学中...）

    Python数据可视化学习(初学中...) 1.使用Matplotlib生成数据图 1.1.安装Matplotlib包 1.2.Matplotlib数据图入门 1.2.1.折线图举例 1.2.2.图表 ...

  4. Python+vtk 实现激光点云数据可视化学习（2021.7.12）

    Python+VTK实现激光点云数据可视化学习 2021.7.12 1.激光点云与VTK简介 2.配置Python环境(Conda+PyCharm+Python3.6+VTK) 3.点云数据(大约60 ...

  5. 27【源码】数据可视化大屏：基于 Echarts + Python Flask 实现的32-9超宽大屏范例 - 监控指挥中心

    目录 效果展示 1. 效果动图 2. 多种主题效果 一. 确定需求方案 1. 屏幕分辨率 2. 部署方式 二. 整体架构设计 三. 编码实现 (基于篇幅及可读性考虑,此处展示部分关键代码) 1. 前端 ...

  6. 29【源码】数据可视化大屏：基于 Echarts + Python Flask 实现的32-9超宽大屏 - 企业综合信息

    我是 YYDataV数据可视化  专注于 数据可视化大屏,工厂扫码装箱系统 等 我的微信 6550523,多多交流 ~ 本案例为32:9超宽分辨率的大屏. 效果展示 1.动态实时更新数据效果图 2.鼠 ...

  7. python数据可视化第三方库有哪些_数据可视化！看看程序员大佬都推荐的几大Python库...

    数据可视化是数据分析中极为重要的部分,而数据可视化图表(如条形图,散点图,折线图,地理图等)也是非常关键的一环.Python作为数据分析中最流行的编程语言之一,有几个库可以创建精美而复杂的数据可视化, ...

  8. 数据可视化（三）基于 Graphviz 实现程序化绘图

    2019独角兽企业重金招聘Python工程师标准>>> 前言 我之前在几篇文章新一代Ntopng网络流量监控-可视化和架构分析. 数据可视化(一)思维利器 OmniGraffle 绘 ...

  9. 数据可视化学习之大屏学习

    一 前言 什么是数据可视化大屏?数据可视化大屏是以大屏为主要展示载体的数据可视化设计.可视化大屏就是一种非常有效的数据可视化工具,它可以将业务的关键指标以可视化的方式展示到一个或多个LED屏幕上,不仅 ...

最新文章

  1. c语言获取指针分配的字节数,c语言指针知识点总结(共6篇).docx
  2. Oracle 12c DG备库Alert报错ORA-01110
  3. Linux修改主机名的方法
  4. 【思维题 状压dp】APC001F - XOR Tree
  5. 2021略阳天津高级中学高考成绩查询,2021年天津高考成绩查询网站查分网址：http://www.zhaokao.net/...
  6. 拥抱时序数据库，构筑IoT时代下智慧康养数据存储底座
  7. 华为MatePad Pro 5G平板正式发布：售价5299元起！
  8. 【github】git 使用命令大全
  9. 视频怎么插入慢动作？
  10. python创建系列_一起学python系列之类（创建和使用类）
  11. 判断对象是否超出屏幕
  12. 深入C++“准”标准库，Boost你的力量
  13. handler更新ui线程的基本用法
  14. HR-Former | 随迟但到，HRNet+Transformer轻装归来（非常值得学习！！！）
  15. 利用sklearn实现adaboost，以单一分类树为例
  16. JDK 和 JRE 有什么区别？
  17. typecho插件仓库集合版，非常方便的使用插件
  18. Android 获取手机号码
  19. 耶斯莫拉~一起来学集合啊~
  20. AM5728+QT的图像采集与处理应用, 中文字库显示

热门文章

  1. 软件测试面试如何写出HR青睐的简历？
  2. 拓扑排序Kahn算法
  3. 15位与18位身份证号码有什么区别和联系
  4. EyouCMS精美简洁作文范文网站模板/易优CMS资讯类企业网站模板
  5. 如何阅读综述文献，以及摘要、结论和引言部分
  6. ThinkPHP5开发(一)实现登录功能
  7. 150行代码写个低配版WPS？：手把手教你实现+附完整源码
  8. cybersource支付对接
  9. JAVA异常处理（三种异常处理机制）
  10. 计算机32位怎么换成64位,怎么把32位CAD换成64位系统