Fix missing tiles in corners for fast gridding by using off-boundary pixels

[joshuacortez]

This PR fixes the missing corner tiles when using FastSquareGridGenerator and FastBingTileGridGenerator

Results

1. In the tutorial, there is no tile mismatch between the fast and regular versions of both square and Bing tile gridding anymore. 🎉
2. After testing for a sample Indonesia ADM2 AOI, the percent match between FastBingTileGridGenerator and BingTileGridGenerator jumps from 99% to 99.9%.

Describing the issue

In the gridding tutorial there is a discrepancy between the regular and the fast version of both square and Bing tile gridding.

The image below shows the mismatch: there is 1 tile available in BingTileGridGenerator but not in FastBingTileGridGenerator (see green tile)

While a 1 tile mismatch might seem innocuous and coincidental, this points to a deeper issue in general. Usage using other AOIs by @tm-jc-nacpil show that FastBingTileGridGenerator can miss even tiles with a large AOI intersection (see bottom left arrow)

Diagnosis

Upon further inspection, I saw this issue is due to 2 things:

1. There is rounding when translating geographic coordinates to pixel coordinates.
2. The way diagonal steps are handled in voxel traversal.

The image below is the same AOI and 1 missing tile (shown in green) in terms of pixel coordinates instead of geographic coordinates. Note that the vertical orientation is flipped in pixel coordinates.

Notice how the 1 missing tile in green has a corner that is tangent to the blue line segment. This points to the underlying issue of off-boundary pixels. These are pixels whose corners are tangential to an AOI polygon edge. Most importantly, these pixels are candidates for adding back to the main tile set for error correction.

You can read more about off-boundary pixels in the 15_polygon_fill.ipynb notebook in this PR. Examples of off-boundary pixels in the image below:

Resolution

Here are the general steps for the fix:

1. Identify the off-boundary pixels. Several functions in 15_polygon_fill are updated to return these.
2. We filter the off-boundary pixels that intersect linestring of the AOI boundary. The intersections should be done in terms of geographic coordinates, not pixel coordinates. This is implemented in the generate_grid method in FastSquareGridGenerator and FastBingTileGridGenerator
   a. Despite this step relying on geometry intersection, the performance hit should be relatively small since we're doing a linestring to polygon intersection (not a polygon to polygon intersection), and the off-boundary pixels should be a small fraction of the total tiles.
3. We then add the filtered off-boundary pixels to the main tile set

Note that even after this fix, there still isn't a 100% guarantee of a perfect match between the regular and fast versions of gridding, which I suspect can still happen due to floating point precision errors. But the percent match should have nonetheless improved.

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[tm-jc-nacpil]

lgtm! keen to merge it so we can test it out in other repos