Interference in gdal_contour at grid lines

[jidanni]

Here we make a perfect wedge,

gdal_create -q -outsize 2 2 -f XYZ /vsistdout/|\
perl -anwle '$F[2]=$F[0]; print "@F"' > 0.xyz
cat 0.xyz
0.5 -0.5 0.5
1.5 -0.5 1.5
0.5 -1.5 0.5
1.5 -1.5 1.5

and drape contours over it,

gdal_contour -3d -q -a Elevation 0.xyz 0.csv -i .05 -lco GEOMETRY=AS_WKT
head -n 1 0.csv && grep -P -C 1 .{99} 0.csv
WKT,ID,Elevation
"LINESTRING Z (0.95 -2 0.95,0.95 -1.5 0.95,0.95 -0.5 0.95,0.95 0.0 0.95)",8,0.95
"LINESTRING Z (0.999999000002 -2 1,0.999999000002 -1.5 1,0.999999000002 -0.5 1,0.999999000002 0.0 1)",9,1
"LINESTRING Z (1.05 -2 1.05,1.05 -1.5 1.05,1.05 -0.5 1.05,1.05 0.0 1.05)",10,1.05
--
"LINESTRING Z (1.45 -2 1.45,1.45 -1.5 1.45,1.45 -0.5 1.45,1.45 0.0 1.45)",18,1.45
"LINESTRING Z (1.499999000002 -2 1.5,1.499999000002 -1.5 1.5,1.499999000002 -0.5 1.5,1.499999000002 0.0 1.5)",19,1.5

Let's take a closer look at the non-exact coordinates. They are composed of

1. x.x99999 +
2. 0.000000000002

and only occur on exact grid lines, the least place we would expect interference to occur. And not one, but two types of interference in fact.

Yes one might say that it is a tolerable amount of interference. But it doesn't happen anywhere else than at exact grid lines so is definitely not simple division residue else it would occur on interpolated contours too.

Perhaps the program has a rounding stage that is used when interpolating, but skipped when on direct grid lines.

Trying again with gdal_contour -3d -i .05 -off 88 has differing Z and Elevation for the first half of the output. Reading the man page about -a, they should be a direct copy.

All it takes is a tiny offset, gdal_contour -3d -i .05 -off .0001, and all the problems are washed away. Same with just gdal_contour -i .03451345: never hits an exact contour, so never encounters the problem.

P.S., when given e.g.,

gdal_contour -i -0.03451345

the program should be enhanced to be sure to bomb out with

ERROR: Intervals must be greater than 0.

Currently the program gets into a loop and must be killed.

[rouault]

Perhaps the program has a rounding stage that is used when interpolating

That seems to be related to the "fudge" mechanism in the marching square algorithm. Applying the below experimental patch removes the unexpected numerical difference on the sample test case of this ticket, but causes failures/crashes in unit tests. So this fudge mechanism seems to be an integral part of the algorithm. Fixing the issue might require lots of work

--- a/alg/marching_squares/square.h
+++ b/alg/marching_squares/square.h
@@ -484,8 +484,8 @@ struct Square
                 y1 = ym;
             }
         }
-        const double fy1 = fudge(level, y1);
-        const double ratio = (level - fy1) / (fudge(level, y2) - fy1);
+        const double fy1 = y1; // fudge(level, y1);
+        const double ratio = (level - fy1) / (y2 /* fudge(level, y2) */ - fy1);
         return x1 * (1. - ratio) + x2 * ratio;
     }

[jidanni]

Also one can tell right away that the program is doing the steps in the wrong order, when I added -off 88 and managed to get different output.

It should have resolved to the same as offset zero, before proceeding with any other calculations.

In other words there is some bias being injected, depending on where we start computing.

For instance if we want 100m contour lines they should be the same no matter if we started at the top of Mount Everest 8800m or at sea level 0m. (8800 & 88 above just coincidence. That 88 was just a random big number I picked.)

Maybe fixing that first would make fixing the rest a lot easier.