grdfilter

贡献者:

周茂

最近更新日期:

2022-07-18


官方文档:

grdfilter

简介:

在空间(或时间)域中对网格滤波

grdfilter 使用卷积、非卷积各向同性或矩形滤波器对网格进行滤波。输出网格可以设置新的区域(-R),间隔(-I)或者配准方式(-T)。这样可以在数据足够大的情况下,可去掉网格边缘以避免滤波的边缘效应。如果滤波器为低通滤波,则输出结果的频率可能低于输入的采样频率。:频域(或称波数域)滤波,见 grdfft

语法

gmt grdfilter ingrid -Ddistance_flag -Fxwidth[/width2][modifiers] -Goutgrid [ -Iincrement ] [ -Ni|p|r ] [ -Rregion ] [ -T ] [ -V[level] ] [ -fflags ] [ -rreg ] [ --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 = p :网格(px,py) 的单位为像素( width 为奇数),使用笛卡尔距离

  • 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 单位为千米,使用球面距离计算

  • flag = 5 :网格(x,y) 的单位为墨卡托 -Jm1 img 单位,width 为千米,使用球面距离计算

-Fxwidth[/width2][modifiers]

设置滤波类型,可从卷积和非卷积滤波中选择。x 为滤波类型代码,后面的 width 为滤波直径,此时进行各向同性滤波;追加 width2 以实现不同方向不同的滤波长度(需要 -Dp-D0 )。默认情况下,执行低通滤波。追加 +h 选项选择高通滤波。对于各向同性滤波,width 可以是一个网格,用来提供权重,这种情况下,网格必须与输出网格具有相同的配准。

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

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

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

  • g :Gaussion,权重通过高斯函数给出,其中宽度是传统高斯标准差的 6 倍

  • f :Custom,权重由网格文件 width 给出,维度必须为奇数,且需 -D0 ,输出网格的间隔必须与输入一致或为其整数倍

  • o :Operator,权重由网格文件 width 给出,维度必须为奇数,且需 -D0 ,输出网格的间隔必须与输入一致或为其整数倍。不同之处在于,假定权重总和为 0,因此不会进行归一化

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

  • m :Median,返回中值;追加 +qquantile 可选择其他的分位数, quantile 范围必须在 0-1 之间,默认为 0.5,即中位数

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

  • h :Histogram mode(另一种众数估计),返回直方图主峰的中心值,追加 /binwidth 来指定直方图每个分区的间隔。如果存在多个众数,则返回其平均数。追加 +l+u 可分别选择多个众数中的最小值和最大值。

  • l :返回最小值

  • L :返回正值中的最小值

  • u :返回最大值

  • U :返回负值中的最大值

L|U 情况下,可能最终的结果不存在数据,则范围 NaN。

-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 覆盖相应的值

-Ni|p|r

确定输入网格中的 NaN 值如何处理。

  • i 计算过程中忽略所有 NaN 值 [默认]

  • ri 相同,但如果输出网格与输入配准相同,则在输入网格的 NaN 值出现的位置将输出网格对应值也设置为 NaN

  • p 如果在滤波搜索范围内出现 NaN 值,则将最终网格点设置为 NaN

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

指定数据范围

-T

转换网格配准方式。或者使用 -r[g|p] 显式指定输出网格的配准方式

-V[level] (more …)

设置 verbose 等级 [w]

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

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

-r[g|p] (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 指定的历元。其中时间变量的单位属性将会默认与上述两个参数一致。

注意事项

  1. 使用 -D5 选项时,输入墨卡托网格必须由 img2grd-C 选项生成,因此 y 值的原点为赤道(即 x = y = 0 对应于 lon = lat = 0)

  2. -I 选项设置的新的 x_incy_inc 如果不是输入数据的间隔的整数倍,则滤波会非常慢

示例

@earth_relief_05m 为 5 分分辨率的远程 DEM 文件,对该文件使用 300 km 为半径的中值滤波;输出范围为 150 E 到 250 E,10 N 到 40 N,输出网格大小设置为 0.5 度;计算过程中使用球面距离

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

如果想执行高通滤波

gmt grdfilter @earth_relief_05m -Gresidual_pacific.nc -Fm600+h -D4 -R150/250/10/40 -I0.5 -V

使用自定义的各向同性滤波器 exp (-0.5*r^2) 对 ripple.nc 滤波,其中距中心的距离 r 为 (2x^2 + y^2 -2xy)/6

gmt grdmath -R-10/10/-10/10 -I1 X 2 POW 2 MUL Y 2 POW ADD X Y MUL 2 MUL \
        SUB 6 DIV NEG 2 DIV EXP DUP SUM DIV = gfilter.nc
gmt grdfilter ripples.nc -Ffgfilter.nc -D0 -Gsmooth.nc -V

相关模块

grdfft, img2grd