Hi
A few lines permit add a new type of interpolation in v.surf.idw, the
exponential interpolation similar to "Geographically Weighted Regression":
weigth = exp( - distance^2 / bandwidth^2)
I attach the diff files, the core is
for (n = 0; n < nsearch; n++)
{
d = exp( - list[n].dist / bw );
sum1 += list[n].z * d;
sum2 += d;
}
P.S. idw = interpolation distance weight, better than inverse distance weight
to future more type of interpolations
--
E. Jorge Tizado
--- /home/jorge/mapas/download/svn/grass6_devel/vector/v.surf.idw/description.html 2009-05-24 23:45:11.000000000 +0200
+++ description.html 2009-05-24 23:40:00.000000000 +0200
@@ -14,6 +14,10 @@
points.</p>
<P>
+When bandwidth is greater 0.0 the weight is Exponential [ exp(-d^2/bandwidth^2) ], and when is 0.0 the
+weight is Inverse Distance Squared [ 1 / d^2^power ].</p>
+
+<P>
This program allows the user to use a GRASS vector point map file,
rather than a raster map layer, as input.</p>
@@ -22,11 +26,11 @@
<p>
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
+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
+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.</p>
<p>
@@ -46,9 +50,9 @@
achieve a similar result.</p>
<p>
-If more than <EM>count</EM> points fall into one target raster cell,
+If more than <EM>count</EM> 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 <EM>count</EM>
+the -n flag is specifed, in which case only the <EM>count</EM>
points closest to the centre of the cell will be interpolated).</p>
<P>
@@ -56,11 +60,11 @@
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 < <EM>p</EM> < 1 and more smoothly for
-larger values. The default value for the power parameter is 2.
+larger values. The default value for the power parameter is 2.
</p>
<P>
-By setting <EM>npoints=1</EM>, the module can be used
+By setting <EM>npoints=1</EM>, the module can be used
to calculate raster Voronoi diagrams (Thiessen polygons).</p>
@@ -78,8 +82,8 @@
<H2>AUTHOR</H2>
-Michael Shapiro,
-U.S. Army Construction Engineering
+Michael Shapiro,
+U.S. Army Construction Engineering
Research Laboratory
<BR>
Improved algorithm (indexes points according to cell and ignores
--- /home/jorge/mapas/download/svn/grass6_devel/vector/v.surf.idw/main.c 2009-05-24 23:45:11.000000000 +0200
+++ main.c 2009-05-24 23:40:15.000000000 +0200
@@ -7,7 +7,7 @@
* 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:
+ * PURPOSE:
* COPYRIGHT: (C) 2003-2006 by the GRASS Development Team
*
* This program is free software under the GNU General Public
@@ -64,10 +64,10 @@
double dist;
double sum1, sum2, interp_value;
int n, field;
- double p;
+ double d, p, bw;
struct
{
- struct Option *input, *npoints, *power, *output, *dfield, *col;
+ struct Option *input, *npoints, *bandw, *power, *output, *dfield, *col;
} parm;
struct cell_list
{
@@ -83,8 +83,7 @@
module = G_define_module();
module->keywords = _("vector, interpolation");
module->description =
- _("Surface interpolation from vector point data by Inverse "
- "Distance Squared Weighting.");
+ _("Interpolation Distance Weighting to generate a surface from vector point data.");
parm.input = G_define_standard_option(G_OPT_V_INPUT);
@@ -98,18 +97,23 @@
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->answer = "2.0";
- parm.power->description =
- _("Power parameter; greater values assign greater influence to closer points");
+ parm.power = G_define_option();
+ parm.power->key = "power";
+ parm.power->type = TYPE_DOUBLE;
+ parm.power->answer = "2.0";
+ parm.power->description = _("Power parameter; weight = 1 / (d^2)^power");
+
+ parm.bandw = G_define_option();
+ parm.bandw->key = "bandwidth";
+ parm.bandw->type = TYPE_DOUBLE;
+ parm.bandw->answer = "0.0";
+ parm.bandw->description = _("Band width; weight = exp( - d^2/bandwidth )");
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.dfield->gisprompt = "old_layer,layer,layer_zero";
parm.col = G_define_option();
parm.col->key = "column";
@@ -144,8 +148,10 @@
list =
(struct list_Point *)G_calloc((size_t) search_points,
sizeof(struct list_Point));
-
- p = atof ( parm.power->answer );
+
+ bw = atof( parm.bandw->answer );
+ bw = bw * bw;
+ p = atof( parm.power->answer );
/* get the window, dimension arrays */
G_get_window(&window);
@@ -383,21 +389,38 @@
/* interpolate */
sum1 = 0.0;
sum2 = 0.0;
- for (n = 0; n < nsearch; n++) {
- if ((dist = list[n].dist)) {
- sum1 += list[n].z / pow ( dist, p );
- sum2 += 1.0 / pow ( dist, p );
- }
- 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;
- }
+ if (bw > 0.0)
+ {
+ for (n = 0; n < nsearch; n++)
+ {
+ d = exp( - list[n].dist / bw );
+ sum1 += list[n].z * d;
+ sum2 += d;
+ }
+ }
+ else
+ {
+ 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;
}
_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev
