r.surf.idw and v.surf.idw results seem seriously wrong

[lrntct]

I've encountered some strange problem with r.surf.idw while trying to interpolate rainfall data. That is the procedure:

• import raster points with r.in.xyz
• convert to integer with r.mapcalc, multiplying value by 10⁶
• run r.surf.idw
• convert back to float with r.mapcalc
• convert original raster points to vector using r.to.vect
• run v.surf.idw
• set same color table to both rasters

Images of the results are attached

[lrntct]

results from r.surf.idw

[wenzeslaus]

Can you please provide a simple shell script with the steps so that this can be easily repeated?

[mlennert]

I can reproduce the issue in the NC data set:

g.region vect=nc_state res=500 -a
v.surf.idw precip_30ynormals col=annual out=vsurfidw --o
v.to.rast precip_30ynormals use=attr attr_col=annual out=precip --o
r.mapcalc "precip_int = int(precip * 100)" --o
r.surf.idw in=precip_int out=rsurfidw_tmp --o
r.mapcalc "rsurfidw = float(rsurfidw_tmp) / 100" --o
r.mapcalc "diff = rsurfidw - vsurfidw"

Analysing the differences, you get:

r.univar map=diff                                                           
total null and non-null cells: 24840200
total null cells: 0
Of the non-null cells:
----------------------
n: 24840200
minimum: -370.28
maximum: 1020.98
range: 1391.26
mean: 48.8595
mean of absolute values: 127.803
standard deviation: 180.027
variance: 32409.7
variation coefficient: 368.458 %
sum: 1213680178.7663

and

r.covar -r map=rsurfidw,vsurfidw
r.covar: complete ...
N = 24840200
1.000000 0.284021 
0.284021 1.000000

Looking at the individual result maps:

r.univar map=rsurfidw
total null and non-null cells: 24840200
total null cells: 0
Of the non-null cells:
----------------------
n: 24840200
minimum: 947.42
maximum: 2329.17
range: 1381.75
mean: 1346.31
mean of absolute values: 1346.31
standard deviation: 183.013
variance: 33493.9
variation coefficient: 13.5937 %
sum: 33442685677.211

and

r.univar map=vsurfidw
total null and non-null cells: 24840200
total null cells: 0
Of the non-null cells:
----------------------
n: 24840200
minimum: 1120.14
maximum: 1546.86
range: 426.72
mean: 1297.45
mean of absolute values: 1297.45
standard deviation: 92.1991
variance: 8500.67
variation coefficient: 7.10616 %
sum: 32229005498.4201

I will also attach some maps, including a map of the original data points. The interpolated maps use the following color table:

947.42 255:255:255
2329.17 0:0:255

v.surf.idw smoothes the map much more...

[mlennert]

maps of NC data test

[neteler]

Looking at the "r.surf.idw.png​" attachment: could it be that some integer rounding happens (just guessing)?

[lrntct]

Here is a way to reproduce the problem:

v.in.ascii input=rainfall.csv output=stations_vect separator=comma z=3
g.region vector=stations_vect res=100
v.surf.idw input=stations_vect output=vsurfidw
v.to.rast input=stations_vect output=stations_rast use=z
r.mapcalc expression='stations_int=int(stations_rast*1000000)'
r.surf.idw input=stations_int output=rsurfidw_int
r.mapcalc expression='rsurfidw=double(rsurfidw_int)/1000000'
r.mapcalc expression='rain_diff=rsurfidw-vsurfidw'
r.colors map=rain_diff color=aspect

The input data:

403368.495212,1866849.929277,0
420183.434764,1853411.873576,0.25
405991.733962,1864563.901667,0.206855792
481344.935675,1946775.612706,0
446846.743957,1940506.802136,0
474982.242862,1851893.486232,0.2413317573
503110.946912,1836822.953657,0.4334121355
333235.275253,1907773.128426,0.0344759653
451811.619719,1933211.337411,0
427242.101414,1874036.5345,0.3940110323

[mlennert]

The more I look at this issue, the more I get the feeling that v.surf.idw is wrong as well. At least, I cannot really explain the output for the NC precipitation data shown in ​idw_tests_NCdata.png.

Here's a test using the NC precipitation data interpolated to maps rsurfidw and vsurfidw according to above instructions (tested in trunk):

echo "673491|55508" | v.in.ascii in=- out=test_interpolation
v.db.addtable test_interpolation
v.db.addcolumn test_interpolation col="idw_vect double precision"
v.db.addcolumn test_interpolation col="idw_rast double precision"
v.db.addcolumn test_interpolation col="owncal double precision"
v.what.rast test_interpolation col=idw_vect rast=vsurfidw
v.what.rast test_interpolation col=idw_rast rast=rsurfidw
v.db.update test_interpolation col=owncal value=$(v.distance -p -a from=test_interpolation@user1 to=precip_30ynormals@PERMANENT dmax=102000 upload=dist,to_attr to_column=annual --quiet | tail -n +2 | awk  -F'|' 'BEGIN{sumvals=0;sumweights=0} {weight=1/$3^2;sumweights+=weight; sumvals+=weight*$4} END{print sumvals/sumweights}')
v.db.select test_interpolation

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.

I would guess the segmentation observed by the OP with r.surf.idw is due to integer rounding issues (although this should be explored), but in the current state, v.surf.idw needs a serious check to see where the error comes from !

I'm renaming this ticket for now, but maybe it needs to be split in two... Bumping it up as a blocker, since this is basic functionality of a GIS which should "just work"...

[mlennert]

Replying to 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...

[mlennert]

Replying to mlennert:

Replying to 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...

[pkelly]

Hi Moritz, Thanks for bringing this to my attention. Looks like the bug with the squared distance was introduced when support for powers other than 2 was added in r32957, and ported to GRASS 7 in r59107.

Regarding the indexing, it's certainly true that it can give slightly different results depending on the layout of the points, and looking back maybe it shouldn't have been made the default mode of operation because of these differences. E.g. one obvious issue is that only points inside the current region are included in the interpolation, whereas without indexing all points are potentially included even if they lie outside the region. And I guess it's still possible there could be a bug in it somewhere as my programming skills weren't quite so well developed back then either.

I haven't been able to test the example in the Trac ticket yet as I don't have GRASS 7 currently installed, only 6.4, and the GRASS6 version of the NC dataset at <https://grass.osgeo.org/sampledata/north_carolina/nc_spm_latest.tar.gz>, when uncompressed, seems to be actually the same as the GRASS7 version.

My original interpolation use case was image pixels from a ground-level camera view that had been perspective-transformed to an overhead view, resulting in highly uneven point density across the region. Without the indexing s.surf.idw (as it was then) could take hours to run. The following extract from my PhD thesis gives some of the background: ​http://www.stjohnspoint.co.uk/gis/idw.pdf (although at the Firefox built-in PDF viewer seems to have trouble displaying it; it's converted from the original PostScript).

One thing I would consider testing first if I had the data, would be to enlarge the region to slightly bigger than the edges of the point cloud (i.e. slightly bigger than the result of g.region vect=stations_vect), just in case there were any rounding issues with points right at the edge falling out of the region and not being included in the indexing.

Paul

[mlennert]

Replying to pkelly:

Hi Moritz, Thanks for bringing this to my attention. Looks like the bug with the squared distance was introduced when support for powers other than 2 was added in r32957, and ported to GRASS 7 in r59107.

Indeed. Before that, no need to take the square root if it was to only square the distance afterwards...

I haven't been able to test the example in the Trac ticket yet as I don't have GRASS 7 currently installed, only 6.4, and the GRASS6 version of the NC dataset at <https://grass.osgeo.org/sampledata/north_carolina/nc_spm_latest.tar.gz>, when uncompressed, seems to be actually the same as the GRASS7 version.

Yes, that's a bug. Running v.build.all in the PERMANENT mapset solves that easily, though.

One thing I would consider testing first if I had the data, would be to enlarge the region to slightly bigger than the edges of the point cloud (i.e. slightly bigger than the result of g.region vect=stations_vect), just in case there were any rounding issues with points right at the edge falling out of the region and not being included in the indexing.

The interpolation was done in a region defined by the nc_state boundary:

g.region vect=nc_state res=500 -ap
projection: 99 (Lambert Conformal Conic)
zone:       0
datum:      nad83
ellipsoid:  a=6378137 es=0.006694380022900787
north:      318500
south:      10500
west:       124000
east:       930500
nsres:      500
ewres:      500
rows:       616
cols:       1613
cells:      993608

Setting the region with station data gives:

g.region vect=precip_30ynormals res=500 -a
projection: 99 (Lambert Conformal Conic)
zone:       0
datum:      nad83
ellipsoid:  a=6378137 es=0.006694380022900787
north:      306500
south:      27500
west:       151500
east:       917500
nsres:      500
ewres:      500
rows:       558
cols:       1532
cells:      854856

so, I don't think that this is the issue, here.

Moritz

[mlennert]

Complete testing script for GRASS6

[mlennert]

Complete testing script for GRASS7

[mlennert]

I've just added complete testing scripts for GRASS6 and GRASS7.

[mlennert]

Replying to mlennert:

Replying to pkelly:

I haven't been able to test the example in the Trac ticket yet as I don't have GRASS 7 currently installed, only 6.4, and the GRASS6 version of the NC dataset at <https://grass.osgeo.org/sampledata/north_carolina/nc_spm_latest.tar.gz>, when uncompressed, seems to be actually the same as the GRASS7 version.

Yes, that's a bug. Running v.build.all in the PERMANENT mapset solves that easily, though.

See #2829.

[neteler]

Replying to mlennert:

Replying to pkelly:

I haven't been able to test the example in the Trac ticket yet as I don't have GRASS 7 currently installed, only 6.4, and the GRASS6 version of the NC dataset at <https://grass.osgeo.org/sampledata/north_carolina/nc_spm_latest.tar.gz>, when uncompressed, seems to be actually the same as the GRASS7 version.

Sorry for the mess. Fixed now (see #2829).

[mlennert]

Investigating further. Adding the following patch for diagnosis (in addition to the above sqrt patch):

Index: main.c			
===================================================================		
--- main.c	(révision 67166)		
+++ main.c	(copie de travail)		
@@ -386,6 +386,7 @@			
 		sum1 = 0.0;	
 		sum2 = 0.0;	
 		for (n = 0; n < nsearch; n++) {	
+		    fprintf(stdout, "%f,%f,%f,%f\n", east, north, list[n].dist, list[n].z);	
 		    if ((dist = list[n].dist)) {	
 			sum1 += list[n].z / pow(dist, p);
 			sum2 += 1.0 / pow(dist, p);

and then running the following commands (i.e. creating one vector point and getting distance and annual precipitation values from the 12 nearest stations to that point, as well as the same values for the equivalent pixel in both runs of v.surf.idw (without and with -n flag):

g.region vect=nc_state res=500 -ap
echo "567250|214750" | v.in.ascii in=- out=center --o
v.distance -ap from=center to=precip_30ynormals up=dist,to_attr to_col=annual dmax=64000 --q | tail -n -12
v.surf.idw precip_30ynormals col=annual out=vsurfidw --o > list.csv
v.surf.idw -n precip_30ynormals col=annual out=vsurfidw_noindex --o > list.csv
grep "567250.000000,214750" list.csv > list_center.csv
grep "567250.000000,214750" list_noindex.csv > list_center_noindex.csv

and putting together the results of the v.distance, v.surf.idw and v.surf.idw -n runs, I get (results ordered by increasing distance):

• v.distance

distance                value
8400,5853273955	        1224,28
33001,5791065737	1183,64
33543,2163702733	1183,64
39603,9217121153	1234,44
41429,8528851821	1158,24
41667,8156183992	1206,5
41677,4928801831	1143
42993,6439600215	1209,04
55288,6854942564	1170,94
58543,7865733219	1216,66
60471,6142897	        1219,2
62763,095429678	        1094,74

• v.surf.idw

distance        value
80047,835648	1191,26
90264,15377	1183,64
91013,285596	1158,24
99493,013889	1219,2
120304,207816	1120,14
171007,225542	1358,9
183102,228906	1440,18
184179,357051	1170,94
194397,039723	1234,44
200403,915046	1163,32
230189,059243	1546,86
231312,997387	1259,84

• v.surf.idw -n

8400,585327	1224,28
33001,579107	1183,64
33543,21637	1183,64
39603,921712	1234,44
41429,852885	1158,24
41667,815618	1206,5
41677,49288	1143
42993,64396	1209,04
55288,685494	1170,94
58543,786573	1216,66
60471,61429	1219,2
62763,09543	1094,74

In other words, v.surf.idw -n gives the exact same values as v.distance, but v.surf.idw without the flag gives a completely different set of distances. So, I have the feeling that the issue is in the ​calculate_distances() function or possibly in the creation of the ​shortlist of cells to look in in main.c, but I haven't been able to pinpoint the culprit, yet.

Any help is more than welcome.

[annakrat]

I hopefully fixed the bug in the indexing part and the missing sqrt in r67211. Please test. Haven't looked at r.surf.idw yet.

[annakrat]

Replying to pkelly:

Hi Moritz, Thanks for bringing this to my attention. Looks like the bug with the squared distance was introduced when support for powers other than 2 was added in r32957, and ported to GRASS 7 in r59107. My original interpolation use case was image pixels from a ground-level camera view that had been perspective-transformed to an overhead view, resulting in highly uneven point density across the region. Without the indexing s.surf.idw (as it was then) could take hours to run. The following extract from my PhD thesis gives some of the background: ​http://www.stjohnspoint.co.uk/gis/idw.pdf (although at the Firefox built-in PDF viewer seems to have trouble displaying it; it's converted from the original PostScript).

Paul

Perhaps you could summarize your algorithm and we can put it in the manual. Currently there is almost nothing about the -n flag.

[mlennert]

Replying to annakrat:

I hopefully fixed the bug in the indexing part and the missing sqrt in r67211. Please test. Haven't looked at r.surf.idw yet.

Beautiful, thank you very much !

v.surf.idw results are now identical with and without the -n flag.

For the NC precipitation example, setting the flag (i.e. not using the indexing method) actually make the module run faster. This merits some more testing. From what I can gather from the PhD thesis and from the code, it seems to me (but I'm really not sure) that the indexing is mostly useful for situations where you have several input points per raster cell. But in this case, why use interpolation and not aggregate statistics per cell (à la r.in.xyz) ?

As for r.surf.idw, the NC precipitation data test gives almost identical results (r.covar -r between v.surf.idw and r.surf.idw results gives 0.999883). There are some very localized differences that might merit a look, but they might just be slight differences in distances calculations and rounding.

However, using the OP's testing data, I still see the issue with the r.surf.idw results.

In my eyes we should backport the fix to v.surf.idw and lower the priority of this bug.

[mlennert]

One remark on v.surf.idw: IIUC, the current distance calculations only work in projected locations. It might be worth introducing geodesic distances...

[neteler]

Replying to mlennert:

One remark on v.surf.idw: IIUC, the current distance calculations only work in projected locations. It might be worth introducing geodesic distances...

FWIW, in v.distance the function Vect_line_geodesic_distance() and related are used for that.

[pkelly]

Replying to annakrat:

Perhaps you could summarize your algorithm and we can put it in the manual. Currently there is almost nothing about the -n flag.

Done in r67230. I guess it's worth considering whether -n should really be the default after all...

And thanks a lot for finding and fixing the bug. I am sorry to have left such a bad bug in GRASS for so many years. Can I maybe just suggest that a slightly neater way of implementing the fix (and more in keeping with the antiquated style of the existing C code!) would be to simply make the variable max static within the function in a similar way as is already done for maxdist, like this:

--- main.c~	2015-12-15 20:06:09.047518721 +0000
+++ main.c	2015-12-18 19:46:01.559351297 +0000
@@ -458,7 +458,8 @@
 void calculate_distances(int row, int column, double north,
 			 double east, int *pointsfound)
 {
-    int j, n, max = 0;
+    int j, n;
+    static int max;
     double dx, dy, dist;
     static double maxdist;

I will say that I haven't tested that, but I think it should work and would be neater. But really, great work finding that as I know the code is hard to follow.

Replying to mlennert:

From what I can gather from the PhD thesis and from the code, it seems to me (but I'm really not sure) that the indexing is mostly useful for situations where you have several input points per raster cell. But in this case, why use interpolation and not aggregate statistics per cell (à la r.in.xyz) ?

Where it is really useful, is where the data is dense in some parts of the region (and should be aggregated) and sparse in others (and should be interpolated). I guess maybe it is a pretty unusual situation, though. The issue with multiple input points falling within one output cell doesn't arise in r.surf.idw, since the distribution of the input points is limited by the raster resolution.

[mlennert]

Replying to pkelly:

I guess it's worth considering whether -n should really be the default after all...

I think we would need a series of tests with very different use cases to decide. At least now, results are identical (but probably won't be in the case of very dense points within a cell), so at least it is mostly only a performance question.

And thanks a lot for finding and fixing the bug. I am sorry to have left such a bad bug in GRASS for so many years. Can I maybe just suggest that a slightly neater way of implementing the fix (and more in keeping with the antiquated style of the existing C code!) would be to simply make the variable max static within the function in a similar way as is already done for maxdist, like this:

Implemented in r67241.

I've also taken the liberty of backporting the fixes to release70 in r67242.

Downgrading priority of the bug back to normal, as now we can see why r.surf.idw gives these weirdly segmented results in the OP's test case.

[wenzeslaus]

Replying to pkelly:

and more in keeping with the antiquated style of the existing C code!

Sorry, but using the same logic we should also introduce some cryptic names and perhaps memory leaks. What we should value in relation to code is readability and reusability. What is means in this particular case might be debatable, but in general we should definitively prefer readable code rather than antiquated style.

[mlennert]

Replying to wenzeslaus:

Replying to pkelly:

and more in keeping with the antiquated style of the existing C code!

Sorry, but using the same logic we should also introduce some cryptic names and perhaps memory leaks. What we should value in relation to code is readability and reusability. What is means in this particular case might be debatable, but in general we should definitively prefer readable code rather than antiquated style.

I agree with you in principle, although when you have a module written in one style, it might be a good idea to not mix too much with another, as that might actually confuse and cause less readability.

And in this particular case, I do find the static declaration within the function more easily understandable and readable than the pointer magic spread over all the code.

[annakrat]

Replying to mlennert:

I agree with you in principle, although when you have a module written in one style, it might be a good idea to not mix too much with another, as that might actually confuse and cause less readability.

Technically, both approaches (explicit passing of a pointer and using static variable) were used in that function, so it's mixing it already.

And in this particular case, I do find the static declaration within the function more easily understandable and readable than the pointer magic spread over all the code.

I spent quite a bit of time on debugging this code and I disagree. I even think that the practice of using static variables caused this serious bug. I am usually not comfortable with using pointers but in this case, it's no magic, it's pretty standard case which is readable. And just because the diff is just one line, it doesn't mean it's neat. I don't plan to change anything now, but that's my opinion.

[mlennert]

Replying to annakrat:

Replying to mlennert:

I agree with you in principle, although when you have a module written in one style, it might be a good idea to not mix too much with another, as that might actually confuse and cause less readability.

Technically, both approaches (explicit passing of a pointer and using static variable) were used in that function, so it's mixing it already.

And in this particular case, I do find the static declaration within the function more easily understandable and readable than the pointer magic spread over all the code.

I spent quite a bit of time on debugging this code and I disagree. I even think that the practice of using static variables caused this serious bug. I am usually not comfortable with using pointers but in this case, it's no magic, it's pretty standard case which is readable. And just because the diff is just one line, it doesn't mean it's neat. I don't plan to change anything now, but that's my opinion.

I apologize if I was a bit rash in deciding to revert your changes and apply Paul's. I should have waited for your reaction first. In my eyes, Paul's solution was more consistent with the rest of the code in the module, but I have no strong opinion on this, and so, if you think that it would be better to reimplement your solution, I can do that.

[annakrat]

Replying to mlennert:

I apologize if I was a bit rash in deciding to revert your changes and apply Paul's. I should have waited for your reaction first. In my eyes, Paul's solution was more consistent with the rest of the code in the module, but I have no strong opinion on this, and so, if you think that it would be better to reimplement your solution, I can do that.

Probably not worth it, I was mainly venting my frustration from the readability of grass code... But I would still argue that consistency is in case of grass not necessarily our main priority but that's for a different discussion.

[neteler]

Ticket retargeted after milestone closed

[neteler]

Ticket retargeted after 7.0.3 milestone closed

[neteler]

See also #2016

[annakrat]

I think this one still needs to be addressed for 6.4.6, right?

[neteler]

Replying to annakrat:

I think this one still needs to be addressed for 6.4.6, right?

done: backport of r67242 fix for severe bug in distance calculations: r68863

Is this ticket now solved?

[annakrat]

I am afraid no, r.surf.idw still seems to give weird results (see the attachment).

[mlennert]

The v.surf.idw related issues are solved, so we should probably either rename this bug once again to make it clear that it is about r.surf.idw, or close this one and open another one for the r.surf.idw issue.

[annakrat]

Opening a new issue would be best.

[hamish]

Hi folks,

ticket renamed. I left the v.surf.idw in the ticket name so as not to be revisionist. Don't think it needs a separate ticket cluttering the trac system.

AFAIU from the ticket/screenshot r.surf.idw results are still pretty wacky.

regards, Hamish