Python读取NC格式数据绘制水汽通量等值线和和流场

这篇具有很好参考价值的文章主要介绍了Python读取NC格式数据绘制水汽通量等值线和和流场。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。


Python读取NC格式数据绘制水汽通量等值线和和流场

一、知识点

计算水汽通量,用到了metpy包,是一个地球科学计算包,内置了很多气象用到的计算函数

小知识点:
1.用湿度计算比湿
2.单位的使用
3.常量的使用,这里涉及了重力加速度g

二、代码拆分

导入包

#处理数据的包
import xarray as xr
import numpy as np
#画图的包
import matplotlib.pyplot as plt
#地图的包
import cartopy.crs as ccrs
import my_class                   #一个自己写的画地图的py文件
#科学计算的包
from metpy.units import units      #里面是单位
import metpy.constants          #里面是常数
import metpy.calc                #里面有各种计算函数

读取数据

ds = xr.open_dataset('download_201306.nc')#用xarray读取NC数据
lat = ds.latitude#读取维度
lon = ds.longitude#读取经度
time = ds.time#读取时间
u = ds['u']#风场U分量
v = ds['v']#风场V分量
tem = ds['t']- 273.15#读取温度,转化为摄氏度
rh = ds['r']#读取湿度

注意:这里读取的数据是全部的格点数据,但是我们画图用不了这么多,所以对数据做分割。只需要一部分

数据筛选分割

#设置经纬度范围
lat_range = lat[(lat>=22) & (lat<=33)]
lon_range = lon[(lon>=106) & (lon<=124)]
#筛选数据(经纬度,时间)
u_region = u.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-28T00:00:00.000000000')
v_region = v.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-28T00:00:00.000000000')
tem_region = tem.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-28T00:00:00.000000000')
rh_region = rh.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-28T00:00:00.000000000')

上面获得了世界时2013年6月28日00时的UV风场、温度、湿度
这里解释一下,作者下载的在分析资料里没有分层的数据,只下载了850hPa一层数据,所以不需要设置高度。

数据计算

#给温度赋予单位摄氏度,正真意义上的摄氏度,metpy有单位体系
tem_region = tem_region * units.degC
#气压赋予hPa
pressure = 850 * units.hPa
#用相对湿度和温度计算露点,函数名字就是字面意思
dewpoint = metpy.calc.dewpoint_from_relative_humidity(tem_region, rh_region)
#用露点和气压算比湿,函数名字就是字面意思,因为计算结果是kg数值,乘以1000变成克
specific_humidity = metpy.calc.specific_humidity_from_dewpoint(pressure, dewpoint)*1000
#计算两个方向上的水汽通量,公式就是  速度*比湿/重力加速度
qu = u_region * specific_humidity / metpy.constants.g
qv = v_region * specific_humidity / metpy.constants.g
#和风速一样,勾股定理就得到了水汽通量的分布数值,为啥这么麻烦呢,因为UV分量一会儿还要画箭头
a = np.sqrt(qu * qu + qv*qv)

这里的知识点:
1.温度和气压赋予单位,因为metpy的计算是有单位体系的
2.利用metpy.calc的内置函数计算物理量
dewpoint_from_relative_humidity:露点_from_相对湿度,参数是温度和湿度
specific_humidity_from_dewpoint:比湿_from_露点,参数是气压和露点。
3.利用水汽通量公式:速度*比湿/重力加速度计算水汽通量的UV分量。

画图设置

#画布
fig = plt.figure(figsize=(9,6))
#子图
ax = fig.add_subplot(111,projection = ccrs.PlateCarree())
ax.set_extent([106, 124, 22, 33])
ax.set_xticks(np.arange(106,125,2),crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(22,34,1),crs=ccrs.PlateCarree())
xticks_str = ['106', '108', '110', '112', '114', '116', '118','120', '122', '124°E']
ax.set_xticklabels(xticks_str,fontsize = 11)
yticks_str = ['22   ','23    ','24    ','25    ','26    ','27    ','28    ','29    ','30    ','31    ','32    ','33°N']
ax.set_yticklabels(yticks_str,fontsize = 11)
my_class.readshapefile('bou2_4l.shp',linewidth=1,ax=ax)

画图设置这里就不多讲了,就是弄个画布、弄个子图,设置一下XY轴,加载个地图。

画图

#画水汽通量等值线
ct=ax.contour(lon_range,lat_range,a,8,colors='k',linewidths=1,levels=range(0,28,2))
ax.clabel(ct,inline=True,fontsize=10,fmt='%.0f')
#下面画箭头
#获得XY矩阵格式坐标,这里对纬度上的格点进行了处理,降尺度,从0.25变成了1,其实就是间隔4取值,目的是因为箭头太密影响效果
h1 = ax.quiver(lon_range[::4],lat_range[::4],qu[::4,::4],qv[::4,::4],     #X,Y,U,V 确定位置和对应的风速
                width = 0.002, #箭杆箭身宽度
                scale = 700, # 箭杆长度,参数scale越小箭头越长
                pivot = 'tail'#箭头的其实位置,这里表示从点起,还有点在中心的‘mid’
                )

#画出风场,和箭头箭轴后,得说明 箭轴长度与风速的对应关系
#调用quiver可以生成 参考箭头 + label。
ax.quiverkey(h1,                      #传入quiver句柄
              X=0.8, Y = -0.07,       #确定 label 所在位置,都限制在[0,1]之间
              U = 20,                    #参考箭头长度 表示20。
              angle = 0,            #参考箭头摆放角度。默认为0,即水平摆放
             label='20',        #箭头的补充:label的内容  +
             labelpos='S',          #label在参考箭头的哪个方向; S表示南边
             color = 'k',labelcolor = 'k', #箭头颜色 + label的颜色
             )

这里画了一个等值线,画了一个水汽通量的方向场,其实就是画风场的方法。

出图

plt.show()

Python读取NC格式数据绘制水汽通量等值线和和流场

三、完整代码

这里的完整代码是两张图,包含了一个循环,就是把两个时次的图放在一起。文章来源地址https://www.toymoban.com/news/detail-400431.html

'''
知识点:计算水汽通量,用到了metpy包,是一个地球科学计算包,内置了很多气象用到的计算函数
小知识点:1.用湿度计算比湿
          2.单位的使用
          3.常量的使用,这里涉及了重力加速度g
'''
#处理数据的包
import xarray as xr
import numpy as np
#画图的包
import matplotlib.pyplot as plt
#地图的包
import cartopy.crs as ccrs
import my_class
#科学计算的包
from metpy.units import units      #里面是单位
import metpy.constants          #里面是常数
import metpy.calc                #里面有各种计算函数



#循环的东西这里就不讲解了
ds = xr.open_dataset('download_201306.nc')
lat = ds.latitude
lon = ds.longitude
time = ds.time
u = ds['u']
v = ds['v']
tem = ds['t']- 273.15#温度变成摄氏度
rh = ds['r']#湿度


lat_range = lat[(lat>=22) & (lat<=33)]
lon_range = lon[(lon>=106) & (lon<=124)]
fig = plt.figure(figsize=(18,9))

#设置2个子图,并且放到列表里面,方便循环的时候用
plt.subplots_adjust(left=0.07, right=0.90, top=0.95, bottom=0.05,wspace=0.2,hspace=0.1)
ax_a = fig.add_subplot(1,2,1,projection = ccrs.PlateCarree())
ax_b = fig.add_subplot(1,2,2,projection = ccrs.PlateCarree())


ax_list = [ax_a,ax_b]
ab = ['(a)','(b)']#图右上角的ab
i = 0#循环变量

for time_i in time:
    #风U
    u_region = u.sel(longitude=lon_range, latitude=lat_range,time = time_i)
    #风V
    v_region = v.sel(longitude=lon_range, latitude=lat_range,time = time_i)
    #温度
    tem_region = tem.sel(longitude=lon_range, latitude=lat_range,time = time_i)
    #湿度
    rh_region = rh.sel(longitude=lon_range, latitude=lat_range,time = time_i)
    #给温度赋予单位摄氏度,正真意义上的摄氏度,metpy有单位体系
    tem_region = tem_region * units.degC
    #气压赋予hPa
    pressure = 850 * units.hPa
    #用相对湿度和温度计算露点,函数名字就是字面意思
    dewpoint = metpy.calc.dewpoint_from_relative_humidity(tem_region, rh_region)
    #用露点和气压算比湿,函数名字就是字面意思,因为计算结果是kg数值,乘以1000变成克
    specific_humidity = metpy.calc.specific_humidity_from_dewpoint(pressure, dewpoint)*1000
    #计算两个方向上的水汽通量,公式就是  速度*比湿/重力加速度
    qu = u_region * specific_humidity / metpy.constants.g
    qv = v_region * specific_humidity / metpy.constants.g
    #和风速一样,勾股定理就得到了水汽通量的分布数值,为啥这么麻烦呢,因为UV分量一会儿还要画箭头
    a = np.sqrt(qu * qu + qv*qv)



    ax_list[i].set_extent([106, 124, 22, 33])
    ax_list[i].set_xticks(np.arange(106,125,2),crs=ccrs.PlateCarree())
    ax_list[i].set_yticks(np.arange(22,34,1),crs=ccrs.PlateCarree())
    xticks_str = ['106', '108', '110', '112', '114', '116', '118','120', '122', '124°E']
    ax_list[i].set_xticklabels(xticks_str,fontsize = 11)
    yticks_str = ['22    ','23    ','24    ','25    ','26    ','27    ','28    ','29    ','30    ','31    ','32    ','33°N']
    ax_list[i].set_yticklabels(yticks_str,fontsize = 11)
    my_class.readshapefile('bou2_4l.shp',linewidth=1,ax=ax_list[i])


    #画等值线
    ct=ax_list[i].contour(lon_range,lat_range,a,8,colors='k',linewidths=1,levels=range(0,28,2))
    ax_list[i].clabel(ct,inline=True,fontsize=10,fmt='%.0f')


    #下面画箭头
    #获得XY矩阵格式坐标,这里对纬度上的格点进行了处理,降尺度,从0.25变成了1,其实就是间隔4取值,目的是因为箭头太密影响效果
    h1 = ax_list[i].quiver(lon_range[::4],lat_range[::4],qu[::4,::4],qv[::4,::4],     #X,Y,U,V 确定位置和对应的风速
                    width = 0.002, #箭杆箭身宽度
                    scale = 700, # 箭杆长度,参数scale越小箭头越长
                    pivot = 'tail'#箭头的其实位置,这里表示从点起,还有点在中心的‘mid’
                    )

    #画出风场,和箭头箭轴后,得说明 箭轴长度与风速的对应关系
    #调用quiver可以生成 参考箭头 + label。
    ax_list[i].quiverkey(h1,                      #传入quiver句柄
                  X=0.8, Y = -0.07,       #确定 label 所在位置,都限制在[0,1]之间
                  U = 20,                    #参考箭头长度 表示风速为5m/s。
                  angle = 0,            #参考箭头摆放角度。默认为0,即水平摆放
                 label='20',        #箭头的补充:label的内容  +
                 labelpos='S',          #label在参考箭头的哪个方向; S表示南边
                 color = 'k',labelcolor = 'k', #箭头颜色 + label的颜色
                 )

    #写abcd
    ax_list[i].text(0.05, 0.95, ab[i],transform=ax_list[i].transAxes, fontsize=11)

    #一次循环结束 i+1
    i += 1


plt.show()

到了这里,关于Python读取NC格式数据绘制水汽通量等值线和和流场的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • 水汽稳定度修正函数\Psi_q对潜热通量影响--模式验证工作

    我之前提出了一个水汽通量廓线关系,这项工作偏理论,如果对下面说的背景不了解的话可以看下 https://agupubs.onlinelibrary.wiley.com/share/YNSG74MV8B8BAAUMCHN3?target=10.1029/2022JD036708 那会没把提出的水汽稳定度修正函数加到CAS-ESM,当时对CAS-ESM模式还没这么熟悉,也想着师兄能帮我,但是师

    2024年02月02日
    浏览(36)
  • 利用MATLAB读取.nc文件单像元数值并转为Excel格式(以中国日降雨量月均数据为例)

     以中国日降雨量月均数据(nc文件包含12月)为例,提取某经纬度下的多月份像元值。 (【数据分享】1960-2020年中国1公里分辨率月降水数据集) 一、确定经纬度所在行列 号 以 92.18E,30.475N 为例,首先在Matlab中输入以下代码: 工作区获取到lat、lon和pre的信息,打开lat和lon文件

    2024年02月07日
    浏览(46)
  • Python读取.nc数据并提取指定时间、经纬度维度对应的变量数值

      本文介绍基于 Python 语言的 netCDF4 库,读取 .nc 格式的数据文件,并提取指定维(时间、经度与纬度)下的变量数据的方法。   我们之前介绍过 .nc 格式的数据,其是 NetCDF (Network Common Data Form)文件的扩展名,是一种常用的科学数据存储格式,多用于存储科学和工程领

    2024年03月08日
    浏览(45)
  • NCDC气象数据的提取与处理(四):python批量读取、写入nc数据经纬度格点数值

    1.问题描述: 2.思路: 3.实现过程: 3.1格点位置匹配 3.2写入表格 4.运行效果 4.1打包站点信息 4.2读取nc文件列表 4.3提取对应格点的nc数据 4.4数据写入 NCDC的站点数据处理在之前三节里已经介绍过了,但是NCDC的就那么几种数据可能不能满足日常使用,比如说辐射数据他就没有。

    2024年02月05日
    浏览(78)
  • Python读取.nc文件的方法与技术详解

    目录 一、引言 二、使用netCDF4库读取.nc文件 安装netCDF4库 导入netCDF4库 打开.nc文件 获取变量 读取变量数据 案例与代码 三、使用xarray库读取.nc文件 安装xarray库 导入xarray库 打开.nc文件 访问变量数据 案例与代码 四、性能与优化 分块读取 使用Dask进行并行计算 减少不必要的变量

    2024年04月23日
    浏览(29)
  • 【GIS】Python多线程转换NC格式文件为TIFF

    【GIS】使用cdsapi下载ERA5和ERA5_land逐小时数据 NC文件读取使用netCDF4,NC文件转换为TIF使用rasterio或者GDAL。 格点数据转换为TIFF文件时候,计算六参数时候,应该要考虑,格点数据存储的坐标属于栅格中心点的位置,转换为TIFF时候,可能需要考虑栅格大小。

    2024年02月12日
    浏览(38)
  • python读取excel数据并用双y轴绘制柱状图和折线图,柱子用渐变颜色填充

    往期python绘图合集: python绘制简单的折线图 python读取excel中数据并绘制多子图多组图在一张画布上 python绘制带误差棒的柱状图 python绘制多子图并单独显示 python读取excel数据并绘制多y轴图像 python绘制柱状图并美化|不同颜色填充柱子 python随机生成数据并用双y轴绘制两条带误差

    2024年02月10日
    浏览(47)
  • 第二章 python-pcl、open3d读取、显示pcd、bin等格式点云数据

    点云数据实际上就是许多组点的集合,每个点由{x,y,z}组成。当然理论上的只包含有3D坐标。 实际激光雷达获取的点云数据还会包含强度、反射率等等。但我们一般只用提取{x,y,z}来处理即可。 点云数据相比于其他传感器数据的核心优势就是在于 精准的深度信息。可惜获取具体

    2024年01月16日
    浏览(62)
  • python读取Excel绘制饼图的两种方式

    matplotlib 简单方便,适合数据作图或科学作图(论文发表) pyecharts 流程略复杂,但功能强大,图形具有交互式,适合项目开发或商业分析报告,但是 它是一个非常新的库,开发不稳定 本文介绍用 pandas库 读取Excel (csv)数据,分别用 matplotlib库 和 pyecharts库 绘制饼图。 注: 实

    2024年02月11日
    浏览(43)
  • Python读取文件相对路径理解以及文件读取路径格式

    绝对路径 :指的是是当前文件在 计算机磁盘中存放的具体位置 就是死的,物理上面的 相对路径 :指的是文件 相对于当前的py文件所处的位置 。就是参照py文件来说明路径,参照物嘛 读取文件路径方式(path.xls文件为例子)   执行命令的py文件同path.xls文件在同一个目录an

    2024年02月06日
    浏览(48)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包