[GRASS-dev] [GRASS GIS] #2820: v.surf.idw results seem seriously wrong and don't match r.surf.idw results

[GRASS GIS]

#2820: v.surf.idw results seem seriously wrong and don't match r.surf.idw results
----------------------+-------------------------------------------------
  Reporter:  lrntct   |      Owner:  grass-dev@…
      Type:  defect   |     Status:  new
  Priority:  blocker  |  Milestone:  7.0.3
 Component:  Vector   |    Version:  7.0.1
Resolution:           |   Keywords:  v.surf.idw r.surf.idw interpolation
       CPU:  x86-64   |   Platform:  Linux
----------------------+-------------------------------------------------

Comment (by mlennert):

 Replying to [comment:6 mlennert]:
 > Replying to [comment:5 mlennert]:
 > > The result is quite unambiguous:
 > >
 > >
 > > {{{
 > > cat|idw_vect|idw_rast|owncal
 > > 1|1545.09361138162|1388.6|1389.27
 > > }}}
 > >
 > > i.e. the calculated value is almost exactly identical to the
 r.surf.idw result, while the v.surf.idw result is very far off.
 > >
 >
 > Using the '-n' flag in v.surf.index seems to solve the issue at least
 partly: the value for the above point is now 1421, i.e. much closer to the
 r.surf.idw and the calculated result (although still significantly
 different).
 >
 > More exploration needed...

 Ok, found part of the problem: distance was always left as dy*dy + dx*dx
 without ever taking the square root. This obviously changes the weights as
 1/10 is closer to 1/100 than 1/100 is to 1/10000...

 The following patch solves this issue. Now I get the same result using the
 -n flag as I do by calculation and with r.surf.idw. Without the -n flag I
 now have a value of 1484 instead of the 1545 above. Indexing thus still
 causes some issues. I think that these should be explored before applying
 below patch. As an intermediate solution we could apply this immediately
 and either reverse the meaning of the -n flag, making no indexing the
 default, or we just completely take away the -n flag and always work
 without indexing.


 {{{
 Index: vector/v.surf.idw/main.c
 ===================================================================
 --- vector/v.surf.idw/main.c    (révision 67143)
 +++ vector/v.surf.idw/main.c    (copie de travail)
 @@ -468,7 +468,7 @@
         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].dist = sqrt(dy * dy + dx * dx);
             list[i].z = points[row][column][j].z;
             i++;

 @@ -486,7 +486,7 @@
             /* 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;
 +           dist = sqrt(dy * dy + dx * dx);

             if (dist < maxdist) {
                 /* replace the largest dist */
 @@ -514,7 +514,7 @@
      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].dist = sqrt(dy * dy + dx * dx);
         list[i].z = noidxpoints[i].z;
      }
      /* find the maximum distance */
 @@ -527,7 +527,7 @@
      for (; i < npoints; i++) {
         dy = noidxpoints[i].north - north;
         dx = noidxpoints[i].east - east;
 -       dist = dy * dy + dx * dx;
 +       dist = sqrt(dy * dy + dx * dx);

         if (dist < maxdist) {
             /* replace the largest dist */
 }}}

 As Paul is the one who introduced indexing, maybe he wants to comment...

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/2820#comment:7>
GRASS GIS <https://grass.osgeo.org>