rotated_map

Put the north pole wherever you like, and draw the map that falls out.

The "Bay Area at the top of the world" poster isn't a new projection — it's an ordinary Mercator applied to rotated coordinates. Spin the globe until your chosen point sits where the north pole normally does, then project as usual. Everything Mercator normally does to Greenland it now does to your home town.

python rotated_map.py --lat 37.7749 --lon -122.4194 -o map.png

That poster is a reproduction of a map going round the internet; the exact parameters were recovered by fitting (see Matching an existing map).

python -m venv .venv && .venv/bin/pip install -r requirements.txt

Natural Earth vector data is downloaded on first use and cached in ./data (gitignored — it comes back on its own).

verify.py checks the geometry against the source data: it samples a grid in rotated coordinates, un-rotates each sample, and asks the original polygons whether that point should be land. Any disagreement is a bug in the rotation and clipping, not a matter of taste.

python verify.py --res 10m

--bearing is the one worth playing with. The chosen point becomes the entire top edge of a Mercator, so the map has to be cut open somewhere; the cut runs along bearing + 180. For San Francisco, --bearing 90 puts the tear out in the Pacific and keeps the continents whole.

--top-km is the other one. Mercator sends the pole to infinity, so the top of the page has to be cut at some radius around your point — and that cut throws away everything inside it. Cut at a comfortable-looking 89.3° of rotated latitude and you have quietly deleted everything within 78 km, which is exactly the part you picked the point for. The default of 3 km puts San Francisco Bay across the top of the page; drop to 1 km and the top edge lands inside the city. Going below a kilometre or so mostly magnifies the limits of the source data rather than revealing anything new.

A world map is usually landscape, inside an ellipse — Mollweide or Hammer. You can rotate one of those too, but not with --place pole: on an ellipse the pole is the single point at the top, so sending your city there shrinks it to nothing. Use --place centre and it lands in the middle of the sheet instead, with true north pointing up, and the whole globe re-centred around it:

python rotated_map.py --lat 37.71 --lon -122.64 --place centre \
    --projection mollweide --layers land,countries,lakes \
    --graticule original -o world.png

Nothing is magnified here — Mollweide is equal-area, that is the point of it — but the Pacific stops being a thing at the edges. The starburst above North America is the real north pole, where the meridians still meet.

--place centre works with every projection: it gives an oblique aspect of whatever you ask for. --place pole is the poster mode.

The azimuthal family puts your point in the middle of a disc even under --place pole. --projection azimuthal is an azimuthal equidistant: every straight line from the centre is a true bearing and distance from home.

python rotated_map.py --lat 39.9042 --lon 116.4074 --projection azimuthal \
    --layers land,countries --graticule original -o beijing.png

Anything in a GeoJSON of lon/lat coordinates can ride along. The poster above uses US county boundaries:

python rotated_map.py --lat 37.71 --lon -122.64 --bearing 38.4 \
  --top-km 11 --bottom-lat -76.9 --res 10m --layers land,countries,lakes \
  --extra "https://raw.githubusercontent.com/plotly/datasets/master/geojson-counties-fips.json@#c3bdb2:0.12" \
  -o poster.png

B=https://geo.datav.aliyun.com/areas_v3/bound
NE=https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson
python rotated_map.py --lat 31.7359 --lon 116.5218 --bearing 270 \
  --top-km 11 --bottom-lat -87 --res 10m --layers land,countries,lakes \
  --fill  "$B/340000_full.json@#f6ead2" \
  --fill  "$B/341500_full.json@#e7c79b" \
  --extra "$B/100000_full.json@#8f8f8f:0.45" \
  --extra "$B/340000_full.json@#b5b0a6:0.22" \
  --extra "$B/341500_full.json@#b09a76:0.20" \
  --label "$B/341500_full.json@name:#5f4c2e" \
  --label "$B/100000_full.json@name:#57606b" \
  --label "$NE/ne_10m_admin_0_countries.geojson@NAME_ZH:#57606b" \
  --label-skip 中华人民共和国 --label-skip 台湾 \
  -o luan.png

Three tints nest: Lu'an in ochre, the rest of Anhui in sand, everywhere else the ordinary land colour. --fill goes under the outlines, so a tinted region still shows its own county lines.

--label writes the names on. Type is sized by how big each feature came out on the page rather than how big it is on the ground, which is the whole point of this projection: Lu'an's seven districts and counties get the largest type on the sheet, and Xinjiang gets some of the smallest. Labels that would collide are nudged around inside their own shape first and dropped only if that fails, and earlier --label sources win — so listing the counties before the provinces keeps 舒城县 from being written over by 安徽省. Countries outside China come from Natural Earth's own NAME_ZH field, so the whole sheet reads in Chinese; --label-skip 中华人民共和国 stops it writing the country name over the provinces that already cover it.

An inland point is harder to frame than a coastal one, and it shows up in --bearing. The tear has to run somewhere, and around Lu'an every direction is land, so the top of the page is land whichever way you cut it. Scoring each bearing by how much of the tear's height on the page falls on land, the best you can do is 44% — due east, out through the East China Sea and then nothing but Pacific — against 91% for the worst. The price is that Shanghai, being almost exactly due east, ends up split across the two edges.

--bottom-lat matters more here too. Lu'an's antipode is in central Argentina rather than open ocean, so cutting the bottom at the usual -84 would take a bite out of South America; -87 keeps it.

The Chinese boundaries come from Aliyun DataV (province, prefecture and county level, nested). They sit about 1.8 km from Natural Earth's coastline, which is generalisation rather than a datum problem, and invisible at the scales where the two datasets meet.

Those numbers are not a guess. Given a picture of a rotated map you can recover the parameters that made it, because each one moves the picture in its own way: --bearing only slides the map horizontally (a roll about the pole is a cyclic shift in rotated longitude), --top-km only restretches it vertically, and the pole position changes the shape itself.

So: classify the target image into land and ocean, do the same to a candidate, and score them by intersection-over-union. Build the candidate by inverse mapping — for each output pixel, un-rotate to true lon/lat and look the answer up in a global land bitmap — and a candidate costs milliseconds rather than a render. Sweep the pole over a grid, and for each one try every horizontal shift with np.roll; the aspect ratio of the target pins --bottom-lat once --top-km is chosen.

For the poster above that converges on 37.71 N, 122.64 W — a point in the Pacific about 10 km off Ocean Beach, not downtown San Francisco — with the top edge 11 km out and a bearing of 38 degrees, at an IoU of 0.92. Downtown San Francisco only manages 0.78, which is what tipped me off that the original was not centred on the city.

The rotation itself is three lines of linear algebra: send the unit vector of (lat, lon) to +Z, then read longitude and latitude back off the rotated vector. Applying it to geometry rather than to points is the whole job.

• Straight lines stop being straight. A border that is a straight segment in lon/lat — the 49th parallel, Egypt's desert borders — is a curve once rotated, so every segment is densified before rotating and again after clipping.
• Longitudes need unwrapping, adaptively. Right beside the new pole a one-degree hop covers an enormous span of rotated longitude. Once two neighbouring vertices land more than 180° apart, unwrapping picks the wrong branch and smears the polygon into a band across the map. _ring_to_xy rotates, finds segments that swing too far, splits them in the source coordinates, and goes round again.
• Polygons that swallow the pole don't close. A continent containing your point lifts to a curve running from L to L ± 360, and the land is everything between it and the pole. Capping that at the top and folding it back into one period is not enough, because real coastlines wander tens of degrees past their own endpoints — fold naively and you either double-cover the overshoot or clip off real land. _cap tiles the lift over three periods (consecutive copies join end to end) and keeps the middle one, so the artificial edges are a whole turn away from anything you see.
• Holes have to be tiled along with the shell. _cap copies the outer ring over three periods; leave the interior rings in one copy and the other two tiles come out solid. That is how the Caspian Sea ends up painted as land — 43 grid samples' worth, which is small enough to miss by eye and obvious the moment you check.
• Clipping can hand back a GeometryCollection — the polygon you asked for plus a stray edge where the cut grazed the outline. Nothing downstream recognised that type, so at 110m the whole of Eurasia silently disappeared from a map centred on Lu'an. parts_of_kind flattens and keeps the dimension actually wanted.
• The map has an edge, and it must not be drawn twice. Every projection's outline comes from walking the boundary of the rotated domain, which gives a rectangle, an ellipse or a disc from one piece of code. But on an azimuthal the two meridians land on top of each other, and the resulting radius traced out and back has zero area — so it survives buffer(0), and stroking it paints a line from the centre to the rim. Layers drawn as outlines have the same problem from the other side: clip a country polygon at the antimeridian and the cut becomes part of its boundary, which then gets stroked as a hard seam across the map. They travel as lines instead.

There is a correctness check for all of this: sample a grid in rotated coordinates, un-rotate each sample, and ask the source polygons whether it should be land. At 10m resolution the rotated geometry and the source data disagree on 0 of 62,280 samples for most pole points.