OSGeo
Interference in gdal_contour at grid lines
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
- x.x99999 +
- 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.
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;
}
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.
No big deal, but we still see odd wave effects in GDAL 3.10. There even is a secondary wave that flips 0002 and 0017.
$ gdal_contour --version
GDAL 3.10.1, released 2025/01/08
$ cat caloocan_xy.1.xyz
-12 12 12
11 12 12
-12 0 0
11 0 0
$ gdal_contour -i 1 -q caloocan_xy.1.xyz /vsistdout/ -f CSV -lco GEOMETRY=AS_WKT
WKT,ID
"LINESTRING (22.5 0.0,11 0,-12 0,-23.5 0.0)",0
"LINESTRING (22.5 1.0,11.0 1.0,-12 1.0,-23.5 1.0)",1
"LINESTRING (22.5 2.0,11.0 2.0,-12 2.0,-23.5 2.0)",2
"LINESTRING (22.5 3.0,11 3,-12 3,-23.5 3.0)",3
"LINESTRING (22.5 4.0,11.0 4.0,-12 4.0,-23.5 4.0)",4
"LINESTRING (22.5 5.0,11.0 5.0,-12 5.0,-23.5 5.0)",5
"LINESTRING (22.5 5.99999900000017,11.0 5.99999900000017,-12 5.99999900000017,-23.5 5.99999900000017)",6
"LINESTRING (22.5 7.0,11.0 7.0,-12 7.0,-23.5 7.0)",7
"LINESTRING (22.5 8.0,11 8,-12 8,-23.5 8.0)",8
"LINESTRING (22.5 9.0,11 9,-12 9,-23.5 9.0)",9
"LINESTRING (22.5 10.0,11 10,-12 10,-23.5 10.0)",10
"LINESTRING (22.5 11.0,11 11,-12 11,-23.5 11.0)",11
"LINESTRING (22.5 11.9999990000002,11.0 11.9999990000002,-12 11.9999990000002,-23.5 11.9999990000002)",12
Please explain what is the wave effect. That the seventh line is at y=5.99999900000017 while you await that it is promptly at 6.0, or what?
All I am saying is way back at comment 0, https://github.com/OSGeo/gdal/issues/10103#issue-2329076437 there were only 0002's. (True, input was different though.) Now there are 0017's swapping back and forth with the 0002's. In an ideal world there should be no long lines of 0000000000's at all. All lines should be about the same length with no weird long lines every few lines.
OK, I just tried all the steps from comment 0, and got the exact same results as back in comment 0.
I.e., whatever got merged above didn't affect this bug. (But the bug was not closed. So maybe that is fine. I see. All that was added was a test. So it is a work in progress.)
Gdal_contour interpolates the locations of lines with fixed height or whatever measure the z coordinate is expressing from the raster data. In this special case it is easy to calculate by other means that the correct y location would be exactly at 6.0 but now the estimation is about one micrometer off. I guess that it is not the 0.000001 m inaccuracy in the interpolation that irritates you, but the number of decimal places in the CSV on some lines.
I consider that the accuracy of the interpolation is the main thing, and it is good enough.