爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 5596|回复: 19

[经验总结] Barnes滤波

[复制链接]

新浪微博达人勋

发表于 2023-4-8 18:04:36 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

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)






barnes_filter.jpg

CodesAndDatas.zip

5.36 MB, 下载次数: 58, 下载积分: 金钱 -5

评分

参与人数 1金钱 +20 贡献 +1 收起 理由
付亚男 + 20 + 1 很给力!

查看全部评分

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2023-4-8 18:08:42 | 显示全部楼层
强!以前还用过Barnes滤波发表过文章
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2023-4-8 20:41:13 | 显示全部楼层
学到啦学到啦
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-4-25 20:47:25 | 显示全部楼层
提示: 作者被禁止或删除 内容自动屏蔽
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-10-27 07:45:42 | 显示全部楼层
ayzqs 发表于 2023-4-8 18:08
强!以前还用过Barnes滤波发表过文章

请问是哪篇?能推荐不,想引用一下!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-10-29 10:42:37 来自手机 | 显示全部楼层
对时间维做滤波的有推荐吗?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2023-10-29 23:51:09 | 显示全部楼层
雪域小丹 发表于 2023-10-29 10:42
对时间维做滤波的有推荐吗?

时间滤波我没做过就是,但其实网上的教程也挺多的,lancozos以及butterworth滤波都挺不错的,可以搜得到。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-10-31 10:53:54 | 显示全部楼层
好的,感谢大佬!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2024-5-25 21:15:56 | 显示全部楼层
楼主,想请教一下,这个函数对WRF模式数据做滤波的时候,为什么一大部分数据会变成0呢
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2024-5-26 10:22:24 | 显示全部楼层
镜灵17 发表于 2024-5-25 21:15
楼主,想请教一下,这个函数对WRF模式数据做滤波的时候,为什么一大部分数据会变成0呢

可能是输入没输入对,你可以详细介绍一下你是怎么做的吗?比如pybarnes是什么版本?以及,你是如何输入的?输入里面有缺测吗?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表