#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@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev