[GRASS-SVN] r69864 - grass-addons/grass7/raster/r.gwr
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Nov 21 12:59:01 PST 2016
Author: mmetz
Date: 2016-11-21 12:59:01 -0800 (Mon, 21 Nov 2016)
New Revision: 69864
Added:
grass-addons/grass7/raster/r.gwr/flag.c
grass-addons/grass7/raster/r.gwr/flag.h
grass-addons/grass7/raster/r.gwr/gwr.h
grass-addons/grass7/raster/r.gwr/gwra.c
grass-addons/grass7/raster/r.gwr/rclist.c
grass-addons/grass7/raster/r.gwr/rclist.h
Modified:
grass-addons/grass7/raster/r.gwr/Makefile
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
Log:
r.gwr: +adaptive bandwidth, +mask option
Modified: grass-addons/grass7/raster/r.gwr/Makefile
===================================================================
--- grass-addons/grass7/raster/r.gwr/Makefile 2016-11-21 20:49:11 UTC (rev 69863)
+++ grass-addons/grass7/raster/r.gwr/Makefile 2016-11-21 20:59:01 UTC (rev 69864)
@@ -2,8 +2,8 @@
PGM = r.gwr
-LIBES = $(RASTERLIB) $(GISLIB) $(GMATHLIB) $(MATHLIB)
-DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+LIBES = $(RASTERLIB) $(BTREE2LIB) $(SEGMENTLIB) $(GISLIB) $(GMATHLIB) $(MATHLIB)
+DEPENDENCIES = $(RASTERDEP) $(BTREE2DEP) $(SEGMENTDEP) $(GISDEP)
EXTRA_LIBS = $(OMPLIB)
EXTRA_CFLAGS = $(OMPCFLAGS)
Modified: grass-addons/grass7/raster/r.gwr/estimate_bw.c
===================================================================
--- grass-addons/grass7/raster/r.gwr/estimate_bw.c 2016-11-21 20:49:11 UTC (rev 69863)
+++ grass-addons/grass7/raster/r.gwr/estimate_bw.c 2016-11-21 20:59:01 UTC (rev 69864)
@@ -140,9 +140,9 @@
readrast(&ybuf3, nrows, ncols);
}
- w1 = w_fn(prevbw);
- w2 = w_fn(bw);
- w3 = w_fn(nextbw);
+ w1 = calc_weights(prevbw);
+ w2 = calc_weights(bw);
+ w3 = calc_weights(nextbw);
/* leave one out: the center cell */
w1[prevbw][prevbw] = 0.;
w2[bw][bw] = 0.;
Added: grass-addons/grass7/raster/r.gwr/flag.c
===================================================================
--- grass-addons/grass7/raster/r.gwr/flag.c (rev 0)
+++ grass-addons/grass7/raster/r.gwr/flag.c 2016-11-21 20:59:01 UTC (rev 69864)
@@ -0,0 +1,60 @@
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "flag.h"
+
+FLAG *flag_create(int nrows, int ncols)
+{
+ unsigned char *temp;
+
+ FLAG *new_flag;
+
+ register int i;
+
+ new_flag = (FLAG *) G_malloc(sizeof(FLAG));
+ new_flag->nrows = nrows;
+ new_flag->ncols = ncols;
+ new_flag->leng = (ncols + 7) / 8;
+ new_flag->array =
+ (unsigned char **)G_malloc(nrows * sizeof(unsigned char *));
+
+ if (!new_flag->array)
+ G_fatal_error(_("Out of memory!"));
+
+ temp =
+ (unsigned char *)G_malloc(nrows * new_flag->leng *
+ sizeof(unsigned char));
+
+ if (!temp)
+ G_fatal_error(_("Out of memory!"));
+
+ for (i = 0; i < nrows; i++) {
+ new_flag->array[i] = temp;
+ temp += new_flag->leng;
+ }
+
+ flag_clear_all(new_flag);
+
+ return (new_flag);
+}
+
+int flag_destroy(FLAG * flags)
+{
+ G_free(flags->array[0]);
+ G_free(flags->array);
+ G_free(flags);
+
+ return 0;
+}
+
+int flag_clear_all(FLAG * flags)
+{
+ register int r, c;
+
+ for (r = 0; r < flags->nrows; r++) {
+ for (c = 0; c < flags->leng; c++) {
+ flags->array[r][c] = 0;
+ }
+ }
+
+ return 0;
+}
Added: grass-addons/grass7/raster/r.gwr/flag.h
===================================================================
--- grass-addons/grass7/raster/r.gwr/flag.h (rev 0)
+++ grass-addons/grass7/raster/r.gwr/flag.h 2016-11-21 20:59:01 UTC (rev 69864)
@@ -0,0 +1,64 @@
+#ifndef __FLAG_H__
+#define __FLAG_H__
+
+
+/* flag.[ch] is a set of routines which will set up an array of bits
+ ** that allow the programmer to "flag" cells in a raster map.
+ **
+ ** FLAG *
+ ** flag_create(nrows,ncols)
+ ** int nrows, ncols;
+ ** opens the structure flag.
+ ** The flag structure will be a two dimensional array of bits the
+ ** size of nrows by ncols. Will initalize flags to zero (unset).
+ **
+ ** flag_destroy(flags)
+ ** FLAG *flags;
+ ** closes flags and gives the memory back to the system.
+ **
+ ** flag_clear_all(flags)
+ ** FLAG *flags;
+ ** sets all values in flags to zero.
+ **
+ ** flag_unset(flags, row, col)
+ ** FLAG *flags;
+ ** int row, col;
+ ** sets the value of (row, col) in flags to zero.
+ **
+ ** flag_set(flags, row, col)
+ ** FLAG *flags;
+ ** int row, col;
+ ** will set the value of (row, col) in flags to one.
+ **
+ ** int
+ ** flag_get(flags, row, col)
+ ** FLAG *flags;
+ ** int row, col;
+ ** returns the value in flags that is at (row, col).
+ **
+ ** idea by Michael Shapiro
+ ** code by Chuck Ehlschlaeger
+ ** April 03, 1989
+ */
+#define FLAG struct _flagsss_
+FLAG {
+ int nrows, ncols, leng;
+ unsigned char **array;
+};
+
+#define FLAG_UNSET(flags,row,col) \
+ (flags)->array[(row)][(col)>>3] &= ~(1<<((col) & 7))
+
+#define FLAG_SET(flags,row,col) \
+ (flags)->array[(row)][(col)>>3] |= (1<<((col) & 7))
+
+#define FLAG_GET(flags,row,col) \
+ (flags)->array[(row)][(col)>>3] & (1<<((col) & 7))
+
+/* flag.c */
+FLAG *flag_create(int, int);
+int flag_destroy(FLAG *);
+int flag_clear_all(FLAG *);
+
+
+#endif /* __FLAG_H__ */
Modified: grass-addons/grass7/raster/r.gwr/gwr.c
===================================================================
--- grass-addons/grass7/raster/r.gwr/gwr.c 2016-11-21 20:49:11 UTC (rev 69863)
+++ grass-addons/grass7/raster/r.gwr/gwr.c 2016-11-21 20:59:01 UTC (rev 69864)
@@ -6,20 +6,13 @@
#include <grass/glocale.h>
#include <grass/raster.h>
#include "local_proto.h"
+#include "gwr.h"
/* geographically weighted regression:
* estimate coefficients for given cell */
-struct MATRIX
+int solvemat(struct MATRIX *m, double a[], double B[])
{
- 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 */
@@ -43,7 +36,7 @@
/* co-linear points result in an undefined matrix, and nearly */
/* co-linear points results in a solution with rounding error */
- if (pivot == 0.0) {
+ if (fabs(pivot) < GRASS_EPSILON) {
G_debug(1, "Matrix is unsolvable");
return 0;
}
@@ -99,7 +92,7 @@
int nsize;
static DCELL *xval = NULL;
DCELL yval;
- int count, isnull;;
+ int count, isnull, solved;
/* OLS */
static double **a = NULL;
static double **B = NULL;
@@ -210,17 +203,17 @@
}
}
}
-
count++;
}
}
-
- /* estimate coefficients */
+
if (count < ninx + 1) {
- G_debug(4, "Not enough valid cells available");
+ G_message(_("Unable to determine coefficients. Consider increasing the bandwidth."));
return 0;
}
+ /* estimate coefficients */
+ solved = ninx + 1;
for (k = 0; k <= ninx; k++) {
m = &(m_all[k]);
@@ -235,17 +228,16 @@
fprintf(stdout, "b%d=0.0\n", i);
}
*/
- G_debug(1, "Try a larger bandwidth");
- count = 0;
+ G_debug(1, "Solving matrix %d failed", k);
+ solved--;
}
}
- if (count < ninx + 1) {
+ if (solved < ninx + 1) {
+ G_debug(3, "%d of %d equation systems could not be solved", ninx + 1 - solved, ninx + 1);
return 0;
}
/* second pass: calculate estimates */
- c = cc;
-
isnull = 0;
for (i = 0; i < ninx; i++) {
Added: grass-addons/grass7/raster/r.gwr/gwr.h
===================================================================
--- grass-addons/grass7/raster/r.gwr/gwr.h (rev 0)
+++ grass-addons/grass7/raster/r.gwr/gwr.h 2016-11-21 20:59:01 UTC (rev 69864)
@@ -0,0 +1,10 @@
+struct MATRIX
+{
+ int n; /* SIZE OF THIS MATRIX (N x N) */
+ double **v;
+};
+
+#define M(m,row,col) (m)->v[(row)][(col)]
+
+
+int solvemat(struct MATRIX *m, double a[], double B[]);
Added: grass-addons/grass7/raster/r.gwr/gwra.c
===================================================================
--- grass-addons/grass7/raster/r.gwr/gwra.c (rev 0)
+++ grass-addons/grass7/raster/r.gwr/gwra.c 2016-11-21 20:59:01 UTC (rev 69864)
@@ -0,0 +1,327 @@
+#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 <grass/rbtree.h>
+#include "rclist.h"
+#include "local_proto.h"
+#include "gwr.h"
+
+/* geographically weighted regression:
+ * estimate coefficients for given cell */
+
+
+struct gwr_pnt
+{
+ int r, c; /* row, col in target window */
+};
+
+static int cmp_rc(const void *first, const void *second)
+{
+ struct rc *a = (struct rc *)first, *b = (struct rc *)second;
+
+ if (a->row == b->row)
+ return (a->col - b->col);
+
+ return (a->row - b->row);
+}
+
+static int cmp_pnts(const void *first, const void *second)
+{
+ struct gwr_pnt *a = (struct gwr_pnt *)first;
+ struct gwr_pnt *b = (struct gwr_pnt *)second;
+
+ if (a->r == b->r)
+ return (a->c - b->c);
+
+ return (a->r - b->r);
+}
+
+static int bfs_search(FLAG *null_flag, struct gwr_pnt *cur_pnts,
+ int max_points, int row, int col, double *distmax)
+{
+ int nrows, ncols, rown, coln;
+ int nextr[8] = {0, -1, 0, 1, -1, -1, 1, 1};
+ int nextc[8] = {1, 0, -1, 0, 1, -1, -1, 1};
+ int n, found;
+ struct rc next, ngbr_rc;
+ struct rclist rilist;
+ struct RB_TREE *visited;
+ double dx, dy, dist;
+
+ visited = rbtree_create(cmp_rc, sizeof(struct rc));
+ ngbr_rc.row = row;
+ ngbr_rc.col = col;
+ rbtree_insert(visited, &ngbr_rc);
+
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ found = 0;
+ *distmax = 0.0;
+
+ if (FLAG_GET(null_flag, row, col)) {
+ /* add to cur_pnts */
+ cur_pnts[found].r = row;
+ cur_pnts[found].c = col;
+ found++;
+ }
+
+ /* breadth-first search */
+ next.row = row;
+ next.col = col;
+ rclist_init(&rilist);
+
+ do {
+ for (n = 0; n < 8; n++) {
+ rown = next.row + nextr[n];
+ coln = next.col + nextc[n];
+
+ if (rown < 0 || rown >= nrows || coln < 0 || coln >= ncols)
+ continue;
+
+ ngbr_rc.row = rown;
+ ngbr_rc.col = coln;
+
+ if (rbtree_find(visited, &ngbr_rc))
+ continue;
+
+ rbtree_insert(visited, &ngbr_rc);
+ rclist_add(&rilist, rown, coln);
+
+ if (FLAG_GET(null_flag, rown, coln)) {
+
+ dx = coln - col;
+ dy = rown - row;
+
+ dist = dx * dx + dy * dy;
+
+ if (*distmax < dist)
+ *distmax = dist;
+
+ cur_pnts[found].r = rown;
+ cur_pnts[found].c = coln;
+ found++;
+ }
+
+ if (found == max_points)
+ break;
+ }
+ if (found == max_points)
+ break;
+ } while (rclist_drop(&rilist, &next)); /* while there are cells to check */
+
+
+ rclist_destroy(&rilist);
+ rbtree_destroy(visited);
+
+ return found;
+}
+
+int gwra(SEGMENT *in_seg, FLAG *null_flag, int ninx, int rr, int cc,
+ int npnts, DCELL *est, double **B0)
+{
+ int r, c, n, nfound;
+ int i, j, k;
+ double bw, maxdist2, dist2, w;
+ static DCELL *xval = NULL, *seg_val = NULL;
+ DCELL yval;
+ int count, isnull, solved;
+ static int npnts_alloc = 0;
+ static struct gwr_pnt *cur_pnts = NULL;
+ /* 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.;
+ seg_val = G_malloc((ninx + 1) * sizeof(DCELL));
+ }
+
+ 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));
+ }
+ }
+
+ if (npnts_alloc < npnts) {
+ if (cur_pnts)
+ G_free(cur_pnts);
+ npnts_alloc = npnts;
+ cur_pnts = G_malloc(sizeof(struct gwr_pnt) * npnts_alloc);
+ }
+
+ 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;
+
+ /* first pass: collect values */
+ nfound = bfs_search(null_flag, cur_pnts, npnts, rr, cc, &maxdist2);
+
+ qsort(cur_pnts, nfound, sizeof(struct gwr_pnt), cmp_pnts);
+
+ /* bandwidth */
+ bw = sqrt(maxdist2);
+
+ /* second pass: load values */
+ count = 0;
+ for (n = 0; n < nfound; n++) {
+
+ r = cur_pnts[n].r;
+ c = cur_pnts[n].c;
+
+ Segment_get(in_seg, (void *)seg_val, r, c);
+
+ isnull = 0;
+ for (i = 0; i < ninx; i++) {
+ xval[i + 1] = seg_val[i];
+ if (Rast_is_d_null_value(&(xval[i + 1]))) {
+ isnull = 1;
+ break;
+ }
+ }
+ if (isnull)
+ continue;
+
+ yval = seg_val[ninx];
+ if (Rast_is_d_null_value(&yval))
+ continue;
+
+ dist2 = (r - rr) * (r - rr) + (c - cc) * (c - cc);
+ w = w_fn(dist2, bw);
+
+ 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;
+
+ /* 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;
+ }
+ }
+ }
+
+ a[0][i] += yval * val1 * w;
+ for (k = 1; k <= ninx; k++) {
+ if (k != i) {
+ int i2 = k > i ? i : i - 1;
+
+ a[k][i2] += yval * val1 * w;
+ }
+ }
+ }
+ count++;
+ }
+ if (count < ninx + 1) {
+ G_message(_("Unable to determine coefficients. Consider increasing the number of points."));
+ return 0;
+ }
+
+ /* estimate coefficients */
+ solved = ninx + 1;
+ 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);
+ }
+ */
+ solved--;
+ }
+ }
+ if (solved < ninx + 1) {
+ G_debug(3, "%d of %d equation systems could not be solved", ninx + 1 - solved, ninx + 1);
+ return 0;
+ }
+
+ /* third pass: calculate estimates */
+ Segment_get(in_seg, (void *)seg_val, rr, cc);
+ isnull = 0;
+ for (i = 0; i < ninx; i++) {
+
+ xval[i + 1] = seg_val[i];
+ 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;
+}
Modified: grass-addons/grass7/raster/r.gwr/local_proto.h
===================================================================
--- grass-addons/grass7/raster/r.gwr/local_proto.h 2016-11-21 20:49:11 UTC (rev 69863)
+++ grass-addons/grass7/raster/r.gwr/local_proto.h 2016-11-21 20:59:01 UTC (rev 69864)
@@ -1,11 +1,18 @@
+#include <grass/segment.h>
+#include "flag.h"
#include "bufs.h"
-extern double ** (*w_fn) (int);
-
void set_wfn(char *name, int vfu);
+extern double (*w_fn) (double, double);
+
+double **calc_weights(int bw);
+
int gwr(struct rb *xbuf, int ninx, struct rb *ybuf, int cc,
int bw, double **w, DCELL *est, double **B0);
+int gwra(SEGMENT *in_seg, FLAG *yflag, int ninx, int rr, int cc,
+ int npnts, DCELL *est, double **B0);
+
int estimate_bandwidth(int *inx, int ninx, int iny, int nrows, int ncols,
DCELL *est, int bw);
Modified: grass-addons/grass7/raster/r.gwr/main.c
===================================================================
--- grass-addons/grass7/raster/r.gwr/main.c 2016-11-21 20:49:11 UTC (rev 69863)
+++ grass-addons/grass7/raster/r.gwr/main.c 2016-11-21 20:59:01 UTC (rev 69864)
@@ -10,7 +10,7 @@
* y = b0 + b1*x1 + b2*x2 + ... + bn*xn + e
* with localized b coefficients
*
- * COPYRIGHT: (C) 2011 by the GRASS Development Team
+ * COPYRIGHT: (C) 2011-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
@@ -31,7 +31,7 @@
int main(int argc, char *argv[])
{
unsigned int r, c, rows, cols, count;
- int *mapx_fd, mapy_fd, mapres_fd, mapest_fd;
+ int *mapx_fd, mapy_fd, mapres_fd, mapest_fd, mask_fd;
int i, j, k, n_predictors;
double yres;
double *B, *Bmin, *Bmax, *Bsum, *Bsumsq, *Bmean, Bstddev;
@@ -41,20 +41,25 @@
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;
+ CELL *mask_buf;
struct rb *xbuf, ybuf;
+ SEGMENT in_seg;
+ DCELL **segx_buf, *seg_val;
+ int segsize, nseg;
+ double mem_mb;
+ FLAG *null_flag;
struct outputb {
int fd;
char name[GNAME_MAX];
DCELL *buf;
} *outb, *outbp;
double **weights;
- int bw;
+ int bw, npnts;
char *name;
- struct Option *input_mapx, *input_mapy,
+ struct Option *input_mapx, *input_mapy, *mask_opt,
*output_res, *output_est, *output_b, *output_opt,
- *bw_opt, *kernel_opt, *vf_opt;
+ *kernel_opt, *vf_opt, *bw_opt, *pnts_opt, *mem_opt;
struct Flag *shell_style, *estimate;
- /*struct Flag *adaptive;*/
struct Cell_head region;
struct GModule *module;
@@ -75,6 +80,12 @@
input_mapy->key = "mapy";
input_mapy->description = (_("Map with Y variable"));
+ mask_opt = G_define_standard_option(G_OPT_R_INPUT);
+ mask_opt->key = "mask";
+ mask_opt->label = _("Raster map to use for masking");
+ mask_opt->description = _("Only cells that are not NULL and not zero are processed");
+ mask_opt->required = NO;
+
output_res = G_define_standard_option(G_OPT_R_OUTPUT);
output_res->key = "residuals";
output_res->required = NO;
@@ -110,7 +121,7 @@
bw_opt->key = "bandwidth";
bw_opt->type = TYPE_INTEGER;
bw_opt->answer = "10";
- bw_opt->required = YES;
+ bw_opt->required = NO;
bw_opt->description =
(_("Bandwidth of the weighing kernel."));
@@ -123,6 +134,21 @@
vf_opt->description =
(_("Variance factor for Gaussian kernel: variance = bandwith / factor."));
+ pnts_opt = G_define_option();
+ pnts_opt->key = "npoints";
+ pnts_opt->type = TYPE_INTEGER;
+ pnts_opt->required = NO;
+ pnts_opt->answer = "0";
+ pnts_opt->label = _("Number of points for adaptive bandwidth");
+ pnts_opt->description = _("If 0, fixed bandwidth is used");
+
+ mem_opt = G_define_option();
+ mem_opt->key = "memory";
+ mem_opt->type = TYPE_INTEGER;
+ mem_opt->required = NO;
+ mem_opt->answer = "300";
+ mem_opt->description = _("Memory in MB for adaptive bandwidth");
+
shell_style = G_define_flag();
shell_style->key = 'g';
shell_style->description = _("Print in shell script style");
@@ -131,10 +157,6 @@
estimate->key = 'e';
estimate->description = _("Estimate optimal bandwidth");
- /*adaptive = G_define_flag();
- adaptive->key = 'a';
- adaptive->description = _("Adaptive kernel size");*/
-
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -163,9 +185,6 @@
if (bw < 2)
G_fatal_error(_("Option <%s> must be > 1"), bw_opt->key);
- /*if(adaptive->answer)*/
- /*IN MAKING*/
-
set_wfn(kernel_opt->answer, atoi(vf_opt->answer));
/* open maps */
@@ -176,6 +195,13 @@
}
mapy_fd = Rast_open_old(input_mapy->answer, "");
+ mask_fd = -1;
+ mask_buf = NULL;
+ if (mask_opt->answer) {
+ mask_fd = Rast_open_old(mask_opt->answer, "");
+ mask_buf = Rast_allocate_c_buf();
+ }
+
for (i = 0; i < n_predictors; i++) {
SSerr_without[i] = 0.0;
}
@@ -192,12 +218,21 @@
exit(EXIT_SUCCESS);
}
+ npnts = 0;
+ if (pnts_opt->answer) {
+ npnts = atoi(pnts_opt->answer);
+ if (npnts != 0 && npnts < n_predictors + 1)
+ G_fatal_error(_("Option <%s> must be > %d"), pnts_opt->key, n_predictors + 1);
+ }
+
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]);
+ if (npnts == 0) {
+ for (i = 0; i < n_predictors; i++) {
+ allocate_bufs(&xbuf[i], cols, bw, mapx_fd[i]);
+ }
+ allocate_bufs(&ybuf, cols, bw, mapy_fd);
}
- allocate_bufs(&ybuf, cols, bw, mapy_fd);
meanY = sumY = 0.0;
@@ -225,32 +260,29 @@
if (count < (unsigned int) (n_predictors + 1))
G_fatal_error(_("Not enough valid cells available"));
+ if (npnts > 0 && (unsigned int) npnts > count)
+ G_fatal_error(_("The number of points for adaptive bandwidth is larger than the number of valid cells"));
+
meanY = sumY / count;
/* residuals output */
+ mapres_fd = -1;
+ mapres_buf = NULL;
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 */
+ mapest_fd = -1;
+ mapest_buf = NULL;
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;
@@ -258,13 +290,82 @@
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);
+ segsize = 0;
+ seg_val = NULL;
+ null_flag = NULL;
+ weights = NULL;
+ if (npnts == 0) {
+ weights = calc_weights(bw);
+
+ /* 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);
+ }
}
+ else {
+ G_message(_("Loading input maps to temporary file..."));
+ segsize = sizeof(DCELL) * (n_predictors + 1);
+ seg_val = G_malloc(segsize);
+
+ segx_buf = G_malloc(sizeof(DCELL *) * n_predictors);
+ for (i = 0; i < n_predictors; i++)
+ segx_buf[i] = Rast_allocate_d_buf();
+
+ mem_mb = 300;
+ if (mem_opt->answer)
+ mem_mb = atoi(mem_opt->answer);
+ if (mem_mb < 10)
+ mem_mb = 10;
+
+ nseg = mem_mb * 1024.0 / (64 * 64 * segsize);
+
+ Segment_open(&in_seg, G_tempfile(), rows, cols, 64, 64,
+ segsize, nseg);
+
+ null_flag = flag_create(rows, cols);
+
+ for (r = 0; r < rows; r++) {
+ G_percent(r, rows, 2);
+
+ for (i = 0; i < n_predictors; i++) {
+ Rast_get_d_row(mapx_fd[i], segx_buf[i], r);
+ }
+ Rast_get_d_row(mapy_fd, mapy_buf, r);
+
+ for (c = 0; c < cols; c++) {
+ int x_null;
+
+ x_null = 0;
+
+ for (i = 0; i < n_predictors; i++) {
+ seg_val[i] = segx_buf[i][c];
+ if (Rast_is_d_null_value(&seg_val[i])) {
+ Rast_set_d_null_value(seg_val, n_predictors + 1);
+ x_null = 1;
+ break;
+ }
+ }
+ if (!x_null) {
+ seg_val[n_predictors] = mapy_buf[c];
+ }
+ if (Segment_put(&in_seg, (void *)seg_val, r, c) != 1)
+ G_fatal_error(_("Unable to write to temporary file"));
+
+ if (!x_null && !Rast_is_d_null_value(&mapy_buf[c]))
+ FLAG_SET(null_flag, r, c);
+ }
+ }
+ G_percent(1, 1, 1);
+
+ for (i = 0; i < n_predictors; i++)
+ G_free(segx_buf[i]);
+ G_free(segx_buf);
+ }
+ G_free(mapy_buf);
+
/* initialize coefficient statistics */
Bmin = G_malloc((n_predictors + 1) * sizeof(double));
Bmax = G_malloc((n_predictors + 1) * sizeof(double));
@@ -293,15 +394,20 @@
}
}
+ G_message(_("Geographically weighted regression..."));
for (r = 0; r < rows; r++) {
G_percent(r, rows, 2);
- for (i = 0; i < n_predictors; i++) {
- readrast(&(xbuf[i]), rows, cols);
+ if (npnts == 0) {
+ for (i = 0; i < n_predictors; i++) {
+ readrast(&(xbuf[i]), rows, cols);
+ }
+ readrast(&ybuf, rows, cols);
}
+
+ if (mask_buf)
+ Rast_get_c_row(mask_fd, mask_buf, r);
- readrast(&ybuf, rows, cols);
-
if (mapres_buf)
Rast_set_d_null_value(mapres_buf, cols);
if (mapest_buf)
@@ -316,21 +422,43 @@
for (c = 0; c < cols; c++) {
int isnull = 0;
+
+ if (mask_buf) {
+ if (Rast_is_c_null_value(&mask_buf[c]) || mask_buf[c] == 0)
+ continue;
+ }
- 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]))) {
+ if (npnts == 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;
+ }
+ }
+ mapy_val = ybuf.buf[bw][c + bw];
+ }
+ else {
+ Segment_get(&in_seg, (void *)seg_val, r, c);
+ if (Rast_is_d_null_value(&(seg_val[0]))) {
isnull = 1;
- break;
}
+ mapy_val = seg_val[n_predictors];
}
+
if (isnull)
continue;
- if (!gwr(xbuf, n_predictors, &ybuf, c, bw, weights, yest, &B)) {
- G_warning(_("Unable to determine coefficients. Consider increasing the bandwidth."));
- continue;
+ if (npnts == 0) {
+ if (!gwr(xbuf, n_predictors, &ybuf, c, bw, weights, yest, &B)) {
+ continue;
+ }
}
+ else {
+ if (!gwra(&in_seg, null_flag, n_predictors, r, c, npnts, yest, &B)) {
+ continue;
+ }
+ }
/* coefficient stats */
for (i = 0; i <= n_predictors; i++) {
@@ -353,7 +481,6 @@
if (mapest_buf)
mapest_buf[c] = yest[0];
- mapy_val = ybuf.buf[bw][c + bw];
if (Rast_is_d_null_value(&mapy_val))
continue;
@@ -507,8 +634,13 @@
Rast_close(mapx_fd[i]);
}
Rast_close(mapy_fd);
- G_free(mapy_buf);
-
+
+ if (mask_buf)
+ Rast_close(mask_fd);
+
+ if (npnts > 0)
+ Segment_close(&in_seg);
+
if (mapres_fd > -1) {
struct History history;
Modified: grass-addons/grass7/raster/r.gwr/r.gwr.html
===================================================================
--- grass-addons/grass7/raster/r.gwr/r.gwr.html 2016-11-21 20:49:11 UTC (rev 69863)
+++ grass-addons/grass7/raster/r.gwr/r.gwr.html 2016-11-21 20:59:01 UTC (rev 69864)
@@ -77,9 +77,18 @@
produce results similar to a global regression, and spatial
non-stationarity can not be explored.
+<h4>Adaptive bandwidth</h4>
+Instead of using a fixed bandwidth (search radius for each cell), an
+adaptive bandwidth can be used by specifying the number of points to be
+used for each local regression with the <em>npoints</em> option. The
+module will find the nearest <em>npoints</em> points for each cell,
+adapt the bandwith accordingly and then calculate a local weighted
+regression.
+
<h4>Kernel functions</h4>
-The kernel function has little influence on the result, much more important is the
-bandwidth. Available kernel funtions to calculate weights are
+The kernel function has little influence on the result, much more
+important is the bandwidth. Available kernel funtions to calculate
+weights are
<dl>
<dt><b>Epanechnikov</b></dt>
@@ -97,6 +106,10 @@
d = distance to the current cell<br>
bw = bandwidth
+<h4>Masking</h4>
+A <em>mask</em> map can be provided to restrict LWR to those cells
+where the mask map is not NULL and not 0 (zero).
+
<h2>REFERENCES</h2>
Brunsdon, C., Fotheringham, A.S., and Charlton, M.E., 1996,
Added: grass-addons/grass7/raster/r.gwr/rclist.c
===================================================================
--- grass-addons/grass7/raster/r.gwr/rclist.c (rev 0)
+++ grass-addons/grass7/raster/r.gwr/rclist.c 2016-11-21 20:59:01 UTC (rev 69864)
@@ -0,0 +1,68 @@
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "rclist.h"
+
+void rclist_init(struct rclist *list)
+{
+ list->head = list->tail = NULL;
+
+ return;
+}
+
+void rclist_add(struct rclist *list, int row, int col)
+{
+ struct rc *new = G_malloc(sizeof(struct rc));
+
+ if (!new)
+ G_fatal_error(_("rclist out of memory"));
+
+ new->next = NULL;
+ new->row = row;
+ new->col = col;
+
+ if (list->head) {
+ list->head->next = new;
+ list->head = new;
+ }
+ else {
+ list->head = list->tail = new;
+ }
+
+ return;
+}
+
+/* return 1 if an element was dropped
+ * return 0 if list is empty
+ */
+int rclist_drop(struct rclist *list, struct rc *rc)
+{
+ if (list->tail) {
+ struct rc *next = list->tail->next;
+
+ rc->row = list->tail->row;
+ rc->col = list->tail->col;
+ G_free(list->tail);
+ list->tail = next;
+ if (!list->tail)
+ list->head = NULL;
+
+ return 1;
+ }
+
+ return 0;
+}
+
+void rclist_destroy(struct rclist *list)
+{
+ struct rc *next = list->tail;
+
+ while (next) {
+ next = next->next;
+ G_free(list->tail);
+ list->tail = next;
+ }
+ list->head = NULL;
+
+ return;
+}
+
Added: grass-addons/grass7/raster/r.gwr/rclist.h
===================================================================
--- grass-addons/grass7/raster/r.gwr/rclist.h (rev 0)
+++ grass-addons/grass7/raster/r.gwr/rclist.h 2016-11-21 20:59:01 UTC (rev 69864)
@@ -0,0 +1,20 @@
+/* row/col list */
+
+struct rc
+{
+ struct rc *next;
+ int row;
+ int col;
+};
+
+struct rclist
+{
+ struct rc *tail, *head;
+};
+
+/* rclist.c */
+void rclist_init(struct rclist *);
+void rclist_add(struct rclist *, int, int);
+int rclist_drop(struct rclist *, struct rc *);
+void rclist_destroy(struct rclist *);
+
Modified: grass-addons/grass7/raster/r.gwr/weights.c
===================================================================
--- grass-addons/grass7/raster/r.gwr/weights.c 2016-11-21 20:49:11 UTC (rev 69863)
+++ grass-addons/grass7/raster/r.gwr/weights.c 2016-11-21 20:59:01 UTC (rev 69864)
@@ -10,99 +10,72 @@
static double vf = -0.5;
-double ** (*w_fn) (int);
+double (*w_fn) (double, double);
-double **epanechnikov(int bw)
+double epanechnikov(double d2, double bw)
{
- int r, c;
- double bw2, bw2d, d, **w;
+ double bw2, bw2d, w;
- w = G_alloc_matrix(bw * 2 + 1, bw * 2 + 1);
bw2 = (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;
- }
- }
+ w = 0;
+
+ if (d2 <= bw2d) {
+ w = 1 - d2 / bw2;
}
-
+
return w;
}
-double **bisquare(int bw)
+double bisquare(double d2, double bw)
{
- int r, c;
- double bw2, bw2d, d, **w, t;
+ double bw2, bw2d, 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;
- }
- }
+ w = 0;
+
+ if (d2 <= bw2d) {
+ t = 1 - d2 / bw2;
+ w = t * t;
}
return w;
}
-double **tricubic(int bw)
+double tricubic(double d2, double bw)
{
- int r, c;
- double bw3, bw3d, d, **w, t;
+ double bw3, bw3d, w, d3, 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;
- }
- }
+ d3 = sqrt(d2);
+ d3 = d3 * d3 * d3;
+ w = 0;
+
+ if (d3 <= bw3d) {
+ t = 1 - d3 / bw3;
+ w = t * t * t;
}
return w;
}
-double **gauss(int bw)
+double gauss(double d2, double bw)
{
- int r, c;
- double bw2, d, **w;
+ double bw2, 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;
+ w = 0;
- if (d <= bw2) {
- w[r + bw][c + bw] = exp(vf * d / bw2);
- }
- }
+ if (d2 <= bw2) {
+ w = exp(vf * d2 / bw2);
}
return w;
@@ -124,3 +97,24 @@
else
G_fatal_error(_("Invalid kernel option '%s'"), name);
}
+
+double **calc_weights(int bw)
+{
+ int r, c, count;
+ double d2, **w;
+
+ w = G_alloc_matrix(bw * 2 + 1, bw * 2 + 1);
+
+ count = 0;
+ for (r = -bw; r <= bw; r++) {
+ for (c = -bw; c <= bw; c++) {
+ d2 = r * r + c * c;
+ w[r + bw][c + bw] = w_fn(d2, bw);
+ if (w[r + bw][c + bw] > 0)
+ count++;
+ }
+ }
+ G_verbose_message(_("%d cells for bandwidth %d"), count, bw);
+
+ return w;
+}
More information about the grass-commit
mailing list