[GRASS-SVN] r61005 - grass-addons/grass7/raster/r.gwr
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jun 26 23:22:32 PDT 2014
Author: mmetz
Date: 2014-06-26 23:22:32 -0700 (Thu, 26 Jun 2014)
New Revision: 61005
Added:
grass-addons/grass7/raster/r.gwr/Makefile
grass-addons/grass7/raster/r.gwr/bufs.c
grass-addons/grass7/raster/r.gwr/bufs.h
grass-addons/grass7/raster/r.gwr/estimate_bw.c
grass-addons/grass7/raster/r.gwr/gwr.c
grass-addons/grass7/raster/r.gwr/local_proto.h
grass-addons/grass7/raster/r.gwr/main.c
grass-addons/grass7/raster/r.gwr/r.gwr.html
grass-addons/grass7/raster/r.gwr/weights.c
Modified:
grass-addons/grass7/raster/r.gwr/
Log:
r.gwr: geographically weighted regression
Property changes on: grass-addons/grass7/raster/r.gwr
___________________________________________________________________
Added: svn:ignore
+ OBJ.* *.tmp.html
Added: grass-addons/grass7/raster/r.gwr/Makefile
===================================================================
--- grass-addons/grass7/raster/r.gwr/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.gwr/Makefile 2014-06-27 06:22:32 UTC (rev 61005)
@@ -0,0 +1,13 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.gwr
+
+LIBES = $(RASTERLIB) $(GISLIB) $(GMATHLIB) $(MATHLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+EXTRA_LIBS = $(OMPLIB)
+EXTRA_CFLAGS = $(OMPCFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass-addons/grass7/raster/r.gwr/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.gwr/bufs.c
===================================================================
--- grass-addons/grass7/raster/r.gwr/bufs.c (rev 0)
+++ grass-addons/grass7/raster/r.gwr/bufs.c 2014-06-27 06:22:32 UTC (rev 61005)
@@ -0,0 +1,83 @@
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include "bufs.h"
+
+/*
+ allocate the i/o bufs
+
+ the i/o bufs will be rotated by the read operation so that the
+ last row read will be in the last i/o buf
+
+ */
+
+int allocate_bufs(struct rb *rbuf, int ncols, int bw, int fd)
+{
+ int i;
+ int ncolsbw;
+ size_t bufsize;
+
+ ncolsbw = ncols + 2 * bw;
+ bufsize = ncolsbw * sizeof(DCELL);
+
+ rbuf->bw = bw;
+ rbuf->nsize = bw * 2 + 1;
+ rbuf->fd = fd;
+ rbuf->row = 0;
+
+ rbuf->buf = (DCELL **) G_malloc(rbuf->nsize * sizeof(DCELL *));
+
+ for (i = 0; i < rbuf->nsize; i++) {
+ rbuf->buf[i] = (DCELL *) G_malloc(bufsize);
+ Rast_set_d_null_value(rbuf->buf[i], ncolsbw);
+ }
+
+ return 0;
+}
+
+int release_bufs(struct rb *rbuf)
+{
+ int i;
+
+ for (i = 0; i < rbuf->nsize; i++) {
+ G_free(rbuf->buf[i]);
+ }
+
+ G_free(rbuf->buf);
+
+ rbuf->bw = 0;
+ rbuf->nsize = 0;
+ rbuf->row = 0;
+
+ return 0;
+}
+
+static int rotate_bufs(struct rb *rbuf)
+{
+ DCELL *temp;
+ int i;
+
+ temp = rbuf->buf[0];
+
+ for (i = 1; i < rbuf->nsize; i++)
+ rbuf->buf[i - 1] = rbuf->buf[i];
+
+ rbuf->buf[rbuf->nsize - 1] = temp;
+
+ return 0;
+}
+
+int readrast(struct rb *rbuf, int nrows, int ncols)
+{
+ rotate_bufs(rbuf);
+
+ if (rbuf->row < nrows) {
+ Rast_get_d_row(rbuf->fd, rbuf->buf[rbuf->nsize - 1] + rbuf->bw, rbuf->row);
+ }
+ else {
+ Rast_set_d_null_value(rbuf->buf[rbuf->nsize - 1] + rbuf->bw, ncols);
+ }
+
+ rbuf->row++;
+
+ return 0;
+}
Property changes on: grass-addons/grass7/raster/r.gwr/bufs.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.gwr/bufs.h
===================================================================
--- grass-addons/grass7/raster/r.gwr/bufs.h (rev 0)
+++ grass-addons/grass7/raster/r.gwr/bufs.h 2014-06-27 06:22:32 UTC (rev 61005)
@@ -0,0 +1,12 @@
+struct rb
+{
+ int fd;
+ int bw; /* bandwidth */
+ int nsize; /* bw * 2 + 1 */
+ int row; /* next row to read */
+ DCELL **buf; /* for reading raster map */
+};
+
+int allocate_bufs(struct rb *rbuf, int ncols, int bw, int fd);
+int release_bufs(struct rb *rbuf);
+int readrast(struct rb *rbuf, int nrows, int ncols);
Property changes on: grass-addons/grass7/raster/r.gwr/bufs.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.gwr/estimate_bw.c
===================================================================
--- grass-addons/grass7/raster/r.gwr/estimate_bw.c (rev 0)
+++ grass-addons/grass7/raster/r.gwr/estimate_bw.c 2014-06-27 06:22:32 UTC (rev 61005)
@@ -0,0 +1,292 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <unistd.h>
+#include <math.h>
+#include <string.h>
+#include <sys/types.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/raster.h>
+#include "local_proto.h"
+
+#ifndef USE_RAND
+
+#ifndef HAVE_DRAND48
+#define lrand48() (((long) rand() ^ ((long) rand() << 16)) & 0x7FFFFFFF)
+#define srand48(sv) (srand((unsigned)(sv)))
+#endif
+
+
+long make_rand(void)
+{
+ return lrand48();
+}
+
+void init_rand(void)
+{
+ srand48((long)time((time_t *) 0));
+}
+
+#else
+
+static long labs(int n)
+{
+ return n < 0 ? (-n) : n;
+}
+
+long make_rand(void)
+{
+ return (labs(rand() + (rand() << 16)));
+}
+
+void init_rand(void)
+{
+ srand(getpid());
+}
+
+#endif
+
+/* estimate bandwidth
+ * start with a small bandwidth */
+
+int estimate_bandwidth(int *inx, int ninx, int iny, int nrows, int ncols,
+ DCELL *est, int bw)
+{
+ int i;
+ int r, c;
+ struct rb *xbuf1, *xbuf2, *xbuf3, ybuf1, ybuf2, ybuf3;
+ int prevbw, nextbw, newbw;
+ int n1, n2, n3;
+ double f11, f12, f2;
+ double ss1, ss2, ss3;
+ double **w1, **w2, **w3;
+ DCELL *ybuf, yval;
+ int nr, nc, nrt, nct, count;
+
+ ybuf = Rast_allocate_d_buf();
+
+ G_message(_("Estimating optimal bandwidth..."));
+
+ /* count non-NULL cells of the dependent variable */
+ nc = 0;
+ for (r = 0; r < nrows; r++) {
+
+ Rast_get_d_row(iny, ybuf, r);
+
+ for (c = 0; c < ncols; c++) {
+ yval = ybuf[c];
+
+ if (!Rast_is_d_null_value(&yval))
+ nc++;
+ }
+ }
+ if (nc == 0)
+ G_fatal_error(_("No non-NULL cells in input map"));
+
+ /* number of cells to use for bandwidth estimation */
+ nr = 10000;
+ if (nr > nc)
+ nr = nc;
+
+ if (bw < 2) {
+ G_warning(_("Initial bandwidth must be > 1"));
+ bw = 2;
+ }
+
+ prevbw = bw - 1;
+ nextbw = bw + 1;
+
+ init_rand();
+
+ xbuf1 = G_malloc(ninx * sizeof(struct rb));
+ xbuf2 = G_malloc(ninx * sizeof(struct rb));
+ xbuf3 = G_malloc(ninx * sizeof(struct rb));
+
+ while (1) {
+
+ G_message(_("Testing bandwidth %d"), bw);
+
+ for (i = 0; i < ninx; i++) {
+ allocate_bufs(&(xbuf1[i]), ncols, prevbw, inx[i]);
+ allocate_bufs(&(xbuf2[i]), ncols, bw, inx[i]);
+ allocate_bufs(&(xbuf3[i]), ncols, nextbw, inx[i]);
+ }
+ allocate_bufs(&ybuf1, ncols, prevbw, iny);
+ allocate_bufs(&ybuf2, ncols, bw, iny);
+ allocate_bufs(&ybuf3, ncols, nextbw, iny);
+
+ /* initialize the raster buffers with 'bw' rows */
+ for (r = 0; r < prevbw; r++) {
+ for (i = 0; i < ninx; i++)
+ readrast(&(xbuf1[i]), nrows, ncols);
+ readrast(&ybuf1, nrows, ncols);
+ }
+ for (r = 0; r < bw; r++) {
+ for (i = 0; i < ninx; i++)
+ readrast(&(xbuf2[i]), nrows, ncols);
+ readrast(&ybuf2, nrows, ncols);
+ }
+ for (r = 0; r < nextbw; r++) {
+ for (i = 0; i < ninx; i++)
+ readrast(&(xbuf3[i]), nrows, ncols);
+ readrast(&ybuf3, nrows, ncols);
+ }
+
+ w1 = w_fn(prevbw);
+ w2 = w_fn(bw);
+ w3 = w_fn(nextbw);
+ /* leave one out: the center cell */
+ w1[prevbw][prevbw] = 0.;
+ w2[bw][bw] = 0.;
+ w3[nextbw][nextbw] = 0.;
+
+ ss1 = ss2 = ss3 = 0;
+ n1 = n2 = n3 = 0;
+
+ nrt = nr;
+ nct = nc;
+ count = 0;
+
+ for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 4);
+
+ for (i = 0; i < ninx; i++) {
+ readrast(&(xbuf1[i]), nrows, ncols);
+ readrast(&(xbuf2[i]), nrows, ncols);
+ readrast(&(xbuf3[i]), nrows, ncols);
+ }
+ readrast(&ybuf1, nrows, ncols);
+ readrast(&ybuf2, nrows, ncols);
+ readrast(&ybuf3, nrows, ncols);
+
+ for (c = 0; c < ncols; c++) {
+ yval = ybuf2.buf[bw][c + bw];
+
+ if (Rast_is_d_null_value(&yval))
+ continue;
+
+ if (make_rand() % nct < nrt) {
+ nrt--;
+ count++;
+
+ if (gwr(xbuf1, ninx, &ybuf1, c, prevbw, w1, est, NULL)) {
+ ss1 += (est[0] - yval) * (est[0] - yval);
+ n1++;
+ }
+ if (gwr(xbuf2, ninx, &ybuf2, c, bw, w2, est, NULL)) {
+ ss2 += (est[0] - yval) * (est[0] - yval);
+ n2++;
+ }
+ if (gwr(xbuf3, ninx, &ybuf3, c, nextbw, w3, est, NULL)) {
+ ss3 += (est[0] - yval) * (est[0] - yval);
+ n3++;
+ }
+ }
+ nct--;
+ }
+ }
+ G_percent(nrows, nrows, 4);
+
+ for (i = 0; i < ninx; i++) {
+ release_bufs(&(xbuf1[i]));
+ release_bufs(&(xbuf2[i]));
+ release_bufs(&(xbuf3[i]));
+ }
+ release_bufs(&ybuf1);
+ release_bufs(&ybuf2);
+ release_bufs(&ybuf3);
+
+ G_debug(3, "count: %d", count);
+
+ if (n1 == 0 || n2 == 0 || n3 == 0)
+ G_fatal_error(_("Unable to calculate sum of squares"));
+
+ /* correcting factor to favour larger bandwiths: bw * 2 + 1 */
+ ss1 /= n1 * prevbw * 2 + 1;
+ ss2 /= n2 * bw * 2 + 1;
+ ss3 /= n3 * nextbw * 2 + 1;
+
+ G_debug(1, "ss1: %g", ss1);
+ G_debug(1, "ss2: %g", ss2);
+ G_debug(1, "ss3: %g", ss3);
+
+#if 0
+ /* activate for debugging */
+ printf("%d|%g\n", prevbw, ss1);
+ printf("%d|%g\n", bw, ss2);
+ printf("%d|%g\n", nextbw, ss3);
+#endif
+
+ /* first gradient */
+ f11 = ss2 - ss1;
+ /* second gradient */
+ f12 = ss3 - ss2;
+ /* gradient of gradient */
+ f2 = f12 - f11;
+
+ G_free(w1[0]);
+ G_free(w2[0]);
+ G_free(w3[0]);
+ G_free(w1);
+ G_free(w2);
+ G_free(w3);
+
+ if (ss2 < ss1 && ss2 < ss3)
+ break;
+
+ if (ss1 < ss3) {
+ /* decrease bandwidth */
+ if (f2 > 0) {
+ /* Newton's method to find the extremum */
+ newbw = bw - (f11 + f12) / (2. * f2) + 0.5;
+ }
+ else {
+ newbw = bw - bw / 2;
+ }
+ }
+ else if (ss1 > ss3) {
+ /* increase bandwidth */
+ if (f2 > 0) {
+ /* Newton's method to find the extremum */
+ newbw = bw - (f11 + f12) / (2. * f2) + 0.5;
+ }
+ else {
+ newbw = bw + bw / 2;
+ }
+ }
+ else {
+ /* ss1 == ss3 && ss2 >= ss1|ss3) */
+ break;
+ }
+
+ if (ss1 < ss2 && newbw > bw) {
+ /* moving bw to wrong direction: stop here? */
+ G_debug(0, "f2: %g", f2);
+ newbw = bw - 1;
+ }
+ if (ss3 < ss2 && newbw < bw) {
+ /* moving bw to wrong direction: stop here? */
+ G_debug(0, "f2: %g", f2);
+ newbw = bw + 1;
+ }
+
+ if (newbw == bw)
+ break;
+
+ if (newbw < 2) {
+ G_debug(1, "Bandwidth is getting too small: %d", newbw);
+ if (bw > 2)
+ bw--;
+ else
+ break;
+ }
+ else
+ bw = newbw;
+
+ prevbw = bw - 1;
+ nextbw = bw + 1;
+ }
+
+ return bw;
+}
Property changes on: grass-addons/grass7/raster/r.gwr/estimate_bw.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.gwr/gwr.c
===================================================================
--- grass-addons/grass7/raster/r.gwr/gwr.c (rev 0)
+++ grass-addons/grass7/raster/r.gwr/gwr.c 2014-06-27 06:22:32 UTC (rev 61005)
@@ -0,0 +1,278 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/raster.h>
+#include "local_proto.h"
+
+/* geographically weighted regression:
+ * estimate coefficients for given cell */
+
+struct MATRIX
+{
+ int n; /* SIZE OF THIS MATRIX (N x N) */
+ double **v;
+};
+
+#define M(m,row,col) (m)->v[(row)][(col)]
+
+static int solvemat(struct MATRIX *m, double a[], double B[])
+{
+ int i, j, i2, j2, imark;
+ double factor, temp, *tempp;
+ double pivot; /* ACTUAL VALUE OF THE LARGEST PIVOT CANDIDATE */
+
+ for (i = 0; i < m->n; i++) {
+ j = i;
+
+ /* find row with largest magnitude value for pivot value */
+
+ pivot = M(m, i, j);
+ imark = i;
+ for (i2 = i + 1; i2 < m->n; i2++) {
+ temp = fabs(M(m, i2, j));
+ if (temp > fabs(pivot)) {
+ pivot = M(m, i2, j);
+ imark = i2;
+ }
+ }
+
+ /* if the pivot is very small then the points are nearly co-linear */
+ /* co-linear points result in an undefined matrix, and nearly */
+ /* co-linear points results in a solution with rounding error */
+
+ if (pivot == 0.0) {
+ G_debug(1, "Matrix is unsolvable");
+ return 0;
+ }
+
+ /* if row with highest pivot is not the current row, switch them */
+
+ if (imark != i) {
+ /*
+ for (j2 = 0; j2 < m->n; j2++) {
+ temp = M(m, imark, j2);
+ M(m, imark, j2) = M(m, i, j2);
+ M(m, i, j2) = temp;
+ }
+ */
+ tempp = m->v[imark];
+ m->v[imark] = m->v[i];
+ m->v[i] = tempp;
+
+ temp = a[imark];
+ a[imark] = a[i];
+ a[i] = temp;
+ }
+
+ /* compute zeros above and below the pivot, and compute
+ values for the rest of the row as well */
+
+ for (i2 = 0; i2 < m->n; i2++) {
+ if (i2 != i) {
+ factor = M(m, i2, j) / pivot;
+ for (j2 = j; j2 < m->n; j2++)
+ M(m, i2, j2) -= factor * M(m, i, j2);
+ a[i2] -= factor * a[i];
+ }
+ }
+ }
+
+ /* SINCE ALL OTHER VALUES IN THE MATRIX ARE ZERO NOW, CALCULATE THE
+ COEFFICIENTS BY DIVIDING THE COLUMN VECTORS BY THE DIAGONAL VALUES. */
+
+ for (i = 0; i < m->n; i++) {
+ B[i] = a[i] / M(m, i, i);
+ }
+
+ return 1;
+}
+
+
+int gwr(struct rb *xbuf, int ninx, struct rb *ybuf, int cc,
+ int bw, double **w, DCELL *est, double **B0)
+{
+ int r, c;
+ int i, j, k;
+ int nsize;
+ static DCELL *xval = NULL;
+ DCELL yval;
+ int count, isnull;;
+ /* OLS */
+ static double **a = NULL;
+ static double **B = NULL;
+ static struct MATRIX *m_all = NULL;
+ struct MATRIX *m;
+
+ if (!xval) {
+ xval = G_malloc((ninx + 1) * sizeof(DCELL));
+ xval[0] = 1.;
+ }
+
+ if (!m_all) {
+ m_all = (struct MATRIX *)G_malloc((ninx + 1) * sizeof(struct MATRIX));
+ a = (double **)G_malloc((ninx + 1) * sizeof(double *));
+ B = (double **)G_malloc((ninx + 1) * sizeof(double *));
+
+ m = &(m_all[0]);
+ m->n = ninx + 1;
+ m->v = (double **)G_malloc(m->n * sizeof(double));
+ m->v[0] = (double *)G_malloc(m->n * m->n * sizeof(double));
+ for (i = 1; i < m->n; i++) {
+ m->v[i] = m->v[i - 1] + m->n;
+ }
+
+ a[0] = (double *)G_malloc(m->n * sizeof(double));
+ B[0] = (double *)G_malloc(m->n * sizeof(double));
+
+ for (k = 1; k <= ninx; k++) {
+ m = &(m_all[k]);
+ m->n = ninx;
+ m->v = (double **)G_malloc(m->n * sizeof(double*));
+ m->v[0] = (double *)G_malloc(m->n * m->n * sizeof(double));
+ for (i = 1; i < m->n; i++) {
+ m->v[i] = m->v[i - 1] + m->n;
+ }
+ a[k] = (double *)G_malloc(m->n * sizeof(double));
+ B[k] = (double *)G_malloc(m->n * sizeof(double));
+ }
+ }
+
+ for (k = 0; k <= ninx; k++) {
+ m = &(m_all[k]);
+ for (i = 0; i < m->n; i++) {
+ for (j = i; j < m->n; j++)
+ M(m, i, j) = 0.0;
+ a[k][i] = 0.0;
+ B[k][i] = 0.0;
+ }
+ }
+
+ Rast_set_d_null_value(est, ninx + 1);
+ if (B0)
+ *B0 = NULL;
+
+ nsize = bw * 2 + 1;
+
+ /* first pass: collect values */
+ count = 0;
+ for (r = 0; r < nsize; r++) {
+
+ for (c = 0; c < nsize; c++) {
+
+ isnull = 0;
+ for (i = 0; i < ninx; i++) {
+ xval[i + 1] = xbuf[i].buf[r][c + cc];
+ if (Rast_is_d_null_value(&(xval[i + 1]))) {
+ isnull = 1;
+ break;
+ }
+ }
+ if (isnull)
+ continue;
+
+ yval = ybuf->buf[r][c + cc];
+ if (Rast_is_d_null_value(&yval))
+ continue;
+
+ for (i = 0; i <= ninx; i++) {
+ double val1 = xval[i];
+
+ for (j = i; j <= ninx; j++) {
+ double val2 = xval[j];
+
+ m = &(m_all[0]);
+ M(m, i, j) += val1 * val2 * w[r][c];
+
+ /* linear model without predictor k */
+ for (k = 1; k <= ninx; k++) {
+ if (k != i && k != j) {
+ int i2 = k > i ? i : i - 1;
+ int j2 = k > j ? j : j - 1;
+
+ m = &(m_all[k]);
+ M(m, i2, j2) += val1 * val2 * w[r][c];
+ }
+ }
+ }
+
+ a[0][i] += yval * val1 * w[r][c];
+ for (k = 1; k <= ninx; k++) {
+ if (k != i) {
+ int i2 = k > i ? i : i - 1;
+
+ a[k][i2] += yval * val1 * w[r][c];
+ }
+ }
+ }
+
+ count++;
+ }
+ }
+
+ /* estimate coefficients */
+ if (count < ninx + 1) {
+ G_debug(4, "Not enough valid cells available");
+ return 0;
+ }
+
+ for (k = 0; k <= ninx; k++) {
+ m = &(m_all[k]);
+
+ /* TRANSPOSE VALUES IN UPPER HALF OF M TO OTHER HALF */
+ for (i = 1; i < m->n; i++)
+ for (j = 0; j < i; j++)
+ M(m, i, j) = M(m, j, i);
+
+ if (!solvemat(m, a[k], B[k])) {
+ /*
+ for (i = 0; i <= ninx; i++) {
+ fprintf(stdout, "b%d=0.0\n", i);
+ }
+ */
+ G_debug(0, "Try a larger bandwidth");
+ count = 0;
+ }
+ }
+ if (count < ninx + 1) {
+ return 0;
+ }
+
+ /* second pass: calculate estimates */
+ c = cc;
+
+ isnull = 0;
+ for (i = 0; i < ninx; i++) {
+
+ xval[i + 1] = xbuf[i].buf[bw][cc + bw];
+ if (Rast_is_d_null_value(&(xval[i + 1]))) {
+ isnull = 1;
+ break;
+ }
+ }
+ if (isnull)
+ return 0;
+
+ est[0] = 0.0;
+ for (k = 0; k <= ninx; k++) {
+ est[0] += B[0][k] * xval[k];
+
+ if (k > 0) {
+ est[k] = 0.0;
+
+ /* linear model without predictor k */
+ for (i = 0; i <= ninx; i++) {
+ if (i != k) {
+ j = k > i ? i : i - 1;
+ est[k] += B[k][j] * xval[i];
+ }
+ }
+ }
+ }
+ if (B0)
+ *B0 = B[0];
+
+ return count;
+}
Property changes on: grass-addons/grass7/raster/r.gwr/gwr.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.gwr/local_proto.h
===================================================================
--- grass-addons/grass7/raster/r.gwr/local_proto.h (rev 0)
+++ grass-addons/grass7/raster/r.gwr/local_proto.h 2014-06-27 06:22:32 UTC (rev 61005)
@@ -0,0 +1,16 @@
+#include "bufs.h"
+
+double **gauss(int);
+double **epanechnikov(int);
+double **quartic(int);
+double **tricubic(int);
+
+extern double ** (*w_fn) (int);
+
+void set_wfn(char *name, int vfu);
+
+int gwr(struct rb *xbuf, int ninx, struct rb *ybuf, int cc,
+ int bw, double **w, DCELL *est, double **B0);
+
+int estimate_bandwidth(int *inx, int ninx, int iny, int nrows, int ncols,
+ DCELL *est, int bw);
Property changes on: grass-addons/grass7/raster/r.gwr/local_proto.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.gwr/main.c
===================================================================
--- grass-addons/grass7/raster/r.gwr/main.c (rev 0)
+++ grass-addons/grass7/raster/r.gwr/main.c 2014-06-27 06:22:32 UTC (rev 61005)
@@ -0,0 +1,479 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.gwr
+ *
+ * AUTHOR(S): Markus Metz
+ *
+ * PURPOSE: Geographically weighted regression
+ * Calculates multiple linear regression from raster maps:
+ * y = b0 + b1*x1 + b2*x2 + ... + bn*xn + e
+ * with localized b coefficients
+ *
+ * COPYRIGHT: (C) 2011 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
+ * for details.
+ *
+ *****************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/raster.h>
+#include "local_proto.h"
+
+
+int main(int argc, char *argv[])
+{
+ unsigned int r, c, rows, cols, count;
+ int *mapx_fd, mapy_fd, mapres_fd, mapest_fd;
+ int i, j, k, n_predictors;
+ double yres;
+ double *B, *Bmin, *Bmax, *Bsum, *Bsumsq, *Bmean, Bstddev;
+ int bcount;
+ DCELL *yest;
+ double sumY, meanY;
+ double SStot, SSerr, SSreg, *SSerr_without;
+ double Rsq, Rsqadj, SE, F, t, AIC, AICc, BIC;
+ DCELL *mapx_val, mapy_val, *mapy_buf, *mapres_buf, *mapest_buf;
+ struct rb *xbuf, ybuf;
+ double **weights;
+ int bw;
+ char *name;
+ struct Option *input_mapx, *input_mapy,
+ *output_res, *output_est, *output_opt,
+ *bw_opt, *kernel_opt, *vf_opt;
+ struct Flag *shell_style, *estimate;
+ struct Cell_head region;
+ struct GModule *module;
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("statistics"));
+ module->description =
+ _("Calculates geographically weighted regression from raster maps.");
+
+ /* Define the different options */
+ input_mapx = G_define_standard_option(G_OPT_R_INPUTS);
+ input_mapx->key = "mapx";
+ input_mapx->description = (_("Map for x coefficient"));
+
+ input_mapy = G_define_standard_option(G_OPT_R_INPUT);
+ input_mapy->key = "mapy";
+ input_mapy->description = (_("Map for y coefficient"));
+
+ output_res = G_define_standard_option(G_OPT_R_OUTPUT);
+ output_res->key = "residuals";
+ output_res->required = NO;
+ output_res->description = (_("Map to store residuals"));
+
+ output_est = G_define_standard_option(G_OPT_R_OUTPUT);
+ output_est->key = "estimates";
+ output_est->required = NO;
+ output_est->description = (_("Map to store estimates"));
+
+ output_opt = G_define_standard_option(G_OPT_F_OUTPUT);
+ output_opt->key = "output";
+ output_opt->required = NO;
+ output_opt->description =
+ (_("ASCII file for storing regression coefficients (output to screen if file not specified)."));
+
+ kernel_opt = G_define_option();
+ kernel_opt->key = "kernel";
+ kernel_opt->type = TYPE_STRING;
+ kernel_opt->options = "gauss,epanechnikov,quartic,tricubic";
+ kernel_opt->answer = "gauss";
+ kernel_opt->required = NO;
+ kernel_opt->description =
+ (_("Weighing kernel function."));
+
+ bw_opt = G_define_option();
+ bw_opt->key = "bandwidth";
+ bw_opt->type = TYPE_INTEGER;
+ bw_opt->answer = "10";
+ bw_opt->required = NO;
+ bw_opt->description =
+ (_("Bandwidth of the weighing kernel."));
+
+ vf_opt = G_define_option();
+ vf_opt->key = "vf";
+ vf_opt->type = TYPE_INTEGER;
+ vf_opt->options = "1,2,4,8";
+ vf_opt->answer = "1";
+ vf_opt->required = NO;
+ vf_opt->description =
+ (_("Variance factor for Gaussian kernel: variance = bandwith / factor."));
+
+ shell_style = G_define_flag();
+ shell_style->key = 'g';
+ shell_style->description = _("Print in shell script style");
+
+ estimate = G_define_flag();
+ estimate->key = 'e';
+ estimate->description = _("Estimate optimal bandwidth");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ name = output_opt->answer;
+ if (name != NULL && strcmp(name, "-") != 0) {
+ if (NULL == freopen(name, "w", stdout)) {
+ G_fatal_error(_("Unable to open file <%s> for writing"), name);
+ }
+ }
+
+ G_get_window(®ion);
+ rows = region.rows;
+ cols = region.cols;
+
+ /* count x maps */
+ for (i = 0; input_mapx->answers[i]; i++);
+ n_predictors = i;
+
+ /* allocate memory for x maps */
+ mapx_fd = (int *)G_malloc(n_predictors * sizeof(int));
+ SSerr_without = (double *)G_malloc(n_predictors * sizeof(double));
+ mapx_val = (DCELL *)G_malloc((n_predictors + 1) * sizeof(DCELL));
+ yest = G_malloc(sizeof(DCELL) * (n_predictors + 1));
+
+ bw = atoi(bw_opt->answer);
+ if (bw < 2)
+ G_fatal_error(_("Option <%s> must be > 1"), bw_opt->key);
+
+ set_wfn(kernel_opt->answer, atoi(vf_opt->answer));
+
+ /* open maps */
+ G_debug(1, "open maps");
+
+ for (i = 0; i < n_predictors; i++) {
+ mapx_fd[i] = Rast_open_old(input_mapx->answers[i], "");
+ }
+ mapy_fd = Rast_open_old(input_mapy->answer, "");
+
+ for (i = 0; i < n_predictors; i++) {
+ SSerr_without[i] = 0.0;
+ }
+ meanY = sumY = 0.0;
+
+ if (estimate->answer) {
+ bw = estimate_bandwidth(mapx_fd, n_predictors, mapy_fd, rows, cols,
+ yest, bw);
+ G_message("Estimated optimal bandwidth: %d", bw);
+
+ exit(EXIT_SUCCESS);
+ }
+
+ xbuf = (struct rb *)G_malloc(n_predictors * sizeof(struct rb));
+
+ for (i = 0; i < n_predictors; i++) {
+ allocate_bufs(&xbuf[i], cols, bw, mapx_fd[i]);
+ }
+ allocate_bufs(&ybuf, cols, bw, mapy_fd);
+
+ meanY = sumY = 0.0;
+
+ G_message(_("Calculating average..."));
+ count = 0;
+ mapy_buf = Rast_allocate_d_buf();
+
+ for (r = 0; r < rows; r++) {
+ G_percent(r, rows, 2);
+
+ Rast_get_d_row(mapy_fd, mapy_buf, r);
+
+ for (c = 0; c < cols; c++) {
+
+ mapy_val = mapy_buf[c];
+
+ if (!Rast_is_d_null_value(&mapy_val)) {
+ sumY += mapy_val;
+ count++;
+ }
+ }
+ }
+ G_percent(rows, rows, 2);
+
+ if (count < (unsigned int) (n_predictors + 1))
+ G_fatal_error(_("Not enough valid cells available"));
+
+ meanY = sumY / count;
+
+ /* residuals output */
+ if (output_res->answer) {
+ mapres_fd = Rast_open_new(output_res->answer, DCELL_TYPE);
+ mapres_buf = Rast_allocate_d_buf();
+ }
+ else {
+ mapres_fd = -1;
+ mapres_buf = NULL;
+ }
+
+ /* estimates output */
+ if (output_est->answer) {
+ mapest_fd = Rast_open_new(output_est->answer, DCELL_TYPE);
+ mapest_buf = Rast_allocate_d_buf();
+ }
+ else {
+ mapest_fd = -1;
+ mapest_buf = NULL;
+ }
+
+ /* gwr for each cell: get estimate */
+ G_message(_("Geographically weighted regression..."));
+
+ weights = w_fn(bw);
+ count = 0;
+ mapx_val[0] = 1.0;
+ SStot = SSerr = SSreg = 0.0;
+ for (i = 0; i < n_predictors; i++) {
+ SSerr_without[i] = 0.0;
+ }
+
+ /* initialize the raster buffers with 'bw' rows */
+ for (i = 0; i < bw; i++) {
+ for (j = 0; j < n_predictors; j++)
+ readrast(&(xbuf[j]), rows, cols);
+ readrast(&ybuf, rows, cols);
+ }
+
+ /* initialize coefficient statistics */
+ Bmin = G_malloc((n_predictors + 1) * sizeof(double));
+ Bmax = G_malloc((n_predictors + 1) * sizeof(double));
+ Bsum = G_malloc((n_predictors + 1) * sizeof(double));
+ Bsumsq = G_malloc((n_predictors + 1) * sizeof(double));
+ Bmean = G_malloc((n_predictors + 1) * sizeof(double));
+
+ for (i = 0; i <= n_predictors; i++) {
+ Bmin[i] = 1.0 / 0.0; /* inf */
+ Bmax[i] = -1.0 / 0.0; /* -inf */
+ Bsum[i] = Bsumsq[i] = Bmean[i] = 0;
+ }
+ bcount = 0;
+
+ for (r = 0; r < rows; r++) {
+ G_percent(r, rows, 2);
+
+ for (i = 0; i < n_predictors; i++) {
+ readrast(&(xbuf[i]), rows, cols);
+ }
+
+ readrast(&ybuf, rows, cols);
+
+ if (mapres_buf)
+ Rast_set_d_null_value(mapres_buf, cols);
+ if (mapest_buf)
+ Rast_set_d_null_value(mapest_buf, cols);
+
+ for (c = 0; c < cols; c++) {
+ int isnull = 0;
+
+ for (i = 0; i < n_predictors; i++) {
+ mapx_val[i + 1] = xbuf[i].buf[bw][c + bw];
+ if (Rast_is_d_null_value(&(mapx_val[i + 1]))) {
+ isnull = 1;
+ break;
+ }
+ }
+ if (isnull)
+ continue;
+
+ if (!gwr(xbuf, n_predictors, &ybuf, c, bw, weights, yest, &B)) {
+ continue;
+ }
+
+ /* coefficient stats */
+ for (k = 0; k <= n_predictors; k++) {
+ if (Bmin[k] > B[k])
+ Bmin[k] = B[k];
+ if (Bmax[k] < B[k])
+ Bmax[k] = B[k];
+ Bsum[k] += B[k];
+ Bsumsq[k] += B[k] * B[k];
+ }
+ bcount++;
+
+ /* set estimate */
+ if (mapest_buf)
+ mapest_buf[c] = yest[0];
+
+ mapy_val = ybuf.buf[bw][c + bw];
+ if (Rast_is_d_null_value(&mapy_val))
+ continue;
+
+ /* set residual */
+ yres = mapy_val - yest[0];
+ if (mapres_buf)
+ mapres_buf[c] = yres;
+
+ SStot += (mapy_val - meanY) * (mapy_val - meanY);
+ SSreg += (yest[0] - meanY) * (yest[0] - meanY);
+ SSerr += yres * yres;
+
+ for (k = 1; k <= n_predictors; k++) {
+
+ /* linear model without predictor k */
+ yres = mapy_val - yest[k];
+
+ /* linear model without predictor k */
+ SSerr_without[k - 1] += yres * yres;
+ }
+ count++;
+ }
+ if (mapres_buf)
+ Rast_put_d_row(mapres_fd, mapres_buf);
+ if (mapest_buf)
+ Rast_put_d_row(mapest_fd, mapest_buf);
+ }
+ G_percent(rows, rows, 2);
+
+ fprintf(stdout, "n=%d\n", count);
+ /* coefficient of determination aka R squared */
+ Rsq = 1 - (SSerr / SStot);
+ fprintf(stdout, "Rsq=%g\n", Rsq);
+ /* adjusted coefficient of determination */
+ Rsqadj = 1 - ((SSerr * (count - 1)) / (SStot * (count - n_predictors - 1)));
+ fprintf(stdout, "Rsqadj=%g\n", Rsqadj);
+ /* F statistic */
+ /* F = ((SStot - SSerr) / (n_predictors)) / (SSerr / (count - n_predictors));
+ * , or: */
+ F = ((SStot - SSerr) * (count - n_predictors - 1)) / (SSerr * (n_predictors));
+ fprintf(stdout, "F=%g\n", F);
+
+ i = 0;
+ /* constant aka estimate for intercept in R */
+ Bmean[i] = Bsum[i] / bcount;
+ fprintf(stdout, "bmean%d=%g\n", i, Bmean[i]);
+ Bstddev = sqrt(Bsumsq[i] / bcount - (Bmean[i] * Bmean[i]));
+ fprintf(stdout, "bstddev%d=%g\n", i, Bstddev);
+ fprintf(stdout, "bmin%d=%g\n", i, Bmin[i]);
+ fprintf(stdout, "bmax%d=%g\n", i, Bmax[i]);
+ /* t score for R squared of the full model, unused */
+ t = sqrt(Rsq) * sqrt((count - 2) / (1 - Rsq));
+ /*
+ fprintf(stdout, "t%d=%f\n", i, t);
+ */
+
+ /* AIC, corrected AIC, and BIC information criteria for the full model */
+ AIC = count * log(SSerr / count) + 2 * (n_predictors + 1);
+ fprintf(stdout, "AIC=%g\n", AIC);
+ AICc = AIC + (2 * n_predictors * (n_predictors + 1)) / (count - n_predictors - 1);
+ fprintf(stdout, "AICc=%g\n", AICc);
+ BIC = count * log(SSerr / count) + log(count) * (n_predictors + 1);
+ fprintf(stdout, "BIC=%g\n", BIC);
+
+ /* error variance of the model, identical to R */
+ SE = SSerr / (count - n_predictors - 1);
+ /*
+ fprintf(stdout, "SE=%f\n", SE);
+ fprintf(stdout, "SSerr=%f\n", SSerr);
+ */
+
+ /* variance = sumsq - (sum / n)^2 */
+
+ for (i = 0; i < n_predictors; i++) {
+
+ fprintf(stdout, "\npredictor%d=%s\n", i + 1, input_mapx->answers[i]);
+ Bmean[i + 1] = Bsum[i + 1] / bcount;
+ fprintf(stdout, "bmean%d=%g\n", i, Bmean[i + 1]);
+ Bstddev = sqrt(Bsumsq[i + 1] / bcount - (Bmean[i + 1] * Bmean[i + 1]));
+ fprintf(stdout, "bstddev%d=%g\n", i + 1, Bstddev);
+ fprintf(stdout, "bmin%d=%g\n", i + 1, Bmin[i + 1]);
+ fprintf(stdout, "bmax%d=%g\n", i + 1, Bmax[i + 1]);
+
+ if (n_predictors > 1) {
+ double Rsqi;
+
+ /* corrected sum of squares for predictor [i] */
+ /* sumsqX_corr = sumsqX[i] - sumX[i] * sumX[i] / (count - n_predictors - 1); */
+
+ /* standard error SE for predictor [i] */
+
+ /* SE[i] with only one predictor: sqrt(SE / sumsqX_corr)
+ * this does not work with more than one predictor */
+ /* in R, SEi is sqrt(diag(R) * resvar) with
+ * R = ???
+ * resvar = rss / rdf = SE global
+ * rss = sum of squares of the residuals
+ * rdf = residual degrees of freedom = count - n_predictors - 1 */
+ /* SEi = sqrt(SE / (Rsq * sumsqX_corr)); */
+ /*
+ fprintf(stdout, "SE%d=%f\n", i + 1, SEi);
+ */
+
+ /* Sum of squares for predictor [i] */
+ /*
+ fprintf(stdout, "SSerr%d=%f\n", i + 1, SSerr_without[i] - SSerr);
+ */
+
+ /* R squared of the model without predictor [i] */
+ /* Rsqi = 1 - SSerr_without[i] / SStot; */
+ /* the additional amount of variance explained
+ * when including predictor [i] :
+ * Rsq - Rsqi */
+ Rsqi = (SSerr_without[i] - SSerr) / SStot;
+ fprintf(stdout, "Rsq%d=%g\n", i + 1, Rsqi);
+
+ /*
+ fprintf(stdout, "t%d=%f\n", i + 1, t);
+ */
+
+ /* F score for Fisher's F distribution
+ * here: F score to test if including predictor [i]
+ * yields a significant improvement
+ * after Lothar Sachs, Angewandte Statistik:
+ * F = (Rsq - Rsqi) * (count - n_predictors - 1) / (1 - Rsq) */
+ /* same like Sumsq / SE */
+ /* same like (SSerr_without[i] / SSerr - 1) * (count - n_predictors - 1) */
+ /* same like R-stats when entered in R-stats as last predictor */
+ F = (SSerr_without[i] / SSerr - 1) * (count - n_predictors - 1);
+ fprintf(stdout, "F%d=%g\n", i + 1, F);
+
+ /* AIC, corrected AIC, and BIC information criteria for
+ * the model without predictor [i] */
+ AIC = count * log(SSerr_without[i] / count) + 2 * (n_predictors);
+ fprintf(stdout, "AIC%d=%g\n", i + 1, AIC);
+ AICc = AIC + (2 * (n_predictors - 1) * n_predictors) / (count - n_predictors - 2);
+ fprintf(stdout, "AICc%d=%g\n", i + 1, AICc);
+ BIC = count * log(SSerr_without[i] / count) + (n_predictors - 1) * log(count);
+ fprintf(stdout, "BIC%d=%g\n", i + 1, BIC);
+ }
+ }
+
+
+ for (i = 0; i < n_predictors; i++) {
+ Rast_close(mapx_fd[i]);
+ }
+ Rast_close(mapy_fd);
+ G_free(mapy_buf);
+
+ if (mapres_fd > -1) {
+ struct History history;
+
+ Rast_close(mapres_fd);
+ G_free(mapres_buf);
+
+ Rast_short_history(output_res->answer, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(output_res->answer, &history);
+ }
+
+ if (mapest_fd > -1) {
+ struct History history;
+
+ Rast_close(mapest_fd);
+ G_free(mapest_buf);
+
+ Rast_short_history(output_est->answer, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(output_est->answer, &history);
+ }
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass-addons/grass7/raster/r.gwr/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.gwr/r.gwr.html
===================================================================
--- grass-addons/grass7/raster/r.gwr/r.gwr.html (rev 0)
+++ grass-addons/grass7/raster/r.gwr/r.gwr.html 2014-06-27 06:22:32 UTC (rev 61005)
@@ -0,0 +1,77 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.gwr</em> calculates a geographically weighted multiple linear
+regression from raster maps, according to the formula
+<div class="code"><pre>
+Y = b0 + sum(bi*Xi) + E
+</pre></div>
+where
+<div class="code"><pre>
+X = {X1, X2, ..., Xm}
+m = number of explaining variables
+Y = {y1, y2, ..., yn}
+Xi = {xi1, xi2, ..., xin}
+E = {e1, e2, ..., en}
+n = number of observations (cases)
+</pre></div>
+
+In R notation:
+<div class="code"><pre>
+Y ~ sum(bi*Xi)
+b0 is the intercept, X0 is set to 1
+</pre></div>
+
+<p>
+The b coefficients are localized, i.e. determined for each cell
+individually.
+
+<p>
+Geographically weighted regressions should be used as a diagnostic tool
+and not as an interpolation method. If a geographically weighted
+regression provides a higher R<sup>2</sup> than the corresponding
+global regression, a crucial predictor is missing in the model. If that
+missing predictor can not be estimated or is supposed to behave
+randomly, a geographically weighted regression might be used for
+interpolation, but the result, in particular the variation of the
+β coefficients needs to be judged according to prior assumptions.
+
+<p>
+<em>r.gwr</em> is designed for large datasets that can not be processed
+in R. A p value is therefore not provided, because even very small,
+meaningless effects will become significant with a large number of
+cells. Instead it is recommended to judge by the amount of variance
+explained (R squared for a given variable) and the gain in AIC (AIC
+without a given variable minus AIC global must be positive) whether the
+inclusion of a given explaining variable in the model is justified.
+
+<h4>The explaining variables</h4>
+R squared for each explaining variable represents the additional amount
+of explained variance when including this variable compared to when
+excluding this variable, that is, this amount of variance is explained
+by the current explaining variable after taking into consideration all
+the other explaining variables.
+<p>
+The F score for each explaining variable allows to test if the inclusion
+of this variable significantly increases the explaining power of the
+model, relative to the global model excluding this explaining variable.
+That means that the F value for a given explaining variable is only
+identical to the F value of the R-function <em>summary.aov</em> if the
+given explaining variable is the last variable in the R-formula. While
+R successively includes one variable after another in the order
+specified by the formula and at each step calculates the F value
+expressing the gain by including the current variable in addition to the
+previous variables, <em>r.regression.multi</em> calculates the F-value
+expressing the gain by including the current variable in addition to all
+other variables, not only the previous variables.
+
+
+<h2>BUGS</h2>
+
+The estimated bandwidth is too small. A correcting factor has been
+introduced, but this is just a workaround.
+
+<h2>AUTHOR</h2>
+
+Markus Metz
+
+<p><i>Last changed: $Date$</i>
Property changes on: grass-addons/grass7/raster/r.gwr/r.gwr.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.gwr/weights.c
===================================================================
--- grass-addons/grass7/raster/r.gwr/weights.c (rev 0)
+++ grass-addons/grass7/raster/r.gwr/weights.c 2014-06-27 06:22:32 UTC (rev 61005)
@@ -0,0 +1,127 @@
+#include <math.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/gmath.h>
+
+/* bandwidth bw
+ * each cell within bw (distance from center <= bw)
+ * gets a weight > 0 */
+
+static double vf = 2;
+
+double ** (*w_fn) (int);
+
+double **epanechnikov(int bw)
+{
+ int r, c;
+ double bw2, bw2d, d, **w;
+
+ w = G_alloc_matrix(bw * 2 + 1, bw * 2 + 1);
+ bw2 = (bw + 1) * (bw + 1);
+ bw2d = bw * bw;
+
+ for (r = -bw; r <= bw; r++) {
+ for (c = -bw; c <= bw; c++) {
+ d = r * r + c * c;
+ w[r + bw][c + bw] = 0;
+
+ if (d <= bw2d) {
+ w[r + bw][c + bw] = 1 - d / bw2;
+ }
+ }
+ }
+
+ return w;
+}
+
+double **quartic(int bw)
+{
+ int r, c;
+ double bw2, bw2d, d, **w, t;
+
+ w = G_alloc_matrix(bw * 2 + 1, bw * 2 + 1);
+ bw2 = (bw + 1) * (bw + 1);
+ bw2d = bw * bw;
+
+ for (r = -bw; r <= bw; r++) {
+ for (c = -bw; c <= bw; c++) {
+ d = r * r + c * c;
+ w[r + bw][c + bw] = 0;
+
+ if (d <= bw2d) {
+ t = 1 - d / bw2;
+ w[r + bw][c + bw] = t * t;
+ }
+ }
+ }
+
+ return w;
+}
+
+double **tricubic(int bw)
+{
+ int r, c;
+ double bw3, bw3d, d, **w, t;
+
+ w = G_alloc_matrix(bw * 2 + 1, bw * 2 + 1);
+ bw3 = (bw + 1) * (bw + 1) * (bw + 1);
+ bw3d = bw * bw * bw;
+
+ for (r = -bw; r <= bw; r++) {
+ for (c = -bw; c <= bw; c++) {
+ d = sqrt(r * r + c * c);
+ d = d * d * d;
+ w[r + bw][c + bw] = 0;
+
+ if (d <= bw3d) {
+ t = 1 - d / bw3;
+ w[r + bw][c + bw] = t * t * t;
+ }
+ }
+ }
+
+ return w;
+}
+
+double **gauss(int bw)
+{
+ int r, c;
+ double bw2, d, **w;
+
+ w = G_alloc_matrix(bw * 2 + 1, bw * 2 + 1);
+ bw2 = bw * bw;
+
+ /* Gaussian function: exp(-x^2 / ( 2 * variance) */
+
+ for (r = -bw; r <= bw; r++) {
+ for (c = -bw; c <= bw; c++) {
+ d = r * r + c * c;
+ w[r + bw][c + bw] = 0;
+
+ if (d <= bw2) {
+ w[r + bw][c + bw] = exp(vf * d / bw2);
+ w[r + bw][c + bw] = 1;
+ }
+ }
+ }
+
+ return w;
+}
+
+/* set weighing kernel function and variance factor */
+void set_wfn(char *name, int vfu)
+{
+ vf = vfu / 2.;
+
+ if (*name == 'g')
+ w_fn = gauss;
+ else if (*name == 'e')
+ w_fn = epanechnikov;
+ else if (*name == 'q')
+ w_fn = quartic;
+ else if (*name == 't')
+ w_fn = tricubic;
+ else
+ G_fatal_error(_("Invalid kernel option '%s'"), name);
+}
Property changes on: grass-addons/grass7/raster/r.gwr/weights.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list