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方向的网格间隔

xinc 和 yinc 为 X 和 Y 方向的网格间隔。对于地理坐标，可以指定网格间隔单位 [默认单位为度]
+e 微调X和Y方向范围的最大值，使得其是网格间隔的整数倍（默认会微调网格间隔以适应给定的数据范围）
+n 表明 xinc 和 yinc 不是网格间隔，而是X和Y方向的节点数。此时会根据节点数、网格区域范围以及网格配准方式重新计算网格间隔。

注意：

若 yinc 设置为0，则表示其与 xinc 相同
若使用 -Rgrdfile 选项，则网格间隔和配准方式已经根据网格文件自动初始化，此时依然可以使用 -I 和 -r 覆盖相应的值

-Ni|p|r

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

i 计算过程中忽略所有 NaN 值 [默认]
r 与 i 相同，但如果输出网格与输入配准相同，则在输入网格的 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_EPOCH 和 TIME_UNIT 指定的历元。其中时间变量的单位属性将会默认与上述两个参数一致。

注意事项

使用 -D5 选项时，输入墨卡托网格必须由 img2grd 的 -C 选项生成，因此 y 值的原点为赤道（即 x = y = 0 对应于 lon = lat = 0）
-I 选项设置的新的 x_inc ， y_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

grdfilter

语法

必选选项

可选选项

地理坐标和时间坐标

注意事项

示例

相关模块