2
votes

I have a set of n geographical coordinates on the surface of Earth, for which I want to calculate a bounding box (find the easternmost, westernmost, northernmost and southernmost location) without falling back to user input (there's no UI to the program). The naive approach is "take max and min of latitude, max and min of longitude, done" - but this returns obviously suboptimal results when the set straddles the 180th meridian (see e.g. Fiji for a similar situation: https://www.openstreetmap.org/relation/571747#map=2/-16.6/0.0 should not zoom out to the whole planet, as the two halves are actually adjacent).

This is acknowledged in multiple solutions, e.g. in Algorithm for determining minimum bounding rectangle for collection of latitude/longitude coordinates , but not solved.

What didn't work:

  • check how others do this (Leaflet has the same problem, see above)
  • enumerate pitfalls (e.g. "if one of the cooordinates gets close to the 180th meridian") - won't work for points whose bounding box is across, but not close to the IDL (e.g. Japan and Hawaii)
  • go through the coordinates and mark if 180 is crossed (depends on the ordering)
2
I don't really need a solution in any specific programming language: as this is a spatial algorithm, it should be possible to calculate without much difficulty regardless of language specifics. - Piskvor left the building

2 Answers

3
votes

I think there is a simply path here. I will only consider longitude, because latitude do not cause any issue. This solution provides an O(NlogN) (due to sorting) solution, that means it is slower than only check for minimum and maximum, but running it on the majority of machines and the majority of programming languages with less then 10^5 points should take not more than a couple of seconds

It is obvious that there are several solutions for boundaries, from a mathematical point of view, since the longitude value can be taken modulo 360.

Several solutions

As we can see in the above example image. The best choice is the green box because it has the smallest size and contains all the points.

Find the smallest box that contains all the points is the same problem as find the biggest box that contains no points

enter image description here

So in this simply image we need to find the two points at the extremes of D

The algorithm to find D (and points associated) require to sort the points for longitude, so between two consecutive sorted point there will certainly be no other points; so we can just check the distances between them (and last-first, too). Here some pseudo-code

Let C = your set of points
Let N = length( C )
Let S = sort( C )
Let maximumDistance = 0
Let easternLongitudeForBoundingBox = undefined
Let westernLongitudeForBoundingBox = undefined
for i = 0 to N-1 :
   Let j = (i + 1) modulo N    # The index of the point "after" point(i)
   Let D = (S[j].longitude - S[i].longitude ) modulo 360
   if (D > maximumDistance) :
      maximumDistance = D
      easternLongitudeForBoundingBox = S[i].longitude
      westernLongitudeForBoundingBox = S[j].longitude

This is totally untested, but if I didn't make mistake, it should work.

Warning: when I use "modulo" it is the real mathematical operator. In some computer languages, (a - b) modulo 360 should be written (360 + a - b) % 360 to give correct results because they think that -1 % 360 == -1 which will give incorrect results here.)

-1
votes

The Leaflet documentation at https://leafletjs.com/reference-1.5.0.html#latlngbounds states:

Caution: if the area crosses the antimeridian (often confused with the International Date Line), you must specify corners outside the [-180, 180] degrees longitude range.

So, for example, if you wanted an area that spanned from 40° to 60° latitude and across the antimeridian from -160° to 160° longitude, it seems you'd need to specify the area's bounds with corners like [40,200], [60,160] or [40,-160], [60,-200].