polygonToCells fails for some simple large polygons

[mutability]

(Originally: uber/h3-py#202)

This polygon (nb: in degrees, not radians) causes polygonToCells to return E_FAILED (at resolution 5, probably true for other resolutions):

{
  {-74, -120},
  {-74, 0},
  {-74, 120},
  {-84, 120},
  {-84, 0},
  {-84, -120}
}

I poked at this a bit and the problem AIUI is:

• polygonToCells ends up generating more than the "maximum" number of cells; because
• maxPolygonToCellsSize underestimates the number of cells needed for this polygon; because
• bboxHexEstimate underestimates the number of cells needed

For this poly bboxHexEstimate returns 16037 cells at resolution 5. By messing with POLYGON_TO_CELLS_BUFFER I got polygonToCells to run to completion, and it generated 20337 cells.

My intuitive guess at what's going wrong is that bboxHexEstimate naively measures the distance between the two corners of the bounding box using distanceKm. But this always returns the shorter great-circle arc, which for large bounding boxes is wrong (for this polygon, the shorter arc is not the diagonal of the bounding box)

[mutability]

Test case showing the problem:

#include <h3api.h>
#include <inttypes.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>

static LatLng goodVerts[] = { {-74, -120},
                             {-74, 0},
                             {-74, 120},
                             {-82, 120},
                             {-82, 0},
                             {-82, -120} };
static GeoLoop goodLoop = { .numVerts = 6, .verts = goodVerts };
static GeoPolygon goodPoly;

static LatLng badVerts[] = { {-74, -120},
                             {-74, 0},
                             {-74, 120},
                             {-84, 120},
                             {-84, 0},
                             {-84, -120} };
static GeoLoop badLoop = { .numVerts = 6, .verts = badVerts };
static GeoPolygon badPoly;

static void initPoly(GeoPolygon *poly, GeoLoop *loop)
{
    for (int i = 0; i < loop->numVerts; ++i) {
        loop->verts[i].lat = loop->verts[i].lat * M_PI / 180.0;
        loop->verts[i].lng = loop->verts[i].lng * M_PI / 180.0;
    }
    poly->numHoles = 0;
    poly->geoloop = *loop;
}

static void fill(const char *what, GeoPolygon *poly)
{
    H3Error error;

    int64_t maxHexagons;
    error = maxPolygonToCellsSize(poly, 5, &maxHexagons);
    fprintf(stderr, "%s: maxPolygonToCellsSize() returns %d\n", what, error);
    if (error != E_SUCCESS)
        return;
    fprintf(stderr, "%s: %" PRIi64 " max hexagons\n", what, maxHexagons);

    H3Index *hexagons = calloc(maxHexagons, sizeof(H3Index));
    assert(hexagons);

    error = polygonToCells(poly, 5, hexagons);
    fprintf(stderr, "%s: polygonToCells() returns %d\n", what, error);

    if (error == E_SUCCESS) {
        int64_t count = 0;
        for (int64_t i = 0; i < maxHexagons; ++i)
            if (hexagons[i] != H3_NULL)
                ++count;

        fprintf(stderr, "%s: %" PRIi64 " non-null hexagons filled\n", what, count);
    }
    free(hexagons);
}

int main(int argc, char **argv)
{
    initPoly(&goodPoly, &goodLoop);
    initPoly(&badPoly, &badLoop);
    
    fill("good", &goodPoly);
    fill("bad", &badPoly);

    return 0;
}

results:

good: maxPolygonToCellsSize() returns 0
good: 18511 max hexagons
good: polygonToCells() returns 0
good: 17728 non-null hexagons filled
bad: maxPolygonToCellsSize() returns 0
bad: 16049 max hexagons
bad: polygonToCells() returns 1

[nrabinowitz]

Thanks for the bug report and test case! I think your diagnosis here is sound, and we can take a look at potential fixes.