polyfill not filling the complete polygon

[nmandery]

Hello,

I just stumbled on an issue where the polyfill function is not able to fully fill polygons. The problem affects long, diagonal polygons, like this:

The blue part is the polygon, the green shape are the generated hexagons. Here is a closeup of the lower end of the polygon:

The polyfill was done at resolution 13. Interestingly the result at resolution 12 looks much better:

The GeoJSON for my test polygon is

{"type":"Polygon","coordinates":[[[5.76910652153309,51.1122595579185],[5.74843667987113,51.0844285870957],[5.74868565718131,51.0841156106137],[5.76939661923247,51.1120877871771],[5.76910652153309,51.1122595579185]]]}

To further investigate this issue I went ahead and removed the line of code which does the removal of heaxgons located outside the polygon in this file:

h3/src/h3lib/lib/algos.c

Line 710 in 96a0e86

out[i] = H3_INVALID_INDEX;

I took this approach because I wondered about the clean cut at the end of the filled part and wanted to see in which radius H3 is generating it indexes.

Running polyfill on my test polygon again gave the following results:

Resolution 12:

Resolution 13:

So the filled minK hex radius seems to be a bit too small for this case. I gave it a try and changed the 1.5 multiplicator in

h3/src/h3lib/lib/bbox.c

Line 120 in 96a0e86

return (int)ceil(bboxRadiusKm / (1.5 * centerHexRadiusKm));

to the - wildly guessed value - 1.35:

return (int)ceil(bboxRadiusKm / (1.35 * centerHexRadiusKm));

This value is certainly not mathematically correct and needs to be validated. It also increases the number of indexes being generated, but gave me a polygon completely covered in indexes:

So I assume the issue is in the bboxHexRadius function in combination with filling shapes with a large number of indexes. This function seems to have some deficits dealing with long, thin, diagonal shapes.
I did not notice any issues when using polyfill with rounder geometries.

[nrabinowitz]

First of all, thanks! This is a great bug report - really thorough. I think the issue here is that the 1.5 constant comes from a regular polyhedron model, and doesn't work when the projection distorts the shape of the hexagon - I think this is less a function of the shape you're trying to polyfill, and more a function of whether it's aligned with the long or short axis of the hexagons in its particular area of the grid. My guess is that the same shape, rotated 90 degrees, would be polyfilled correctly with the current settings (and in fact that's what happens at res 12 - suggesting that res 11 would also potentially fail to fill the shape).

We've been considering an alternative approach for polyfill that would solve this problem and be much more efficient for this kind of shape, using h3ToIJ (#102). I'd rather pursue that approach than try to fiddle with the scaling factor here, though I'd like thoughts from @isaacbrodsky and @dfellis on whether that's a reasonable stop-gap. One option would be to calculate the scale factor from a hexagon at res - 1, which is the approximate shape of the k-ring.

[dfellis]

Definitely an amazing bug report.

I disagree on leaving the issue alone until h3ToIJ is ready for use, though. I think we should alter the scaling factor to guarantee correctness and deal with the polyfill slowdown for the short term and then get the performance back up.

We don't have to do the math to determine the maximum extent of hex distortion, if that's what you're worried about wasting time on; just set the constant to something we know could never be reached even with distortion, such as 1.1, and move on.

[nrabinowitz]

just set the constant to something we know could never be reached even with distortion, such as 1.1, and move on.

I'm fine with this, though given the perf implications it seems worthwhile calculating the worst-case scenario offline and using that as a constant.

[nmandery]

Thank you both!

You are correct @nrabinowitz , resolution 11 also has this problem:

I had a brief look at #102 and will keep an eye on this functionality in the future. For now I will wait and use the released version of H3.

In the meantime I am thinking of additionally tackling this issue outside of H3 by splitting the polygon into smaller parts and running polyfill on each of them separately. I hope to compensate the increased radius and keep performance up by having to generate and check less indexes.

[nmandery]

I just noticed something embarrassing - seems I had a quite stupid bug in my code I simply missed.

I will re-check this issue tomorrow, but there is a good chance this issue is invalid and H3 works as expected. So please do not spend any more time on this right now.

[nmandery]

Please close this issue as invalid. The problem was caused on my side where I had a systematic mixup of coordinate positions (x,y to lon,lat and back). I am working on an Postgresql binding for H3 and the bug happened during the translation from native postgresql geometry types to the H3 types. Because this error was a systematic mistake, it took a while for me to notice it. It can be frustrating when such simple errors can be the cause for bugs.

Thank you both for your support. I am sorry for the effort this caused on your side.

[isaacbrodsky]

Glad to hear it's working for you. I'll close this issue since it seems the current solutions is working well enough.

[nrabinowitz]

I think it might still be worth looking into this, though - the original bug report seemed legit, and even if the OP's current issue is solved, there may still be areas of the world/H3 resolutions where this is an issue. Seeing if we can get a repro case would allow us to say whether we fixed it with an h3ToIJ implementation.

[nrabinowitz]

This issue is legitimate, even if it doesn't affect the OP. Reversing the lat/long in the test case is still a valid polygon, and does hit this issue.

Repro case:

static GeoCoord testVerts[] = {{0.10068990369902957, 0.8920772174196191},
                               {0.10032914690616246, 0.8915914753447348},
                               {0.10033349237998787, 0.8915860128746426},
                               {0.10069496685903621, 0.8920742194546231}};
static Geofence testGeofence = {.numVerts = 4, .verts = testVerts};
static GeoPolygon testPolygon;

SUITE(polyfill) {

    testPolygon.geofence = testGeofence;
    testPolygon.numHoles = 0;

    TEST(weirdShape) {
        int res = 13;
        int numHexagons = H3_EXPORT(maxPolyfillSize)(&testPolygon, res);
        H3Index* hexagons = calloc(numHexagons, sizeof(H3Index));

        H3_EXPORT(polyfill)(&testPolygon, res, hexagons);
        int actualNumHexagons = countActualHexagons(hexagons, numHexagons);

        printf("Hex count: %d\n", actualNumHexagons);

        free(hexagons);
    }
}

Output with the current 1.5 coeffecient in bboxHexRadius:

TEST polyfill
Hex count: 3945

Output with the coefficient set to 1.0:

TEST polyfill
Hex count: 4353

[greatathrun]

Sorry for my Englist!

I try to creat a KML file that contain the global grid in res 2. I find indexes within the given geofence by Polyfill and get the boundary by GeoToBoundary. I found an issue where the grid which the Function's return out of the geofence.You can see the geofence vertices and the grids in the image:

[nrabinowitz]

@greatathrun There's currently a known issue where you can't polyfill a polygon containing a pole. For this particular purpose, though, you ought to be able to use the h3ToGeoBoundaryHier utility, passing in each of the base cells.