grdfilter
- 贡献者:
- 最近更新日期:
2022-07-18
- 官方文档:
- 简介:
在空间(或时间)域中对网格滤波
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 指定的历元。其中时间变量的单位属性将会默认与上述两个参数一致。
注意事项
示例
@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