[GRASS-SVN] r69763 - grass-addons/grass7/raster/r.series.lwr
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Nov 2 14:50:43 PDT 2016
Author: mmetz
Date: 2016-11-02 14:50:43 -0700 (Wed, 02 Nov 2016)
New Revision: 69763
Added:
grass-addons/grass7/raster/r.series.lwr/r.series.lwr.html
Removed:
grass-addons/grass7/raster/r.series.lwr/r.hants.html
Modified:
grass-addons/grass7/raster/r.series.lwr/Makefile
grass-addons/grass7/raster/r.series.lwr/main.c
Log:
new grass addon r.series.lwr
Modified: grass-addons/grass7/raster/r.series.lwr/Makefile
===================================================================
--- grass-addons/grass7/raster/r.series.lwr/Makefile 2016-11-02 21:47:26 UTC (rev 69762)
+++ grass-addons/grass7/raster/r.series.lwr/Makefile 2016-11-02 21:50:43 UTC (rev 69763)
@@ -1,6 +1,6 @@
MODULE_TOPDIR = ../..
-PGM = r.hants
+PGM = r.series.lwr
LIBES = $(GMATHLIB) $(RASTERLIB) $(GISLIB)
DEPENDENCIES = $(GMATHDEP) $(RASTERDEP) $(GISDEP)
Modified: grass-addons/grass7/raster/r.series.lwr/main.c
===================================================================
--- grass-addons/grass7/raster/r.series.lwr/main.c 2016-11-02 21:47:26 UTC (rev 69762)
+++ grass-addons/grass7/raster/r.series.lwr/main.c 2016-11-02 21:50:43 UTC (rev 69763)
@@ -1,13 +1,11 @@
/****************************************************************************
*
- * MODULE: r.hants
- * AUTHOR(S): Wout Verhoef, NLR, Remote Sensing Dept., June 1998
- * Mohammad Abouali, converted to MATLAB, 2011
- * Markus Metz, converted to C, 2013
- * based on r.series
+ * MODULE: r.series.lwr
+ * AUTHOR(S): Markus Metz
+ * based on r.hants
* PURPOSE: time series correction
- * COPYRIGHT: (C) 2013 by the GRASS Development Team
+ * COPYRIGHT: (C) 2016 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
@@ -26,8 +24,6 @@
#include <grass/glocale.h>
#include <grass/gmath.h>
-/* TODO: add option for amplitude and phase output */
-
struct input
{
const char *name;
@@ -42,6 +38,87 @@
DCELL *buf;
};
+static double uniform(double ref, double x, double max)
+{
+ double dist = fabs(x - ref) / max;
+
+ if (dist < 1)
+ return 1;
+
+ return 0;
+}
+
+static double triangular(double ref, double x, double max)
+{
+ double dist = fabs(x - ref) / max;
+
+ if (dist == 0)
+ return 1;
+
+ if (dist >= 1)
+ return 0;
+
+ return (1.0 - dist);
+}
+
+static double epanechnikov(double ref, double x, double max)
+{
+ double dist = fabs(x - ref) / max;
+
+ if (dist == 0)
+ return 1;
+
+ if (dist >= 1)
+ return 0;
+
+ return (1.0 - dist * dist);
+}
+
+static double quartic(double ref, double x, double max)
+{
+ double dist = fabs(x - ref) / max;
+ double tmp;
+
+ if (dist == 0)
+ return 1;
+
+ if (dist >= 1)
+ return 0;
+
+ tmp = 1.0 - dist * dist;
+
+ return (tmp * tmp);
+}
+
+static double tricube(double ref, double x, double max)
+{
+ double dist = fabs(x - ref) / max;
+ double tmp;
+
+ if (dist == 0)
+ return 1;
+
+ if (dist >= 1)
+ return 0;
+
+ tmp = 1.0 - dist * dist * dist;
+
+ return (tmp * tmp * tmp);
+}
+
+static double cosine(double ref, double x, double max)
+{
+ double dist = fabs(x - ref) / max;
+
+ if (dist == 0)
+ return 1;
+
+ if (dist >= 1)
+ return 0;
+
+ return (cos(M_PI_2 * dist));
+}
+
static int solvemat(double **m, double a[], double B[], int n)
{
int i, j, i2, j2, imark;
@@ -68,7 +145,7 @@
/* co-linear points results in a solution with rounding error */
if (pivot == 0.0) {
- G_warning(_("Matrix is unsolvable"));
+ G_debug(4, "Matrix is unsolvable");
return 0;
}
@@ -115,46 +192,65 @@
return 1;
}
+static double term(int term, double x)
+{
+ switch (term) {
+ case 0:
+ return 1.0;
+ case 1:
+ return x;
+ case 2:
+ return x * x;
+ case 3:
+ return x * x * x;
+ }
+ return 0.0;
+}
+
int main(int argc, char *argv[])
{
struct GModule *module;
struct
{
struct Option *input, *file, *suffix,
- *amp, /* prefix for amplitude output */
- *phase, /* prefix for phase output */
- *nf, /* number of harmonics */
+ *order, /* polynomial order */
+ *weight, /* weighing function */
*fet, /* fit error tolerance */
+ *ts, /* time steps */
*dod, /* degree of over-determination */
- *range, /* low/high threshold */
- *ts, /* time steps*/
- *bl, /* length of base period */
+ *range, /* range of valid values */
*delta; /* threshold for high amplitudes */
} parm;
struct
{
struct Flag *lo, *hi, *lazy;
} flag;
- int i, j, k;
- int num_inputs;
+ int i, j, k, n;
+ int num_inputs, in_lo, in_hi;
struct input *inputs = NULL;
int num_outputs;
struct output *outputs = NULL;
- struct output *out_amp = NULL;
- struct output *out_phase = NULL;
char *suffix;
struct History history;
- DCELL *values = NULL, *rc = NULL;
+ struct Colors colors;
+ DCELL *values = NULL, *values2 = NULL;
int nrows, ncols;
int row, col;
- double lo, hi, fet, *cs, *sn, *ts, delta;
- int bl;
- double **mat, **mat_t, **A, *Av, *Azero, *za, *zr, maxerrlo, maxerrhi;
- int asize;
- int dod, nf, nr, nout, noutmax;
- int rejlo, rejhi, *useval;
- int do_amp, do_phase;
+ int order;
+ double fet, maxerrlo, maxerrhi, lo, hi;
+ DCELL *resultn;
+ int msize;
+ double **m, **m2, *a, *a2, *B;
+ double *ts, weight;
+ double (*weight_func)(double, double, double);
+ int dod;
+ int *isnull, n_nulls;
+ int min_points, n_points;
+ int this_margin;
+ double max_ts, tsdiff1, tsdiff2;
+ double delta;
+ int rejlo, rejhi;
G_gisinit(argv[0]);
@@ -178,28 +274,25 @@
parm.suffix->key = "suffix";
parm.suffix->type = TYPE_STRING;
parm.suffix->required = NO;
- parm.suffix->answer = "_hants";
+ parm.suffix->answer = "lwr";
parm.suffix->label = _("Suffix for output maps");
parm.suffix->description = _("The suffix will be appended to input map names");
- parm.amp = G_define_option();
- parm.amp->key = "amplitude";
- parm.amp->type = TYPE_STRING;
- parm.amp->required = NO;
- parm.amp->label = _("Prefix for output amplitude maps");
+ parm.order = G_define_option();
+ parm.order->key = "order";
+ parm.order->type = TYPE_INTEGER;
+ parm.order->required = NO;
+ parm.order->options = "1,2,3";
+ parm.order->answer = "2";
+ parm.order->description = _("order number");
- parm.phase = G_define_option();
- parm.phase->key = "phase";
- parm.phase->type = TYPE_STRING;
- parm.phase->required = NO;
- parm.phase->label = _("Prefix for output phase maps");
+ parm.weight = G_define_option();
+ parm.weight->key = "weight";
+ parm.weight->type = TYPE_STRING;
+ parm.weight->required = NO;
+ parm.weight->options = "uniform,triangular,epanechnikov,quartic,tricube,cosine";
+ parm.weight->answer = "tricube";
- parm.nf = G_define_option();
- parm.nf->key = "nf";
- parm.nf->type = TYPE_INTEGER;
- parm.nf->required = YES;
- parm.nf->description = _("Number of frequencies");
-
parm.fet = G_define_option();
parm.fet->key = "fet";
parm.fet->type = TYPE_DOUBLE;
@@ -211,7 +304,7 @@
parm.dod->key = "dod";
parm.dod->type = TYPE_INTEGER;
parm.dod->required = NO;
- parm.dod->answer = "0";
+ parm.dod->answer = "1";
parm.dod->description = _("Degree of over-determination");
parm.range = G_define_option();
@@ -226,11 +319,6 @@
parm.ts->multiple = YES;
parm.ts->description = _("Time steps of the input maps");
- parm.bl = G_define_option();
- parm.bl->key = "base_period";
- parm.bl->type = TYPE_INTEGER;
- parm.bl->description = _("Length of the base period");
-
parm.delta = G_define_option();
parm.delta->key = "delta";
parm.delta->type = TYPE_DOUBLE;
@@ -260,18 +348,12 @@
if (!parm.input->answer && !parm.file->answer)
G_fatal_error(_("Please specify input= or file="));
- if (parm.amp->answer && parm.phase->answer) {
- if (strcmp(parm.amp->answer, parm.phase->answer) == 0)
- G_fatal_error(_("'%s' and '%s' options must be different"),
- parm.amp->key, parm.phase->key);
- }
- do_amp = parm.amp->answer != NULL;
- do_phase = parm.phase->answer != NULL;
+ order = atoi(parm.order->answer);
+ if (order < 1)
+ G_fatal_error(_("The order number must be > 0"));
+ if (order > 3)
+ G_fatal_error(_("The order number must be <= 3"));
- nf = atoi(parm.nf->answer);
- if (nf < 1)
- G_fatal_error(_("The number of frequencies must be > 0"));
-
fet = DBL_MAX;
if (parm.fet->answer) {
fet = atof(parm.fet->answer);
@@ -302,6 +384,23 @@
if ((rejlo || rejhi) && !parm.fet->answer)
G_fatal_error(_("Fit error tolerance is required when outliers should be rejected"));
+ weight_func = tricube;
+ if (strcmp(parm.weight->answer, "uniform") == 0)
+ weight_func = uniform;
+ else if (strcmp(parm.weight->answer, "triangular") == 0)
+ weight_func = triangular;
+ else if (strcmp(parm.weight->answer, "epanechnikov") == 0)
+ weight_func = epanechnikov;
+ else if (strcmp(parm.weight->answer, "quartic") == 0)
+ weight_func = quartic;
+ else if (strcmp(parm.weight->answer, "tricube") == 0)
+ weight_func = tricube;
+ else if (strcmp(parm.weight->answer, "cosine") == 0)
+ weight_func = cosine;
+ else
+ G_fatal_error(_("Unknown weighing function '%s'"),
+ parm.weight->answer);
+
/* process the input maps from the file */
if (parm.file->answer) {
FILE *in;
@@ -369,12 +468,24 @@
if (num_inputs < 3)
G_fatal_error(_("At least 3 input maps are required"));
- /* length of base period */
- if (parm.bl->answer)
- bl = atoi(parm.bl->answer);
- else
- bl = num_inputs;
+ ts = G_malloc(num_inputs * sizeof(double));
+ if (parm.ts->answer) {
+ for (i = 0; parm.ts->answers[i]; i++);
+ if (i != num_inputs)
+ G_fatal_error(_("Number of time steps does not match number of input maps"));
+ for (i = 0; parm.ts->answers[i]; i++) {
+ ts[i] = atof(parm.ts->answers[i]);
+ if (i > 0 && ts[i - 1] >= ts[i])
+ G_fatal_error(_("Time steps must increase"));
+ }
+ }
+ else {
+ for (i = 0; i < num_inputs; i++) {
+ ts[i] = i;
+ }
+ }
+
/* open output maps */
num_outputs = num_inputs;
@@ -382,7 +493,7 @@
suffix = parm.suffix->answer;
if (!suffix)
- suffix = "_hants";
+ suffix = "_lwr";
for (i = 0; i < num_outputs; i++) {
struct output *out = &outputs[i];
@@ -395,128 +506,35 @@
out->fd = Rast_open_new(output_name, DCELL_TYPE);
}
- nr = 2 * nf + 1;
+ min_points = 1 + order + dod;
- if (nr > num_inputs)
- G_fatal_error(_("The maximum number of frequencies for %d input maps is %d."),
- num_inputs, (num_inputs - 1) / 2);
+ if (min_points > num_inputs)
+ G_fatal_error(_("At least %d input maps are required for "
+ "order number %d and "
+ "degree of over-determination %d."),
+ min_points, order, dod);
- noutmax = num_inputs - nr - dod;
-
- if (noutmax < 0)
- G_fatal_error(_("The degree of overdetermination can not be larger than %d "
- "for %d input maps and %d frequencies."),
- dod + noutmax, num_inputs, nf);
-
- if (noutmax == 0)
- G_warning(_("Missing values can not be reconstructed, "
- "please reduce either '%s' or '%s'"),
- parm.nf->key, parm.dod->key);
-
- if (do_amp) {
- /* open output amplitude */
- out_amp = G_calloc(nf, sizeof(struct output));
-
- for (i = 0; i < nf; i++) {
- struct output *out = &out_amp[i];
- char output_name[GNAME_MAX];
-
- sprintf(output_name, "%s.%d", parm.amp->answer, i);
-
- out->name = G_store(output_name);
- out->buf = Rast_allocate_d_buf();
- out->fd = Rast_open_new(output_name, DCELL_TYPE);
- }
- }
-
- if (do_phase) {
- /* open output phase */
- out_phase = G_calloc(nf, sizeof(struct output));
-
- for (i = 0; i < nf; i++) {
- struct output *out = &out_phase[i];
- char output_name[GNAME_MAX];
-
- sprintf(output_name, "%s.%d", parm.phase->answer, i);
-
- out->name = G_store(output_name);
- out->buf = Rast_allocate_d_buf();
- out->fd = Rast_open_new(output_name, DCELL_TYPE);
- }
- }
-
/* initialise variables */
values = G_malloc(num_inputs * sizeof(DCELL));
+ values2 = G_malloc(num_inputs * sizeof(DCELL));
+ resultn = G_malloc(num_inputs * sizeof(DCELL));
nrows = Rast_window_rows();
ncols = Rast_window_cols();
- cs = G_alloc_vector(bl);
- sn = G_alloc_vector(bl);
- ts = G_alloc_vector(num_inputs);
- rc = G_malloc(num_inputs * sizeof(DCELL));
- useval = G_malloc(num_inputs * sizeof(int));
+ isnull = G_malloc(num_inputs * sizeof(int));
- if (parm.ts->answer) {
- for (i = 0; parm.ts->answers[i]; i++);
- if (i != num_inputs)
- G_fatal_error(_("Number of time steps does not match number of input maps"));
- for (i = 0; parm.ts->answers[i]; i++)
- ts[i] = atof(parm.ts->answers[i]);
- }
- else {
- for (i = 0; i < num_inputs; i++) {
- ts[i] = i;
- }
- }
+ msize = 1 + order;
+ m = G_alloc_matrix(msize, msize);
+ m2 = G_alloc_matrix(msize, msize);
+ a = G_alloc_vector(msize);
+ a2 = G_alloc_vector(msize);
+ B = G_alloc_vector(msize);
- mat = G_alloc_matrix(nr, num_inputs);
- mat_t = G_alloc_matrix(num_inputs, nr);
- A = G_alloc_matrix(nr, nr);
- Av = *A;
- Azero = G_alloc_vector(nr * nr);
- asize = nr * nr * sizeof(double);
- za = G_alloc_vector(nr);
- zr = G_alloc_vector(nr);
- rc = G_alloc_vector(num_inputs);
+ /* calculate weights for different margins */
- for (i = 0; i < nr * nr; i++)
- Azero[i] = 0;
-
- for (i = 0; i < bl; i++) {
- double ang = 2.0 * M_PI * i / bl;
-
- cs[i] = cos(ang);
- sn[i] = sin(ang);
- }
- for (j = 0; j < num_inputs; j++) {
- mat[0][j] = 1.;
- mat_t[j][0] = 1.;
- }
-
- for (i = 0; i < nr / 2; i++) {
- int i1 = 2 * i + 1;
- int i2 = 2 * i + 2;
-
- for (j = 0; j < num_inputs; j++) {
-
- int idx = (int) ((i + 1) * (ts[j])) % bl;
-
- if (idx >= bl)
- G_fatal_error("cs/sn index out of range: %d, %d", idx, bl);
-
- if (i2 >= nr)
- G_fatal_error("mat index out of range: %d, %d", 2 * i + 2, nr);
-
- mat[i1][j] = cs[idx];
- mat[i2][j] = sn[idx];
- mat_t[j][i1] = mat[i1][j];
- mat_t[j][i2] = mat[i2][j];
- }
- }
-
/* process the data */
- G_message(_("Harmonic analysis of %d input maps..."), num_inputs);
+ G_message(_("Local weighted regression of %d input maps..."), num_inputs);
for (row = 0; row < nrows; row++) {
G_percent(row, nrows, 4);
@@ -535,146 +553,205 @@
}
for (col = 0; col < ncols; col++) {
- int null = 0, non_null = 0;
+ n_nulls = 0;
for (i = 0; i < num_inputs; i++) {
DCELL v = inputs[i].buf[col];
- useval[i] = 0;
+ isnull[i] = 0;
if (Rast_is_d_null_value(&v)) {
- null++;
+ isnull[i] = 1;
+ n_nulls++;
}
else if (parm.range->answer && (v < lo || v > hi)) {
Rast_set_d_null_value(&v, 1);
- null++;
+ isnull[i] = 1;
+ n_nulls++;
}
- else {
- non_null++;
- useval[i] = 1;
- }
-
values[i] = v;
}
- nout = null;
- /* HANTS */
- if (nout <= noutmax) {
- int n = 0, done = 0;
+ /* LWR */
+ if (num_inputs - n_nulls >= min_points) {
- while (!done) {
+ for (i = 0; i < num_inputs; i++) {
+ DCELL result;
+
+ /* margin around i */
+ n_points = 0;
+ in_lo = in_hi = i;
+ this_margin = 0;
+ if (!isnull[i])
+ n_points++;
+ for (j = 1; j < num_inputs; j++) {
+ if (i - j >= 0) {
+ if (!isnull[i - j]) {
+ n_points++;
+ in_lo = i - j;
+ }
+ }
+ if (i + j < num_inputs) {
+ if (!isnull[i + j]) {
+ n_points++;
+ in_hi = i + j;
+ }
+ }
+ if (n_points >= min_points) {
+ this_margin = j;
+ break;
+ }
+ }
+
+ tsdiff1 = ts[i] - ts[in_lo];
+ tsdiff2 = ts[in_hi] - ts[i];
+
+ max_ts = tsdiff1;
+ if (max_ts < tsdiff2)
+ max_ts = tsdiff2;
- /* za = mat * y */
+ max_ts *= (1.0 + 1.0 / this_margin);
- /* A = mat * diag(p) * mat' */
+ /* initialize matrix and vectors */
+ for (j = 0; j <= order; j++) {
+ a[j] = 0;
+ B[j] = 0;
+ m[j][j] = 0;
+ for (k = 0; k < j; k++) {
+ m[j][k] = m[k][j] = 0;
+ }
+ }
+
+ /* load points */
+ for (n = in_lo; n <= in_hi; n++) {
+ if (isnull[n])
+ continue;
- /* mat: nr, num_inputs
- * diag(p): num_inputs, num_inputs
- * mat_t: num_inputs, nr
- * A temp: nr, num_inputs
- * A: nr, nr */
+ weight = weight_func(ts[i], ts[n], max_ts);
+ for (j = 0; j <= order; j++) {
+ double val1 = term(j, ts[n]);
- memcpy(Av, Azero, asize);
- for (i = 0; i < nr; i++) {
- za[i] = 0;
- for (j = 0; j < num_inputs; j++) {
- if (useval[j]) {
- za[i] += mat[i][j] * values[j];
+ for (k = j; k <= order; k++) {
+ double val2 = term(k, ts[n]);
- for (k = 0; k < nr; k++)
- A[i][k] += mat[i][j] * mat_t[j][k];
+ m[j][k] += val1 * val2 * weight;
}
+ a[j] += values[n] * val1 * weight;
}
+ }
- if (i > 0) {
- A[i][i] += delta;
+ /* TRANSPOSE VALUES IN UPPER HALF OF M TO OTHER HALF */
+ m2[0][0] = m[0][0];
+ a2[0] = a[0];
+ for (j = 1; j <= order; j++) {
+ for (k = 0; k < j; k++) {
+ m[j][k] = m[k][j];
+ m2[j][k] = m2[k][j] = m[k][j];
}
+ m[j][j] *= (1 + delta);
+ m2[j][j] = m[j][j];
+ a2[j] = a[j];
}
-
- /* zr = A \ za
- * solve A * zr = za */
- if (!solvemat(A, za, zr, nr)) {
- done = -1;
- Rast_set_d_null_value(rc, num_outputs);
- break;
- }
- /* G_math_solver_gauss(A, zr, za, nr) is much slower */
- /* rc = mat' * zr */
- maxerrlo = maxerrhi = 0;
- for (i = 0; i < num_inputs; i++) {
- rc[i] = 0;
- for (j = 0; j < nr; j++) {
- rc[i] += mat_t[i][j] * zr[j];
+ if (solvemat(m2, a2, B, order + 1) != 0) {
+ /* get estimate */
+ result = 0.0;
+ for (j = 0; j <= order; j++) {
+ result += B[j] * term(j, ts[i]);
}
- if (useval[i]) {
- if (maxerrlo < rc[i] - values[i])
- maxerrlo = rc[i] - values[i];
- if (maxerrhi < values[i] - rc[i])
- maxerrhi = values[i] - rc[i];
- }
- }
- if (rejlo || rejhi) {
- done = 1;
- if (rejlo && maxerrlo > fet)
- done = 0;
- if (rejhi && maxerrhi > fet)
- done = 0;
-
- if (!done) {
- /* filter outliers */
- for (i = 0; i < num_inputs; i++) {
+
+ if (rejlo || rejhi) {
+ int done = 0;
- if (useval[i]) {
- if (rejlo && rc[i] - values[i] > maxerrlo * 0.5) {
- useval[i] = 0;
- nout++;
+ for (n = in_lo; n <= in_hi; n++) {
+ if (isnull[n])
+ continue;
+
+ values2[n] = values[n];
+ }
+
+ while (!done) {
+ done = 1;
+
+ maxerrlo = maxerrhi = 0;
+ for (n = in_lo; n <= in_hi; n++) {
+ if (isnull[n])
+ continue;
+
+ resultn[n] = 0.0;
+ for (j = 0; j <= order; j++) {
+ resultn[n] += B[j] * term(j, ts[n]);
}
- if (rejhi && values[i] - rc[i] > maxerrhi * 0.5) {
- useval[i] = 0;
- nout++;
- }
+ if (maxerrlo < resultn[n] - values2[n])
+ maxerrlo = resultn[n] - values2[n];
+ if (maxerrhi < values2[n] - resultn[n])
+ maxerrhi = values2[n] - resultn[n];
}
- }
- }
- }
- n++;
- if (n >= num_inputs)
- done = 1;
- if (nout > noutmax)
- done = 1;
- }
+ if (rejlo && maxerrlo > fet)
+ done = 0;
+ if (rejhi && maxerrhi > fet)
+ done = 0;
- for (i = 0; i < num_outputs; i++) {
- struct output *out = &outputs[i];
+ if (!done) {
- out->buf[col] = rc[i];
- if (rc[i] < lo)
- out->buf[col] = lo;
- else if (rc[i] > hi)
- out->buf[col] = hi;
- }
+ a2[0] = 0;
+ m2[0][0] = m[0][0];
+ for (j = 1; j <= order; j++) {
+ for (k = 0; k < j; k++) {
+ m2[j][k] = m2[k][j] = m[k][j];
+ }
+ m2[j][j] = m[j][j];
+ a2[j] = 0;
+ }
- if (do_amp || do_phase) {
- /* amplitude and phase */
- /* skip constant */
+ /* replace outliers */
+ for (n = in_lo; n <= in_hi; n++) {
+ if (isnull[n])
+ continue;
- for (i = 1; i < nr; i += 2) {
- int ifr = i >> 1;
+ weight = weight_func(ts[i], ts[n], max_ts);
+ if (rejlo && resultn[n] - values2[n] > maxerrlo * 0.5) {
+ values2[n] = (resultn[n] + values2[n]) * 0.5;
+ }
+ if (rejhi && values2[n] - resultn[n] > maxerrhi * 0.5) {
+ values2[n] = (values2[n] + resultn[n]) * 0.5;
+ }
+ for (j = 0; j <= order; j++) {
+ double val1 = term(j, ts[n]);
- if (do_amp) {
- out_amp[ifr].buf[col] = sqrt(zr[i] * zr[i] +
- zr[i + 1] * zr[i + 1]);
+ a2[j] += values2[n] * val1 * weight;
+ }
+ }
+ done = 1;
+ if (solvemat(m2, a2, B, order + 1) != 0) {
+ /* update estimate */
+ result = 0.0;
+ for (j = 0; j <= order; j++) {
+ result += B[j] * term(j, ts[i]);
+ }
+ done = 0;
+ }
+ }
+ }
}
+ }
+ else {
+ double wsum = 0.0;
- if (do_phase) {
- double angle = atan2(zr[i + 1], zr[i]) * 180 / M_PI;
+ G_warning(_("Points are (nearly) co-linear, using weighted average"));
- if (angle < 0)
- angle += 360;
- out_phase[ifr].buf[col] = angle;
+ result = 0.0;
+ for (n = in_lo; n <= in_hi; n++) {
+ if (isnull[n])
+ continue;
+
+ weight = weight_func(ts[i], ts[n], max_ts);
+ result += values[n] * weight;
+ wsum += weight;
}
+ result /= wsum;
}
+ outputs[i].buf[col] = result;
}
}
else {
@@ -683,74 +760,32 @@
Rast_set_d_null_value(&out->buf[col], 1);
}
- if (do_amp || do_phase) {
- for (i = 0; i < nf; i++) {
- if (do_amp)
- Rast_set_d_null_value(&out_amp[i].buf[col], 1);
- if (do_phase)
- Rast_set_d_null_value(&out_phase[i].buf[col], 1);
- }
- }
}
}
for (i = 0; i < num_outputs; i++)
Rast_put_d_row(outputs[i].fd, outputs[i].buf);
-
- if (do_amp || do_phase) {
- for (i = 0; i < nf; i++) {
- if (do_amp)
- Rast_put_d_row(out_amp[i].fd, out_amp[i].buf);
- if (do_phase)
- Rast_put_d_row(out_phase[i].fd, out_phase[i].buf);
- }
- }
}
G_percent(row, nrows, 2);
- /* Close input maps */
- if (!flag.lazy->answer) {
- for (i = 0; i < num_inputs; i++)
- Rast_close(inputs[i].fd);
- }
-
- /* close output maps */
+ /* close input and output maps */
for (i = 0; i < num_outputs; i++) {
+ struct input *in = &inputs[i];
struct output *out = &outputs[i];
+ if (!flag.lazy->answer)
+ Rast_close(in->fd);
+
Rast_close(out->fd);
Rast_short_history(out->name, "raster", &history);
Rast_command_history(&history);
Rast_write_history(out->name, &history);
- }
- if (do_amp) {
- /* close output amplitudes */
- for (i = 0; i < nf; i++) {
- struct output *out = &out_amp[i];
-
- Rast_close(out->fd);
-
- Rast_short_history(out->name, "raster", &history);
- Rast_command_history(&history);
- Rast_write_history(out->name, &history);
- }
+ if (Rast_read_colors(in->name, "", &colors) == 1)
+ Rast_write_colors(out->name, G_mapset(), &colors);
}
- if (do_phase) {
- /* close output phases */
- for (i = 0; i < nf; i++) {
- struct output *out = &out_phase[i];
-
- Rast_close(out->fd);
-
- Rast_short_history(out->name, "raster", &history);
- Rast_command_history(&history);
- Rast_write_history(out->name, &history);
- }
- }
-
exit(EXIT_SUCCESS);
}
Deleted: grass-addons/grass7/raster/r.series.lwr/r.hants.html
===================================================================
--- grass-addons/grass7/raster/r.series.lwr/r.hants.html 2016-11-02 21:47:26 UTC (rev 69762)
+++ grass-addons/grass7/raster/r.series.lwr/r.hants.html 2016-11-02 21:50:43 UTC (rev 69763)
@@ -1,204 +0,0 @@
-<h2>DESCRIPTION</h2>
-
-<em>r.hants</em> performs a harmonic analysis of time series in order to
-estimate missing values and identify outliers. For each input map, an
-output map with the suffix <em>suffix</em> (default: _hants) is created.
-
-<p>
-The option <b>nf</b>, number of frequencies, should be carefully chosen.
-Different numbers of frequencies should be tested first on a small test
-region before running the module on the full region. As a rule of thumb,
-the number of frequencies should be at least
-<em>estimated periodicity + 3</em>, e.g. for NDVI with an annual cycle
-(one peak per year), the number of frequencies should be at least 4 when
-analysing one year. If two peaks are assumed per year, the number of
-frequencies should be at least 5 when analysing one year.
-
-<p>The number of frequencies should also not be too large. With a large
-number of frequencies, outliers can no longer be identified because the
-fit is "too good", outliers can be represented by the estimate.
-The number of frequencies should also be smaller than <em>n input maps / 2</em>
-if missing values should be reconstructed.
-
-
-<h2>NOTES</h2>
-
-The optional <em>amplitude</em> and <em>phase</em> output maps contain
-the amplitude and phase for each frequency. The amplitude maps can be
-used to identify the dominant frequency with <em>r.series
-method=max_raster</em>. The baseline frequeny (base period) has the
-suffix <em>.0</em>, its first harmonic has the suffix <em>.1</em>, its
-second harmonic has the suffix <em>.2</em>, etc. The value of the
-output of <em>r.series method=max_raster</em> is identical to the
-number of the suffix. With the <em>amplitude</em> output maps for NDVI
-input, this can be used to determine the number of peaks in vegetation
-growth within the base period, where 0 (zero) means that the dominant
-frequency is the base period, i.e. one peak per base period.
-
-<p>
-If the <em>range</em> option is given, any values which fall outside
-that range will be treated as if they were NULL.
-The <em>range</em> parameter can be set to <em>low,high</em> thresholds:
-values outside of this range are treated as NULL.
-
-<p>
-The <em>low,high</em> thresholds are floating point, so use <em>-inf</em> or
-<em>inf</em> for a single threshold (e.g., <em>range=0,inf</em> to ignore
-negative values, or <em>range=-inf,-200.4</em> to ignore values above -200.4).
-
-<p>
-The maximum number of raster maps that can be processed is given by the
-per-user limit of the operating system. For example, the soft limits
-for users are typically 1024. The soft limit can be changed with e.g.
-<tt>ulimit -n 4096</tt> (UNIX-based operating systems) but not higher
-than the hard limit. If it is too low, you can as superuser add an
-entry in
-
-<div class="code"><pre>
-/etc/security/limits.conf
-# <domain> <type> <item> <value>
-your_username hard nofile 4096
-</pre></div>
-
-This would raise the hard limit to 4096 files. Also have a look at the
-overall limit of the operating system
-<div class="code"><pre>
-cat /proc/sys/fs/file-max
-</pre></div>
-which is on modern Linux systems several 100,000 files.
-
-
-<p>Use the <em>file</em> option to analyze large amount of raster maps
-without hitting open files limit and the size limit of command line
-arguments. The computation is slower than the <em>input</em> option
-method. For every single row in the output map(s) all input maps are
-opened and closed. The amount of RAM will rise linear with the number of
-specified input maps. The input and file options are mutually exclusive.
-Input is a text file with a new line separated list of raster map names.
-
-<p>For the usage of the option <b>fet</b> see Roerink et al. (2000). Its
-value is relative to the value range of the variable being considered.
-
-<h2>EXAMPLES</h2>
-
-<h3>Precipitation data example</h3>
-This small example is based on a climatic dataset for North Carolina which
-was from publicly available data (monthly temperature averages and monthly
-precipitation sums from 2000 to 2012, downloadable as
-<a href="http://courses.ncsu.edu/mea592/common/media/02/nc_climate_spm_2000_2012.zip">GRASS GIS 7 location</a>):
-
-
-<div class="code"><pre>
-# set computational region to one of the maps
-g.region raster=2004_03_tempmean -p
-</pre></div>
-
-Visualize the time series as animation:
-<div class="code"><pre>
-# note: color table is different from standard "celsius"
-g.gui.animation rast=`g.list type=raster pattern="*tempmean" sep=comma`
-</pre></div>
-
-<p>
-Since HANTS is CPU intensive, we test for now at lower resolution:
-<div class="code"><pre>
-g.region -p res=5000
-
-# HANTS: Harmonic analysis of the 156 input maps...
-# just wildly guessing the parameters for a test run:
-
-# generate and check list of input maps (the order matters!)
-g.list type=raster pattern="20??_??_tempmean" output=tempmean.csv
-
-r.hants file=tempmean.csv \
- nf=6 fet=0.1 dod=5 delta=0.1 base_period=12
-
-# assign reasonable color tables for temperature
-for map in `g.list type=raster pattern="*tempmean_hants"` ; do
- r.colors $map color=celsius
-done
-
-# assign degree Celsius color table
-r.colors 2000_06_tempmean_hants color=celsius
-
-# verify with one of the 156 results (still at reduced resolution):
-r.mapcalc "2000_06_tempmean_diff = 2000_06_tempmean - 2000_06_tempmean_hants"
-
-r.colors 2000_06_tempmean_diff color=differences
-d.mon wx0
-d.rast 2000_06_tempmean_hants
-d.rast 2000_06_tempmean_diff
-
-
-r.univar 2000_06_tempmean_diff -g
-n=5066
-null_cells=5040
-cells=10106
-min=-0.0899336115228095
-max=0.359362050140941
-range=0.449295661663751
-mean=0.188579838052468
-...
-
-# see HANTS time series as animation
-g.gui.animation rast=`g.list type=raster pattern="*tempmean_hants" sep=comma`
-
-
-# Check HANTS behaviour in a given point
-east=740830
-north=168832
-
-for map in `g.list rast pat="20??_??_tempmean"` ; do
- r.what map=$map coordinates=$east,$north >> time_series_orig.csv
-done
-
-for map in `g.list rast pat="*tempmean_hants"` ; do
- r.what map=$map coordinates=$east,$north >> time_series_hants.csv
-done
-
-# merge files:
-echo "east|north|temp_orig|temp_hants" > time_series_final.csv
-
-paste -d'|' time_series_orig.csv time_series_hants.csv | \
- cut -d'|' -f1,2,4,8 >> time_series_final.csv
-
-# Resulting CSV file: 'time_series_final.csv'
-</pre></div>
-
-<h3>Using t.* modules to check for basic statistics</h3>
-
-The temporal framework (t.* modules) can be used to assess basic statistics:
-
-<div class="code"><pre>
-# create spatio temporal data set with hants output maps
-t.create type=strds temporaltype=absolute output=tempmean_hants \
- title="Mean Temperature HANTS" description="Mean Temperature reconstructed with HANTS"
-
-# register maps in the strds
-t.register -i type=raster input=tempmean_hants \
- maps=`g.list raster pattern=*tempmean_hants sep=,` start="2000-01-01" \
- increment="1 months"
-
-# getting general info of the strds (including max and min of the whole series)
-t.info type=strds input=tempmean_hants
-
-# getting statistics for each map in the series
-t.rast.univar -h tempmean_hants > stats_hants.txt
-</pre></div>
-
-<h2>SEE ALSO</h2>
-
-<em><a href="r.series.html">r.series</a></em>
-
-<h2>REFERENCES</h2>
-<p>
-Roerink, G. J., Menenti, M. and Verhoef, W., 2000. Reconstructing
-cloudfree NDVI composites using Fourier analysis of time series.
-International Journal of Remote Sensing, 21 (9), 1911-1917.
-DOI: <a href=http://dx.doi.org/10.1080/014311600209814>10.1080/014311600209814</a>
-
-<h2>AUTHOR</h2>
-
-Markus Metz
-
-<p><i>Last changed: $Date$</i>
Copied: grass-addons/grass7/raster/r.series.lwr/r.series.lwr.html (from rev 69762, grass-addons/grass7/raster/r.series.lwr/r.hants.html)
===================================================================
--- grass-addons/grass7/raster/r.series.lwr/r.series.lwr.html (rev 0)
+++ grass-addons/grass7/raster/r.series.lwr/r.series.lwr.html 2016-11-02 21:50:43 UTC (rev 69763)
@@ -0,0 +1,88 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.series.lwr</em> performs a local weighted regression (LWR) of
+time series in order to estimate missing values and identify outliers.
+For each input map, an output map with the suffix <em>suffix</em>
+(default: _lwr) is created.
+
+<p>
+For each point in the time series, the (in time) neighboring points are
+considered to estimate a polynomial function that best fits the
+observations. The points are weighted according to their distance (in
+time) to the current point with points farther away getting a lower
+weight.
+
+<p>
+The option <b>order</b> determines the order of the polynomial function
+used to fit the observations. An order of 1 is a linear regression,
+recommended and default is <em>order=2</em>. An order of 3 might cause
+over-fitting, resulting in no smoothing, and should only be used with
+data that show high fluctuations.
+
+<p>
+The <em>delta</em> option is experimental. It smoothes the polynomial
+function: with a high <em>delta</em> value, a higher order polynomial
+function fitting becomes similar to a moving average.
+
+<p>
+Optionally, outliers can be removed with the <em>-l</em> and
+<em>-h</em> flags when fitting a polynomial. In this case, a fit error
+tolerance (option <em>fit</em>) must be provided. These flags are
+slowing down the module and are only needed when the time series
+contains many outliers. A few outliers are handled well with the
+default settings.
+
+<h2>NOTES</h2>
+
+If the <em>range</em> option is given, any values which fall outside
+that range will be treated as if they were NULL.
+The <em>range</em> parameter can be set to <em>low,high</em> thresholds:
+values outside of this range are treated as NULL.
+
+<p>
+The <em>low,high</em> thresholds are floating point, so use
+<em>-inf</em> or <em>inf</em> for a single threshold (e.g.,
+<em>range=0,inf</em> to ignore negative values, or
+<em>range=-inf,-200.4</em> to ignore values above -200.4).
+
+<p>
+The maximum number of raster maps that can be processed is given by the
+per-user limit of the operating system. For example, the soft limits
+for users are typically 1024. The soft limit can be changed with e.g.
+<tt>ulimit -n 4096</tt> (UNIX-based operating systems) but not higher
+than the hard limit. If it is too low, you can as superuser add an
+entry in
+
+<div class="code"><pre>
+/etc/security/limits.conf
+# <domain> <type> <item> <value>
+your_username hard nofile 4096
+</pre></div>
+
+This would raise the hard limit to 4096 files. Also have a look at the
+overall limit of the operating system
+<div class="code"><pre>
+cat /proc/sys/fs/file-max
+</pre></div>
+which is on modern Linux systems several 100,000 files.
+
+<p>Use the <em>-z</em> flag to analyze large amount of raster maps
+without hitting open files limit and the size limit of command line
+arguments. This will however increase the processing time. For every
+single row in the output map(s) all input maps are opened and closed.
+The amount of RAM will rise linear with the number of specified input
+maps.
+<p>
+The input and file options are mutually exclusive. Input is a
+text file with a new line separated list of raster map names.
+
+<h2>SEE ALSO</h2>
+
+<em><a href="r.series.html">r.series</a></em>,
+<em><a href="r.hants.html">r.hants</a></em>
+
+<h2>AUTHOR</h2>
+
+Markus Metz
+
+<p><i>Last changed: $Date$</i>
More information about the grass-commit
mailing list