- 积分
- 10278
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2019-7-29
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 雨落森林 于 2024-7-15 17:34 编辑
Barnes是一种常用的空间滤波方法,主要使用两个常数g和c来计算高斯权重,并对每个网格点进行空间插值,从而成为滤除高频波动的低通滤波器。当使用两种不同的常数g和c方案时,两者是分别保留了不同尺度的低频波动。而两种方法的滤波结果做差就能得到中尺度的波动了。
在家园能找到南信大王咏青老师的Barnes带通滤波的fortran程序。但是由于fortran不方便用,因此我改写成了python代码并且发布了pybarnes库,可以直接用pip install pybarnes来安装pybarnes库(需要安装一个metpy作为前置条件),关于这个库,也可以去见github主页:https://github.com/LinOuyang/pybarnes。
现在我们用Python来实现,从库里面导入一个BarnesFilter,将输入原始场(比如500hpa位势高度z)放入其中,得到滤波器f。可以选择调用其中的lowpass设置常数g和c进行低通滤波。或者直接用其中的bandpass,需要设置g1、c1、g2、c2,然后才能得到中尺度的波动。
接下来直接看代码。
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib.colors import TwoSlopeNorm
from pybarnes import BarnesFilter
# 读取数据,选区东亚区域,并且分辨率0.5度,确保数据不要太大
ds = xr.open_dataset("./era5uvz500.nc").sortby("latitude")
ds = ds.sel(longitude=slice(70, 140), latitude=slice(0, 60))
ds = ds.isel(longitude=slice(None, None, 2), latitude=slice(None, None, 2))
z1 = ds.z[0]/98
# 创建Barnes滤波器,当输入的数据是xr的array或者dataset,可以不指定lon, lat,否则要指定。
f = BarnesFilter(z1)
# 用g=0.3, c=150000做低通滤波
z2 = f.lowpass(g=0.3, c=150000)
# 用两套g和c方案做差得到中尺度波动
z3 = f.bandpass(g1=0.3, c1=30000, g2=0.3, c2=150000)
# 画图看看效果
plt.figure(figsize=(16, 4.8))
plt.subplot(131)
plt.colorbar(plt.contourf(z1, levels=np.arange(560, 590), extend="both", cmap=plt.cm.RdYlBu_r))
plt.title("Raw data field")
plt.subplot(132)
plt.title("The filtered field (g=0.3,c=150000)")
plt.colorbar(plt.contourf(z2, levels=np.arange(560, 590), extend="both", cmap=plt.cm.RdYlBu_r))
plt.subplot(133)
plt.title("Mesoscale fluctuations field")
plt.colorbar(plt.contourf(z3, cmap=plt.cm.bwr, norm=TwoSlopeNorm(0)))
plt.savefig("./barnes_filter.jpg", dpi=300, bbox_inches="tight")
plt.show()
## 2024年4月更新 ##
推出了0.2版本:
1. 新增对于一维的站点数据的支持(BarnesFilter1d)
2. 修复0.1版本对于部分形状的数组不支持的重大bug
3. 优化了内存使用机制,在不超过内存的情况下尽可能加快运算
* 使用最多的是对wrfout数据做barnes滤波,使用方法如下:
比如对slp做滤波:
slp = getwar(ncfile, "slp")
f = BarnesFilter(slp.data, lon=slp.XLONG.data, lat=slp.XLAT.data)
z2 = f.lowpass(g=0.3, c=150000)
z3 = f.bandpass(g1=0.3, c1=30000, g2=0.3, c2=150000)
|
评分
-
查看全部评分
|