Hengefinder: Finding when the sun aligns with your street

Next week in Manhattan, the sunset will align perfectly with the east-west streets of the city grid. It’s beautiful, and people know it. Twice a year, crowds gather to see the brief moment when the sun sits perfectly on the horizon, framed by skyscrapers on either side. It’s called Manhattanhenge, after Stonehenge.

I wanted to know how astronomers figure out when Manhattanhenge happens. And if I could figure that out, why limit it to Manhattan?

This was one of my first projects at the Recurse Center: Hengefinder, a tool that lets you find a henge pretty much anywhere the sun sets.

• Source code
• Hengefinder website
• Hengefinder mobile app (A follow up to the website, created by fellow Recurser John Pribyl)

The basic steps (and some terminology I would soon learn):

• find the angle of a road (its bearing, relative to true north)
• find the angle of the sun at sunset each day (its azimuth)
• find the dates when those angles match.

It was appealing to me that this was mostly made up of many sub problems I could either think through thoroughly, or choose to mostly skip (by handing off to a library, a brute-force solution, an approximation, etc.). It’s like a bunch of closed boxes. I left plenty of those boxes closed — I didn’t build my own astronomical model, for instance. Other boxes I opened. And then I found many things I thought would be straightforward that turned out to be somewhat less trivial in practice as my assumptions broke; Like, in reality, you can’t treat roads like flat lines, latitude and longitude don’t behave like a Cartesian grid, and the word “sunset” is more ambiguous than I initially thought. I’ve enjoyed closing the gap between my assumptions and reality.

So, rather than just talking about the project as a whole, I’m going to walk through a couple challenges I ran into in building this, and how I went about solving them. One challenge for each of those supposedly “simple” steps above.

Challenge #1: Finding the road bearing (and rediscovering the Earth is not flat)

The first challenge was calculating the bearing of a street (its angle relative to true north). If you take the latitude and longitude of one address, and the latitude and longitude of another address down the street, you end up with the coordinates of two points on Earth. We can then use some trig to get the angle. My initial (incorrect) guess for getting the road bearing was to take the difference in latitude and the difference in longitude, then just get the angle with atan2 (the inverse tangent).

The problem is, that would only work if the Earth were flat.

Latitude and longitude lines aren’t arranged the same way on the Earth. Latitude lines (which run east-west) are evenly spaced: one degree of latitude is basically the same distance everywhere on Earth. Longitude lines (which run north-south) aren’t. They converge as you move toward the poles. For example, one degree of longitude spans roughly 52 miles (84 km) in New York City, but way less at higher latitudes (about 33 miles / 53 km in Anchorage).

How do you manage that?

TL;DR, if you just want the practical takeaway

Scale longitude by cos(latitude), which puts longitude and latitude in the same “units.” Only then do atan2.

”Ok, but why?”

(And if you don’t really care about the math, feel free to skip to the “corrected approach” section below.)

When you do tan = opposite / adjacent, you’re assuming the dimensions of the triangle are in consistent units. They aren’t here. Like I said, degrees of longitude shrink as you move toward the poles, while latitude stays consistent. That means the horizontal and vertical sides of the triangle are in different units.

To fix this, we need a little spherical geometry.

Imagine the Earth as a unit sphere (radius = 1), and take a vertical cross-section through it. A point at latitude φ forms a right triangle:

• the hypotenuse is the Earth’s radius (1)
• the vertical leg is sin(φ)
• the horizontal leg is cos(φ)

That horizontal leg is the distance from the Earth’s axis to the point of our “observer.” It’s also the radius of the circle you trace when you move east-west at latitude φ.

This means that at the equator (φ = 0°), the radius is cos(0) = 1, and at the poles (φ = 90°), the radius is effectively cos(90°) = 0.

So, east-west distances shrink as you move away from the equator. You’re walking around smaller and smaller circles. To compare longitude meaningfully with latitude, you have to scale it by the radius of that circle, which is cos(latitude).

Once longitude is scaled this way, the two axes live in comparable units, and you can use the inverse tangent.

Corrected approach

# assuming lat/lon are in degrees
delta_y = lat_2 - lat_1
# scale longitude so it's in the same units as latitude
mean_lat = math.radians((lat_1 + lat_2) / 2)
delta_x = (lon_2 - lon_1) * math.cos(mean_lat)
# now we can safely use atan2
bearing_rad = math.atan2(delta_x, delta_y)
bearing_deg = math.degrees(bearing_rad)

(Note that this is still an approximation, but for short street segments, it’s accurate enough.)

Challenge #2: Finding the sun’s azimuth (or, what do we mean by sunset?)

For a perfect “henge” moment, the sun needs to be sitting directly on the horizon.

I used the Python library Astral to compute the timing, altitude, and azimuth (angle) of solar events at a given location. But Astral’s definition of a sunset is slightly different from what I want here. Astral uses the standard astronomical definition of sunset, i.e. when the sun has fully dipped below the horizon. For my purpose, that’s too late. I want to see when the sun is sitting right on top of the horizon. Astral’s version is still useful to me, however, since time-wise, it gets me in the ballpark on a given day.

For a moment to be a “henge,” I want the sun’s disk to be fully visible and just touching the horizon. I have a target “altitude” I want the sun to be at, above the horizon, based on the apparent size of the sun.

Why this is a boundary search, not a value search

So I’m trying to find the last moment when the sun’s altitude is still above my target threshold altitude.

You could just do a linear search minute-by-minute in some time window before Astral’s sunset. It’d probably be fine, but each evaluation is an Astral API call. Plus, that’s less fun. This is one of the “boxes” I wanted to open. Because the sun’s altitude changes monotonically as sunset approaches, this indicates a binary search.

The textbook approach to binary search is this:

def basic_binary_search(arr, target):
    left = 0
    right = len(arr) - 1
    while left <= right:
        mid = (left + right) // 2
        if arr[mid] == target:
            return mid
        elif arr[mid] < target:
            left = mid + 1
        else:
            right = mid - 1
    # Not found
    return -1

But that formulation makes several assumptions that don’t really apply here. For one, we’re not searching for an exact target match. Because with minute-level resolution, we may never get the exact moment. Instead, in looking for a henge, I’m really evaluating a boolean each minute: “Is the sun still above the target altitude at this minute?” That gives a pattern like:

True, True, True, True, False, False, False

Rather than looking for a specific value, now we’re just looking for the last True before it flips.

A boundary (last-true) binary search

So, I used a different approach to the binary search. Here’s the core loop, searching within a window of time before Astral’s sunset:

while left < right:
    mid = (left + right + 1) // 2 # upper-biased midpoint for a "last true" search
    if altitude > target: # still valid, keep it
        left = mid # this minute is valid, could be the last one
    else: # sun at or below horizon
        right = mid - 1 # this minute is invalid, look earlier

This turns it into a “last true” binary search, which lets me find the final minute when the sun is still above the target altitude.