dimfilter

贡献者:

周茂

最近更新日期:

2022-07-16


官方文档:

dimfilter

简介:

在空间域对网格进行方向中值(DIrectional Median, DIM)滤波

dimfilter 对空间域(或时间域)中的网格滤波。该模块将滤波范围(圆)划分 为多个扇区,对每个扇区进行一级或者二级滤波,并根据二级滤波选择最终的结果。 这种划分多个扇区的做法与普通中值滤波的区别在于可以避免在具有一定倾斜趋势的 区域中某些特征(以对海底地形滤波为例,特征可能为海山)导致的中值偏离的情况。 输出网格可以使用 -R 进行裁剪,或者使用 -I 设置新的网格间隔;裁剪可以避 免边缘效应。-Q 选项用于误差分析模式并且输入文件中包含滤波深度。 dimfilter 不会向其他滤波一样产生平滑的输出,因为其返回 N 个扇区中的 N 个中值中的最小值。除非数据数据中无噪声,否则输出结果可能很不平滑。因此,通 常建议对网格进行一次额外的滤波(例如,使用 grdfilter )。

语法

gmt dimfilter ingrid -D0-4 -Fxwidth[modifier] -Goutgrid -Nxsectors[modifier] [ -L ] [ -Q ] [ -Iincrement ] [ -Rregion ] [ -T ] [ -V[level] ] [ -fflags ] [ -hheaders ] [ --PAR=value ]

必选选项

ingrid[=ID|?varname][+bband][+ddivisor][+ninvalid][+ooffset][+sscale]

输入网格名。通过追加 =ID 可指定 网格格式 [默认为 =nf]; 追加 ?varname 可指定 NetCDF 变量 [默认为 GMT 找到的第一个 2-D 网格]。 对网格追加子选项可进行额外设置:

  • +b 选取一个波段 band(仅用于图片)[默认为 0]

  • +d 将网格除以一个数 divisor [默认为 1]

  • +n 将网格中的 invalid 替换为 NaN

  • +o 将网格中的值进行一定的偏移,即加一个数 offset [默认为 0]

  • +s 将网格中的值缩放,即乘以 scale [默认为 1]

offset 操作位于 scale 之前。

-Dflag

距离 flag 用来设置滤波相关的宽度的单位和类型等信息:

  • flag = 0 :网格(x,y) 的单位和 width 相同,使用笛卡尔距离

  • flag = 1 :网格(x,y) 的单位为度,width 单位为千米,使用笛卡尔距离

  • flag = 2 :网格(x,y) 的单位为度,width 单位为千米,dx 乘以 cos(lat), lat 为所有纬度中值,使用笛卡尔距离

上述计算都狠快,因为只需计算一次权重矩阵。下面的三个选项相对较慢,因为对于 每个纬度都需重新计算权重矩阵

  • flag = 3 :网格(x,y) 的单位为度,width 单位为千米,dx 乘以 cos(lat), lat 为对应的纬度,使用笛卡尔距离

  • flag = 4 :网格(x,y) 的单位为度,width 单位为千米,使用球面距离计

-Fxwidth[+l|u]

设置一级滤波类型,可从卷积和非卷积滤波中选择。x 为滤波类型代码,后面的 width 为滤波直径。

其中卷积滤波对应的代码包括:

  • b :Boxcar,即滤波窗口内所有点等权

  • c :Cosine Arch,滤波窗口内的权中为 Cosin 曲线

  • g :Gaussion,权重通过高斯函数给出

非卷积滤波对应的代码包括:

  • m :Median,返回中值

  • p :Maximum likelihood probability (一种众数估计),返回众数; 如果存在多个众数,则返回其平均数。追加 +l+u 可分别选择 多个众数中的最小值和最大值

-Nxsectors[+l|u]

设置二级滤波类型,以及扇区的个数 sectors ,扇区个数必须为整数并大于 0;当扇区格式为 0 时,则不使用二级滤波。

其中滤波对应的代码包括:

  • l :返回最小值

  • u :返回最大值

  • a :返回平均值

  • m :Median,返回中值;

  • p :Maximum likelihood probability (一种众数估计),返回众数; 如果存在多个众数,则返回其平均数。追加 +l+u 可分别选择 多个众数中的最小值和最大值

-Goutgrid[=ID][+ddivisor][+ninvalid][+ooffset|a][+sscale|a][:driver[dataType][+coptions]]

输出网格文件名。通过追加 =ID 可指定 网格格式。 对网格追加子选项可进行额外设置:

  • +d 将网格除以一个数 divisor [默认为 1]

  • +n 将网格中的 invalid 替换为 NaN

  • +o 将网格中的值进行一定的偏移,即加一个数 offset,或使用 a 自动对值 进行调整以保证整数网格的精度 [默认为 0]

  • +s 将网格中的值缩放,即乘以 scale,或使用 a 自动对网格缩放以保证整数 网格的精度 [默认为 1]

offset 操作位于 scale 之前;+sa 将会同时设置 +oa。如果需 使用 GDAL 指定网格格式,则 ID 应设置为 gd,并指定 :driver 和可选的数据 类型 dataType,以及 +coptions 选项传递给 GDAL 的 -co 选项。

可选选项

-Ixinc[+e|n][/yinc[+e|n]]

指定X和Y方向的网格间隔

  • xincyinc 为 X 和 Y 方向的网格间隔。对于地理坐标,可以指定网格间隔单位 [默认单位为度]

  • +e 微调X和Y方向范围的最大值,使得其是网格间隔的整数倍(默认会微调 网格间隔以适应给定的数据范围)

  • +n 表明 xincyinc 不是网格间隔,而是X和Y方向的节点数。 此时会根据节点数、网格区域范围以及网格配准方式重新计算网格间隔。

注意:

  • yinc 设置为0,则表示其与 xinc 相同

  • 若使用 -Rgrdfile 选项,则网格间隔和配准方式已经根据网格文件自动初始化, 此时依然可以使用 -I-r 覆盖相应的值

-L

使用此选项将 dim.template.sh 脚本内容写到标准输出,不与其他选项组合使用

-Rxmin/xmax/ymin/ymax[+r][+uunit] (more …)

指定数据范围

-T

转换输出网格的配准方式,使其与输入网格相反

-Q

该模式输入数据不是常见的网格格式,而是 grd2xyz-Z 选项输出的形式作为输入。使用该模式输出的结果为均值,MAD 以及中值, 可用于进一步的误差分析,可见 脚本模版

-V[level] (more …)

设置 verbose 等级 [w]

-f[i|o]colinfo (more …)

指定输入或输出列的数据类型

-h[i|o][n][+c][+d][+msegheader][+rremark][+ttitle] (more …)

跳过或生成指定数目的头段记录

-^-

显示简短的帮助信息,包括模块简介和基本语法信息(Windows下只能使用 -

-++

显示帮助信息,包括模块简介、基本语法以及模块特有选项的说明

-? 或无参数

显示完整的帮助信息,包括模块简介、基本语法以及所有选项的说明

--PAR=value

临时修改GMT参数的值,可重复多次使用。参数列表见 配置参数

地理坐标和时间坐标

当输出文件为 netCDF 格式时,根据输入数据、网格或者 -f 以及 -R 选项,会自动将输出坐标命名为“longitute”,“lattitude” 或者 “time”。例如: -f0x -f1t-R90w/90e/0t/3t 都会生成 longitude/time 网格。 当 x,y 或者 z 坐标为时间时,在网格中将会存储为相对时间,其相对于 TIME_EPOCHTIME_UNIT 指定的历元。其中时间变量的单位 属性将会默认与上述两个参数一致。

示例

远程文件 @earth_relief_05m 是一个 5 分分辨率的测深网格,对其使用 300 km 半径的中值滤波,设定范围为 150E ~ 250E,10N ~ 40N,输出结果分辨率 为 0.5 度。为了防止中值被地形坡度影响,将滤波窗口划分为 6 个扇形,从 6 个扇形中的中值中选择最小值,使用球面距离计算

gmt dimfilter @earth_relief_05m -Gfiltered_pacific.nc -Fm600 -D4 -Nl6 -R150/250/10/40 -I0.5 -V

假定文件 cape_verde.nc 是一个范围为 32W ~ 15W,8N ~ 25N,分辨率为 0.5 分才到测深文件,如果要去除该文件中的小尺度特征以观察 swell,为防止边缘 效应的影响,将输出结果限定为 27.5W ~ 20.5W,12.5N ~ 19.5N,分辨率为 2 分。使用笛卡尔距离计算

gmt dimfilter cape_verde.nc -Gt.nc -Fm220 -Nl8 -D2 -R-27.5/-20.5/12.5/19.5 -I2m -V
gmt grdfilter t.nc -Gcape_swell.nc -Fg50 -D2 -V

假定对某一区域的测深文件进行一系列的滤波,结果为 fxxx.nc,其中 xxx 表示 滤波半径,可对这一系列文件进行趋势分析,并得到每个点的 MAD(median absolute devuation)

gmt grd2xyz f100.nc -Z > f100.txt
gmt grd2xyz f110.nc -Z > f110.txt
gmt grd2xyz f120.nc -Z > f120.txt
gmt grd2xyz f130.nc -Z > f130.txt
paste f100.txt f110.txt f120.txt f130.txt > depths.txt
gmt dimfilter depths.txt -Q > output.z

上述 paste 为 Unix/Linux 命令,用于在水平方向合并多个文件。

注意事项

当输入为地理坐标网格时,卷积滤波都会对滤波权重进行归一化,以使用滤波 窗口大小的随纬度的变化。并能在跨越 360 度或者包含极点时进行正确的处理。 但是非卷积滤波器,不使用权重,因此只能在笛卡尔坐标网格使用。如果要使用 这些滤波,则应先进行投影后使用。

脚本模版

dim.template.sh 是一个脚本框架,可在其中进行完整的 dim 分析,包括 MAD 分析, 通过 -L 选项可以获取。

参考文献

Kim, S.-S., and Wessel, P. (2008), Directional Median Filtering for Regional-Residual Separation of Bathymetry, Geochem. Geophys. Geosyst., 9, Q03005, https://doi.org/10.1029/2007GC001850.

相关模块

grdfilter