前往小程序，Get更优阅读体验！
立即前往
文档建议反馈控制台
首页
学习
活动
专区
工具
TVP
最新优惠活动
文章/答案/技术大牛
发布
首页
学习
活动
专区
工具
TVP最新优惠活动
返回腾讯云官网
社区首页 >专栏 >回旋镖！meteva也能绘制wrfout气象要素分布

回旋镖！meteva也能绘制wrfout气象要素分布

作者头像
用户11172986
发布于 2024-06-20 19:09:41
820
发布于 2024-06-20 19:09:41
举报
文章被收录于专栏：气python风雨气python风雨

前言

博主在早期对meteva的使用写了一个笔记，就是meteva，这可能是气象萌新最需要的python库

在使用中发现它不能对有兰伯特投影的wrfout数据直接绘图，所以使用了其他库进行重新网格插值再绘图

今天在逛meteva的showdoc时刷新出了一个官方教程，大体是将wrfout数据转为pandas格式

然后使用idw进行插值绘图

下面让我们开始实践吧

温馨提示

由于可视化代码过长隐藏，可点击回旋镖！meteva也能绘制wrfout气象要素分布运行Fork查看 🔜🔜若没有成功加载可视化图，点击运行可以查看 ps:隐藏代码在【代码已被隐藏】所在行，点击所在行，可以看到该行的最右角，会出现个三角形，点击查看即可

导入库与读取数据

In [23]:

代码语言：javascript
复制
代码语言：javascript
复制
import xarray as xr
import matplotlib.pyplot as plt
#由于meteva函数调用的是宋体，当前镜像的matplotlib字体库无宋体，先设置现有的tff
plt.rcParams['font.sans-serif'] = ['Source Han Sans CN']  
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

path = "/home/mw/input/wrfout3385/wrfout_d02_2022-07-14_1100.nc"   
grd0 = xr.open_dataset(path)   #通过xarray读入网格数据
print(grd0)
代码语言：javascript
复制
代码语言：javascript
复制
<xarray.Dataset> Size: 51MB
Dimensions:                (Time: 1, south_north: 90, west_east: 90,
                            bottom_top: 49, bottom_top_stag: 50,
                            soil_layers_stag: 4, west_east_stag: 91,
                            south_north_stag: 91, seed_dim_stag: 8)
Coordinates:
    XLAT                   (Time, south_north, west_east) float32 32kB ...
    XLONG                  (Time, south_north, west_east) float32 32kB ...
    XTIME                  (Time) datetime64[ns] 8B ...
    XLAT_U                 (Time, south_north, west_east_stag) float32 33kB ...
    XLONG_U                (Time, south_north, west_east_stag) float32 33kB ...
    XLAT_V                 (Time, south_north_stag, west_east) float32 33kB ...
    XLONG_V                (Time, south_north_stag, west_east) float32 33kB ...
Dimensions without coordinates: Time, south_north, west_east, bottom_top,
                                bottom_top_stag, soil_layers_stag,
                                west_east_stag, south_north_stag, seed_dim_stag
Data variables: (12/213)
    Times                  (Time) |S19 19B ...
    LU_INDEX               (Time, south_north, west_east) float32 32kB ...
    ZNU                    (Time, bottom_top) float32 196B ...
    ZNW                    (Time, bottom_top_stag) float32 200B ...
    ZS                     (Time, soil_layers_stag) float32 16B ...
    DZS                    (Time, soil_layers_stag) float32 16B ...
    ...                     ...
    PCB                    (Time, south_north, west_east) float32 32kB ...
    PC                     (Time, south_north, west_east) float32 32kB ...
    LANDMASK               (Time, south_north, west_east) float32 32kB ...
    LAKEMASK               (Time, south_north, west_east) float32 32kB ...
    SST                    (Time, south_north, west_east) float32 32kB ...
    SST_INPUT              (Time, south_north, west_east) float32 32kB ...
Attributes: (12/134)
    TITLE:                            OUTPUT FROM WRF V4.4.2 MODEL
    START_DATE:                      2022-07-14_00:00:00
    SIMULATION_START_DATE:           2022-07-14_00:00:00
    WEST-EAST_GRID_DIMENSION:        91
    SOUTH-NORTH_GRID_DIMENSION:      91
    BOTTOM-TOP_GRID_DIMENSION:       50
    ...                              ...
    ISLAKE:                          21
    ISICE:                           15
    ISURBAN:                         13
    ISOILWATER:                      14
    HYBRID_OPT:                      2
    ETAC:                            0.2

转为pandas

In [25]:

代码语言：javascript
复制
代码语言：javascript
复制
import pandas as pd
import numpy as np
sta = pd.DataFrame({"level":0, "time":grd0.coords["XTIME"].values[0],"dtime":0,
                     "id":np.arange(grd0.coords["XLAT"].values.size),
                    "lon":grd0.coords["XLONG"].values.flatten(),
                   "lat":grd0.coords["XLAT"].values.flatten(),
                  "RAIN":grd0.variables["RAINC"].values.flatten()+grd0.variables["RAINNC"].values.flatten()
                 })
代码语言：javascript
复制
print(sta)
代码语言：javascript
复制
      level                time  dtime    id         lon        lat       RAIN
0         0 2022-07-14 11:00:00      0     0  100.901947  27.242496   0.323906
1         0 2022-07-14 11:00:00      0     1  101.002350  27.238125   0.238994
2         0 2022-07-14 11:00:00      0     2  101.102722  27.233673   0.315282
3         0 2022-07-14 11:00:00      0     3  101.203064  27.229111   0.267089
4         0 2022-07-14 11:00:00      0     4  101.303406  27.224480   0.215387
...     ...                 ...    ...   ...         ...        ...        ...
8095      0 2022-07-14 11:00:00      0  8095  110.649567  34.472286  19.087238
8096      0 2022-07-14 11:00:00      0  8096  110.757538  34.459393  17.735456
8097      0 2022-07-14 11:00:00      0  8097  110.865448  34.446400  15.124860
8098      0 2022-07-14 11:00:00      0  8098  110.973358  34.433327  13.774797
8099      0 2022-07-14 11:00:00      0  8099  111.081207  34.420151  13.101636

[8100 rows x 7 columns]

IDW插值

In [26]:

代码语言：javascript
复制
代码语言：javascript
复制
import meteva.base as meb
grid = meb.grid([100,112,0.05],[28,36,0.05])  # 设置一个等经纬度的网格信息变量
grd1 = meb.interp_sg_idw(sta,grid = grid,nearNum=4,effectR=20)  # 将站点数据插值到网格上
print(grd1)
代码语言：javascript
复制
代码语言：javascript
复制
<xarray.DataArray 'data0' (member: 1, level: 1, time: 1, dtime: 1, lat: 161,
                           lon: 241)> Size: 310kB
array([[[[[[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]]]]]])
Coordinates:
  * member   (member) <U4 16B 'RAIN'
  * level    (level) int64 8B 0
  * time     (time) datetime64[ns] 8B 2022-07-14T11:00:00
  * dtime    (dtime) int64 8B 0
  * lat      (lat) float64 1kB 28.0 28.05 28.1 28.15 ... 35.85 35.9 35.95 36.0
  * lon      (lon) float64 2kB 100.0 100.0 100.1 100.2 ... 111.9 112.0 112.0
Attributes:
    data_start_columns:  6

绘图

In [27]:

代码语言：javascript
复制
代码语言：javascript
复制
import cmaps
cmaps0=cmaps.prcp_1
meb.contourf_2d_grid(grd1,cmap=cmaps0) #绘制插值结果
代码语言：javascript
复制
代码语言：javascript
复制

高度层绘制

按照官方的思路，我们对高度层的混合比数据进行绘制

In [28]:

代码语言：javascript
复制
代码语言：javascript
复制
from netCDF4 import Dataset
import numpy as np
import pandas as pd
from wrf import (to_np, getvar, ALL_TIMES, CoordPair, interplevel)
#由于meteva函数调用的是宋体，当前镜像的matplotlib字体库无宋体，先设置现有的tff
plt.rcParams['font.sans-serif'] = ['Source Han Sans CN']  
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

wrf_file = '/home/mw/input/wrfout3385/wrfout_d02_2022-07-14_0600.nc'

ncfile = Dataset(wrf_file)

# 使用wrf-python的getvar函数获取整个时间序列的湿度数据
qvapor = getvar(ncfile, 'QVAPOR', timeidx=ALL_TIMES, method='cat')

# 定义需要插值的目标气压值，这里为500hPa
target_plev = 500.0  # in hPa

# 获取模型中的气压数据，这是进行插值所必需的
pressure = getvar(ncfile, 'P')

# 使用interplevel函数进行垂直插值，得到500hpa高度的QVAPOR数据
qvapor_500hpa = interplevel(qvapor, pressure, target_plev)

# 获取必要的坐标信息
time_coord = ncfile.START_DATE
lat = getvar(ncfile, 'XLAT').values
lon = getvar(ncfile, 'XLONG').values

# 将数据整理到Pandas DataFrame
sta1 = pd.DataFrame({
    "level": target_plev,  # 直接指定为500hPa
    "time": [time_coord][0],  "dtime":0,
    "id": np.arange(lat.size).flatten(),  # 网格点ID，可选
    "lon": lon.flatten(),
    "lat": lat.flatten(),
    "QVAPOR": qvapor_500hpa.values.flatten(),  # 使用插值后的QVAPOR作为湿度数据
})

print(sta1.head())  # 打印DataFrame的前几行以检查结果

# 关闭NetCDF文件
代码语言：javascript
复制
ncfile.close()
代码语言：javascript
复制
   level                 time  dtime  id         lon        lat    QVAPOR
0  500.0  2022-07-14_00:00:00      0   0  100.901947  27.242496  0.001734
1  500.0  2022-07-14_00:00:00      0   1  101.002350  27.238125  0.001666
2  500.0  2022-07-14_00:00:00      0   2  101.102722  27.233673  0.001562
3  500.0  2022-07-14_00:00:00      0   3  101.203064  27.229111  0.001444
4  500.0  2022-07-14_00:00:00      0   4  101.303406  27.224480  0.001443

In [29]:

代码语言：javascript
复制
代码语言：javascript
复制
grid = meb.grid([100,112,0.05],[28,35,0.05])  # 设置一个等经纬度的网格信息变量
grd1 = meb.interp_sg_idw(sta1,grid = grid,nearNum=4,effectR=20)  # 将站点数据插值到网格上
代码语言：javascript
复制
print(grd1)
代码语言：javascript
复制
<xarray.DataArray 'data0' (member: 1, level: 1, time: 1, dtime: 1, lat: 141,
                           lon: 241)> Size: 272kB
array([[[[[[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]]]]]])
Coordinates:
  * member   (member) <U6 24B 'QVAPOR'
  * level    (level) int64 8B 500
  * time     (time) datetime64[ns] 8B 2022-07-14
  * dtime    (dtime) int64 8B 0
  * lat      (lat) float64 1kB 28.0 28.05 28.1 28.15 ... 34.85 34.9 34.95 35.0
  * lon      (lon) float64 2kB 100.0 100.0 100.1 100.2 ... 111.9 112.0 112.0
Attributes:
    data_start_columns:  6
代码语言：javascript
复制

下面我绘制了基于cartopy的两种版本的500hPa混合比分布，可作对比

等经纬度网格版本

兰伯特版本

小结

能摸索出新的方法对数据处理绘图是愉快的

图片空白部分不是数据有误或者作图有错，而是地形的原因

总的来说meteva绘图一句话，但是前期的数据转换较为复杂 ，还是会劝退萌新

本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2024-05-03，如有侵权请联系 cloudcommunity@tencent.com 删除
函数
可视化
数据
time
变量

本文分享自 气python风雨 微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ，欢迎热爱写作的你一起参与！

函数
可视化
数据
time
变量
评论
登录后参与评论
0 条评论
热度
最新
登录 后参与评论
推荐阅读
LV.
文章
0
获赞
0
目录
  • 前言
  • 温馨提示
  • 导入库与读取数据
  • 转为pandas
  • IDW插值
  • 绘图
  • 高度层绘制
  • 等经纬度网格版本
  • 兰伯特版本
  • 小结
领券

腾讯云开发者

扫码关注腾讯云开发者

扫码关注腾讯云开发者

领取腾讯云代金券

热门产品

热门推荐

更多推荐

Copyright © 2013 - 2024 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有 

深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 深公网安备号 44030502008569

腾讯云计算（北京）有限责任公司 京ICP证150476号 |  京ICP备11018762号 | 京公网安备号11010802020287

问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档

Copyright © 2013 - 2024 Tencent Cloud.

All Rights Reserved. 腾讯云 版权所有

登录 后参与评论