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

[Moritz Lennert]

Paul,

Do you have any opinion on this ?


-------- Forwarded Message --------
Subject: Re: [GRASS-dev] [GRASS GIS] #2820: v.surf.idw results seem 
seriously wrong and don't match r.surf.idw results
Date: Tue, 15 Dec 2015 16:08:07 -0000
From: GRASS GIS <trac at osgeo.org>
Reply-To: grass-dev at lists.osgeo.org

#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>

_______________________________________________
grass-dev mailing list
grass-dev at lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev