[GRASS-SVN] r66996 - grass-addons/grass7/raster/r.series.filter
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Dec 1 07:57:51 PST 2015
Author: DmitryKolesov
Date: 2015-12-01 07:57:51 -0800 (Tue, 01 Dec 2015)
New Revision: 66996
Modified:
grass-addons/grass7/raster/r.series.filter/r.series.filter.html
grass-addons/grass7/raster/r.series.filter/r.series.filter.py
Log:
r.series.filter: add median filter
Modified: grass-addons/grass7/raster/r.series.filter/r.series.filter.html
===================================================================
--- grass-addons/grass7/raster/r.series.filter/r.series.filter.html 2015-12-01 10:11:10 UTC (rev 66995)
+++ grass-addons/grass7/raster/r.series.filter/r.series.filter.html 2015-12-01 15:57:51 UTC (rev 66996)
@@ -10,7 +10,8 @@
<p>
<em>result_prefix</em> Prefix for raster names of filterd <em>X</em>.
<p>
-<em>method</em> filtering method. Savitzky-Golay filter is only implemented now.
+<em>method</em> filtering method. Implemented filters are Savitzky-Golay filter <em>savgol</em>
+ and median filter <em>median</em>.
<p>
<em>winsize</em> The length of the filter window. <em>winsize</em> must be a positive odd integer.
<p>
@@ -32,20 +33,20 @@
<p>
There is a procedure for searching for good filtering parameters: it uses <em>opt_points</em>
random points and perfoms filtering in that points. The result of the filtering can be tested
-for quality. The quality function is a tradeof two features: accuracy and smoothing.
+for quality. The quality function is a trade of two features: accuracy and smoothing.
Accuracy can be estimated as the (abs) difference between original and filtered data,
quality of smoothing can be estimated as absalute values of the derivates. So there are two
-parameters <em>diff_penalty</em> and <em>deriv_penalty</em> that can ajust the tradeof.
+parameters <em>diff_penalty</em> and <em>deriv_penalty</em> that can ajust the trade-of.
<p>
So the optimizing procedure performs loop over filtering parameters and calculates
- thequality function
+the next penalty function:
<div class="code"><pre>
-penalty =quality diff_penalty * sum(abs(Xi-Fi)) + sum(abs(dFi))
+penalty = diff_penalty * sum(abs(Xi-Fi)) + sum(abs(dFi))
</pre></div>
where <em>Xi</em> are original signals in the samplig points, <em>Fi</em>
are filtered signals in the sampling points.
<p>
-Optimal parameters are used for signal filtering in the whole region.
+The optimal parameters are used for signal filtering in the whole region.
<h2>EXAMPLES</h2>
Modified: grass-addons/grass7/raster/r.series.filter/r.series.filter.py
===================================================================
--- grass-addons/grass7/raster/r.series.filter/r.series.filter.py 2015-12-01 10:11:10 UTC (rev 66995)
+++ grass-addons/grass7/raster/r.series.filter/r.series.filter.py 2015-12-01 15:57:51 UTC (rev 66996)
@@ -49,7 +49,7 @@
#% multiple: no
#% answer: savgol
#% description: Used method
-#% descriptions: savgol; Savitzky–Golay filter
+#% descriptions: savgol; Savitzky–Golay filter; median; Median filter
#%end
#%option
#% key: winsize
@@ -101,7 +101,9 @@
import numpy as np
from scipy.signal import savgol_filter
+from scipy.signal import medfilt
+
if "GISBASE" not in os.environ:
sys.stderr.write("You must be in GRASS GIS to run this program.\n")
sys.exit(1)
@@ -147,12 +149,17 @@
r.close()
-def _filter(row_data, winsize, order):
+def _filter(method, row_data, winsize, order):
result = np.empty(row_data.shape)
_, cols = row_data.shape
for i in range(cols):
arr = _fill_nulls(row_data[:, i])
- arr = savgol_filter(arr, winsize, order, mode='nearest')
+ if method == 'savgol':
+ arr = savgol_filter(arr, winsize, order, mode='nearest')
+ elif method == 'median':
+ arr = medfilt(arr, kernel_size=winsize)
+ else:
+ grass.fatal('The method is not implemented')
result[:, i] = arr
return result
@@ -183,7 +190,7 @@
return diff_penalty * difference + deriv_penalty * deriv_diff
-def optimize_params(names, npoints, diff_penalty, deriv_penalty):
+def optimize_params(method, names, npoints, diff_penalty, deriv_penalty):
"""Perform crossvalidation:
take 'npoints' random points,
find winsize and order that minimize the quality function
@@ -216,23 +223,58 @@
close_rasters(inputs)
# Find the optima
+ best_winsize = best_order = None
+ if method == 'savgol':
+ best_winsize, best_order = _optimize_savgol(input_data)
+ elif method == 'median':
+ best_winsize = _optimize_median(input_data)
+ else:
+ grass.fatal('The method is not implemented')
+
+ return best_winsize, best_order
+
+
+def _optimize_savgol(input_data, diff_penalty, deriv_penalty):
+ """Find optimal params for savgol_filter.
+
+ Returns winsize and order
+ """
+ map_count, npoints = input_data.shape
best = np.inf
best_winsize = best_order = None
for winsize in range(5, map_count/2, 2):
for order in range(2, min(winsize - 2, 10)): # 10 is a 'magic' number: we don't want very hight polynomyal fitting usually
- for i in range(npoints):
- test_data = np.copy(input_data)
- test_data = _filter(test_data, winsize, order)
- penalty = fitting_quality(input_data, test_data,
- diff_penalty, deriv_penalty)
- if penalty < best:
- best = penalty
- best_winsize, best_order = winsize, order
+ test_data = np.copy(input_data)
+ test_data = _filter('savgol', test_data, winsize, order)
+ penalty = fitting_quality(input_data, test_data,
+ diff_penalty, deriv_penalty)
+ if penalty < best:
+ best = penalty
+ best_winsize, best_order = winsize, order
return best_winsize, best_order
+def _optimize_median(input_data, diff_penalty, deriv_penalty):
+ """Find optimal params for median filter.
-def filter(names, winsize, order, prefix):
+ Returns winsize
+ """
+ map_count, npoints = input_data.shape
+ best = np.inf
+ best_winsize = order = None
+ for winsize in range(3, map_count/2, 2):
+ test_data = np.copy(input_data)
+ test_data = _filter('median', test_data, winsize, order)
+ penalty = fitting_quality(input_data, test_data,
+ diff_penalty, deriv_penalty)
+ if penalty < best:
+ best = penalty
+ best_winsize = winsize
+
+ return best_winsize
+
+
+def filter(method, names, winsize, order, prefix):
inputs = init_rasters(names)
output_names = [prefix + name for name in names]
outputs = init_rasters(output_names)
@@ -244,7 +286,7 @@
for i in range(reg.rows):
# import ipdb; ipdb.set_trace()
row_data = np.array([_get_row_or_nan(r, i) for r in inputs])
- filtered_rows = _filter(row_data, winsize, order)
+ filtered_rows = _filter(method, row_data, winsize, order)
for map_num in range(len(outputs)):
map = outputs[map_num]
row = filtered_rows[map_num, :]
@@ -278,6 +320,8 @@
optimize = flags['c']
+ method = options['method']
+
xnames = options['input']
xnames = xnames.split(',')
@@ -313,13 +357,13 @@
grass.error("Order of the filter must be less than window length")
if optimize:
- winsize, order = optimize_params(xnames, opt_points,
+ winsize, order = optimize_params(method, xnames, opt_points,
diff_penalty, deriv_penalty)
if winsize is None:
grass.error("Optimization procedure doesn't convergence.")
sys.exit(1)
- filter(xnames, winsize, order, res_prefix)
+ filter(method, xnames, winsize, order, res_prefix)
if __name__ == "__main__":
options, flags = grass.parser()
More information about the grass-commit
mailing list