Hi, > Thanks, Jorge and Hamish for elaborating on this. > > If I understand correctly, then Jorge's additions could be integrated > seamlessly into v.surf.idw, that is without breaking/changing its > default behaviour?
Yes. I send a the proposed files. The only apparent change of default behaviour is (now) a required value to power parameter (no 2.0 by default) to be more conscient with the real weight applied to squared distance (and it is not default for gaussian). The gaussian function is implemented with the flag -g. I remove the previous 'bandwidth' parameter and use the power parameter to do it (less changes and optimizing resources). > If so, I am in favour of integrating these additions > into SVN. > I agree that we should stick with the established terminology. I think > the best we might be to extend the module description to something > like: > > "Surface interpolation from vector point data by Inverse Distance > Squared Weighting and Radial Basis Functions" ??? I like this description. Set in main.c > Some additional details for the HTML man page would also be good. few details added to description.html Jorge > Ben > > Hamish wrote: > > Benjamin wrote: > >> A bit of web searching for "Interpolation Distance Weight" or > >> "Exponential Interpolation" (i.e. related to GIS and spatial > >> interpolation). did not produce anything meaningful for me. > >> Could you point me to some literature where I can find > >> more information on this? > > > > do a search for "radial basis functions". > > > > e.g. a quick visit google finds: > > http://en.wikipedia.org/wiki/Radial_basis_function > > http://www.cs.ualberta.ca/~sutton/book/8/node7.html > > > > http://www.ncgia.ucsb.edu/conf/SANTA_FE_CD-ROM/sf_papers/fogel_david/sant > >afe.html > > > > > > my knowledge of it is mostly from its application to interpolating 2 or > > 3 coupled-component velocity vectors between sampling sites. e.g. this > > from the Journal of Geophysical Research Sept 2006: > > > > http://www.otago.ac.nz/marinescience/po/pdfs/VennellBeatsonADCPtidalRBF20 > >06preprint.pdf > > > > > > you can play with some simple decay function variations in the v.surf.icw > > and r.surf.volcano scripts in wiki addons. > > > > e.g. 1/d^2 gives better transitions for interpolations, but 1/d^3 better > > deals with minimizing errors due to unconstrained boundaries for > > extrapolations. > > > > > > > > IMHO the well known meaning of the "IDW" acronym shouldn't be changed by > > us as people may try a search for the known meaning. Augmented, well > > have fun, especially if you can cite a validation in the literature. > > > > Personally my main hope for GRASS in this area is to see the segmentation > > code improved in the RST modules so that the "block effect" in areas of > > high point-density gradient can be automatically tuned away. ie. for > > most situations why worry much about simple IDW when you have a much more > > advanced thin plate spline solution already there as a drop-in > > replacement? > > > > oh yeah, and update r.surf.contour from its current GRASS 4-era level. > > (0 elev in input data is treated as null; integer-only output) > > > > > > Hamish > > > > > > > > > > > > _______________________________________________ > > grass-dev mailing list > > [email protected] > > http://lists.osgeo.org/mailman/listinfo/grass-dev -- E. Jorge Tizado
DESCRIPTION
v.surf.idw fills a raster matrix with interpolated values generated from a set of irregularly spaced data points using numerical approximation (weighted averaging) techniques. The interpolated value of a cell is determined by values of nearby data points and the distance of the cell from those input points. In comparison with other methods, numerical approximation allows representation of more complex surfaces (particularly those with anomalous features), restricts the spatial influence of any errors, and generates the interpolated surface from the data points.
This program allows the user to use a GRASS vector point map file, rather than a raster map layer, as input.
NOTES
The amount of memory used by this program is related to the number of vector points in the current region. If the vector point map is very dense (i.e., contains many data points), the program may not be able to get all the memory it needs from the system. The time required to execute is related to the resolution of the current region, after an initial delay determined by the time taken to read the input vector points map.
To read and interpolate from the elevation co-ordinates as 3rd dimension of the vector geometry, use layer=0. In this case no column parameter has to be specified.
If the user has a mask set, then interpolation is only done for those cells that fall within the mask. However, all vector points in the current region are used even if they fall outside the mask. Vector points outside the current region are not used in the interpolation. A larger region may be set and a mask used to limit interpolation to a smaller area if it is desired to use vector points from outside the region in the interpolation. The -n flag may also be used to achieve a similar result.
If more than count points fall into one target raster cell, the mean of all the site values will determine the cell value (unless the -n flag is specifed, in which case only the count points closest to the centre of the cell will be interpolated).
The power= parameter defines an exponential distance weight.
Greater values assign greater influence to values closer to the
point to be interpolated. The interpolation function peaks sharply over
the given data points for 0 < p < 1 and more smoothly for
larger values.
With default mode (Inverse Squared Distance: weight = 1.0/[dist^2^power] );
the recommended value for the power parameter is 2.
With flag -g (Gaussian Distance: weight = exp[-dist^2/power^2] );
the value for the power parameter is a distance.
By setting npoints=1, the module can be used to calculate raster Voronoi diagrams (Thiessen polygons).
SEE ALSO
d.vectg.region
r.surf.contour
r.surf.idw
r.surf.idw2
r.surf.gauss
r.surf.fractal
r.surf.random
v.surf.rst
AUTHOR
Michael Shapiro, U.S. Army Construction Engineering Research LaboratoryImproved algorithm (indexes points according to cell and ignores points outside current region) by Paul Kelly
Last changed: $Date: 2009-04-03 18:40:07 +0200 (Fri, 03 Apr 2009) $
/****************************************************************************
*
* MODULE: v.surf.idw
* AUTHOR(S): Michael Shapiro, U.S. Army Construction Engineering Research Laboratory
* Improved algorithm (indexes points according to cell and ignores
* points outside current region) by Paul Kelly
* further: Radim Blazek <radim.blazek gmail.com>, Huidae Cho <grass4u gmail.com>,
* Glynn Clements <glynn gclements.plus.com>, Markus Neteler <neteler itc.it>
* PURPOSE:
* COPYRIGHT: (C) 2003-2006 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
* for details.
*
*****************************************************************************/
#include <stdlib.h>
#include <math.h>
#include <grass/gis.h>
#include <grass/glocale.h>
int search_points = 12;
long npoints = 0;
long **npoints_currcell;
int nsearch;
static int i;
struct Point
{
double north, east;
double z;
};
struct list_Point
{
double north, east;
double z;
double dist;
};
struct Point ***points;
struct Point *noidxpoints = NULL;
struct list_Point *list;
static struct Cell_head window;
static struct Flag *noindex;
void calculate_distances(int, int, double, double, int *);
void calculate_distances_noindex(double, double);
/* read_sites.c */
void read_sites(char *, int, char *);
int main(int argc, char *argv[])
{
int fd, maskfd;
CELL *mask;
DCELL *dcell;
struct GModule *module;
struct History history;
int row, col;
int searchrow, searchcolumn, pointsfound;
int *shortlistrows = NULL, *shortlistcolumns = NULL;
long ncells = 0;
double north, east;
double dist;
double sum1, sum2, interp_value;
int n, field;
double d, p, bw;
struct
{
struct Option *input, *npoints, *power, *output, *dfield, *col;
} parm;
struct Flag *gauss, *multi, *poly;
struct cell_list
{
int row, column;
struct cell_list *next;
};
struct cell_list **search_list = NULL, **search_list_start = NULL;
int max_radius, radius;
int searchallpoints = 0;
G_gisinit(argv[0]);
module = G_define_module();
module->keywords = _("vector, interpolation");
module->description = _("Surface interpolation from vector point data by Inverse Distance "
"Squared Weighting and Gaussian Radial Basis Functions");
parm.input = G_define_standard_option(G_OPT_V_INPUT);
parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
parm.npoints = G_define_option();
parm.npoints->key = "npoints";
parm.npoints->key_desc = "count";
parm.npoints->type = TYPE_INTEGER;
parm.npoints->required = NO;
parm.npoints->description = _("Number of interpolation points");
parm.npoints->answer = "12";
parm.power = G_define_option();
parm.power->key = "power";
parm.power->type = TYPE_DOUBLE;
parm.power->required = YES;
parm.power->description = _("Power parameter");
parm.dfield = G_define_standard_option(G_OPT_V_FIELD);
parm.dfield->description = _("If set to 0, z coordinates are used (3D vector only)");
parm.dfield->answer = "1";
parm.dfield->gisprompt = "old_layer,layer,layer_zero";
parm.col = G_define_option();
parm.col->key = "column";
parm.col->type = TYPE_STRING;
parm.col->required = NO;
parm.col->label = _("Attribute table column with values to interpolate");
parm.col->description = _("Required if layer > 0");
gauss = G_define_flag();
gauss->key = 'g';
gauss->description = _("Gaussian weights function [ exp(-dist^2/power^2) ]");
gauss->answer = 0;
noindex = G_define_flag();
noindex->key = 'n';
noindex->label = _("Don't index points by raster cell");
noindex->description = _("Slower but uses"
" less memory and includes points from outside region"
" in the interpolation");
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
if (G_legal_filename(parm.output->answer) < 0)
G_fatal_error(_("<%s> is an illegal file name"), parm.output->answer);
if (sscanf(parm.npoints->answer, "%d", &search_points) != 1 ||
search_points < 1)
G_fatal_error(_("%s=%s - illegal number of interpolation points"),
parm.npoints->key, parm.npoints->answer);
sscanf(parm.dfield->answer, "%d", &field);
if (field > 0 && !parm.col->answer)
G_fatal_error(_("No attribute column specified"));
list =
(struct list_Point *)G_calloc((size_t) search_points,
sizeof(struct list_Point));
p = atof( parm.power->answer );
if (gauss->answer)
p *= p;
/* get the window, dimension arrays */
G_get_window(&window);
if (!noindex->answer) {
npoints_currcell = (long **)G_malloc(window.rows * sizeof(long *));
points =
(struct Point ***)G_malloc(window.rows * sizeof(struct Point **));
for (row = 0; row < window.rows; row++) {
npoints_currcell[row] =
(long *)G_malloc(window.cols * sizeof(long));
points[row] =
(struct Point **)G_malloc(window.cols *
sizeof(struct Point *));
for (col = 0; col < window.cols; col++) {
npoints_currcell[row][col] = 0;
points[row][col] = NULL;
}
}
}
/* read the elevation points from the input sites file */
read_sites(parm.input->answer, field, parm.col->answer);
if (npoints == 0)
G_fatal_error(_("No data points found"));
nsearch = npoints < search_points ? npoints : search_points;
if (!noindex->answer) {
/* Arbitrary point to switch between searching algorithms. Could do
* with refinement PK */
if ((window.rows * window.cols) / npoints > 400) {
/* Using old algorithm.... */
searchallpoints = 1;
ncells = 0;
/* Make an array to contain the row and column indices that have
* sites in them; later will just search through all these. */
for (searchrow = 0; searchrow < window.rows; searchrow++)
for (searchcolumn = 0; searchcolumn < window.cols;
searchcolumn++)
if (npoints_currcell[searchrow][searchcolumn] > 0) {
shortlistrows = (int *)G_realloc(shortlistrows,
(1 +
ncells) *
sizeof(int));
shortlistcolumns =
(int *)G_realloc(shortlistcolumns,
(1 + ncells) * sizeof(int));
shortlistrows[ncells] = searchrow;
shortlistcolumns[ncells] = searchcolumn;
ncells++;
}
}
else {
/* Fill look-up table of row and column offsets for
* doing a circular region growing search looking for sites */
/* Use units of column width */
max_radius = (int)(0.5 + sqrt(window.cols * window.cols +
(window.rows * window.ns_res /
window.ew_res) * (window.rows *
window.ns_res /
window.ew_res)));
search_list =
(struct cell_list **)G_malloc(max_radius *
sizeof(struct cell_list *));
search_list_start =
(struct cell_list **)G_malloc(max_radius *
sizeof(struct cell_list *));
for (radius = 0; radius < max_radius; radius++)
search_list[radius] = NULL;
for (row = 0; row < window.rows; row++)
for (col = 0; col < window.cols; col++) {
radius = (int)sqrt(col * col +
(row * window.ns_res / window.ew_res) *
(row * window.ns_res / window.ew_res));
if (search_list[radius] == NULL)
search_list[radius] =
search_list_start[radius] =
G_malloc(sizeof(struct cell_list));
else
search_list[radius] =
search_list[radius]->next =
G_malloc(sizeof(struct cell_list));
search_list[radius]->row = row;
search_list[radius]->column = col;
search_list[radius]->next = NULL;
}
}
}
/* allocate buffers, etc. */
dcell = G_allocate_d_raster_buf();
if ((maskfd = G_maskfd()) >= 0)
mask = G_allocate_cell_buf();
else
mask = NULL;
fd = G_open_raster_new(parm.output->answer, DCELL_TYPE);
if (fd < 0)
G_fatal_error(_("Unable to create raster map <%s>"),
parm.output->answer);
G_important_message(_("Interpolating raster map <%s> (%d rows, %d cols)... "),
parm.output->answer, window.rows, window.cols);
north = window.north + window.ns_res / 2.0;
for (row = 0; row < window.rows; row++) {
G_percent(row, window.rows - 1, 1);
if (mask) {
if (G_get_map_row(maskfd, mask, row) < 0)
exit(1);
}
north -= window.ns_res;
east = window.west - window.ew_res / 2.0;
for (col = 0; col < window.cols; col++) {
east += window.ew_res;
/* don't interpolate outside of the mask */
if (mask && mask[col] == 0) {
G_set_d_null_value(&dcell[col], 1);
continue;
}
/* If current cell contains more than nsearch points just average
* all the points in this cell and don't look in any others */
if (!(noindex->answer) && npoints_currcell[row][col] >= nsearch) {
sum1 = 0.0;
for (i = 0; i < npoints_currcell[row][col]; i++)
sum1 += points[row][col][i].z;
interp_value = sum1 / npoints_currcell[row][col];
}
else {
if (noindex->answer)
calculate_distances_noindex(north, east);
else {
pointsfound = 0;
i = 0;
if (searchallpoints == 1) {
/* If there aren't many sites just check them all to find
* the nearest */
for (n = 0; n < ncells; n++)
calculate_distances(shortlistrows[n],
shortlistcolumns[n], north,
east, &pointsfound);
}
else {
radius = 0;
while (pointsfound < nsearch) {
/* Keep widening the search window until we find
* enough points */
search_list[radius] = search_list_start[radius];
while (search_list[radius] != NULL) {
/* Always */
if (row <
(window.rows - search_list[radius]->row)
&& col <
(window.cols -
search_list[radius]->column)) {
searchrow =
row + search_list[radius]->row;
searchcolumn =
col + search_list[radius]->column;
calculate_distances(searchrow,
searchcolumn, north,
east, &pointsfound);
}
/* Only if at least one offset is not 0 */
if ((search_list[radius]->row > 0 ||
search_list[radius]->column > 0) &&
row >= search_list[radius]->row &&
col >= search_list[radius]->column) {
searchrow =
row - search_list[radius]->row;
searchcolumn =
col - search_list[radius]->column;
calculate_distances(searchrow,
searchcolumn, north,
east, &pointsfound);
}
/* Only if both offsets are not 0 */
if (search_list[radius]->row > 0 &&
search_list[radius]->column > 0) {
if (row <
(window.rows -
search_list[radius]->row) &&
col >= search_list[radius]->column) {
searchrow =
row + search_list[radius]->row;
searchcolumn =
col - search_list[radius]->column;
calculate_distances(searchrow,
searchcolumn,
north, east,
&pointsfound);
}
if (row >= search_list[radius]->row &&
col <
(window.cols -
search_list[radius]->column)) {
searchrow =
row - search_list[radius]->row;
searchcolumn =
col + search_list[radius]->column;
calculate_distances(searchrow,
searchcolumn,
north, east,
&pointsfound);
}
}
search_list[radius] =
search_list[radius]->next;
}
radius++;
}
}
}
/* interpolate */
sum1 = 0.0;
sum2 = 0.0;
if (gauss->answer)
{
for (n = 0; n < nsearch; n++)
{
d = exp( - list[n].dist / p );
sum1 += list[n].z * d;
sum2 += d;
}
}
else /* default: inverse squared distance */
{
for (n = 0; n < nsearch; n++)
{
if ((dist = list[n].dist))
{
d = pow ( dist, p );
sum1 += list[n].z / d;
sum2 += 1.0 / d;
}
else
{
/* If one site is dead on the centre of the cell, ignore
* all the other sites and just use this value.
* (Unlikely when using floating point numbers?) */
sum1 = list[n].z;
sum2 = 1.0;
break;
}
}
}
interp_value = sum1 / sum2;
}
dcell[col] = (DCELL) interp_value;
}
G_put_d_raster_row(fd, dcell);
}
G_close_cell(fd);
/* writing history file */
G_short_history(parm.output->answer, "raster", &history);
G_command_history(&history);
G_write_history(parm.output->answer, &history);
G_done_msg("");
exit(EXIT_SUCCESS);
}
void newpoint(double z, double east, double north)
{
int row, column;
row = (int)((window.north - north) / window.ns_res);
column = (int)((east - window.west) / window.ew_res);
if (!noindex->answer) {
if (row < 0 || row >= window.rows || column < 0 ||
column >= window.cols) ;
else { /* Ignore sites outside current region as can't be indexed */
points[row][column] =
(struct Point *)G_realloc(points[row][column],
(1 +
npoints_currcell[row][column]) *
sizeof(struct Point));
points[row][column][npoints_currcell[row][column]].north = north;
points[row][column][npoints_currcell[row][column]].east = east;
points[row][column][npoints_currcell[row][column]].z = z;
npoints_currcell[row][column]++;
npoints++;
}
}
else {
noidxpoints = (struct Point *)G_realloc(noidxpoints,
(1 +
npoints) *
sizeof(struct Point));
noidxpoints[npoints].north = north;
noidxpoints[npoints].east = east;
noidxpoints[npoints].z = z;
npoints++;
}
}
void calculate_distances(int row, int column, double north,
double east, int *pointsfound)
{
int j, n, max = 0;
double dx, dy, dist;
static double maxdist;
/* Check distances and find the points to use in interpolation */
for (j = 0; j < npoints_currcell[row][column]; j++) {
/* fill list with first nsearch points */
if (i < nsearch) {
dy = points[row][column][j].north - north;
dx = points[row][column][j].east - east;
list[i].dist = dy * dy + dx * dx;
list[i].z = points[row][column][j].z;
i++;
/* find the maximum distance */
if (i == nsearch) {
maxdist = list[max = 0].dist;
for (n = 1; n < nsearch; n++) {
if (maxdist < list[n].dist)
maxdist = list[max = n].dist;
}
}
}
else {
/* go thru rest of the points now */
dy = points[row][column][j].north - north;
dx = points[row][column][j].east - east;
dist = dy * dy + dx * dx;
if (dist < maxdist) {
/* replace the largest dist */
list[max].z = points[row][column][j].z;
list[max].dist = dist;
maxdist = list[max = 0].dist;
for (n = 1; n < nsearch; n++) {
if (maxdist < list[n].dist)
maxdist = list[max = n].dist;
}
}
}
}
*pointsfound += npoints_currcell[row][column];
}
void calculate_distances_noindex(double north, double east)
{
int n, max;
double dx, dy, dist;
double maxdist;
/* fill list with first nsearch points */
for (i = 0; i < nsearch; i++) {
dy = noidxpoints[i].north - north;
dx = noidxpoints[i].east - east;
list[i].dist = dy * dy + dx * dx;
list[i].z = noidxpoints[i].z;
}
/* find the maximum distance */
maxdist = list[max = 0].dist;
for (n = 1; n < nsearch; n++) {
if (maxdist < list[n].dist)
maxdist = list[max = n].dist;
}
/* go thru rest of the points now */
for (; i < npoints; i++) {
dy = noidxpoints[i].north - north;
dx = noidxpoints[i].east - east;
dist = dy * dy + dx * dx;
if (dist < maxdist) {
/* replace the largest dist */
list[max].z = noidxpoints[i].z;
list[max].dist = dist;
maxdist = list[max = 0].dist;
for (n = 1; n < nsearch; n++) {
if (maxdist < list[n].dist)
maxdist = list[max = n].dist;
}
}
}
}
_______________________________________________ grass-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/grass-dev
