Cardinal directions

If you go to Jerusalem, you’ll find a sign saying Jesus was crucified here, and another saying this is the cave where he awoke after three days. And if you go to Troy, you’ll read about Homer sitting on this very rock while (researching?) writing the Iliad.

And if you go to any bit of land that juts impressively enough out to sea, you might find something claiming that it’s the southernmost/southeasternmost/etc part of whichever land mass you happen to be on.

With the four points of the compass, this is generally easy enough (except for Greenland, where bits of gravel make things complicated), but what exactly does southeasternmost mean?

Prompted by an overly confident plaque, I dusted off my trig knowledge and decided to find out. My starting assumption is that the something-most point of land is the one where, if you move a perpendicular line something-ward, that point will be the last land that the line touches.

To see which point wins is then a matter of finding the point that maximises the following equation, which blends the latitude and longitude:

d = lon * sin(a) + lat * cos(a)

Where lon and lat are self-explanatory, and a is the bearing of interest (north is 0°, east is 90° and so on). I’d like to explain this equation but I’m not sure why it works, I just stumbled my way into it and it seems that it does work.

So for north, south, east and west, this simply collapses to maximising/minimising the latitude or longitude. For the other points of the compass, it’s a bit more interesting. SE (for example) is a half way between S and E. ESE is a bit east of that. EbS is just a hair south of east.

I was feeling generous, so I calculated the something-most points not only for the cardinal (N, S etc) and inter-cardinal (SE and co) but also the half-winds (ESE et al) and even the quarter winds (EbS inter alia). I wasn’t feeling generous enough to include the fractional divisions past that, so each landmass gets 32 extremities!

Calculating all these points looks something like this:

for b in bearings:
    for p in points:
        d = p.x*sin(b)+p.y*cos(b)
        if d > biggest:
           save p

My first geometry file (South African coastline) had about 50,000 points in it, and this took about 20 minutes to run. I thought of doing clever things with hash tables but then I remembered: numba. By converting the code from Pythonic -> pure numpy -> numba jitted, the time went from 20 minutes -> 20 seconds -> 0.2 seconds, for a speed up of ~5,000 times, without even having to think too hard! Then I ran the same code for 53,623,594 points (all countries in the world, plus landmasses separately), and it still only takes a second or two.

Go play!

Naturally, I threw the data together with some Vue and Mapbox so everyone can enjoy it, and go around setting records straight.

It’s live here at rdrn.me/cardinality.

The code to generate the data (and for the site) is on GitHub.

There are probably a lot of weird cases, so either blame me and my algorithm, or the input data I used. But probably the latter.