[GRASS-SVN] r57980 - grass-addons/grass7/vector/v.surf.mass
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Oct 11 04:58:06 PDT 2013
Author: mmetz
Date: 2013-10-11 04:58:05 -0700 (Fri, 11 Oct 2013)
New Revision: 57980
Modified:
grass-addons/grass7/vector/v.surf.mass/globals.h
grass-addons/grass7/vector/v.surf.mass/main.c
Log:
v.surf.mass: add weight option
Modified: grass-addons/grass7/vector/v.surf.mass/globals.h
===================================================================
--- grass-addons/grass7/vector/v.surf.mass/globals.h 2013-10-10 22:16:07 UTC (rev 57979)
+++ grass-addons/grass7/vector/v.surf.mass/globals.h 2013-10-11 11:58:05 UTC (rev 57980)
@@ -11,7 +11,7 @@
int area; /* area id */
double interp; /* interpolated value */
double adj; /* adjustment */
-
+ double weight; /* weighing factor */
};
/* mass area */
Modified: grass-addons/grass7/vector/v.surf.mass/main.c
===================================================================
--- grass-addons/grass7/vector/v.surf.mass/main.c 2013-10-10 22:16:07 UTC (rev 57979)
+++ grass-addons/grass7/vector/v.surf.mass/main.c 2013-10-11 11:58:05 UTC (rev 57980)
@@ -36,17 +36,17 @@
int row, col, nrows, ncols, rrow, ccol, nrow, ncol;
SEGMENT out_seg;
- int out_fd;
+ int w_fd, out_fd, have_weights;
DCELL *drastbuf;
double seg_size;
int seg_mb, segments_in_memory;
int doit, iter, maxiter;
double threshold, maxadj;
- int negative = 0, have_neg;
+ int negative = 0;
int layer;
int i, j, nareas;
- struct marea *areas;
+ struct marea *areas, *ap;
struct Map_info In;
struct line_pnts *Points;
@@ -55,10 +55,14 @@
struct History history;
struct lcell thiscell;
- double outside_val, value, interp, avgdiff;
+ double outside_val, value;
+#ifdef TOBLER_STRICT
+ double interp, avgdiff;
+#endif
+
struct GModule *module;
- struct Option *in_opt, *out_opt, *dfield_opt, *col_opt,
+ struct Option *in_opt, *w_opt, *out_opt, *dfield_opt, *col_opt,
*memory_opt, *iter_opt, *threshold_opt;
/* border condition */
struct Flag *withz_flag; /* allow negative */
@@ -87,6 +91,11 @@
in_opt = G_define_standard_option(G_OPT_V_INPUT);
in_opt->label = _("Name of input vector point map");
+ w_opt = G_define_standard_option(G_OPT_R_INPUT);
+ w_opt->label = _("Name of optional weighing raster map");
+ w_opt->key = "weight";
+ w_opt->required = NO;
+
dfield_opt = G_define_standard_option(G_OPT_V_FIELD);
col_opt = G_define_standard_option(G_OPT_DB_COLUMN);
@@ -160,6 +169,21 @@
if (withz_flag->answer)
layer = 0;
+ /* input weights */
+ w_fd = -1;
+ have_weights = 0;
+ drastbuf = NULL;
+ if (w_opt->answer) {
+ char * mapset = (char *) G_find_raster2(w_opt->answer, "");
+
+ if (mapset == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), w_opt->answer);
+
+ have_weights = 1;
+ w_fd = Rast_open_old(w_opt->answer, "");
+ drastbuf = Rast_allocate_d_buf();
+ }
+
/* raster output */
out_fd = Rast_open_new(out_opt->answer, DCELL_TYPE);
@@ -178,14 +202,14 @@
segments_in_memory = seg_mb / seg_size + 0.5;
G_debug(1, "%d %dx%d segments held in memory", segments_in_memory, SEGSIZE, SEGSIZE);
+ /* initialize */
+ G_message(_("Initializing..."));
+
if (segment_open(&out_seg, G_tempfile(),
nrows, ncols, SEGSIZE, SEGSIZE,
sizeof(struct lcell), segments_in_memory) != 1)
G_fatal_error(_("Can not create temporary file"));
- /* initialize */
- G_message(_("Initializing..."));
-
Points = Vect_new_line_struct();
Cats = Vect_new_cats_struct();
db_init_string(&dbstring);
@@ -312,6 +336,9 @@
y = Rast_row_to_northing(row + 0.5, &window);
+ if (have_weights)
+ Rast_get_d_row(w_fd, drastbuf, row);
+
for (col = 0; col < ncols; col++) {
x = Rast_col_to_easting(col + 0.5, &window);
@@ -320,10 +347,28 @@
thiscell.interp = value;
thiscell.adj = 0.0;
thiscell.area = Vect_find_area(&In, x, y);
+ thiscell.weight = 1.;
+ if (have_weights) {
+ if (Rast_is_d_null_value(&drastbuf[col])) {
+ thiscell.weight = 0; /* OK ? */
+ }
+ else {
+ thiscell.weight = drastbuf[col];
+ if (thiscell.weight < 0) {
+ G_warning(_("Weight is negative at row %d, col %d, using 0 instead"),
+ row, col);
+ thiscell.weight = 0;
+ }
+ }
+ }
if (thiscell.area > 0) {
- thiscell.interp = areas[thiscell.area].value;
- areas[thiscell.area].count++;
+ ap = &areas[thiscell.area];
+ thiscell.interp = ap->value;
+ ap->count++;
+ if (have_weights) {
+ ap->interp += thiscell.weight;
+ }
}
segment_put(&out_seg, &thiscell, row, col);
@@ -332,6 +377,11 @@
G_percent(row, nrows, 2);
Vect_close(&In);
+ if (have_weights) {
+ G_free(drastbuf);
+ Rast_close(w_fd);
+ }
+
G_message(_("Adjust cell values"));
for (row = 0; row < nrows; row++) {
G_percent(row, nrows, 2);
@@ -341,7 +391,12 @@
if (areas[thiscell.area].count > 0 &&
!Rast_is_d_null_value(&thiscell.interp)) {
- value = thiscell.interp / areas[thiscell.area].count;
+ ap = &areas[thiscell.area];
+ value = thiscell.interp / ap->count;
+ if (have_weights && ap->interp > 0) {
+ thiscell.weight = thiscell.weight * ap->count / ap->interp;
+ value *= thiscell.weight;
+ }
thiscell.interp = value;
segment_put(&out_seg, &thiscell, row, col);
}
@@ -369,9 +424,9 @@
for (row = 0; row < nrows; row++) {
for (col = 0; col < ncols; col++) {
int count = 0;
+ double sum = 0.;
value = .0;
- interp = .0;
for (rrow = -1; rrow < 2; rrow++) {
nrow = row + rrow;
@@ -385,23 +440,24 @@
segment_get(&out_seg, &thiscell, nrow, ncol);
if (!Rast_is_d_null_value(&thiscell.interp)) {
- value += thiscell.interp;
+ value += thiscell.interp * thiscell.weight;
count++;
+ sum += thiscell.weight;
}
}
}
}
- if (count) {
+ if (sum > 0) {
segment_get(&out_seg, &thiscell, row, col);
if (!Rast_is_d_null_value(&thiscell.interp)) {
- value /= count;
- value -= thiscell.interp;
+ value /= sum;
+ value -= thiscell.interp / thiscell.weight;
/* relax */
/* value /= 8; */
thiscell.adj = value;
- if (thiscell.area && !negative) {
- if (areas[thiscell.area].value == 0)
+ if (!negative) {
+ if (thiscell.area && areas[thiscell.area].value == 0)
thiscell.adj = 0;
if (thiscell.interp + thiscell.adj < 0)
thiscell.adj = -thiscell.interp;
@@ -439,15 +495,15 @@
segment_put(&out_seg, &thiscell, row, col);
value = thiscell.adj + areas[thiscell.area].adj;
+ if (have_weights && thiscell.weight > 1)
+ value /= thiscell.weight;
+
if (maxadj < value * value) {
maxadj = value * value;
l_row = row;
l_col = col;
}
}
- if (!negative && interp < 0)
- have_neg = 1;
-
if (thiscell.area) {
areas[thiscell.area].interp += thiscell.interp;
}
@@ -506,24 +562,22 @@
/* deviation from Tobler 1979: less steps, faster convergence */
/* Step 2 */
- areas[0].adj = 1;
- for (i = 1; i <= nareas; i++) {
- areas[i].adj = 1;
- if (areas[i].interp != 0)
- areas[i].adj = areas[i].value / areas[i].interp;
- }
-
- /* Step 3 */
for (row = 0; row < nrows; row++) {
for (col = 0; col < ncols; col++) {
segment_get(&out_seg, &thiscell, row, col);
if (!Rast_is_d_null_value(&thiscell.interp)) {
/* multiplication with area scale factor */
- thiscell.interp = areas[thiscell.area].adj *
- (thiscell.interp + thiscell.adj);
+ value = 1;
+ ap = &areas[thiscell.area];
+ if (thiscell.area > 0 && ap->interp != 0)
+ value = ap->value / ap->interp;
+ thiscell.interp = value * (thiscell.interp + thiscell.adj);
segment_put(&out_seg, &thiscell, row, col);
- value = thiscell.adj;
+ if (ap->count)
+ value = thiscell.adj - ap->adj / ap->count;
+ if (have_weights && thiscell.weight > 1)
+ value /= thiscell.weight;
if (maxadj < value * value) {
maxadj = value * value;
l_row = row;
More information about the grass-commit
mailing list