[GRASS-SVN] r67885 - grass-addons/grass7/raster/r.change.info
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Feb 18 13:38:31 PST 2016
Author: mmetz
Date: 2016-02-18 13:38:31 -0800 (Thu, 18 Feb 2016)
New Revision: 67885
Added:
grass-addons/grass7/raster/r.change.info/chisq.c
grass-addons/grass7/raster/r.change.info/dist.c
grass-addons/grass7/raster/r.change.info/gain.c
grass-addons/grass7/raster/r.change.info/gini.c
grass-addons/grass7/raster/r.change.info/r.change.info.html
grass-addons/grass7/raster/r.change.info/ratio.c
grass-addons/grass7/raster/r.change.info/window.h
Removed:
grass-addons/grass7/raster/r.change.info/divr_cats.c
grass-addons/grass7/raster/r.change.info/intr_cats.c
grass-addons/grass7/raster/r.change.info/r.neighbors.html
grass-addons/grass7/raster/r.change.info/readweights.c
Modified:
grass-addons/grass7/raster/r.change.info/Makefile
grass-addons/grass7/raster/r.change.info/bufs.c
grass-addons/grass7/raster/r.change.info/gather.c
grass-addons/grass7/raster/r.change.info/local_proto.h
grass-addons/grass7/raster/r.change.info/main.c
grass-addons/grass7/raster/r.change.info/ncb.h
grass-addons/grass7/raster/r.change.info/readcell.c
Log:
add r.change.info
Modified: grass-addons/grass7/raster/r.change.info/Makefile
===================================================================
--- grass-addons/grass7/raster/r.change.info/Makefile 2016-02-18 21:32:17 UTC (rev 67884)
+++ grass-addons/grass7/raster/r.change.info/Makefile 2016-02-18 21:38:31 UTC (rev 67885)
@@ -1,6 +1,6 @@
MODULE_TOPDIR = ../..
-PGM = r.neighbors
+PGM = r.change.info
LIBES = $(STATSLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB)
DEPENDENCIES = $(STATSDEP) $(RASTERDEP) $(GISDEP)
Modified: grass-addons/grass7/raster/r.change.info/bufs.c
===================================================================
--- grass-addons/grass7/raster/r.change.info/bufs.c 2016-02-18 21:32:17 UTC (rev 67884)
+++ grass-addons/grass7/raster/r.change.info/bufs.c 2016-02-18 21:38:31 UTC (rev 67885)
@@ -12,15 +12,19 @@
int allocate_bufs(void)
{
- int i;
- int bufsize;
+ int i, j;
+ int ncols;
+ size_t bufsize;
- bufsize = (Rast_window_cols() + 2 * ncb.dist) * sizeof(DCELL);
+ ncols = Rast_input_window_cols();
+ bufsize = ncols * sizeof(CELL);
- ncb.buf = (DCELL **) G_malloc(ncb.nsize * sizeof(DCELL *));
- for (i = 0; i < ncb.nsize; i++) {
- ncb.buf[i] = (DCELL *) G_malloc(bufsize);
- Rast_set_d_null_value(ncb.buf[i], Rast_window_cols() + 2 * ncb.dist);
+ for (i = 0; i < ncb.nin; i++) {
+ ncb.in[i].buf = (CELL **) G_malloc(ncb.nsize * sizeof(CELL *));
+ for (j = 0; j < ncb.nsize; j++) {
+ ncb.in[i].buf[j] = (CELL *) G_malloc(bufsize);
+ Rast_set_c_null_value(ncb.in[i].buf[j], ncols);
+ }
}
return 0;
@@ -28,15 +32,17 @@
int rotate_bufs(void)
{
- DCELL *temp;
- int i;
+ CELL *temp;
+ int i, j;
- temp = ncb.buf[0];
+ for (i = 0; i < ncb.nin; i++) {
+ temp = ncb.in[i].buf[0];
- for (i = 1; i < ncb.nsize; i++)
- ncb.buf[i - 1] = ncb.buf[i];
+ for (j = 1; j < ncb.nsize; j++)
+ ncb.in[i].buf[j - 1] = ncb.in[i].buf[j];
- ncb.buf[ncb.nsize - 1] = temp;
+ ncb.in[i].buf[ncb.nsize - 1] = temp;
+ }
return 0;
}
Added: grass-addons/grass7/raster/r.change.info/chisq.c
===================================================================
--- grass-addons/grass7/raster/r.change.info/chisq.c (rev 0)
+++ grass-addons/grass7/raster/r.change.info/chisq.c 2016-02-18 21:38:31 UTC (rev 67885)
@@ -0,0 +1,65 @@
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include "ncb.h"
+#include "window.h"
+
+/* CHI-square statistics to assess the relevance of an attribute
+ * from Quinlain 1986
+ * expected frequency for each dist bin and year:
+ * pexp[i] = p * (ncells[i] / ncellstot)
+ * then
+ * sum of all (p[i] - pexp[i])^2 / pexp[i]
+ *
+ * dist: array of distributions
+ * nd: number of cells in each distribution
+ * n: number of distributions
+ * dsize: distributions' size */
+
+DCELL chisq(double **dist, int *nd, int n, int dsize)
+{
+ int i, d, n_comb;
+ double p_sum, p_exp, result;
+
+ result = 0;
+ n_comb = 0;
+ for (d = 0; d < n; d++) {
+ n_comb += nd[d];
+ }
+ for (i = 0; i < dsize; i++) {
+
+ /* create combined distribution */
+ p_sum = 0;
+ for (d = 0; d < n; d++) {
+ if (dist[d][i] > 0 && nd[d] > 0) {
+ p_sum += dist[d][i];
+ }
+ }
+ /* squared and normalized difference to combined distribution */
+ if (p_sum > 0) {
+ for (d = 0; d < n; d++) {
+ if (nd[d] > 0) {
+ p_exp = p_sum * nd[d] / n_comb;
+ result += (dist[d][i] - p_exp) * (dist[d][i] - p_exp) / p_exp;
+ }
+ }
+ }
+ }
+
+ return result;
+}
+
+DCELL chisq1(void)
+{
+ return chisq(ci.dt, ci.n, ncb.nin, ci.ntypes);
+}
+
+DCELL chisq2(void)
+{
+ return chisq(ci.ds, ci.n, ncb.nin, ci.nsizebins);
+}
+
+DCELL chisq3(void)
+{
+ return chisq(ci.dts, ci.n, ncb.nin, ci.dts_size);
+}
Added: grass-addons/grass7/raster/r.change.info/dist.c
===================================================================
--- grass-addons/grass7/raster/r.change.info/dist.c (rev 0)
+++ grass-addons/grass7/raster/r.change.info/dist.c 2016-02-18 21:38:31 UTC (rev 67885)
@@ -0,0 +1,67 @@
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include "ncb.h"
+#include "window.h"
+
+/* average statistical distance between average distribution and
+ * observed distributions:
+ * sum of the absolute differences divided by the number of distributions
+ *
+ * dist: array of distributions
+ * nd: number of cells in each distribution
+ * n: number of distributions
+ * dsize: distributions' size */
+
+DCELL dist(double **dist, int *nd, int n, int dsize)
+{
+ int i, d, n_comb;
+ double dp;
+ double p_avg, d_dist;
+
+ d_dist = 0;
+ n_comb = 0;
+ for (d = 0; d < n; d++) {
+ n_comb += nd[d];
+ }
+ for (i = 0; i < dsize; i++) {
+
+ /* create combined distribution */
+ p_avg = 0;
+ for (d = 0; d < n; d++) {
+ if (dist[d][i] > 0 && nd[d] > 0) {
+ p_avg += dist[d][i];
+ }
+ }
+ /* difference to combined distribution */
+ if (p_avg > 0) {
+ p_avg /= n_comb;
+ for (d = 0; d < n; d++) {
+ if (nd[d] > 0) {
+ dp = p_avg - dist[d][i] / nd[d];
+ d_dist += fabs(dp);
+ }
+ }
+ }
+ }
+
+ if (d_dist == 0)
+ return 0;
+
+ return d_dist / (2 * (n - 1));
+}
+
+DCELL dist1(void)
+{
+ return dist(ci.dt, ci.n, ncb.nin, ci.ntypes);
+}
+
+DCELL dist2(void)
+{
+ return dist(ci.ds, ci.n, ncb.nin, ci.nsizebins);
+}
+
+DCELL dist3(void)
+{
+ return dist(ci.dts, ci.n, ncb.nin, ci.dts_size);
+}
Deleted: grass-addons/grass7/raster/r.change.info/divr_cats.c
===================================================================
--- grass-addons/grass7/raster/r.change.info/divr_cats.c 2016-02-18 21:32:17 UTC (rev 67884)
+++ grass-addons/grass7/raster/r.change.info/divr_cats.c 2016-02-18 21:38:31 UTC (rev 67885)
@@ -1,10 +0,0 @@
-#include <grass/gis.h>
-#include <grass/raster.h>
-#include "ncb.h"
-int divr_cats(void)
-{
- Rast_set_cats_fmt("$1 $?different categories$category$", 1.0, 0.0, 0.0, 0.0,
- &ncb.cats);
-
- return 0;
-}
Added: grass-addons/grass7/raster/r.change.info/gain.c
===================================================================
--- grass-addons/grass7/raster/r.change.info/gain.c (rev 0)
+++ grass-addons/grass7/raster/r.change.info/gain.c 2016-02-18 21:38:31 UTC (rev 67885)
@@ -0,0 +1,90 @@
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include "ncb.h"
+#include "window.h"
+
+
+/* information gain = amount of change
+ * gain1: gain for patch type frequencies
+ * gain2: gain for clump size class frequencies
+ * gain3: gain for patch type x clump size class frequencies
+ * pc: proportion of changes =
+ * number of transitions / ( n unmasked cells x n inputs ) */
+
+
+/* gain for any number of equally sized distributions using log2
+ * thus the upper bound is log2(n)
+ * dist: array of distributions
+ * nd: number of cells in each distribution
+ * n: number of distributions
+ * dsize: distributions' size
+ * h: entropy for each distribution */
+DCELL gain(double **dist, int *nd, int n, int dsize, double *h)
+{
+ int i, d, n_comb;
+ double p, pe, p_comb, H_comb, H_avg;
+
+ H_comb = 0;
+ n_comb = 0;
+ for (d = 0; d < n; d++) {
+ h[d] = 0;
+ n_comb += nd[d];
+ }
+ for (i = 0; i < dsize; i++) {
+ p_comb = 0;
+ for (d = 0; d < n; d++) {
+ if (dist[d][i] > 0 && nd[d] > 0) {
+ p = dist[d][i] / nd[d];
+ p_comb += dist[d][i];
+ pe = entropy_p(p);
+ h[d] += pe;
+ }
+ }
+ if (p_comb > 0) {
+ H_comb += entropy_p(p_comb / n_comb);
+ }
+ }
+ H_avg = 0;
+ for (d = 0; d < n; d++) {
+ H_avg += entropy(h[d]) * nd[d];
+ }
+ H_avg /= n_comb;
+ H_comb = entropy(H_comb);
+
+ if (H_comb < 0)
+ H_comb = 0;
+ if (H_avg < 0)
+ H_avg = 0;
+
+ if (H_comb < H_avg)
+ return 0;
+
+ if (H_comb == 0)
+ return 0;
+
+ return (H_comb - H_avg);
+}
+
+
+DCELL pc(void)
+{
+ /* proportion of changes
+ * theoretical max: ncb.n * (ncb.nin - 1) */
+ return (double) ci.nchanges / (ncb.n * (ncb.nin - 1));
+}
+
+DCELL gain1(void)
+{
+ return gain(ci.dt, ci.n, ncb.nin, ci.ntypes, ci.ht);
+}
+
+DCELL gain2(void)
+{
+ return gain(ci.ds, ci.n, ncb.nin, ci.nsizebins, ci.hs);
+}
+
+DCELL gain3(void)
+{
+ return gain(ci.dts, ci.n, ncb.nin, ci.dts_size, ci.hts);
+}
Modified: grass-addons/grass7/raster/r.change.info/gather.c
===================================================================
--- grass-addons/grass7/raster/r.change.info/gather.c 2016-02-18 21:32:17 UTC (rev 67884)
+++ grass-addons/grass7/raster/r.change.info/gather.c 2016-02-18 21:38:31 UTC (rev 67885)
@@ -1,18 +1,21 @@
+#include <math.h>
#include <grass/gis.h>
#include <grass/raster.h>
#include "ncb.h"
+#include "window.h"
-/*
- given the starting col of the neighborhood,
- copy the cell values from the bufs into the array of values
- and return the number of values copied.
- */
#define sqr(x) ((x) * (x))
+double (*entropy) (double);
+double (*entropy_p) (double);
+
+static double alpha = 1.;
+
void circle_mask(void)
{
int i, j;
+ double dist, distsq;
if (ncb.mask)
return;
@@ -22,67 +25,282 @@
for (i = 0; i < ncb.nsize; i++)
ncb.mask[i] = G_malloc(ncb.nsize);
- for (i = 0; i < ncb.nsize; i++)
- for (j = 0; j < ncb.nsize; j++)
+ dist = (ncb.nsize - 1.0) / 2.0;
+ distsq = sqr(dist);
+ if (ncb.nsize % 2 == 0) {
+ /* nsize is even */
+ distsq += 0.25;
+ }
+
+ ncb.n = 0;
+ for (i = 0; i < ncb.nsize; i++) {
+ for (j = 0; j < ncb.nsize; j++) {
ncb.mask[i][j] =
- sqr(i - ncb.dist) + sqr(j - ncb.dist) <= sqr(ncb.dist);
+ sqr(i - dist) + sqr(j - dist) <= distsq;
+ ncb.n += ncb.mask[i][j];
+ }
+ }
}
-void weights_mask(void)
+/* argument of the sum for general entropy */
+double eai(double p)
{
- int i, j;
+ if (p)
+ return pow(p, alpha);
+ return 0;
+}
- if (ncb.mask)
- return;
+/* argument of the sum for Shannon entropy */
+double shi(double p)
+{
+ if (p)
+ return -p * log2(p);
+ return 0;
+}
- ncb.mask = G_malloc(ncb.nsize * sizeof(char *));
+/* general entropy for alpha */
+double eah(double h)
+{
+ if (h)
+ return log2(h) / (1 - alpha);
+ return 0;
+}
- for (i = 0; i < ncb.nsize; i++)
- ncb.mask[i] = G_malloc(ncb.nsize);
+/* Shannon entropy */
+double shh(double h)
+{
+ return h;
+}
- for (i = 0; i < ncb.nsize; i++)
- for (j = 0; j < ncb.nsize; j++)
- ncb.mask[i][j] = ncb.weights[i][j] != 0;
+int set_alpha(double a)
+{
+ if (a == 1) {
+ entropy_p = shi;
+ entropy = shh;
+ }
+ else {
+ alpha = a;
+ entropy_p = eai;
+ entropy = eah;
+ }
+
+ return (a == 1);
}
-int gather(DCELL *values, int offset)
+
+int types_differ(CELL a, CELL b)
{
+ if (!Rast_is_c_null_value(&a) && !Rast_is_c_null_value(&b))
+ return (a != b);
+
+ if (Rast_is_c_null_value(&a) && Rast_is_c_null_value(&b))
+ return 0;
+
+ /* either a or b is NULL */
+ return 1;
+
+ if ((Rast_is_c_null_value(&a) && !Rast_is_c_null_value(&b)) ||
+ (!Rast_is_c_null_value(&a) && Rast_is_c_null_value(&b)))
+ return 1;
+
+ return 0;
+}
+
+/* patch identification */
+int gather(int offset)
+{
int row, col;
+ int i, j;
+ int idx;
int n = 0;
+ int nchanges = 0;
+ int connected, old_pid, new_pid;
+ int *pid_curr, *pid_prev, *pid_tmp;
+ struct c_h *ch;
+
+ /* reset stats */
+ for (i = 0; i < ncb.nin; i++) {
+ ci.n[i] = 0;
+ ci.ht[i] = 0;
+ ci.hts[i] = 0;
- *values = 0;
+ for (j = 0; j < ci.ntypes; j++) {
+ ci.dt[i][j] = 0;
+ ci.dts[i][j] = 0;
+ }
+ for (j = ci.ntypes; j < ci.dts_size; j++)
+ ci.dts[i][j] = 0;
+ for (j = 0; j < ci.nsizebins; j++) {
+ ci.ds[i][j] = 0;
+ }
+
+ ch = &ci.ch[i];
+ ch->pid = 0;
+
+ Rast_set_c_null_value(&ch->up, 1);
+ Rast_set_c_null_value(&ch->left, 1);
+ Rast_set_c_null_value(&ch->curr, 1);
+
+ for (j = 0; j < ch->palloc; j++) {
+ ch->pst[j].size = 0;
+ ch->pst[j].type = 0;
+ }
+ ch->np = 0;
+
+ for (j = 0; j < ncb.nsize; j++) {
+ ch->pid_curr[j] = 0;
+ ch->pid_prev[j] = 0;
+ }
+ }
+
+ /* collect distributions */
+ n = 0;
for (row = 0; row < ncb.nsize; row++) {
+
+ for (i = 0; i < ncb.nin; i++) {
+ ch = &ci.ch[i];
+
+ Rast_set_c_null_value(&ch->left, 1);
+
+ pid_tmp = ch->pid_prev;
+ ch->pid_prev = ch->pid_curr;
+ ch->pid_curr = pid_tmp;
+ }
+
for (col = 0; col < ncb.nsize; col++) {
if (ncb.mask && !ncb.mask[row][col])
continue;
- values[n] = ncb.buf[row][offset + col];
+ for (i = 0; i < ncb.nin; i++) {
- n++;
- }
- }
+ ch = &ci.ch[i];
+ ch->pid_curr[col] = 0;
- return n;
-}
+ ch->curr = ncb.in[i].buf[row][offset + col];
+ /* number of changes */
+ if (i > 0)
+ nchanges += types_differ(ch->curr, ci.ch[i - 1].curr);
+
+ if (Rast_is_c_null_value(&ch->curr)) {
+ ch->left = ch->curr;
+ continue;
+ }
+ Rast_set_c_null_value(&ch->up, 1);
+ if (row > 0) {
+ if (ncb.mask) {
+ if (ncb.mask[row - 1][col])
+ ch->up = ncb.in[i].buf[row - 1][offset + col];
+ }
+ else {
+ ch->up = ncb.in[i].buf[row - 1][offset + col];
+ }
+ }
+
+ /* type count */
+ ci.dt[i][ch->curr - ci.tmin] += 1;
+
+ /* trace patch clumps */
+ pid_curr = ch->pid_curr;
+ pid_prev = ch->pid_prev;
-int gather_w(DCELL *values, DCELL (*values_w)[2], int offset)
-{
- int row, col;
- int n = 0;
+ connected = 0;
+ if (!Rast_is_c_null_value(&ch->left) && ch->curr == ch->left) {
+ pid_curr[col] = pid_curr[col - 1];
+ connected = 1;
+ ch->pst[pid_curr[col]].size++;
+ }
- values_w[0][0] = 0;
- values_w[0][1] = 1;
+ if (!Rast_is_c_null_value(&ch->up) && ch->curr == ch->up) {
- for (row = 0; row < ncb.nsize; row++) {
- for (col = 0; col < ncb.nsize; col++) {
- values[n] = values_w[n][0] = ncb.buf[row][offset + col];
- values_w[n][1] = ncb.weights[row][col];
+ if (pid_curr[col] != pid_prev[col]) {
- n++;
+ if (connected) {
+ ch->np--;
+
+ if (ch->np == 0) {
+ G_fatal_error("npatch == 0 at row %d, col %d", row, col);
+ }
+ }
+
+ old_pid = pid_curr[col];
+ new_pid = pid_prev[col];
+ pid_curr[col] = new_pid;
+ if (old_pid > 0) {
+ /* merge */
+ /* update left side of the current row */
+ for (j = 0; j < col; j++) {
+ if (pid_curr[j] == old_pid)
+ pid_curr[j] = new_pid;
+ }
+ /* update right side of the previous row */
+ for (j = col + 1; j < ncb.nsize; j++) {
+ if (pid_prev[j] == old_pid)
+ pid_prev[j] = new_pid;
+ }
+ ch->pst[new_pid].size += ch->pst[old_pid].size;
+ ch->pst[old_pid].size = 0;
+
+ if (old_pid == ch->pid)
+ ch->pid--;
+ }
+ else {
+ ch->pst[new_pid].size++;
+ }
+ }
+ connected = 1;
+ }
+ if (!connected) {
+ /* start new patch */
+ ch->np++;
+ ch->pid++;
+ pid_curr[col] = ch->pid;
+
+ if (ch->pid >= ch->palloc) {
+ ch->pst = (struct pst *)G_realloc(ch->pst, (ch->pid + 10) * sizeof(struct pst));
+
+ for (j = ch->palloc; j < ch->pid + 10; j++)
+ ch->pst[j].size = 0;
+
+ ch->palloc = ch->pid + 10;
+ }
+
+ ch->pst[ch->pid].size = 1;
+ ch->pst[ch->pid].type = ch->curr;
+ }
+
+ ch->left = ch->curr;
+ ci.n[i]++;
+ n++;
+ }
}
}
+
+ if (n == 0) {
+ return 0;
+ }
+ /* convert patches to distribution of types and size classes */
+ for (i = 0; i < ncb.nin; i++) {
+ int nsum;
+
+ ch = &ci.ch[i];
+ nsum = 0;
+ for (j = 0; j <= ch->pid; j++) {
+ if (ch->pst[j].size > 0) {
+ frexp(ch->pst[j].size, &idx);
+ ci.ds[i][idx - 1] += ch->pst[j].size;
+ idx = (ch->pst[j].type - ci.tmin) * ci.nsizebins + idx - 1;
+ ci.dts[i][idx] += ch->pst[j].size;
+ }
+ nsum += ch->pst[j].size;
+ }
+ if (nsum != ci.n[i])
+ G_fatal_error("patch sum is %d, should be %d", nsum, ci.n[i]);
+ }
+
+ ci.nchanges = nchanges;
+
return n;
}
Added: grass-addons/grass7/raster/r.change.info/gini.c
===================================================================
--- grass-addons/grass7/raster/r.change.info/gini.c (rev 0)
+++ grass-addons/grass7/raster/r.change.info/gini.c 2016-02-18 21:38:31 UTC (rev 67885)
@@ -0,0 +1,75 @@
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include "ncb.h"
+#include "window.h"
+
+/* Gini index = 1 - Simpson
+ *
+ * dist: array of distributions
+ * nd: number of cells in each distribution
+ * n: number of distributions
+ * dsize: distributions' size */
+
+DCELL gini(double **dist, int *nd, int n, int dsize)
+{
+ int i, d, n_comb;
+ double p, p_avg, gini, gini_avg, *gini_d;
+
+ gini_d = G_malloc(n * sizeof(double));
+ gini = 0;
+ gini_avg = 0;
+ n_comb = 0;
+ for (d = 0; d < n; d++) {
+ n_comb += nd[d];
+ gini_d[d] = 1;
+ }
+ for (i = 0; i < dsize; i++) {
+
+ /* comnined distribution */
+ p_avg = 0;
+ for (d = 0; d < n; d++) {
+ if (dist[d][i] > 0 && nd[d] > 0) {
+ p_avg += dist[d][i];
+ p = dist[d][i] / nd[d];
+ gini_d[d] -= p * p;
+ }
+ }
+
+ /* gini of combined distribution */
+ if (p_avg > 0) {
+ p_avg /= n_comb;
+ gini -= p_avg * p_avg;
+
+ }
+ }
+ /* average gini of single distributions */
+ gini_avg = 0;
+ for (d = 0; d < n; d++) {
+ gini_avg += ((double) nd[d] / n_comb) * gini_d[d];
+ }
+ G_free(gini_d);
+
+ if (gini == 0)
+ return 0;
+
+ return gini - gini_avg;
+
+ /* Gini ratio: the max decrease is (n - 1) / n */
+ return gini_avg * n / (n - 1);
+}
+
+DCELL gini1(void)
+{
+ return gini(ci.dt, ci.n, ncb.nin, ci.ntypes);
+}
+
+DCELL gini2(void)
+{
+ return gini(ci.ds, ci.n, ncb.nin, ci.nsizebins);
+}
+
+DCELL gini3(void)
+{
+ return gini(ci.dts, ci.n, ncb.nin, ci.dts_size);
+}
Deleted: grass-addons/grass7/raster/r.change.info/intr_cats.c
===================================================================
--- grass-addons/grass7/raster/r.change.info/intr_cats.c 2016-02-18 21:32:17 UTC (rev 67884)
+++ grass-addons/grass7/raster/r.change.info/intr_cats.c 2016-02-18 21:38:31 UTC (rev 67885)
@@ -1,10 +0,0 @@
-#include <grass/gis.h>
-#include <grass/raster.h>
-#include "ncb.h"
-
-int intr_cats(void)
-{
- Rast_set_cats_fmt("$1% dispersion", 1.0, -1.0, 0.0, 0.0, &ncb.cats);
-
- return 0;
-}
Modified: grass-addons/grass7/raster/r.change.info/local_proto.h
===================================================================
--- grass-addons/grass7/raster/r.change.info/local_proto.h 2016-02-18 21:32:17 UTC (rev 67884)
+++ grass-addons/grass7/raster/r.change.info/local_proto.h 2016-02-18 21:38:31 UTC (rev 67885)
@@ -2,24 +2,42 @@
extern int allocate_bufs(void);
extern int rotate_bufs(void);
-/* gather */
+/* mask */
extern void circle_mask(void);
-extern void weights_mask(void);
-extern int gather(DCELL *, int);
-extern int gather_w(DCELL *, DCELL(*)[2], int);
+/* gather */
+int gather(int);
+int set_alpha(double);
+double eai(double);
+double eah(double);
+double shi(double);
+double shh(double);
+
/* readcell.c */
-extern int readcell(int, int, int, int);
+extern int readcell(int, int, int);
-/* divr_cats.c */
-extern int divr_cats(void);
+/* gain.c */
+DCELL pc(void);
+DCELL gain1(void);
+DCELL gain2(void);
+DCELL gain3(void);
-/* intr_cats.c */
-extern int intr_cats(void);
+/* ratio.c */
+DCELL ratio1(void);
+DCELL ratio2(void);
+DCELL ratio3(void);
-/* null_cats.c */
-extern int null_cats(const char *);
+/* gini.c */
+DCELL gini1(void);
+DCELL gini2(void);
+DCELL gini3(void);
-/* read_weights.c */
-extern void read_weights(const char *);
-extern void gaussian_weights(double);
+/* dist.c */
+DCELL dist1(void);
+DCELL dist2(void);
+DCELL dist3(void);
+
+/* chisq.c */
+DCELL chisq1(void);
+DCELL chisq2(void);
+DCELL chisq3(void);
Modified: grass-addons/grass7/raster/r.change.info/main.c
===================================================================
--- grass-addons/grass7/raster/r.change.info/main.c 2016-02-18 21:32:17 UTC (rev 67884)
+++ grass-addons/grass7/raster/r.change.info/main.c 2016-02-18 21:38:31 UTC (rev 67885)
@@ -1,17 +1,13 @@
/****************************************************************************
*
- * MODULE: r.neighbors
- * AUTHOR(S): Michael Shapiro, CERL (original contributor)
- * Markus Neteler <neteler itc.it>, Bob Covill <bcovill tekmap.ns.ca>,
- * Brad Douglas <rez touchofmadness.com>, Glynn Clements <glynn gclements.plus.com>,
- * Jachym Cepicky <jachym les-ejk.cz>, Jan-Oliver Wagner <jan intevation.de>,
- * Radim Blazek <radim.blazek gmail.com>
+ * MODULE: r.change.info
+ * AUTHOR(S): Markus Metz
+ * based on r.neighbors
*
- * PURPOSE: Makes each cell category value a function of the category values
- * assigned to the cells around it, and stores new cell values in an
- * output raster map layer
- * COPYRIGHT: (C) 1999-2006 by the GRASS Development Team
+ * PURPOSE: Change assessment for categorical raster series
+ *
+ * COPYRIGHT: (C) 2014 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
@@ -21,74 +17,57 @@
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
+#include <math.h>
#include <grass/gis.h>
#include <grass/raster.h>
#include <grass/glocale.h>
-#include <grass/stats.h>
#include "ncb.h"
+#include "window.h"
#include "local_proto.h"
-typedef int (*ifunc) (void);
+typedef DCELL dfunc(void);
struct menu
{
- stat_func *method; /* routine to compute new value */
- stat_func_w *method_w; /* routine to compute new value (weighted) */
- ifunc cat_names; /* routine to make category names */
- int copycolr; /* flag if color table can be copied */
- int half; /* whether to add 0.5 to result (redundant) */
- int otype; /* output type */
+ dfunc *method; /* routine to compute new value */
char *name; /* method name */
char *text; /* menu display - full description */
};
-enum out_type {
- T_FLOAT = 1,
- T_INT = 2,
- T_COUNT = 3,
- T_COPY = 4,
- T_SUM = 5
-};
-
#define NO_CATS 0
/* modify this table to add new methods */
static struct menu menu[] = {
- {c_ave, w_ave, NO_CATS, 1, 1, T_FLOAT, "average", "average value"},
- {c_median, w_median, NO_CATS, 1, 0, T_FLOAT, "median", "median value"},
- {c_mode, w_mode, NO_CATS, 1, 0, T_COPY, "mode", "most frequently occurring value"},
- {c_min, NULL, NO_CATS, 1, 0, T_COPY, "minimum", "lowest value"},
- {c_max, NULL, NO_CATS, 1, 0, T_COPY, "maximum", "highest value"},
- {c_range, NULL, NO_CATS, 1, 0, T_COPY, "range", "range value"},
- {c_stddev, w_stddev, NO_CATS, 0, 1, T_FLOAT, "stddev", "standard deviation"},
- {c_sum, w_sum, NO_CATS, 1, 0, T_SUM, "sum", "sum of values"},
- {c_count, w_count, NO_CATS, 0, 0, T_COUNT, "count", "count of non-NULL values"},
- {c_var, w_var, NO_CATS, 0, 1, T_FLOAT, "variance", "statistical variance"},
- {c_divr, NULL, divr_cats, 0, 0, T_INT, "diversity",
- "number of different values"},
- {c_intr, NULL, intr_cats, 0, 0, T_INT, "interspersion",
- "number of values different than center value"},
- {c_quart1, w_quart1, NO_CATS, 1, 0, T_FLOAT, "quart1", "first quartile"},
- {c_quart3, w_quart3, NO_CATS, 1, 0, T_FLOAT, "quart3", "third quartile"},
- {c_perc90, w_perc90, NO_CATS, 1, 0, T_FLOAT, "perc90", "ninetieth percentile"},
- {c_quant, w_quant, NO_CATS, 1, 0, T_FLOAT, "quantile", "arbitrary quantile"},
- {0, 0, 0, 0, 0, 0, 0, 0}
+ {pc, "pc", "proportion of changes"},
+ {gain1, "gain1", "Information gain for category distributions"},
+ {gain2, "gain2", "Information gain for size distributions"},
+ {gain3, "gain3", "Information gain for category and size distributions"},
+ {ratio1, "ratio1", "Information gain ratio for category distributions"},
+ {ratio2, "ratio2", "Information gain ratio for size distributions"},
+ {ratio3, "ratio3", "Information gain ratio for category and size distributions"},
+ {gini1, "gini1", "Gini impurity for category distributions"},
+ {gini2, "gini2", "Gini impurity for size distributions"},
+ {gini3, "gini3", "Gini impurity for category and size distributions"},
+ {dist1, "dist1", "Statistical distance for category distributions"},
+ {dist2, "dist2", "Statistical distance for size distributions"},
+ {dist3, "dist3", "Statistical distance for category and size distributions"},
+ {chisq1, "chisq1", "CHI-square for category distributions"},
+ {chisq2, "chisq2", "CHI-square for size distributions"},
+ {chisq3, "chisq3", "CHI-square for category and size distributions"},
+ {NULL, NULL, NULL}
};
struct ncb ncb;
+struct changeinfo ci;
+
struct output
{
const char *name;
char title[1024];
int fd;
DCELL *buf;
- stat_func *method_fn;
- stat_func_w *method_fn_w;
- int copycolr;
- ifunc cat_names;
- int map_type;
- double quantile;
+ dfunc *method_fn;
};
static int find_method(const char *method_name)
@@ -104,96 +83,72 @@
return -1;
}
-static RASTER_MAP_TYPE output_type(RASTER_MAP_TYPE input_type, int weighted, int mode)
-{
- switch (mode) {
- case T_FLOAT:
- return DCELL_TYPE;
- case T_INT:
- return CELL_TYPE;
- case T_COUNT:
- return weighted ? DCELL_TYPE : CELL_TYPE;
- case T_COPY:
- return input_type;
- case T_SUM:
- return weighted ? DCELL_TYPE : input_type;
- default:
- G_fatal_error(_("Invalid out_type enumeration: %d"), mode);
- return -1;
- }
-}
+int make_colors(struct Colors *colr, DCELL min, DCELL max);
+
int main(int argc, char *argv[])
{
- char *p;
- int in_fd;
- int selection_fd;
+ char *p, *pp;
+ size_t dlen;
int num_outputs;
struct output *outputs = NULL;
- int copycolr, weights, have_weights_mask;
- char *selection;
RASTER_MAP_TYPE map_type;
- int row, col;
+ int row, col, ocol, roff, coff;
+ int rspill, cspill;
+ int orows, ocols;
int readrow;
int nrows, ncols;
+ struct Range range;
+ struct FPRange drange;
+ DCELL dmin, dmax;
+ CELL min, max, imin, imax;
int i, n;
+ int step;
+ double alpha;
struct Colors colr;
struct Cell_head cellhd;
- struct Cell_head window;
+ struct Cell_head window, owind;
struct History history;
struct GModule *module;
struct
{
- struct Option *input, *output, *selection;
- struct Option *method, *size;
- struct Option *title;
- struct Option *weight;
- struct Option *gauss;
- struct Option *quantile;
+ struct Option *input, *output;
+ struct Option *method, *wsize, *step, *alpha;
} parm;
struct
{
struct Flag *align, *circle;
} flag;
- DCELL *values; /* list of neighborhood values */
- DCELL *values_tmp; /* list of neighborhood values */
- DCELL(*values_w)[2]; /* list of neighborhood values and weights */
- DCELL(*values_w_tmp)[2]; /* list of neighborhood values and weights */
G_gisinit(argv[0]);
module = G_define_module();
G_add_keyword(_("raster"));
- G_add_keyword(_("algebra"));
G_add_keyword(_("statistics"));
- G_add_keyword(_("aggregation"));
- G_add_keyword(_("neighbor"));
- G_add_keyword(_("focal statistics"));
- G_add_keyword(_("filter"));
+ G_add_keyword(_("change detection"));
+ G_add_keyword(_("landscape structure"));
module->description =
- _("Makes each cell category value a "
- "function of the category values assigned to the cells "
- "around it, and stores new cell values in an output raster "
- "map layer.");
+ _("Landscape change assessment");
- parm.input = G_define_standard_option(G_OPT_R_INPUT);
+ parm.input = G_define_standard_option(G_OPT_R_INPUTS);
- parm.selection = G_define_standard_option(G_OPT_R_INPUT);
- parm.selection->key = "selection";
- parm.selection->required = NO;
- parm.selection->description = _("Name of an input raster map to select the cells which should be processed");
-
parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.output->description = _("Name for output raster map(s)");
parm.output->multiple = YES;
parm.method = G_define_option();
parm.method->key = "method";
parm.method->type = TYPE_STRING;
parm.method->required = NO;
- parm.method->answer = "average";
- p = G_malloc(1024);
+ parm.method->answer = "ratio3";
+ dlen = 0;
for (n = 0; menu[n].name; n++) {
+ dlen += strlen(menu[n].name);
+ }
+ dlen += n;
+ p = G_malloc(dlen);
+ for (n = 0; menu[n].name; n++) {
if (n)
strcat(p, ",");
else
@@ -201,87 +156,169 @@
strcat(p, menu[n].name);
}
parm.method->options = p;
- parm.method->description = _("Neighborhood operation");
+ parm.method->description = _("Change assessment");
+ dlen = 0;
+ for (n = 0; menu[n].name; n++) {
+ dlen += strlen(menu[n].name);
+ dlen += strlen(menu[n].text);
+ }
+ dlen += n * 2;
+ pp = G_malloc(dlen);
+ for (n = 0; menu[n].name; n++) {
+ if (n)
+ strcat(pp, ";");
+ else
+ *pp = 0;
+ strcat(pp, menu[n].name);
+ strcat(pp, ";");
+ strcat(pp, menu[n].text);
+ }
+ parm.method->descriptions = pp;
parm.method->multiple = YES;
- parm.method->guisection = _("Neighborhood");
- parm.size = G_define_option();
- parm.size->key = "size";
- parm.size->type = TYPE_INTEGER;
- parm.size->required = NO;
- parm.size->description = _("Neighborhood size");
- parm.size->answer = "3";
- parm.size->guisection = _("Neighborhood");
+ parm.wsize = G_define_option();
+ parm.wsize->key = "size";
+ parm.wsize->type = TYPE_INTEGER;
+ parm.wsize->required = NO;
+ parm.wsize->description = _("Window size (cells)");
+ parm.wsize->answer = "40";
+ parm.wsize->guisection = _("Moving window");
- parm.title = G_define_option();
- parm.title->key = "title";
- parm.title->key_desc = "phrase";
- parm.title->type = TYPE_STRING;
- parm.title->required = NO;
- parm.title->description = _("Title for output raster map");
+ parm.step = G_define_option();
+ parm.step->key = "step";
+ parm.step->type = TYPE_INTEGER;
+ parm.step->required = NO;
+ parm.step->description = _("Processing step (cells)");
+ parm.step->answer = "40";
+ parm.step->guisection = _("Moving window");
- parm.weight = G_define_standard_option(G_OPT_F_INPUT);
- parm.weight->key = "weight";
- parm.weight->required = NO;
- parm.weight->description = _("Text file containing weights");
+ parm.alpha = G_define_option();
+ parm.alpha->key = "alpha";
+ parm.alpha->type = TYPE_DOUBLE;
+ parm.alpha->required = NO;
+ parm.alpha->label = _("Alpha for general entropy");
+ parm.alpha->description = _("Default = 1 for Shannon Entropy");
+ parm.alpha->answer = "1";
- parm.gauss = G_define_option();
- parm.gauss->key = "gauss";
- parm.gauss->type = TYPE_DOUBLE;
- parm.gauss->required = NO;
- parm.gauss->description = _("Sigma (in cells) for Gaussian filter");
-
- parm.quantile = G_define_option();
- parm.quantile->key = "quantile";
- parm.quantile->type = TYPE_DOUBLE;
- parm.quantile->required = NO;
- parm.quantile->multiple = YES;
- parm.quantile->description = _("Quantile to calculate for method=quantile");
- parm.quantile->options = "0.0-1.0";
-
flag.align = G_define_flag();
flag.align->key = 'a';
- flag.align->description = _("Do not align output with the input");
+ flag.align->description = _("Do not align input region with input maps");
flag.circle = G_define_flag();
flag.circle->key = 'c';
- flag.circle->description = _("Use circular neighborhood");
- flag.circle->guisection = _("Neighborhood");
+ flag.circle->description = _("Use circular mask");
+ flag.circle->guisection = _("Moving window");
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
- sscanf(parm.size->answer, "%d", &ncb.nsize);
- if (ncb.nsize <= 0)
- G_fatal_error(_("Neighborhood size must be positive"));
- if (ncb.nsize % 2 == 0)
- G_fatal_error(_("Neighborhood size must be odd"));
- ncb.dist = ncb.nsize / 2;
+ sscanf(parm.wsize->answer, "%d", &ncb.nsize);
+ if (ncb.nsize <= 1)
+ G_fatal_error(_("Window size must be > 1"));
+ ncb.n = ncb.nsize * ncb.nsize;
- if (parm.weight->answer && flag.circle->answer)
- G_fatal_error(_("-%c and %s= are mutually exclusive"),
- flag.circle->key, parm.weight->key);
+ sscanf(parm.step->answer, "%d", &step);
+ if (step <= 0)
+ G_fatal_error(_("Processing step must be positive"));
+
+ if (step > ncb.nsize)
+ G_fatal_error(_("Processing step can not be larger than window size"));
- if (parm.weight->answer && parm.gauss->answer)
- G_fatal_error(_("%s= and %s= are mutually exclusive"),
- parm.weight->key, parm.gauss->key);
+ sscanf(parm.alpha->answer, "%lf", &alpha);
+ if (alpha <= 0)
+ G_fatal_error(_("Alpha for general entropy must be positive"));
+ set_alpha(alpha);
- ncb.oldcell = parm.input->answer;
+ for (i = 0; parm.input->answers[i]; i++)
+ ;
+ ncb.nin = i;
+ if (ncb.nin < 2)
+ G_fatal_error(_("At least two input maps are required"));
+
+ ncb.in = G_malloc(ncb.nin * sizeof(struct input));
+
+ for (i = 0; i < ncb.nin; i++)
+ ncb.in[i].name = parm.input->answers[i];
+ Rast_get_cellhd(ncb.in[0].name, "", &cellhd);
+ G_get_window(&window);
if (!flag.align->answer) {
- Rast_get_cellhd(ncb.oldcell, "", &cellhd);
- G_get_window(&window);
Rast_align_window(&window, &cellhd);
- Rast_set_window(&window);
}
+ if (cellhd.rows < ncb.nsize)
+ G_fatal_error(_("The current region is too small, it must have at least %d rows"), ncb.nsize);
+ nrows = cellhd.rows;
+ if (cellhd.cols < ncb.nsize)
+ G_fatal_error(_("The current region is too small, it must have at least %d cols"), ncb.nsize);
+ ncols = cellhd.cols;
- nrows = Rast_window_rows();
- ncols = Rast_window_cols();
+ /* adjust the input window */
+ roff = coff = 0;
- /* open raster maps */
- in_fd = Rast_open_old(ncb.oldcell, "");
- map_type = Rast_get_map_type(in_fd);
+ /* number of steps fitting into input rows */
+ orows = (nrows - ncb.nsize) / step + 1;
+ /* number of steps fitting into input cols */
+ ocols = (ncols - ncb.nsize) / step + 1;
+ rspill = nrows - ncb.nsize - (orows - 1) * step;
+ if (rspill) {
+ /* shrink input window */
+ roff = rspill / 2;
+ if (roff) {
+ /* shift input north */
+ window.north -= roff * window.ns_res;
+ }
+ /* shift input south */
+ window.south += (rspill - roff) * window.ns_res;
+ window.rows -= rspill;
+ }
+ cspill = ncols - ncb.nsize - (ocols - 1) * step;
+ if (cspill) {
+ /* shrink input window */
+ coff = cspill / 2;
+ if (coff) {
+ /* shift input west */
+ window.west += coff * window.ew_res;
+ }
+ /* shift input east */
+ window.east -= (cspill - coff) * window.ew_res;
+ window.cols -= cspill;
+ }
+ /* not needed */
+ G_adjust_Cell_head(&window, 0, 0);
+
+ Rast_set_input_window(&window);
+
+ nrows = Rast_input_window_rows();
+ ncols = Rast_input_window_cols();
+
+ /* construct the output window */
+ owind = window;
+ owind.rows = orows;
+ owind.cols = ocols;
+ owind.ns_res = step * window.ns_res;
+ owind.ew_res = step * window.ew_res;
+
+ if (ncb.nsize != step) {
+ owind.north = owind.north - ((ncb.nsize - step) / 2.0) * window.ns_res;
+ owind.south = owind.south + ((ncb.nsize - step) / 2.0) * window.ns_res;
+
+ owind.east = owind.east - ((ncb.nsize - step) / 2.0) * window.ew_res;
+ owind.west = owind.west + ((ncb.nsize - step) / 2.0) * window.ew_res;
+ }
+ /* needed */
+ G_adjust_Cell_head(&owind, 1, 1);
+
+ Rast_set_output_window(&owind);
+
+ /* open input raster maps */
+ for (i = 0; i < ncb.nin; i++) {
+ ncb.in[i].fd = Rast_open_old(ncb.in[i].name, "");
+ map_type = Rast_get_map_type(ncb.in[i].fd);
+ if (map_type != CELL_TYPE)
+ G_warning(_("Input raster <%s> is not of type CELL"), ncb.in[i].name);
+ }
+
/* process the output maps */
for (i = 0; parm.output->answers[i]; i++)
;
@@ -290,159 +327,108 @@
for (i = 0; parm.method->answers[i]; i++)
;
if (num_outputs != i)
- G_fatal_error(_("%s= and %s= must have the same number of values"),
- parm.output->key, parm.method->key);
+ G_fatal_error(_("output= and method= must have the same number of values"));
outputs = G_calloc(num_outputs, sizeof(struct output));
- /* read the weights */
- weights = 0;
- ncb.weights = NULL;
ncb.mask = NULL;
- if (parm.weight->answer) {
- read_weights(parm.weight->answer);
- weights = 1;
- }
- else if (parm.gauss->answer) {
- gaussian_weights(atof(parm.gauss->answer));
- weights = 1;
- }
-
- copycolr = 0;
- have_weights_mask = 0;
for (i = 0; i < num_outputs; i++) {
struct output *out = &outputs[i];
const char *output_name = parm.output->answers[i];
const char *method_name = parm.method->answers[i];
int method = find_method(method_name);
- RASTER_MAP_TYPE otype = output_type(map_type, weights, menu[method].otype);
out->name = output_name;
- if (weights) {
- if (menu[method].method_w) {
- out->method_fn = NULL;
- out->method_fn_w = menu[method].method_w;
- }
- else {
- if (parm.weight->answer) {
- G_warning(_("Method %s not compatible with weighing window, using weight mask instead"),
- method_name);
- if (!have_weights_mask) {
- weights_mask();
- have_weights_mask = 1;
- }
- }
- else if (parm.gauss->answer) {
- G_warning(_("Method %s not compatible with Gaussian filter, using unweighed version instead"),
- method_name);
- }
-
- out->method_fn = menu[method].method;
- out->method_fn_w = NULL;
- }
- }
- else {
- out->method_fn = menu[method].method;
- out->method_fn_w = NULL;
- }
- out->copycolr = menu[method].copycolr;
- out->cat_names = menu[method].cat_names;
- if (out->copycolr)
- copycolr = 1;
- out->quantile = (parm.quantile->answer && parm.quantile->answers[i])
- ? atof(parm.quantile->answers[i])
- : 0;
- out->buf = Rast_allocate_d_buf();
- out->fd = Rast_open_new(output_name, otype);
- /* TODO: method=mode should propagate its type */
+ out->method_fn = menu[method].method;
- /* get title, initialize the category and stat info */
- if (parm.title->answer)
- strcpy(out->title, parm.title->answer);
- else
- sprintf(out->title, "%dx%d neighborhood: %s of %s",
- ncb.nsize, ncb.nsize, menu[method].name, ncb.oldcell);
- }
+ out->buf = Rast_allocate_d_output_buf();
+ out->fd = Rast_open_new(output_name, DCELL_TYPE);
- /* copy color table? */
- if (copycolr) {
- G_suppress_warnings(1);
- copycolr =
- (Rast_read_colors(ncb.oldcell, "", &colr) > 0);
- G_suppress_warnings(0);
+ sprintf(out->title, "%s, %dx%d window, step %d",
+ menu[method].text, ncb.nsize, ncb.nsize, step);
}
- /* allocate the cell buffers */
- allocate_bufs();
+ if (flag.circle->answer)
+ circle_mask();
- /* initialize the cell bufs with 'dist' rows of the old cellfile */
- readrow = 0;
- for (row = 0; row < ncb.dist; row++)
- readcell(in_fd, readrow++, nrows, ncols);
+ /* initialize change info */
+ frexp(ncb.n, &ci.nsizebins);
+ G_debug(1, "n cells: %d, size * size: %d", ncb.n, ncb.nsize * ncb.nsize);
+ G_debug(1, "nsizebins: %d", ci.nsizebins);
- /* open the selection raster map */
- if (parm.selection->answer) {
- G_message(_("Opening selection map <%s>"), parm.selection->answer);
- selection_fd = Rast_open_old(parm.selection->answer, "");
- selection = Rast_allocate_null_buf();
- } else {
- selection_fd = -1;
- selection = NULL;
+ /* max number of different types */
+ Rast_init_range(&range);
+ Rast_read_range(ncb.in[0].name, "", &range);
+ Rast_get_range_min_max(&range, &min, &max);
+ for (i = 1; i < ncb.nin; i++) {
+ Rast_init_range(&range);
+ Rast_read_range(ncb.in[i].name, "", &range);
+ Rast_get_range_min_max(&range, &imin, &imax);
+ if (min > imin)
+ min = imin;
+ if (max < imax)
+ max = imax;
}
+ ci.tmin = min;
+ ci.ntypes = max - min + 1;
+ ci.dts_size = ci.ntypes * ci.nsizebins;
- if (flag.circle->answer)
- circle_mask();
+ ci.n = G_malloc(ncb.nin * sizeof(int));
+ ci.dt = G_malloc(ncb.nin * sizeof(double *));
+ ci.dt[0] = G_malloc(ncb.nin * ci.ntypes * sizeof(double));
+ ci.ds = G_malloc(ncb.nin * sizeof(double *));
+ ci.ds[0] = G_malloc(ncb.nin * ci.nsizebins * sizeof(double));
+ ci.dts = G_malloc(ncb.nin * sizeof(double *));
+ ci.dts[0] = G_malloc(ncb.nin * ci.dts_size * sizeof(double));
+ ci.ht = G_malloc(ncb.nin * sizeof(double));
+ ci.hs = G_malloc(ncb.nin * sizeof(double));
+ ci.hts = G_malloc(ncb.nin * sizeof(double));
- values_w = NULL;
- values_w_tmp = NULL;
- if (weights) {
- values_w =
- (DCELL(*)[2]) G_malloc(ncb.nsize * ncb.nsize * 2 * sizeof(DCELL));
- values_w_tmp =
- (DCELL(*)[2]) G_malloc(ncb.nsize * ncb.nsize * 2 * sizeof(DCELL));
+ ci.ch = G_malloc(ncb.nin * sizeof(struct c_h));
+
+ for (i = 0; i < ncb.nin; i++) {
+ ci.ch[i].palloc = ci.ntypes;
+ ci.ch[i].pst = G_malloc(ci.ntypes * sizeof(struct pst));
+
+ ci.ch[i].pid_curr = G_malloc(ncb.nsize * sizeof(int));
+ ci.ch[i].pid_prev = G_malloc(ncb.nsize * sizeof(int));
+
+ if (i > 0) {
+ ci.dt[i] = ci.dt[i - 1] + ci.ntypes;
+ ci.ds[i] = ci.ds[i - 1] + ci.nsizebins;
+ ci.dts[i] = ci.dts[i - 1] + ci.dts_size;
+ }
}
- values = (DCELL *) G_malloc(ncb.nsize * ncb.nsize * sizeof(DCELL));
- values_tmp = (DCELL *) G_malloc(ncb.nsize * ncb.nsize * sizeof(DCELL));
- for (row = 0; row < nrows; row++) {
- G_percent(row, nrows, 2);
- readcell(in_fd, readrow++, nrows, ncols);
+ /* allocate the cell buffers */
+ allocate_bufs();
- if (selection)
- Rast_get_null_value_row(selection_fd, selection, row);
+ /* initialize the cell bufs with 'nsize - 1' rows of the old cellfile */
+ readrow = 0;
+ for (row = 0; row < ncb.nsize - 1; row++)
+ readcell(readrow++, nrows, ncols);
- for (col = 0; col < ncols; col++) {
+ for (row = 0; row < nrows - ncb.nsize + 1; row++) {
+ G_percent(row, nrows, 5);
+ readcell(readrow++, nrows, ncols);
+
+ if (row % step)
+ continue;
- if (selection && selection[col]) {
- /* ncb.buf length is region row length + 2 * ncb.dist (eq. floor(neighborhood/2))
- * Thus original data start is shifted by ncb.dist! */
- for (i = 0; i < num_outputs; i++)
- outputs[i].buf[col] = ncb.buf[ncb.dist][col + ncb.dist];
- continue;
- }
+ for (col = 0, ocol = 0; col < ncols - ncb.nsize + 1; col += step, ocol++) {
- if (weights)
- n = gather_w(values, values_w, col);
- else
- n = gather(values, col);
+ n = gather(col);
for (i = 0; i < num_outputs; i++) {
struct output *out = &outputs[i];
- DCELL *rp = &out->buf[col];
+ DCELL *rp = &out->buf[ocol];
if (n == 0) {
Rast_set_d_null_value(rp, 1);
}
else {
- if (out->method_fn_w) {
- memcpy(values_w_tmp, values_w, n * 2 * sizeof(DCELL));
- (*out->method_fn_w)(rp, values_w_tmp, n, &out->quantile);
- }
- else {
- memcpy(values_tmp, values, n * sizeof(DCELL));
- (*out->method_fn)(rp, values_tmp, n, &out->quantile);
- }
+ *rp = (*out->method_fn)();
}
}
}
@@ -455,29 +441,87 @@
}
G_percent(row, nrows, 2);
- Rast_close(in_fd);
+ for (i = 0; i < ncb.nin; i++)
+ Rast_close(ncb.in[i].fd);
- if (selection)
- Rast_close(selection_fd);
for (i = 0; i < num_outputs; i++) {
+ char ncs;
+ int nc;
+
Rast_close(outputs[i].fd);
- /* put out category info */
- null_cats(outputs[i].title);
- if (outputs[i].cat_names)
- outputs[i].cat_names();
-
- Rast_write_cats(outputs[i].name, &ncb.cats);
-
- if (copycolr && outputs[i].copycolr)
- Rast_write_colors(outputs[i].name, G_mapset(), &colr);
-
Rast_short_history(outputs[i].name, "raster", &history);
Rast_command_history(&history);
+ ncs = parm.method->answers[i][strlen(parm.method->answers[i]) - 1];
+ nc = 0;
+ if (ncs == '1')
+ nc = ci.ntypes;
+ else if (ncs == '2')
+ nc = ci.nsizebins;
+ else if (ncs == '3')
+ nc = ci.dts_size;
+ Rast_format_history(&history, HIST_DATSRC_1,
+ "Change assessment with %s, %d classes, window size %d, step %d",
+ parm.method->answers[i], nc, ncb.nsize, step);
Rast_write_history(outputs[i].name, &history);
+
+ Rast_put_cell_title(outputs[i].name, outputs[i].title);
+
+ Rast_init_fp_range(&drange);
+ Rast_read_fp_range(outputs[i].name, G_mapset(), &drange);
+ Rast_get_fp_range_min_max(&drange, &dmin, &dmax);
+ make_colors(&colr, dmin, dmax);
+ Rast_write_colors(outputs[i].name, G_mapset(), &colr);
}
exit(EXIT_SUCCESS);
}
+int make_colors(struct Colors *colr, DCELL min, DCELL max)
+{
+ DCELL val1, val2;
+ int r1, g1, b1, r2, g2, b2;
+ DCELL rng;
+
+ if (Rast_is_d_null_value(&min)) {
+ min = 0;
+ max = 1;
+ }
+ if (min > 0)
+ min = 0;
+
+ rng = max - min;
+
+ /* colors: green -> yellow -> orange -> red */
+ Rast_init_colors(colr);
+ val1 = min;
+ val2 = min + rng * 0.2;
+ r1 = 150;
+ g1 = 255;
+ b1 = 150;
+ r2 = 255;
+ g2 = 255;
+ b2 = 0;
+ Rast_add_d_color_rule(&val1, r1, g1, b1, &val2, r2, g2, b2, colr);
+ val1 = val2;
+ r1 = r2;
+ g1 = g2;
+ b1 = b2;
+ val2 = min + rng * 0.5;
+ r2 = 255;
+ g2 = 122;
+ b2 = 0;
+ Rast_add_d_color_rule(&val1, r1, g1, b1, &val2, r2, g2, b2, colr);
+ val1 = val2;
+ r1 = r2;
+ g1 = g2;
+ b1 = b2;
+ val2 = max + 1.0e-12;
+ r2 = 255;
+ g2 = 0;
+ b2 = 0;
+ Rast_add_d_color_rule(&val1, r1, g1, b1, &val2, r2, g2, b2, colr);
+
+ return 1;
+}
Modified: grass-addons/grass7/raster/r.change.info/ncb.h
===================================================================
--- grass-addons/grass7/raster/r.change.info/ncb.h 2016-02-18 21:32:17 UTC (rev 67884)
+++ grass-addons/grass7/raster/r.change.info/ncb.h 2016-02-18 21:38:31 UTC (rev 67885)
@@ -1,15 +1,23 @@
#include <grass/raster.h>
+struct input
+{
+ int fd;
+ const char *name;
+ CELL **buf; /* for reading raster map */
+};
+
struct ncb /* neighborhood control block */
{
- DCELL **buf; /* for reading raster map */
- int *value; /* neighborhood values */
- int nsize; /* size of the neighborhood */
- int dist; /* nsize/2 */
+ int nsize; /* radius * 2 + 1 */
+#if 0
+ double dist; /* radius of the neighborhood */
+#endif
+ int n; /* number of unmasked cells */
+ char **mask;
struct Categories cats;
- char **mask;
- DCELL **weights;
- const char *oldcell;
+ int nin; /* number of input maps */
+ struct input *in;
};
extern struct ncb ncb;
Copied: grass-addons/grass7/raster/r.change.info/r.change.info.html (from rev 67884, grass-addons/grass7/raster/r.change.info/r.neighbors.html)
===================================================================
--- grass-addons/grass7/raster/r.change.info/r.change.info.html (rev 0)
+++ grass-addons/grass7/raster/r.change.info/r.change.info.html 2016-02-18 21:38:31 UTC (rev 67885)
@@ -0,0 +1,229 @@
+<h2>DESCRIPTION</h2>
+
+<em><b>r.change.info</b></em> calculates landscape change assessment
+for a series of categorical maps, e.g. land cover/land use, with
+different measures based on information theory and machine learning.
+More than two <b>input</b> maps can be specified.
+
+<p>
+<em><b>r.change.info</b></em> moves a processing window over the
+<b>input</b> maps. This processing window is the current landscape
+under consideration. The size of the window is defined with
+<b>size</b>. Change assessment is done for each processing window
+(landscape) separately. The centers of the processing windows are
+<b>step</b> cells apart and the <b>output</b> maps will have a
+resolution of <b>step</b> input cells. <b>step</b> should not be larger
+than <b>size</b>, otherwise some cells will be skipped. If <b>step</b>
+is half of <b>size</b> , the moving windows will overlap by 50%. The
+overlap increases when <b>step</b> becomes smaller. A smaller
+<b>step</b> and/or a larger <b>size</b> will require longer processing
+time.
+
+<p>
+The measures <em>information gain</em>, <em>information gain
+ratio</em>, <em>CHI-square</em> and <em>Gini-impurity</em> are commonly
+used in decision tree modelling (Quinlan 1986) to compare
+distributions. These measures as well as the statistical distance are
+based on landscape structure and are calculated for the distributions
+of patch categories and/or patch sizes. A patch is a contiguous block
+of cells with the same category (class), for example a forest fragment.
+The proportion of changes is based on cell changes in the current
+landscape.
+
+<h3>Cell-based change assessment</h3>
+
+The method <b>pc</b> calculates the <em>proportion of changes</em> as
+the actual number of cell changes in the current landscape divided by
+the theoretical maximum number of changes (number of cells in the
+processing window x (number of input maps - 1)).
+
+<h3>Landscape structure change assessment</h3>
+
+<h4>Landscape structure</h4>
+For each processing window, the number of cells per category are
+counted and patches are identified.
+The size and category of each patch is recorded. From these cell and
+patch statistics, three kinds of patterns (distributions) are
+calculated:
+<dl>
+<dt><strong>1. Distributions over categories (e.g land cover class)</strong>
+<dd>This provides information about changes in categories (e.g land
+cover class), e.g. if one category becomes more prominent. This detects
+changes in category composition.
+<dt><strong>2. Distributions over size classes</strong>
+<dd>This provides information about fragmentation, e.g. if a few large
+fragments are broken up into many small fragments. This detects changes
+in fragmentation.
+<dt><strong>3. Distributions over categories and size classes.</strong>
+<dd>This provides information about whether particular combinations of
+category and size class changed between input maps. This detects
+changes in the general landscape structure.
+</dl>
+
+The latter is obtained from the category and size of each patch, i.e.
+each unique category - size class combination becomes a separate class.
+<p>
+The numbers indicate which distribtution will be used for the selected
+method (see below).
+<p>
+A low change in category distributions and a high change in
+size distributions means that the frequency of categories did not
+change much whereas that the size of patches changed.
+
+<h4>Information gain</h4>
+The methods <b>gain1,gain2,gain3</b> calculate the <em>information
+gain</em> after Quinlan (1986). The information gain is the difference
+between the entropy of the combined distribution and the average
+entropy of the observed distributions (conditional entropy). A larger
+value means larger differences between input maps.
+<p>
+Information gain indicates the absolute amount of information gained
+(to be precise, reduced uncertainty) when considering the individual
+input maps instead of their combination. When cells and patches are
+distributed over a large number of categories and a large number of size
+classes, information gain tends to over-estimate changes.
+<p>
+The information gain can be zero even if all cells changed, but the
+distributions (frequencies of occurence) remained identical. The square
+root of the information gain is sometimes used as a distance measure
+and is closely related to Fisher's information metric.
+
+<h4>Information gain ratio</h4>
+The methods <b>ratio1,ratio2,ratio3</b> calculate the <em>information
+gain ratio</em> that changes occured, estimated with the ratio of the
+average entropy of the observed distributions to the entropy of the
+combined distribution. In other words, the ratio is equivalent to the
+ratio of actual change to maximum possible change (in uncertainty). The
+gain ratio is better suited than absolute information gain when the
+cells are distributed over a large number of categories and a large number
+of size classes. The gain ratio here follows the same rationale like
+the gain ratio of Quinlan (1986), but is calculated differently.
+<p>
+The gain ratio is always in the range (0, 1). A larger value means
+larger differences between input maps.
+
+<h4>CHI-square</h4>
+The methods <b>chisq1,chisq2,chisq3</b> calculate <em>CHI square</em>
+after Quinlan (1986) to estimate the relevance of the different input
+maps. If the input maps are identical, the relevance measured as
+CHI-square is zero, no change occured. If the input maps differ from
+each other substantially, major changes occured and the relevance
+measured as CHI-square is large.
+
+<h4>Gini impurity</h4>
+The methods <b>gini1,gini2,gini3</b> calculate the <em>Gini
+impurity</em>, which is 1 - Simpson's index, or 1 - 1 / diversity, or 1
+- 1 / 2^entropy for alpha = 1. The Gini impurity can thus be regarded
+as a modified measure of the diversity of a distribution. Changes
+occured when the diversity of the combined distribution is larger than
+the average diversity of the observed distributions, thus a larger
+value means larger differences between input maps.
+<p>
+The Gini impurity is always in the range (0, 1) and calculated with<br><br>
+G = 1 - <big>∑</big> p<sub>i</sub><sup>2</sup>
+
+<p>
+The methods <em>information gain</em> and <em>CHI square</em> are the
+most sensitive measures, but also most susceptible to noise. The
+<em>information gain ratio</em> is less sensitive, but more robust
+against noise. The <em>Gini impurity</em> is the least sensitive and
+detects only drastic changes.
+
+
+<h4>Distance</h4>
+The methods <b>dist1,dist2,dist3</b> calculate the statistical
+<em>distance</em> from the absolute differences between the average
+distribution and the observed distributions. The distance is always in
+the range (0, 1). A larger value means larger differences between input
+maps.
+
+<p>
+Methods using the category or size class distributions (<em>gain1</em>,
+<em>gain2</em>, <em>ratio1</em>, <em>ratio2</em> <em>gini1</em>,
+<em>gini2</em>, <em>dist1</em>, <em>dist2</em>) are less sensitive than
+methods using the combined category and size class distributions
+(<em>gain3</em>, <em>ratio3</em>, <em>gini3</em>, <em>dist3</em>).
+
+<p>
+For a thorough change assessment it is recommended to calculate
+different change assessment measures (at least information gain and
+information gain ratio) and investigate the differences between these
+change assessments.
+
+<h2>NOTES</h2>
+<h3>Shannon's entropy</h3>
+Entropies for information gain and its ratio are by default Shannon
+entropies <em>H</em>, calculated with<br><br>
+H = <big>∑</big> p<sub>i</sub> * log<sub>2</sub>(p<sub>i</sub>)
+
+<p>
+The entropies are here calculated with base 2 logarithms. The upper
+bound of information gain is thus log<sub>2</sub>(number of classes).
+Classes can be categories, size classes, or unique combinations of
+categories and size classes.
+
+<h3>Rényi's entropy</h3>
+The <b>alpha</b> option can be used to calculate general entropies
+<em>H<sub>α</sub></em> after Rényi (1961) with the formula<br><br>
+H<sub>α</sub> = 1 / (1 - α) * log<sub>2</sub>
+(<big>∑</big> p<sub>i</sub><sup>α</sup>)
+
+<p>
+An <b>alpha</b> < 1 gives higher weight to small frequencies,
+whereas an <b>alpha</b> > 1 gives higher weight to large
+frequencies. This is useful for noisy input data such as the MODIS land
+cover/land use products MCD12*. These data often differ only in
+single-cell patches. These differences can be due to the applied
+classification procedure. Moreover, the probabilities that a cell has
+been assigned to class A or class B are often very similar, i.e.
+different classes are confused by the applied classification procedure.
+In such cases an <b>alpha</b> > 1, e.g. 2, will give less weight to
+small changes and more weight to large changes, to a degree alleviating
+the problem of class confusion.
+
+<h2>EXAMPLES</h2>
+
+Assuming there is a time series of the MODIS land cover/land use
+product MCD12Q1, land cover type 1, available, and the raster maps have
+the names
+<div class="code"><pre>
+MCD12Q1.A2001.Land_Cover_Type_1
+MCD12Q1.A2002.Land_Cover_Type_1
+MCD12Q1.A2003.Land_Cover_Type_1
+...
+</pre></div>
+
+then a change assement can be done with
+<div class="code"><pre>
+r.change.info in=`g.mlist type=rast pat=MCD12Q1.A*.Land_Cover_Type_1 sep=,` \
+ method=pc,gain1,gain2,ratio1,ratio2,dist1,dist2
+ out=MCD12Q1.pc,MCD12Q1.gain1,MCD12Q1.gain2,MCD12Q1.ratio1,MCD12Q1.ratio2,MCD12Q1.dist1,MCD12Q1.dist2 \
+ radius=20 step=40 alpha=2
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em><a href="r.neighbors.html">r.neighbors</a></em><br>
+<em><a href="g.region.html">g.region</a></em><br>
+<em><a href="r.li.shannon.html">r.li.shannon</a></em><br>
+<em><a href="r.li.simpson.html">r.li.simpson</a></em><br>
+<em><a href="r.li.renyi.html">r.li.renyi</a></em><br>
+<em><a href="r.clump.html">r.clump</a></em>
+
+<h2>REFERENCES</h2>
+
+<ul>
+<li>
+Quinlan, J.R. 1986. Induction of decision trees. Machine Learning 1: 81-106.
+<a href="http://dx.doi.org/10.1007/BF00116251">DOI:10.1007/BF00116251</a>
+<li>
+Rényi, A. 1961. <a href="http://digitalassets.lib.berkeley.edu/math/ucb/text/math_s4_v1_article-27.pdf">On measures of information and entropy.</a> Proceedings of the fourth Berkeley Symposium on Mathematics, Statistics and Probability 1960: 547-561.
+<li>
+Shannon, C.E. 1948. A Mathematical Theory of Communication. Bell System Technical Journal 27(3): 379-423. <a href="http://dx.doi.org/10.1002/j.1538-7305.1948.tb01338.x">DOI:10.1002/j.1538-7305.1948.tb01338.x</a>
+</ul>
+
+<h2>AUTHOR</h2>
+
+Markus Metz
+
+<p><i>Last changed: $Date$</i>
Deleted: grass-addons/grass7/raster/r.change.info/r.neighbors.html
===================================================================
--- grass-addons/grass7/raster/r.change.info/r.neighbors.html 2016-02-18 21:32:17 UTC (rev 67884)
+++ grass-addons/grass7/raster/r.change.info/r.neighbors.html 2016-02-18 21:38:31 UTC (rev 67885)
@@ -1,392 +0,0 @@
-<h2>DESCRIPTION</h2>
-
-<em><b>r.neighbors</b></em> looks at each cell in a raster input
-file, and examines the values assigned to the
-cells in some user-defined "neighborhood" around it. It
-outputs a new raster map layer in which each cell is
-assigned a value that is some (user-specified)
-function of the values in that cell's neighborhood. For
-example, each cell in the output layer might be assigned a
-value equal to the average of the values
-appearing in its 3 x 3 cell "neighborhood" in the input
-layer. Note that the centre cell is also included in the calculation.
-
-<h3>OPTIONS</h3>
-
-The user must specify the names of the raster map layers to
-be used for <b>input</b> and <b>output</b>, the
-<b>method</b> used to analyze neighborhood
-values (i.e., the neighborhood function or operation to be
-performed), and the <b>size</b> of the neighborhood.
-<p>The user can optionally
-specify a <b>selection</b> map, to compute new values only where the raster
-cells of the selection map are not NULL. In case of a NULL cells,
-the values from the input map are copied into the output map.
-This may useful to smooth only parts of an elevation map (pits, peaks, ...).
-
-<p><em>Example how to use a selection map with method=average:</em><br>
-input map:
-<div class="code"><pre>
-1 1 1 1 1
-1 1 1 1 1
-1 1 10 1 1
-1 1 1 1 1
-1 1 1 1 1
-</pre></div>
-selection map, NULL values are marked as *:
-<div class="code"><pre>
-* * * * *
-* * 1 * *
-* 1 1 1 *
-* * 1 * *
-* * * * *
-</pre></div>
-The output map:
-<div class="code"><pre>
-1 1 1 1 1
-1 1 2 1 1
-1 2 2 2 1
-1 1 2 1 1
-1 1 1 1 1
-</pre></div>
-Without using the selection map, the output map would look like this:
-<div class="code"><pre>
-1 1 1 1 1
-1 2 2 2 1
-1 2 2 2 1
-1 2 2 2 1
-1 1 1 1 1
-</pre></div>
-
-<p>Optionally, the user can also specify the <b>TITLE</b> to
-be assigned to the raster map layer <b>output</b>, elect
-to not align the resolution of the output with that of the
-input (the <b>-a</b> option), and run <em><b>r.neighbors</b></em>
-with a custom matrix weights with the <em>weight</em> option.
-These options are described further below.
-<p>
-<em>Neighborhood Operation Methods:</em>
-The <b>neighborhood</b> operators determine what new
-value a center cell in a neighborhood will have after examining
-values inside its neighboring cells.
-Each cell in a raster map layer becomes the center cell of a neighborhood
-as the neighborhood window moves from cell to cell throughout the map layer.
-<em><b>r.neighbors</b></em> can perform the following operations:
-
-<p><dl>
-
-<dt><b>average</b>
-
-<dd>The average value within the neighborhood.
-In the following example, the result would be:
-
-<br>
-
-(7*4 + 6 + 5 + 4*3)/9 = 5.66
-
-<br>
-
-The result is rounded to the nearest integer (in this case 6).
-
-<dt><b>median</b>
-
-<dd>The value found half-way through a list of the
-neighborhood's values,
-when these are ranged in numerical order.
-
-<dt><b>mode</b>
-
-<dd>The most frequently occurring value in the neighborhood.
-
-<dt><b>minimum</b>
-
-<dd>The minimum value within the neighborhood.
-
-<dt><b>maximum</b>
-
-<dd>The maximum value within the neighborhood.
-<div class="code"><pre>
- Raw Data Operation New Data
- +---+---+---+ +---+---+---+
- | 7 | 7 | 5 | | | | |
- +---+---+---+ average +---+---+---+
- | 4 | 7 | 4 |--------->| | 6 | |
- +---+---+---+ +---+---+---+
- | 7 | 6 | 4 | | | | |
- +---+---+---+ +---+---+---+
-</pre></div>
-
-<dt><b>range</b>
-
-<dd>The range value within the neighborhood.
-
-<dt><b>stddev</b>
-
-<dd>The statistical standard deviation of values
-within the neighborhood (rounded to the nearest integer).
-
-<dt><b>sum</b>
-
-<dd>The sum of values within the neighborhood.
-
-<dt><b>variance</b>
-
-<dd>The statistical variance of values
-within the neighborhood (rounded to the nearest integer).
-
-<dt><b>diversity</b>
-
-<dd>The number of different values within the neighborhood.
-In the above example, the diversity is 4.
-
-<dt><b>interspersion</b>
-
-<dd>The percentage of cells containing values which differ from the values
-assigned to the center cell in the neighborhood, plus 1.
-In the above example, the interspersion is:
-
-<br>
-
-5/8 * 100 + 1 = 63.5
-
-<br>
-
-The result is rounded to the nearest integer (in this case 64).
-
-</dl>
-<p><br>
-
-<em>Neighborhood Size:</em>
-The neighborhood <b>size</b> specifies which cells surrounding any given
-cell fall into the neighborhood for that cell.
-The <b>size</b> must be an odd integer and represent the length of one of
-moving window edges in cells.
-For example, a size value of 3 will result in
-<div class="code"><pre>
- _ _ _
- |_|_|_|
- 3 x 3 neighborhood ---> |_|_|_|
- |_|_|_|
-
-</pre></div>
-<p>
-<em>Matrix weights:</em>
-A custom matrix can be used if none of the neighborhood operation
-methods are desirable by using the <b>weight</b>. This option must
-be used in conjunction with the <b>size</b> option to specify the
-matrix size. The weights desired are to be entered into a text file.
-For example, to calculate the focal mean with a matrix <b>size</b> of
-3,
-<div class="code"><pre>
-r.neigbors in=input.map out=output.map size=3 weight=weights.txt
-</pre></div>
-
-The contents of the weight.txt file:
-<div class="code"><pre>
-3 3 3
-1 4 8
-9 5 3
-</pre></div>
-
-This corresponds to the following 3x3 matrix:
-<div class="code"><pre>
-+-+-+-+
-|3|3|3|
-+-+-+-+
-|1|4|8|
-+-+-+-+
-|9|5|3|
-+-+-+-+
-</pre></div>
-
-To calculate an annulus shaped neighborhood the contents of weight.txt file
-may be e.g. for size=5:
-<div class="code"><pre>
- 0 1 1 1 0
- 1 0 0 0 1
- 1 0 0 0 1
- 1 0 0 0 1
- 0 1 1 1 0
-</pre></div>
-
-The way that weights are used depends upon the specific aggregate
-(<b>method</b>) being used.
-
-However, most of the aggregates have the property that multiplying all
-of the weights by the same factor won't change the final result (an
-exception is <b>method=count</b>).
-
-Also, most (if not all) of them have the properties that an integer
-weight of N is equivalent to N occurrences of the cell value, and
-having all weights equal to one produces the same result as when
-weights are not used.
-
-When weights are used, the calculation for <b>method=average</b> is:
-
-<div class="code"><pre>
- sum(w[i]*x[i]) / sum(w[i])
-</pre></div>
-
-In the case where all weights are zero, this will end up with both the
-numerator and denominator to zero, which produces a NULL result.
-
-<h3>FLAGS</h3>
-<dl>
-<dt><b>-a</b>
-
-<dd>If specified, <em><b>r.neighbors</b></em> will not align the output
-raster map layer with that of the input raster map layer.
-The <em><b>r.neighbors</b></em> program works in the current geographic region.
-It is recommended, but not required, that the resolution
-of the geographic region be the same as that of the raster map layer.
-By default, if unspecified,
-<em><b>r.neighbors</b></em> will align these geographic region settings.
-<p>
-<dt><b>-c</b>
-<dd>
-This flag will use a circular neighborhood for the moving analysis window,
-centered on the current cell.
-
-<p>The exact masks for the first few neighborhood sizes are as follows:
-<div class="code"><pre>
-3x3 . X . 5x5 . . X . . 7x7 . . . X . . .
- X O X . X X X . . X X X X X .
- . X . X X O X X . X X X X X .
- . X X X . X X X O X X X
- . . X . . . X X X X X .
- . X X X X X .
- . . . X . . .
-
-9x9 . . . . X . . . . 11x11 . . . . . X . . . . .
- . . X X X X X . . . . X X X X X X X . .
- . X X X X X X X . . X X X X X X X X X .
- . X X X X X X X . . X X X X X X X X X .
- X X X X O X X X X . X X X X X X X X X .
- . X X X X X X X . X X X X X O X X X X X
- . X X X X X X X . . X X X X X X X X X .
- . . X X X X X . . . X X X X X X X X X .
- . . . . X . . . . . X X X X X X X X X .
- . . X X X X X X X . .
- . . . . . X . . . . .
-</pre></div>
-
-</dl>
-
-<h2>NOTES</h2>
-
-The <em><b>r.neighbors</b></em> program works in the current geographic region
-with the current mask, if any. It is recommended, but not required,
-that the resolution of the geographic region be the same as that
-of the raster map layer. By default, <em><b>r.neighbors</b></em> will align
-these geographic region settings. However, the user can select to keep
-original input and output resolutions which are not aligned by specifying
-this (e.g., using the <b>-a</b> option).
-<p>
-<em><b>r.neighbors</b></em> doesn't propagate NULLs, but computes the
-aggregate over the non-NULL cells in the neighborhood.
-<p>
-The <b>-c</b> flag and the <b>weights</b> parameter are mutually exclusive. Any
-use of the two together will produce an error. Differently-shaped neighborhood
-analysis windows may be achieved by using the <b>weight=</b> parameter to
-specify a weights file where all values are equal. The user can also vary the
-weights at the edge of the neighborhood according to the proportion of the cell
-that lies inside the neighborhood circle, effectively anti-aliasing the analysis
-mask.
-<p>
-For aggregates where a weighted calculation isn't meaningful
-(specifically: minimum, maximum, diversity and interspersion), the
-weights are used to create a binary mask, where zero causes the cell
-to be ignored and any non-zero value causes the cell to be used.
-<p>
-<em><b>r.neighbors</b></em> copies the GRASS <em>color</em> files associated with
-the input raster map layer for those output map layers that are based
-on the neighborhood average, median, mode, minimum, and maximum.
-Because standard deviation, variance, diversity, and interspersion are indices,
-rather than direct correspondents to input values,
-no <em>color</em> files are copied for these map layers.
-(The user should note that although the <em>color</em> file is copied
-for <em>average</em> neighborhood function output,
-whether or not the color file makes sense for the output
-will be dependent on the input data values.)
-
-<h3>Propagation of output precision</h3>
-
-The following logic has been implemented: For any aggregate, there are
-two factors affecting the output type:
-
-<ol>
-<li> Whether the input map is integer or floating-point.</li>
-<li> Whether the weighted or unweighted version of the aggregate is used.</li>
-</ol>
-
-These combine to create four possibilities:
-<p>
-<table border="1">
- <tr><th>input type</th><th>integer</th><th>integer</th><th>float</th><th>float</th></tr>
- <tr><td>weighted</td><td>no</td><td>yes</td><td>no</td><td>yes</td></tr>
- <tr><td> </td><td> </td><td> </td><td> </td><td> </td></tr>
- <tr><td>average</td><td>float</td><td>float</td><td>float</td><td>float</td></tr>
- <tr><td>median</td><td>[1]</td><td>[1]</td><td>float</td><td>float</td></tr>
- <tr><td>mode</td><td>integer</td><td>integer</td><td>[2]</td><td>[2]</td></tr>
- <tr><td>minimum</td><td>integer</td><td>integer</td><td>float</td><td>float</td></tr>
- <tr><td>maximum</td><td>integer</td><td>integer</td><td>float</td><td>float</td></tr>
- <tr><td>range</td><td>integer</td><td>integer</td><td>float</td><td>float</td></tr>
- <tr><td>stddev</td><td>float</td><td>float</td><td>float</td><td>float</td></tr>
- <tr><td>sum</td><td>integer</td><td>float</td><td>float</td><td>float</td></tr>
- <tr><td>count</td><td>integer</td><td>float</td><td>integer</td><td>float</td></tr>
- <tr><td>variance</td><td>float</td><td>float</td><td>float</td><td>float</td></tr>
- <tr><td>diversity</td><td>integer</td><td>integer</td><td>integer</td><td>integer</td></tr>
- <tr><td>interspersion</td><td>integer</td><td>integer</td><td>integer</td><td>integer</td></tr>
- <tr><td>quart1</td><td>[1]</td><td>[1]</td><td>float</td><td>float</td></tr>
- <tr><td>quart3</td><td>[1]</td><td>[1]</td><td>float</td><td>float</td></tr>
- <tr><td>perc90</td><td>[1]</td><td>[1]</td><td>float</td><td>float</td></tr>
- <tr><td>quantile</td><td>[1]</td><td>[1]</td><td>float</td><td>float</td></tr>
-</table>
-<p>
-[1] For integer input, quantiles may produce float results from
-interpolating between adjacent values.
-<br>
-[2] Calculating the mode of floating-point data is essentially
-meaningless.
-<p>
-
-With the current aggregates, there are 5 cases:
-
-<ol>
-<li> Output is always float: average, variance, stddev, quantiles (with
-interpolation).</li>
-<li> Output is always integer: diversity, interspersion.</li>
-<li> Output is integer if unweighted, float if weighted: count.</li>
-<li> Output matches input: minimum, maximum, range, mode (subject to
-note 2 above), quantiles (without interpolation).</li>
-<li> Output is integer for integer input and unweighted aggregate,
-otherwise float: sum.</li>
-</ol>
-
-<!-- TODO
-<h2>EXAMPLES</h2>
-
-<div class="code"><pre>
-r.neighbors
-</pre></div>
-
--->
-
-<h2>SEE ALSO</h2>
-
-<em><a href="g.region.html">g.region</a></em><br>
-<em><a href="r.clump.html">r.clump</a></em><br>
-<em><a href="r.mapcalc.html">r.mapcalc</a></em><br>
-<em><a href="r.mfilter.html">r.mfilter</a></em><br>
-<em><a href="r.statistics.html">r.statistics</a></em><br>
-<em><a href="r.support.html">r.support</a></em>
-
-
-<h2>AUTHORS</h2>
-
-Original version: Michael Shapiro,
-U.S.Army Construction Engineering Research Laboratory
-<br>
-Updates for GRASS GIS 7 by Glynn Clements and others
-
-<p><i>Last changed: $Date$</i>
Added: grass-addons/grass7/raster/r.change.info/ratio.c
===================================================================
--- grass-addons/grass7/raster/r.change.info/ratio.c (rev 0)
+++ grass-addons/grass7/raster/r.change.info/ratio.c 2016-02-18 21:38:31 UTC (rev 67885)
@@ -0,0 +1,89 @@
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include "ncb.h"
+#include "window.h"
+
+/* information gain ratio = proportion of change
+ * dist: array of distributions
+ * nd: number of cells in each distribution
+ * n: number of distributions
+ * dsize: distributions' size
+ * h: entropy for each distribution */
+
+DCELL ratio(double **dist, int *nd, int n, int dsize, double *h)
+{
+ int i, d, n_comb;
+ double p, pe, p_comb, H_comb, H_avg, IV, igr;
+
+ /* entropy of combined distribution */
+ /* average entropy of all distributions */
+
+ H_comb = 0;
+ H_avg = 0;
+ IV = 0;
+ n_comb = 0;
+ for (d = 0; d < n; d++) {
+ h[d] = 0;
+ n_comb += nd[d];
+ }
+ for (i = 0; i < dsize; i++) {
+ p_comb = 0;
+ for (d = 0; d < n; d++) {
+ if (dist[d][i] > 0 && nd[d] > 0) {
+ p = dist[d][i] / nd[d];
+ p_comb += dist[d][i];
+ pe = entropy_p(p);
+ h[d] += pe;
+ }
+ }
+ if (p_comb > 0) {
+ H_comb += entropy_p(p_comb / n_comb);
+ }
+ }
+ H_avg = 0;
+ IV = 0;
+ for (d = 0; d < n; d++) {
+ H_avg += entropy(h[d]) * nd[d] / n_comb;
+ /* intrinsic value in decision tree modelling
+ * here a constant if there are no NULL cells */
+ IV += entropy_p((double)nd[d] / n_comb);
+ }
+ H_comb = entropy(H_comb);
+ IV = entropy(IV);
+
+ if (H_comb < 0)
+ H_comb = 0;
+ if (H_avg < 0)
+ H_avg = 0;
+ if (IV < 0)
+ IV = 0;
+
+ if (H_comb == 0)
+ return 0;
+
+ if (H_comb < H_avg)
+ return 0;
+
+ /* information gain rate:
+ * actual change divided by maximum possible change
+ * (H_comb - H_avg) / H_comb */
+ igr = 1 - H_avg / H_comb;
+
+ return igr;
+}
+
+DCELL ratio1(void)
+{
+ return ratio(ci.dt, ci.n, ncb.nin, ci.ntypes, ci.ht);
+}
+
+DCELL ratio2(void)
+{
+ return ratio(ci.ds, ci.n, ncb.nin, ci.nsizebins, ci.hs);
+}
+
+DCELL ratio3(void)
+{
+ return ratio(ci.dts, ci.n, ncb.nin, ci.dts_size, ci.hts);
+}
Modified: grass-addons/grass7/raster/r.change.info/readcell.c
===================================================================
--- grass-addons/grass7/raster/r.change.info/readcell.c 2016-02-18 21:32:17 UTC (rev 67884)
+++ grass-addons/grass7/raster/r.change.info/readcell.c 2016-02-18 21:38:31 UTC (rev 67885)
@@ -4,14 +4,20 @@
#include "ncb.h"
#include "local_proto.h"
-int readcell(int fd, int row, int nrows, int ncols)
+int readcell(int row, int nrows, int ncols)
{
+ int i;
+
rotate_bufs();
- if (row < nrows)
- Rast_get_d_row(fd, ncb.buf[ncb.nsize - 1] + ncb.dist, row);
- else
- Rast_set_d_null_value(ncb.buf[ncb.nsize - 1] + ncb.dist, ncols);
+ if (row < nrows) {
+ for (i = 0; i < ncb.nin; i++)
+ Rast_get_c_row(ncb.in[i].fd, ncb.in[i].buf[ncb.nsize - 1], row);
+ }
+ else {
+ for (i = 0; i < ncb.nin; i++)
+ Rast_set_c_null_value(ncb.in[i].buf[ncb.nsize - 1], ncols);
+ }
return 0;
}
Deleted: grass-addons/grass7/raster/r.change.info/readweights.c
===================================================================
--- grass-addons/grass7/raster/r.change.info/readweights.c 2016-02-18 21:32:17 UTC (rev 67884)
+++ grass-addons/grass7/raster/r.change.info/readweights.c 2016-02-18 21:38:31 UTC (rev 67885)
@@ -1,44 +0,0 @@
-#include <stdio.h>
-#include <math.h>
-#include <grass/gis.h>
-#include <grass/glocale.h>
-#include "ncb.h"
-#include "local_proto.h"
-
-void read_weights(const char *filename)
-{
- FILE *fp = fopen(filename, "r");
- int i, j;
-
- ncb.weights = G_malloc(ncb.nsize * sizeof(DCELL *));
- for (i = 0; i < ncb.nsize; i++)
- ncb.weights[i] = G_malloc(ncb.nsize * sizeof(DCELL));
-
- if (!fp)
- G_fatal_error(_("Unable to open weights file %s"), filename);
-
- for (i = 0; i < ncb.nsize; i++)
- for (j = 0; j < ncb.nsize; j++)
- if (fscanf(fp, "%lf", &ncb.weights[i][j]) != 1)
- G_fatal_error(_("Error reading weights file %s"), filename);
-
- fclose(fp);
-}
-
-void gaussian_weights(double sigma)
-{
- double sigma2 = sigma * sigma;
- int i, j;
-
- ncb.weights = G_malloc(ncb.nsize * sizeof(DCELL *));
- for (i = 0; i < ncb.nsize; i++)
- ncb.weights[i] = G_malloc(ncb.nsize * sizeof(DCELL));
-
- for (i = 0; i < ncb.nsize; i++) {
- double y = i - ncb.dist;
- for (j = 0; j < ncb.nsize; j++) {
- double x = j - ncb.dist;
- ncb.weights[i][j] = exp(-(x*x+y*y)/(2*sigma2))/(2*M_PI*sigma2);
- }
- }
-}
Added: grass-addons/grass7/raster/r.change.info/window.h
===================================================================
--- grass-addons/grass7/raster/r.change.info/window.h (rev 0)
+++ grass-addons/grass7/raster/r.change.info/window.h 2016-02-18 21:38:31 UTC (rev 67885)
@@ -0,0 +1,45 @@
+
+/* type and size of each patch */
+struct pst {
+ CELL type; /* patch type */
+ int size; /* patch size (number of cells) */
+};
+
+/* clumping helper */
+struct c_h {
+ struct pst *pst;
+ int palloc; /* alloc'd patches */
+ int np; /* number of patches */
+
+ CELL up, left, curr;
+ int *pid_curr, *pid_prev;
+ int pid;
+};
+
+/* size bins: 1 - < 2, 2 - < 4, 4 - < 8, 8 - < 16, etc */
+struct changeinfo {
+ /* globals */
+ int nsizebins; /* number of bins for patch sizes */
+ int ntypes; /* number of different patch types */
+ int dts_size; /* dts size */
+ int tmin; /* smallest patch type number */
+
+ int nchanges; /* number of type changes */
+
+ /* for each input map */
+ int *n; /* number of valid cells */
+ double **dt; /* distribution over types */
+ double **ds; /* distribution over size bins */
+ double **dts; /* distribution over types and size bins */
+ double *ht; /* entropy for dt */
+ double *hs; /* entropy for ds */
+ double *hts; /* entropy for dts */
+
+ struct c_h *ch; /* clumping helper */
+};
+
+extern struct changeinfo ci;
+
+
+extern double (*entropy) (double);
+extern double (*entropy_p) (double);
More information about the grass-commit
mailing list