[GRASS-SVN] r58215 - grass/trunk/vector/v.what.rast
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Nov 13 22:52:26 PST 2013
Author: hamish
Date: 2013-11-13 22:52:25 -0800 (Wed, 13 Nov 2013)
New Revision: 58215
Modified:
grass/trunk/vector/v.what.rast/local_proto.h
grass/trunk/vector/v.what.rast/main.c
grass/trunk/vector/v.what.rast/v.what.rast.html
Log:
add flag to interpolate values from the nearest four cells instead of using nearest neighbor
Modified: grass/trunk/vector/v.what.rast/local_proto.h
===================================================================
--- grass/trunk/vector/v.what.rast/local_proto.h 2013-11-14 03:16:22 UTC (rev 58214)
+++ grass/trunk/vector/v.what.rast/local_proto.h 2013-11-14 06:52:25 UTC (rev 58215)
@@ -6,6 +6,7 @@
int count; /* nuber of points with category 'cat' */
int row;
int col;
+ double x, y; /* used with interp flag */
CELL value;
DCELL dvalue; /* used for FCELL and DCELL */
};
Modified: grass/trunk/vector/v.what.rast/main.c
===================================================================
--- grass/trunk/vector/v.what.rast/main.c 2013-11-14 03:16:22 UTC (rev 58214)
+++ grass/trunk/vector/v.what.rast/main.c 2013-11-14 06:52:25 UTC (rev 58215)
@@ -20,12 +20,10 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
-
#include <grass/raster.h>
#include <grass/dbmi.h>
#include <grass/vector.h>
#include <grass/glocale.h>
-
#include "local_proto.h"
int main(int argc, char *argv[])
@@ -33,11 +31,11 @@
int i, j, nlines, type, field, cat;
int fd;
- /* struct Categories RCats; *//* TODO */
+ /* struct Categories RCats; */ /* TODO */
struct Cell_head window;
RASTER_MAP_TYPE out_type;
- CELL *cell;
- DCELL *dcell;
+ CELL *cell_row, *prev_c_row, *next_c_row;
+ DCELL *dcell_row, *prev_d_row, *next_d_row;
int width;
int row, col;
char buf[2000];
@@ -45,7 +43,7 @@
{
struct Option *vect, *rast, *field, *col, *where;
} opt;
- struct Flag *print_flag;
+ struct Flag *interp_flag, *print_flag;
int Cache_size;
struct order *cache;
int cur_row;
@@ -96,6 +94,11 @@
opt.where = G_define_standard_option(G_OPT_DB_WHERE);
+ interp_flag = G_define_flag();
+ interp_flag->key = 'i';
+ interp_flag->description =
+ _("Interpolate values from the nearest four cells");
+
print_flag = G_define_flag();
print_flag->key = 'p';
print_flag->description =
@@ -209,6 +212,10 @@
cache[point_cnt].row = row;
cache[point_cnt].col = col;
+ if (interp_flag->answer) {
+ cache[point_cnt].x = Points->x[0];
+ cache[point_cnt].y = Points->y[0];
+ }
cache[point_cnt].cat = cat;
cache[point_cnt].count = 1;
point_cnt++;
@@ -250,10 +257,23 @@
/* Allocate space for raster row */
if (out_type == CELL_TYPE)
- cell = Rast_allocate_c_buf();
+ cell_row = Rast_allocate_c_buf();
else
- dcell = Rast_allocate_d_buf();
+ dcell_row = Rast_allocate_d_buf();
+ if (interp_flag->answer) {
+ if (out_type == CELL_TYPE) {
+ prev_c_row = Rast_allocate_c_buf();
+ next_c_row = Rast_allocate_c_buf();
+ }
+ else {
+ prev_d_row = Rast_allocate_d_buf();
+ next_d_row = Rast_allocate_d_buf();
+ }
+
+ G_begin_distance_calculations();
+ }
+
/* Extract raster values from file and store in cache */
G_debug(1, "Extracting raster values");
@@ -264,18 +284,228 @@
continue; /* duplicate cats */
if (cur_row != cache[point].row) {
- if (out_type == CELL_TYPE)
- Rast_get_c_row(fd, cell, cache[point].row);
- else
- Rast_get_d_row(fd, dcell, cache[point].row);
+ if (out_type == CELL_TYPE) {
+ Rast_get_c_row(fd, cell_row, cache[point].row);
+
+ if (interp_flag->answer) {
+ if (cache[point].row <= 0)
+ Rast_set_null_value(prev_c_row, window.cols,
+ out_type);
+ else
+ Rast_get_c_row(fd, prev_c_row, cache[point].row - 1);
+
+ if (cache[point].row + 1 > window.rows - 1)
+ Rast_set_null_value(next_c_row, window.cols,
+ out_type);
+ else
+ Rast_get_c_row(fd, next_c_row, cache[point].row + 1);
+ }
+ }
+ else {
+ Rast_get_d_row(fd, dcell_row, cache[point].row);
+
+ if (interp_flag->answer) {
+ if (cache[point].row <= 0)
+ Rast_set_null_value(prev_d_row, window.cols,
+ out_type);
+ else
+ Rast_get_d_row(fd, prev_d_row, cache[point].row - 1);
+
+ if (cache[point].row + 1 > window.rows - 1)
+ Rast_set_null_value(next_d_row, window.cols,
+ out_type);
+ else
+ Rast_get_d_row(fd, next_d_row, cache[point].row + 1);
+ }
+ }
}
cur_row = cache[point].row;
- if (out_type == CELL_TYPE) {
- cache[point].value = cell[cache[point].col];
+ if (!interp_flag->answer) {
+ if (out_type == CELL_TYPE)
+ cache[point].value = cell_row[cache[point].col];
+ else
+ cache[point].dvalue = dcell_row[cache[point].col];
}
else {
- cache[point].dvalue = dcell[cache[point].col];
+ /* do four-way IDW */
+ /* TODO: optimize, parallelize, and function-ize! */
+ double distance[4], weight[4];
+ double weightsum, valweight;
+ int col_offset, row_offset;
+
+ weightsum = valweight = 0;
+
+
+ if (cache[point].x <
+ Rast_col_to_easting(cache[point].col,
+ &window) + window.ew_res / 2)
+ col_offset = -1;
+ else
+ col_offset = +1;
+
+ if (cache[point].y >
+ Rast_row_to_northing(cache[point].row,
+ &window) - window.ns_res / 2)
+ row_offset = -1;
+ else
+ row_offset = +1;
+
+ distance[0] = G_distance(cache[point].x, cache[point].y,
+ Rast_col_to_easting(cache[point].col,
+ &window) +
+ window.ew_res / 2,
+ Rast_row_to_northing(cache[point].row,
+ &window) -
+ window.ns_res / 2);
+
+ distance[1] = G_distance(cache[point].x, cache[point].y,
+ Rast_col_to_easting(cache[point].col +
+ col_offset,
+ &window) +
+ window.ew_res / 2,
+ Rast_row_to_northing(cache[point].row,
+ &window) -
+ window.ns_res / 2);
+
+ if (row_offset == -1) {
+ distance[2] = G_distance(cache[point].x, cache[point].y,
+ Rast_col_to_easting(cache[point].
+ col + col_offset,
+ &window) +
+ window.ew_res / 2,
+ Rast_row_to_northing(cache[point].
+ row - 1,
+ &window) -
+ window.ns_res / 2);
+ distance[3] =
+ G_distance(cache[point].x, cache[point].y,
+ Rast_col_to_easting(cache[point].col,
+ &window) +
+ window.ew_res / 2,
+ Rast_row_to_northing(cache[point].row - 1,
+ &window) -
+ window.ns_res / 2);
+ }
+ else {
+ distance[2] = G_distance(cache[point].x, cache[point].y,
+ Rast_col_to_easting(cache[point].
+ col + col_offset,
+ &window) +
+ window.ew_res / 2,
+ Rast_row_to_northing(cache[point].
+ row + 1,
+ &window) -
+ window.ns_res / 2);
+ distance[3] =
+ G_distance(cache[point].x, cache[point].y,
+ Rast_col_to_easting(cache[point].col,
+ &window) +
+ window.ew_res / 2,
+ Rast_row_to_northing(cache[point].row + 1,
+ &window) -
+ window.ns_res / 2);
+ }
+
+
+ if (out_type == CELL_TYPE) {
+
+ CELL nearby_c_val[4];
+
+ /* avoid infinite weights */
+ if (distance[0] < GRASS_EPSILON) {
+ cache[point].value = cell_row[cache[point].col];
+ continue;
+ }
+
+ nearby_c_val[0] = cell_row[cache[point].col];
+
+ if (cache[point].col + col_offset < 0 ||
+ cache[point].col + col_offset >= window.cols) /* UNTESTED */
+ Rast_set_null_value(&nearby_c_val[1], 1, out_type); /* UNTESTED */
+ else
+ nearby_c_val[1] = cell_row[cache[point].col + col_offset];
+
+ if (row_offset == -1) {
+ if (cache[point].col + col_offset < 0 ||
+ cache[point].col + col_offset >= window.cols)
+ Rast_set_null_value(&nearby_c_val[2], 1, out_type);
+ else
+ nearby_c_val[2] =
+ prev_c_row[cache[point].col + col_offset];
+
+ nearby_c_val[3] = prev_c_row[cache[point].col];
+ }
+ else {
+ if (cache[point].col + col_offset < 0 ||
+ cache[point].col + col_offset >= window.cols)
+ Rast_set_null_value(&nearby_c_val[2], 1, out_type);
+ else
+ nearby_c_val[2] =
+ next_c_row[cache[point].col + col_offset];
+
+ nearby_c_val[3] = next_c_row[cache[point].col];
+ }
+
+ for (i = 0; i < 4; i++) {
+ if (!Rast_is_null_value(&nearby_c_val[i], out_type)) {
+ weight[i] = 1.0 / (distance[i] * distance[i]);
+ weightsum += weight[i];
+ valweight += weight[i] * nearby_c_val[i];
+ }
+ }
+
+ cache[point].value = valweight / weightsum;
+ }
+ else {
+ DCELL nearby_d_val[4];
+
+ /* avoid infinite weights */
+ if (distance[0] < GRASS_EPSILON) {
+ cache[point].dvalue = dcell_row[cache[point].col];
+ continue;
+ }
+
+ nearby_d_val[0] = dcell_row[cache[point].col];
+
+ if (cache[point].col + col_offset < 0 ||
+ cache[point].col + col_offset >= window.cols)
+ Rast_set_null_value(&nearby_d_val[1], 1, out_type);
+ else
+ nearby_d_val[1] =
+ dcell_row[cache[point].col + col_offset];
+
+ if (row_offset == -1) {
+ if (cache[point].col + col_offset < 0 ||
+ cache[point].col + col_offset >= window.cols)
+ Rast_set_null_value(&nearby_d_val[2], 1, out_type);
+ else
+ nearby_d_val[2] =
+ prev_d_row[cache[point].col + col_offset];
+
+ nearby_d_val[3] = prev_d_row[cache[point].col];
+ }
+ else {
+ if (cache[point].col + col_offset < 0 ||
+ cache[point].col + col_offset >= window.cols)
+ Rast_set_null_value(&nearby_d_val[2], 1, out_type);
+ else
+ nearby_d_val[2] =
+ next_d_row[cache[point].col + col_offset];
+
+ nearby_d_val[3] = next_d_row[cache[point].col];
+ }
+
+ for (i = 0; i < 4; i++) {
+ if (!Rast_is_null_value(&nearby_d_val[i], out_type)) {
+ weight[i] = 1.0 / (distance[i] * distance[i]);
+ weightsum += weight[i];
+ valweight += weight[i] * nearby_d_val[i];
+ }
+ }
+
+ cache[point].dvalue = valweight / weightsum;
+ }
}
} /* point loop */
Rast_close(fd);
Modified: grass/trunk/vector/v.what.rast/v.what.rast.html
===================================================================
--- grass/trunk/vector/v.what.rast/v.what.rast.html 2013-11-14 03:16:22 UTC (rev 58214)
+++ grass/trunk/vector/v.what.rast/v.what.rast.html 2013-11-14 06:52:25 UTC (rev 58215)
@@ -3,16 +3,21 @@
<em>v.what.rast</em> reads the raster value for each point in the vector map
and updates <b>column</b> in the vector attribute table by this value. The
column should be type number (integer, float, double, ... ).
-<br>
-If multiple points have the same category, the attribute value is set to NULL.
-If the raster value is NULL, then attribute value is set to NULL.
<p>
If the <b>-p</b> flag is used, then the attribute table is not updated
and the results are printed to <tt>stdout</tt>.
+<p>
+If the <b>-i</b> flag is used, then the value to be uploaded to the database
+is interpolated from the four nearest raster cells values using an inverse
+distance weighting method (IDW). This is useful for cases when the vector
+point density is much higher than the raster cell size.
<h2>NOTES</h2>
+If multiple points have the same category, the attribute value is set to NULL.
+If the raster value is NULL, then attribute value is set to NULL.
+<p>
<em>v.what.rast</em> operates on the attribute table. To modify the vector
geometry instead, use <em>v.drape</em>.
<p>
@@ -20,23 +25,35 @@
pipe the output of this module into the UNIX <tt>sort</tt> tool
(<tt>sort -n</tt>). If you need coordinates, after sorting use
<em>v.out.ascii</em> and the UNIX <tt>paste</tt> tool
-(<tt>paste -d'|'</tt>).
+(<tt>paste -d'|'</tt>). In the case of a NULL result, a "<tt>*</tt>"
+will be printed in lieu of the value.
+<p>
+The interpolation flag is only useful for continuous value raster maps,
+if a categorical raster is given as input the results will be nonsense.
+Since the search window is limited to four raster cells there may still
+be raster cell-edge artifacts visible in the results, this compromise
+has been made for processing speed. If one or more of the nearest four
+raster cells is NULL, then only the raster cells containing values will
+be used in the weighted average.
<h2>EXAMPLES</h2>
-A) Reading values from raster map at position of vector points, writing these values
- into a column of the attribute table connected to the vector map:
+A) Reading values from raster map at position of vector points,
+ writing these values into a column of the attribute table
+ connected to the vector map:
<br>
+
<div class="code"><pre>
v.what.rast map=pnts raster=elevation column=heights
</pre></div>
<p>
B) In case of a vector map without attached attribute table, first add
-a new attribute table. This table is then populated with values
-queried from the raster map:
+ a new attribute table. This table is then populated with values
+ queried from the raster map:
<br>
+
<div class="code"><pre>
# create new random vector points map
v.random pnts n=100
@@ -64,12 +81,14 @@
<a href="v.drape.html">v.drape</a>,
<a href="v.univar.html">v.univar</a>,
<a href="v.rast.stats.html">v.rast.stats</a>,
+v.what.rast.buffer (AddOn module),
<a href="v.what.vect.html">v.what.vect</a>
</em>
-<h2>AUTHOR</h2>
-Radim Blazek
+<h2>AUTHORS</h2>
+Radim Blazek<br>
+Hamish Bowman (interpolation)
<p>
<i>Last changed: $Date$</i>
More information about the grass-commit
mailing list