Support for getBounds() with conic and thematic projections

[SnailBones]

Fixes https://github.com/mapbox/gl-js-team/issues/494

In some projections, the straight edges of the screen become curved when transformed into Mercator space. This issue is the cause of bugs with bounds-related APIs with projections and is most noticeable in conic projections like Albers. While getBounds() previously only checked the corner points, this PR checks additional points to ensure that lng/lat extremums are correctly included in the returned bounds.

This PR introduces a new function _getBoundsProjection() which checks 19 points along each edge in addition to the corners (20 * 4 = 80 points total), the extremums of all of them forming the bounds. This is the same approach used in the debug page introduced in #12199, and is distinct from the recursive error-checking approach introduced for globe in #12286.

EDIT: The new behavior is implemented in _getBoundsNonRectangular(), which replaces _getGlobeBounds(). Logic remains the same except we now check points in in LngLat space instead of Mercator coordinates. The new function is used for globe, conic and thematic projections (i.e. all besides Mercator and Equirectangular) and supports terrain.

When in non-globe projections, we check a minimum of 16 points per side. We check more points outside of globe since projecting a screen point (called in pointLocation) in non-globe projections is a faster operation (in globe or with terrain it requires raytracing), and missing areas are more visible in projections (the curved edges of globe hide them) .

Albers:

Before:

After:

Equal Earth:

Before:

After:

Launch Checklist

• briefly describe the changes in this PR
• include before/after visuals or gifs if this PR includes visual changes
• write tests for all new functionality
• document any changes to public APIs
• manually test the debug page
• apply changelog label ('bug', 'feature', 'docs', etc) or use the label 'skip changelog'
• add an entry inside this element for inclusion in the mapbox-gl-js changelog: <changelog>Fix map.getBounds returning incorrect bounds with adaptive projections.</changelog>

[mourner]

This is the same approach used in the debug page introduced in #12199, and is distinct from the recursive error-checking approach introduced for globe in #12199.

Just wondering — do we need to use a distinct approach here? Looking at the globe bounds approach, there's nothing specific to globe aside from pole handling — in theory, could we reuse the same code and thus make the PR smaller?

[SnailBones]

This is the same approach used in the debug page introduced in #12199, and is distinct from the recursive error-checking approach introduced for globe in #12199.

Just wondering — do we need to use a distinct approach here? Looking at the globe bounds approach, there's nothing specific to globe aside from pole handling — in theory, could we reuse the same code and thus make the PR smaller?

I tried that first and it returns incorrect bounds, like so:

I think the reason is that in non-globe and non-Mercator projections, pointCoordinate() returns a Mercator coordinate that's not yet projected (It's in the grid of the map plane but doesn't correspond to geographic space).

To work around this, the new function interpolates builds the bounds directly in projected lngLat values instead of Mercator coordinates. (i.e. replacing pointCoordinate() with pointLocation())

But I think that the lngLat approach should work in both Globe and Mercator, so I'll see if I can combine the two functions!

[SnailBones]

Just pushed a change combining the two functions.

One issue is that now your concern @mourner becomes shared with projections:

We should also be careful of the case where midpoint is in the middle but the line gets curvy on either side of the midpoint

This issue is likely not noticeable in globe, since regions near the edge of the earth are very compressed.
In conic projections, it's less subtle:

But given it's unlikely for developers to be using conic projections at low zoom (especially with high pitch, which brings out the worst behavior), maybe these buggy edge cases are an acceptable cost given the improvement to code simplicity.

[mourner]

@SnailBones hmm, this doesn't look quite right... The adaptive sampling should only be valid in screen space — this guarantees the midpoint we take between two points lies on a straight line between them, so that it also belongs to the screen bounds line, but that no longer holds with lng/lat space which is non-linear. The latter also makes it less likely to converge on a settled result for the screen, which likely leads to artifacts above. Maybe let's revert to the previous approach of simple uniform sampling, even if it's more code? We could revisit this later — I'd be interest in digging into this.

[SnailBones]

The adaptive sampling should only be valid in screen space — this guarantees the midpoint we take between two points lies on a straight line between them, so that it also belongs to the screen bounds line, but that no longer holds with lng/lat space which is non-linear.

Not sure if I follow this. To clarify, my modified adaptive sampling function still follows the same conversion: Screenspace -> Map plane/Mercator space -> Lng/Lat space. The difference is which step the max/min comparisons take place to construct the bounds: previously Mercator space, now Lng/Lat space.

Besides performance, the relevant change on Globe is the different space used for comparing error values to maxErr. One way to think of this is that previously maximum error was represented as distance on a Mercator map, now it's on an Equirectangular projection. I don't think there's any inaccuracy introduced here: a straight line in screen space in Globe is a curve in both in Mercator and Equirectangular. In fact, representing error in Equirectangular Lng/Lat space should be consistent across latitude, though it shares longitude exaggeration toward the poles with Mercator.

While I don't think the validity of the algorithm's behavior is affected, that's not to say that this is necessarily an improvement or that it justifies the extra conversions to Lng/Lat which call Math.atan!

The latter also makes it less likely to converge on a settled result for the screen, which likely leads to artifacts above.

It's true that the result can be very jumpy, especially at high pitch, where a small amount of screen space can cover the majority of map space due to perspective. In cases where the screen midpoint does not expand bounds but an extremum is in the upper 1/2 of the screen (as seen in my second screenshot above), adaptive sampling will stop traversing and miss the extremum.

Similar gaps can of course be missed with uniform sampling, but with something like 20 samples on each side any missed gaps will be much smaller. We could keep the combined function and solve this with a minimum recursion depth by projection type: i.e. in Albers, always check a minimum of 8 or 16 points per side.

[SnailBones]

I added a minimum of 15 samples per side on projections (16 including corners).

Bounds now look correct in the above cases, and manual testing in all projections look good. There's still edge cases where the bounds look strange, but the missing area takes up a small amount of screenspace (maximum 1/16 of screen height.)

There are some missing edges on Globe, but these are consistent between this branch and main and occur close enough to the horizon that this is probably not noticeable in most cases.

main:

this branch:

main:

this branch:

There's also a strange issue where during the globe transition at high pitch, the bounds covers an area behind the camera. Likely the top two corners are being placed on the map incorrectly. This is the same between main and this PR.

main:

this branch:

[mourner]

Not sure if I follow this. To clarify, my modified adaptive sampling function still follows the same conversion: Screenspace -> Map plane/Mercator space -> Lng/Lat space. The difference is which step the max/min comparisons take place to construct the bounds: previously Mercator space, now Lng/Lat space.

@SnailBones Apologies, I misread the code — sorry for derailing you here! The PR is a great improvement — I'm sure we'll be able to track down the mentioned edge cases later when they matter.

[karimnaaji]

LGTM as this is clearly an improvement for alternate projections. Regarding the cases mentioned above that could still be improved with globe or high pitch, it'd be worth capturing the issue with potential ideas of improvements on how to make it more robust:

[SnailBones]

Thanks @karimnaaji !

I'll go ahead and merge since this covers most cases, but here's notes for future improvements:

As seen in your image above, the biggest issue still remaining is highly pitched maps revealing the sky. In Mercator and globe, pointCoordinate() snaps points in the sky to the closest point on the map edge, but this behavior is missing a full implementation with other projections, resulting in points snapping to strange locations. Fixing this should lead to much more accurate bounds at high pitch.

[gagecarto]

Are there solutions to this issue? It would be great to initialize a map in Albers and specify bounds or to fitBounds after the map is initialized.