oxedyne/fe2o3 — fe2o3_geom/src/proj.rs

72.3 KiB, 47 runs

created by r1870400018:19776, which is this file's identity for as long as the history lasts, whatever it is later renamed to

download · who wrote it · its history

1	//! Map projections: turning a position on a sphere into a position on a page, and back.
2	//!
3	//! The projection functions work on the plane every projection reference uses: right-handed,
4	//! `y` increasing northward, scaled by a sphere radius the caller passes. Equal Earth draws
5	//! the whole world at once, Web Mercator draws a flat map that street maps have made familiar,
6	//! and orthographic draws the half of the world a viewer is looking at, as a globe. Each is
7	//! checked against PROJ 9.7.1.
8	//!
9	//! A [`Viewport`] puts the last two on a screen: one camera -- a centre, metres per pixel,
10	//! a heading -- read by either projection, with `y` downward as a canvas has it. It projects
11	//! whole rings of unit vectors at a time, clipped and ready to paint, into a [`ScreenPaths`].
12	//! Further projections belong beside these in this module rather than in a module of their
13	//! own.
14
15	use crate::planar::Pt;
16
17	use oxedyne_fe2o3_core::prelude::*;
18
19	use std::f64::consts::{
20	FRAC_PI_2,
21	PI,
22	TAU,
23	};
24
25	// ---------------------------------------------------------------------------------------------
26	// Equal Earth
27	// ---------------------------------------------------------------------------------------------
28
29	/// First polynomial coefficient of the Equal Earth projection.
30	///
31	/// The four coefficients are those published by Šavrič, Patterson and Jenny, *The Equal
32	/// Earth map projection*, International Journal of Geographical Information Science 33(3),
33	/// 2019, and are the same values PROJ carries for `+proj=eqearth`.
34	pub const EQUAL_EARTH_A1: f64 = 1.340264;
35
36	/// Second polynomial coefficient of the Equal Earth projection.
37	pub const EQUAL_EARTH_A2: f64 = -0.081106;
38
39	/// Third polynomial coefficient of the Equal Earth projection.
40	pub const EQUAL_EARTH_A3: f64 = 0.000893;
41
42	/// Fourth polynomial coefficient of the Equal Earth projection.
43	pub const EQUAL_EARTH_A4: f64 = 0.003796;
44
45	/// Half the width of the projected world on a unit sphere.
46	///
47	/// The value of `x` at the equator on the antimeridian, which is where the map is widest.
48	pub const EQUAL_EARTH_HALF_WIDTH: f64 = 2.706_629_983_696_074_3;
49
50	/// Half the height of the projected world on a unit sphere.
51	///
52	/// The value of `y` at a pole. The world is therefore 2.0546 times as wide as it is tall.
53	pub const EQUAL_EARTH_HALF_HEIGHT: f64 = 1.317_362_759_157_413;
54
55	/// Projects a position onto the Equal Earth plane.
56	///
57	/// Equal Earth is pseudocylindrical and equal-area: every square metre of ground occupies
58	/// the same area on the page wherever it is, which is what makes it honest about how much
59	/// of the world a continent is. Parallels are straight and evenly spaced enough to read,
60	/// meridians are curved, and the shape of the land is close to what a person expects.
61	///
62	/// # Arguments
63	/// * `lat` - Latitude in signed decimal degrees, positive north.
64	/// * `lng` - Longitude in signed decimal degrees, positive east.
65	/// * `lon_0` - The central meridian, in degrees. Zero for Greenwich down the middle.
66	/// * `radius` - The radius of the sphere, in whatever units the result should come out in.
67	///
68	/// # Returns
69	/// The projected point, `x` eastward and `y` northward from the map's centre.
70	///
71	/// Latitude is clamped to the poles and longitude wrapped into a half turn either side of
72	/// `lon_0`, so a fix arriving one part in a billion outside its range projects rather than
73	/// producing a not-a-number. A non-finite argument yields a non-finite result: the caller
74	/// is the one that knows whether that is a hole in its data or an error.
75	pub fn equal_earth(lat: f64, lng: f64, lon_0: f64, radius: f64) -> Pt {
76	let phi = lat.clamp(-90.0, 90.0).to_radians();
77	let lam = wrap_half_turn(lng - lon_0).to_radians();
78
79	// The parametric latitude, on which the whole projection is a polynomial.
80	let theta = ((3.0f64).sqrt() / 2.0 * phi.sin()).clamp(-1.0, 1.0).asin();
81	let t2 = theta * theta;
82	let t3 = t2 * theta;
83	let t6 = t3 * t3;
84	let t7 = t6 * theta;
85	let t8 = t6 * t2;
86	let t9 = t8 * theta;
87
88	// The denominator of x is dy/dθ, which is what makes the projection equal-area.
89	let dy = 9.0 * EQUAL_EARTH_A4 * t8
90	+ 7.0 * EQUAL_EARTH_A3 * t6
91	+ 3.0 * EQUAL_EARTH_A2 * t2
92	+ EQUAL_EARTH_A1;
93
94	let x = 2.0 * (3.0f64).sqrt() * lam * theta.cos() / (3.0 * dy);
95	let y = EQUAL_EARTH_A4 * t9
96	+ EQUAL_EARTH_A3 * t7
97	+ EQUAL_EARTH_A2 * t3
98	+ EQUAL_EARTH_A1 * theta;
99
100	Pt::new(radius * x, radius * y)
101	}
102
103	/// Projects a position onto the Equal Earth plane of a unit sphere, Greenwich centred.
104	///
105	/// The bare form, for a caller that will scale the result itself.
106	pub fn equal_earth_unit(lat: f64, lng: f64) -> Pt {
107	equal_earth(lat, lng, 0.0, 1.0)
108	}
109
110	/// The sphere radius that makes an Equal Earth map `2 * half_width` across.
111	///
112	/// A map is usually specified by the box it has to fill rather than by the size of the
113	/// world it draws, and this converts one into the other. The height that follows is
114	/// `half_width * EQUAL_EARTH_HALF_HEIGHT / EQUAL_EARTH_HALF_WIDTH`, or very nearly half
115	/// the width.
116	pub fn equal_earth_radius_for_half_width(half_width: f64) -> f64 {
117	half_width / EQUAL_EARTH_HALF_WIDTH
118	}
119
120	// ---------------------------------------------------------------------------------------------
121	// Orthographic
122	// ---------------------------------------------------------------------------------------------
123
124	/// Half the width of the orthographic projection on a unit sphere.
125	///
126	/// The projection is the sphere seen from infinitely far away, so the map is the sphere's own
127	/// disc and its half width is the radius. A caller wanting a globe `2 * half_width` across
128	/// therefore passes `half_width` as the radius; there is no conversion to do.
129	pub const ORTHOGRAPHIC_HALF_WIDTH: f64 = 1.0;
130
131	/// Projects a position onto the orthographic plane, as seen from above `lat_0`, `lon_0`.
132	///
133	/// The orthographic projection is what a sphere looks like from far enough away that the rays
134	/// arrive parallel: the near hemisphere fills a disc, foreshortened towards the rim, and the
135	/// far hemisphere is behind it. It is neither equal-area nor conformal, and it is the only
136	/// projection that does not have to choose, because it is not flattening the world at all —
137	/// it is drawing the object. That makes it the honest one for a globe a viewer turns.
138	///
139	/// # Arguments
140	/// * `lat` - Latitude in signed decimal degrees, positive north.
141	/// * `lng` - Longitude in signed decimal degrees, positive east.
142	/// * `lat_0` - The latitude the viewer is above, in degrees.
143	/// * `lon_0` - The longitude the viewer is above, in degrees.
144	/// * `radius` - The radius of the sphere, in whatever units the result should come out in.
145	///
146	/// # Returns
147	/// The projected point, `x` rightward and `y` upward from the centre of the disc.
148	///
149	/// A position on the far hemisphere still projects, onto the point of the near hemisphere
150	/// directly in front of it, because the formula cannot tell them apart. Ask
151	/// [`orthographic_cos_c`] which side of the globe a position is on; it is separate so that a
152	/// caller clipping a coastline can interpolate along an edge to the horizon, where that
153	/// cosine is zero, rather than being handed a hole.
154	pub fn orthographic(lat: f64, lng: f64, lat_0: f64, lon_0: f64, radius: f64) -> Pt {
155	let phi = lat.clamp(-90.0, 90.0).to_radians();
156	let phi_0 = lat_0.clamp(-90.0, 90.0).to_radians();
157	let lam = wrap_half_turn(lng - lon_0).to_radians();
158
159	let x = phi.cos() * lam.sin();
160	let y = phi_0.cos() * phi.sin() - phi_0.sin() * phi.cos() * lam.cos();
161
162	Pt::new(radius * x, radius * y)
163	}
164
165	/// The cosine of the angle between a position and the centre of an orthographic projection.
166	///
167	/// One where the position is under the viewer, zero on the horizon, and negative on the far
168	/// side of the globe. A caller drawing a coastline keeps the vertices where this is positive
169	/// and finds the horizon crossing by interpolating an edge to where it is zero.
170	///
171	/// # Arguments
172	/// * `lat` - Latitude in signed decimal degrees, positive north.
173	/// * `lng` - Longitude in signed decimal degrees, positive east.
174	/// * `lat_0` - The latitude the viewer is above, in degrees.
175	/// * `lon_0` - The longitude the viewer is above, in degrees.
176	pub fn orthographic_cos_c(lat: f64, lng: f64, lat_0: f64, lon_0: f64) -> f64 {
177	let phi = lat.clamp(-90.0, 90.0).to_radians();
178	let phi_0 = lat_0.clamp(-90.0, 90.0).to_radians();
179	let lam = wrap_half_turn(lng - lon_0).to_radians();
180	phi_0.sin() * phi.sin() + phi_0.cos() * phi.cos() * lam.cos()
181	}
182
183	/// Brings an angle in degrees into the half turn either side of zero.
184	///
185	/// A longitude difference of 190° east is 170° west, and a map has to draw it there.
186	fn wrap_half_turn(deg: f64) -> f64 {
187	if !deg.is_finite() {
188	return deg;
189	}
190	let full = 360.0;
191	let mut d = deg % full;
192	if d > 180.0 {
193	d -= full;
194	} else if d < -180.0 {
195	d += full;
196	}
197	d
198	}
199
200	// ---------------------------------------------------------------------------------------------
201	// Web Mercator
202	// ---------------------------------------------------------------------------------------------
203
204	pub const EARTH_RADIUS_M: f64 = 6_371_008.8; // IUGG mean radius
205	pub const WEB_MERCATOR_MAX_LAT: f64 = 85.051_128_779_806_59; // atan(sinh(pi)): the map is square
206
207	/// Projects a position onto the spherical Web Mercator plane.
208	///
209	/// Web Mercator is conformal, so a small square of ground stays square on the page at every
210	/// latitude, and it is the projection street maps have taught everyone to read. Its price is
211	/// the poles, which are infinitely far away: latitude is clamped to
212	/// [`WEB_MERCATOR_MAX_LAT`], where the map is exactly as tall as it is wide, `2 * pi * radius`.
213	/// Longitude is wrapped into a half turn either side of `lon_0`.
214	pub fn web_mercator(lat: f64, lng: f64, lon_0: f64, radius: f64) -> Pt {
215	let phi = lat.clamp(-WEB_MERCATOR_MAX_LAT, WEB_MERCATOR_MAX_LAT).to_radians();
216	let lam = wrap_half_turn(lng - lon_0).to_radians();
217	// PROJ writes the ordinate as asinh(tan(phi)); it equals ln(tan(pi/4 + phi/2)).
218	Pt::new(radius * lam, radius * phi.tan().asinh())
219	}
220
221	/// The position a Web Mercator point came from, as `(lat, lng)` in degrees.
222	///
223	/// Longitude comes back wrapped into `[-180, 180]`, so a point on a copy of the world to either
224	/// side of the central one answers with the same position as its twin. A point above or below
225	/// the square map still answers, with a latitude past [`WEB_MERCATOR_MAX_LAT`]; a caller that
226	/// draws only the square asks [`Viewport::inverse`] instead.
227	pub fn web_mercator_inverse(pt: Pt, lon_0: f64, radius: f64) -> (f64, f64) {
228	let lat = (pt.y / radius).sinh().atan().to_degrees();
229	let lng = wrap_half_turn(lon_0 + (pt.x / radius).to_degrees());
230	(lat, lng)
231	}
232
233	/// The position an orthographic point came from, as `(lat, lng)` in degrees.
234	///
235	/// The inverse of [`orthographic`] for the near hemisphere, which is the only one a point on
236	/// the disc can mean. A point off the disc is off the globe and answers `None`; one within a
237	/// part in ten billion of the rim is taken to be on it.
238	pub fn orthographic_inverse(pt: Pt, lat_0: f64, lon_0: f64, radius: f64) -> Option<(f64, f64)> {
239	let rho = pt.x.hypot(pt.y) / radius;
240	if !rho.is_finite() \|\| rho > 1.0 + 1.0e-10 {
241	return None;
242	}
243	let phi_0 = lat_0.clamp(-90.0, 90.0).to_radians();
244	if rho < 1.0e-15 {
245	return Some((lat_0.clamp(-90.0, 90.0), wrap_half_turn(lon_0)));
246	}
247	let sin_c = rho.min(1.0);
248	let cos_c = (1.0 - sin_c * sin_c).max(0.0).sqrt();
249	let (x, y) = (pt.x / radius, pt.y / radius);
250	// Snyder, Map Projections: A Working Manual, equations 20-14 and 20-15.
251	let phi = (cos_c * phi_0.sin() + y * sin_c * phi_0.cos() / rho).clamp(-1.0, 1.0).asin();
252	let lam = (x * sin_c).atan2(rho * cos_c * phi_0.cos() - y * sin_c * phi_0.sin());
253	Some((phi.to_degrees(), wrap_half_turn(lon_0 + lam.to_degrees())))
254	}
255
256	// ---------------------------------------------------------------------------------------------
257	// Unit vectors
258	// ---------------------------------------------------------------------------------------------
259
260	/// A position as a unit vector: `x` towards Greenwich on the equator, `z` towards the north
261	/// pole. The same frame [`crate::cell`] uses.
262	pub fn unit_vec(lat: f64, lng: f64) -> [f64; 3] {
263	let (a, b) = (lat.to_radians(), lng.to_radians());
264	let c = a.cos();
265	[c * b.cos(), c * b.sin(), a.sin()]
266	}
267
268	/// The position of a vector of any non-zero length, as `(lat, lng)` in degrees.
269	pub fn vec_lat_lng(v: &[f64; 3]) -> (f64, f64) {
270	let len = (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]).sqrt();
271	let z = (v[2] / len).clamp(-1.0, 1.0);
272	(z.asin().to_degrees(), v[1].atan2(v[0]).to_degrees())
273	}
274
275	fn dot3(a: &[f64; 3], b: &[f64; 3]) -> f64 { a[0] * b[0] + a[1] * b[1] + a[2] * b[2] }
276
277	// ---------------------------------------------------------------------------------------------
278	// Viewport
279	// ---------------------------------------------------------------------------------------------
280
281	/// The projection a [`Viewport`] draws with.
282	#[derive(Clone, Copy, Debug, PartialEq, Eq)]
283	pub enum Projection {
284	WebMercator, // a flat map, north up unless turned, wrapping east-west without end
285	Orthographic, // a globe seen from above the centre
286	}
287
288	/// How [`Viewport::project_rings`] treats each ring it is given.
289	#[derive(Clone, Copy, Debug, PartialEq, Eq)]
290	pub enum RingMode {
291	Fill, // closed ring, painted: clipped pieces are stitched along the clip edge
292	Outline, // closed ring, stroked: its closing edge is drawn, and it is cut where clipped
293	Line, // open line, stroked, cut where clipped
294	}
295
296	/// One camera over the sphere, shared by the flat map and the globe.
297	///
298	/// A viewport is the centre being looked at, how much ground a pixel covers there, which way
299	/// is up, and the size of the screen. The two projections read the same numbers, so a caller
300	/// switching between them keeps the place and the scale and changes only the picture.
301	///
302	/// Screen coordinates are pixels from the top left corner, `y` downward, as a canvas has them.
303	/// `heading` is the compass bearing, in degrees clockwise from north, that points up the
304	/// screen; zero is north up. `m_per_px` is ground metres per pixel at the centre, on a sphere
305	/// of [`EARTH_RADIUS_M`], so a Mercator map at 60 degrees draws the world twice as wide as the
306	/// same `m_per_px` does at the equator, and a cell the same size on the ground stays the same
307	/// size on the screen either way.
308	#[derive(Clone, Copy, Debug, PartialEq)]
309	pub struct Viewport {
310	pub kind: Projection,
311	pub lat_0: f64, // degrees, the centre of the screen
312	pub lon_0: f64, // degrees
313	pub heading: f64, // degrees clockwise from north, pointing up the screen
314	pub m_per_px: f64, // ground metres per pixel at the centre
315	pub w: f64, // screen width in pixels
316	pub h: f64, // screen height in pixels
317	}
318
319	/// The clip and the margin around the screen, in pixels, inside which geometry is kept.
320	///
321	/// Wide enough that a stroke's own width never shows the cut, and small enough that nothing
322	/// far off the screen is carried.
323	const CLIP_MARGIN_PX: f64 = 8.0;
324
325	/// A viewport's numbers turned into the arithmetic that draws with them, once per call.
326	pub(crate) struct Frame {
327	pub(crate) cx: f64, // screen centre
328	pub(crate) cy: f64,
329	pub(crate) ch: f64, // cosine and sine of the heading
330	pub(crate) sh: f64,
331	pub(crate) r_px: f64, // the sphere's radius in pixels (Mercator: at the equator)
332	pub(crate) c: [f64; 3], // orthographic basis: centre, east, north
333	pub(crate) e: [f64; 3],
334	pub(crate) n: [f64; 3],
335	pub(crate) cos_clip: f64, // orthographic clip cap, as the cosine of its radius
336	pub(crate) rho_clip: f64, // and as a circle on the screen, in pixels
337	pub(crate) lam_0: f64, // Mercator central meridian, radians
338	pub(crate) y_0: f64, // Mercator ordinate of the centre, pixels
339	pub(crate) x_lo: f64, // Mercator: the screen and its margin on the unturned plane
340	pub(crate) x_hi: f64,
341	pub(crate) y_lo: f64,
342	pub(crate) y_hi: f64,
343	}
344
345	impl Frame {
346	/// A point on the unturned plane, `yd` already downward, to the screen.
347	pub(crate) fn to_screen(&self, x: f64, yd: f64) -> (f64, f64) {
348	(self.cx + x * self.ch + yd * self.sh, self.cy - x * self.sh + yd * self.ch)
349	}
350
351	/// A screen point back to the unturned plane, `yd` downward.
352	pub(crate) fn from_screen(&self, sx: f64, sy: f64) -> (f64, f64) {
353	let (dx, dy) = (sx - self.cx, sy - self.cy);
354	(dx * self.ch - dy * self.sh, dx * self.sh + dy * self.ch)
355	}
356
357	/// An orthographic vertex: where it lands on the unturned plane, and how far in front of
358	/// the clip it is. The last is linear in the vector, which is what lets a crossing be
359	/// found by one division along the straight edge drawn between two vertices.
360	fn ortho(&self, p: &[f64; 3]) -> (f64, f64, f64) {
361	(self.r_px * dot3(p, &self.e), -self.r_px * dot3(p, &self.n), dot3(p, &self.c) - self.cos_clip)
362	}
363
364	/// A Mercator vertex's longitude from the central meridian, in radians, and its ordinate,
365	/// both unscaled. `None` for a pole, whose longitude is undefined.
366	fn merc(&self, p: &[f64; 3]) -> Option<(f64, f64)> {
367	if p[0] * p[0] + p[1] * p[1] < 1.0e-24 {
368	return None;
369	}
370	let len = (p[0] * p[0] + p[1] * p[1] + p[2] * p[2]).sqrt();
371	let y = (p[2] / len).clamp(-1.0, 1.0).atanh().clamp(-PI, PI);
372	Some((wrap_pi(p[1].atan2(p[0]) - self.lam_0), y))
373	}
374	}
375
376	/// Brings an angle in radians into the half turn either side of zero.
377	fn wrap_pi(a: f64) -> f64 {
378	let mut d = a % TAU;
379	if d > PI {
380	d -= TAU;
381	} else if d < -PI {
382	d += TAU;
383	}
384	d
385	}
386
387	impl Viewport {
388	/// A viewport, refusing numbers that could not draw anything.
389	pub fn new(
390	kind: Projection,
391	lat_0: f64,
392	lon_0: f64,
393	heading: f64,
394	m_per_px: f64,
395	w: f64,
396	h: f64,
397	)
398	-> Outcome<Self>
399	{
400	let v = Self { kind, lat_0, lon_0, heading, m_per_px, w, h };
401	res!(v.check());
402	Ok(v)
403	}
404
405	pub(crate) fn check(&self) -> Outcome<()> {
406	let finite = self.lat_0.is_finite() && self.lon_0.is_finite() && self.heading.is_finite();
407	if !finite {
408	return Err(err!("A viewport centred on {}, {} with heading {} is not a place.",
409	self.lat_0, self.lon_0, self.heading; Invalid, Input));
410	}
411	if !(self.m_per_px > 0.0 && self.m_per_px.is_finite()) {
412	return Err(err!("A viewport of {} m per pixel draws nothing.", self.m_per_px;
413	Invalid, Input, Range));
414	}
415	if !(self.w > 0.0 && self.h > 0.0 && self.w.is_finite() && self.h.is_finite()) {
416	return Err(err!("A viewport {} by {} pixels has no screen.", self.w, self.h;
417	Invalid, Input, Range));
418	}
419	Ok(())
420	}
421
422	pub(crate) fn frame(&self) -> Frame {
423	let lat_0 = self.lat_0.clamp(-90.0, 90.0);
424	let (phi, lam) = (lat_0.to_radians(), self.lon_0.to_radians());
425	let hd = self.heading.to_radians();
426	let (cx, cy) = (self.w / 2.0, self.h / 2.0);
427	let reach = cx.hypot(cy) + CLIP_MARGIN_PX;
428	let mut f = Frame {
429	cx, cy, ch: hd.cos(), sh: hd.sin(), r_px: 0.0,
430	c: [0.0; 3], e: [0.0; 3], n: [0.0; 3], cos_clip: 0.0, rho_clip: 0.0,
431	lam_0: lam, y_0: 0.0, x_lo: 0.0, x_hi: 0.0, y_lo: 0.0, y_hi: 0.0,
432	};
433	match self.kind {
434	Projection::Orthographic => {
435	f.r_px = EARTH_RADIUS_M / self.m_per_px;
436	f.c = unit_vec(lat_0, self.lon_0);
437	f.e = [-lam.sin(), lam.cos(), 0.0];
438	f.n = [-phi.sin() * lam.cos(), -phi.sin() * lam.sin(), phi.cos()];
439	// Clip to the smaller of the horizon and the circle round the screen, so that
440	// close in, a ring is cut just off the screen rather than at a horizon
441	// thousands of pixels away.
442	if reach >= f.r_px {
443	f.cos_clip = 0.0;
444	f.rho_clip = f.r_px;
445	} else {
446	let s = reach / f.r_px;
447	f.cos_clip = (1.0 - s * s).sqrt();
448	f.rho_clip = reach;
449	}
450	},
451	Projection::WebMercator => {
452	let phi_c = lat_0.clamp(-WEB_MERCATOR_MAX_LAT, WEB_MERCATOR_MAX_LAT).to_radians();
453	// The Mercator scale factor at the centre is sec(phi), so the plane is drawn
454	// cos(phi) as large as the ground resolution alone would say.
455	f.r_px = EARTH_RADIUS_M * phi_c.cos() / self.m_per_px;
456	f.y_0 = f.r_px * phi_c.tan().asinh();
457	let m = CLIP_MARGIN_PX;
458	let corners = [(-m, -m), (self.w + m, -m), (-m, self.h + m), (self.w + m, self.h + m)];
459	let (mut xl, mut xh, mut yl, mut yh) = (f64::MAX, f64::MIN, f64::MAX, f64::MIN);
460	for (sx, sy) in corners {
461	let (x, yd) = f.from_screen(sx, sy);
462	xl = xl.min(x);
463	xh = xh.max(x);
464	yl = yl.min(-yd);
465	yh = yh.max(-yd);
466	}
467	f.x_lo = xl;
468	f.x_hi = xh;
469	f.y_lo = yl;
470	f.y_hi = yh;
471	},
472	}
473	f
474	}
475
476	/// Where a position lands on the screen, or `None` behind the globe.
477	///
478	/// On the flat map every position lands somewhere, on the copy of the world nearest the
479	/// centre; whether that is within the screen is the caller's question.
480	pub fn forward(&self, lat: f64, lng: f64) -> Option<Pt> {
481	if self.check().is_err() {
482	return None;
483	}
484	let f = self.frame();
485	self.forward_in(&f, lat, lng)
486	}
487
488	/// [`Viewport::forward`] with the frame already made, for a caller projecting many points.
489	pub(crate) fn forward_in(&self, f: &Frame, lat: f64, lng: f64) -> Option<Pt> {
490	match self.kind {
491	Projection::Orthographic => {
492	let p = unit_vec(lat, lng);
493	if dot3(&p, &f.c) < 0.0 {
494	return None;
495	}
496	let (x, yd, _) = f.ortho(&p);
497	let (sx, sy) = f.to_screen(x, yd);
498	Some(Pt::new(sx, sy))
499	},
500	Projection::WebMercator => {
501	let q = web_mercator(lat, lng, self.lon_0, f.r_px);
502	let (sx, sy) = f.to_screen(q.x, -(q.y - f.y_0));
503	Some(Pt::new(sx, sy))
504	},
505	}
506	}
507
508	/// The position under a screen point, as `(lat, lng)` in degrees, or `None` where the
509	/// point is off the globe or above or below the square flat map.
510	pub fn inverse(&self, pt: Pt) -> Option<(f64, f64)> {
511	if self.check().is_err() {
512	return None;
513	}
514	let f = self.frame();
515	let (x, yd) = f.from_screen(pt.x, pt.y);
516	match self.kind {
517	Projection::Orthographic =>
518	orthographic_inverse(Pt::new(x, -yd), self.lat_0, self.lon_0, f.r_px),
519	Projection::WebMercator => {
520	let y = -yd + f.y_0;
521	if y.abs() > PI * f.r_px {
522	return None;
523	}
524	Some(web_mercator_inverse(Pt::new(x, y), self.lon_0, f.r_px))
525	},
526	}
527	}
528
529	/// A spherical cap holding everything the screen shows, as its centre and its angular
530	/// radius in radians.
531	///
532	/// It is what a caller asks for the cells in view ([`crate::cell::cover_cap`]) and what
533	/// culls a ring before it is projected. On the globe it is exact: the circle round the
534	/// screen, or the horizon if that is nearer. On the flat map it is the furthest point of
535	/// the screen's edge from the centre, with two per cent to spare, and the whole sphere
536	/// once the screen spans half the world either side.
537	pub fn bounding_cap(&self) -> ([f64; 3], f64) {
538	let lat_0 = self.lat_0.clamp(-90.0, 90.0);
539	let centre = unit_vec(lat_0, self.lon_0);
540	if self.check().is_err() {
541	return (centre, PI);
542	}
543	let f = self.frame();
544	match self.kind {
545	Projection::Orthographic => {
546	let reach = f.cx.hypot(f.cy) + CLIP_MARGIN_PX;
547	if reach >= f.r_px {
548	(centre, FRAC_PI_2)
549	} else {
550	(centre, (reach / f.r_px).asin())
551	}
552	},
553	Projection::WebMercator => {
554	let half = PI * f.r_px;
555	if f.x_lo <= -half \|\| f.x_hi >= half {
556	return (centre, PI);
557	}
558	// Walk the screen's edge, clamped to the square map, for the furthest point.
559	const STEPS: usize = 64;
560	let m = CLIP_MARGIN_PX;
561	let mut far: f64 = 0.0;
562	for k in 0..(4 * STEPS) {
563	let t = (k % STEPS) as f64 / STEPS as f64;
564	let (sx, sy) = match k / STEPS {
565	0 => (-m + t * (self.w + 2.0 * m), -m),
566	1 => (self.w + m, -m + t * (self.h + 2.0 * m)),
567	2 => (self.w + m - t * (self.w + 2.0 * m), self.h + m),
568	_ => (-m, self.h + m - t * (self.h + 2.0 * m)),
569	};
570	let (x, yd) = f.from_screen(sx, sy);
571	let y = (-yd + f.y_0).clamp(-half, half);
572	let (lat, lng) = web_mercator_inverse(Pt::new(x, y), self.lon_0, f.r_px);
573	let p = unit_vec(lat, lng);
574	far = far.max(dot3(&p, &centre).clamp(-1.0, 1.0).acos());
575	}
576	(centre, (far * 1.02 + 1.0e-9).min(PI))
577	},
578	}
579	}
580
581	/// Projects rings of unit vectors onto the screen, clipped, into `out`.
582	///
583	/// This is the one routine that draws a coastline, a border or a cell outline, on either
584	/// projection, so that a caller painting a frame hands over its rings and gets back
585	/// screen paths with nothing left to work out.
586	///
587	/// On the globe a ring is clipped against a cap: the horizon, or the circle round the
588	/// screen when that is nearer. A filled ring's visible stretches are stitched together
589	/// along that circle's rim in the ring's own winding, so a continent cut by the horizon
590	/// is still a closed shape, and a ring with nothing showing still fills the screen when the
591	/// view is inside it. On the flat map a ring is unwrapped across the antimeridian, closed
592	/// round the pole it encircles if it encircles one, repeated for every copy of the world
593	/// the screen shows, and clipped to the screen.
594	///
595	/// # Arguments
596	/// * `tol_m` - How far the rings may already stray from the truth, in metres of ground: a
597	/// simplified coastline's own tolerance, or zero for exact rings such as cell outlines.
598	/// A visible stretch of a filled ring no larger than this is noise from simplification
599	/// and is dropped, so that a ring which crosses itself cannot break the stitching.
600	/// * `eps_px` - A vertex this close to the one before it on the screen is dropped.
601	pub fn project_rings<R: AsRef<[[f64; 3]]>>(
602	&self,
603	rings: &[R],
604	mode: RingMode,
605	tol_m: f64,
606	eps_px: f64,
607	out: &mut ScreenPaths,
608	)
609	-> Outcome<()>
610	{
611	res!(self.check());
612	let f = self.frame();
613	let eps = if eps_px.is_finite() { eps_px.max(0.0) } else { 0.0 };
614	let mut pen = Pen { out, eps2: eps * eps, open: false };
615	match self.kind {
616	Projection::Orthographic => {
617	let slack = if tol_m.is_finite() { f.r_px * tol_m.max(0.0) / EARTH_RADIUS_M } else { 0.0 };
618	// The rim is walked in steps whose sagitta is under the vertex tolerance.
619	let sag = eps.max(0.25).min(f.rho_clip);
620	let step = (2.0 * (1.0 - sag / f.rho_clip).acos()).clamp(1.0e-4, 0.1);
621	for ring in rings {
622	let ring = ring.as_ref();
623	match mode {
624	RingMode::Fill => ortho_fill(&f, ring, slack, step, &mut pen),
625	_ => ortho_line(&f, ring, mode == RingMode::Outline, &mut pen),
626	}
627	}
628	},
629	Projection::WebMercator => {
630	for ring in rings {
631	merc_ring(&f, ring.as_ref(), mode, &mut pen);
632	}
633	},
634	}
635	Ok(())
636	}
637	}
638
639	// ---------------------------------------------------------------------------------------------
640	// Screen paths
641	// ---------------------------------------------------------------------------------------------
642
643	/// One path in a [`ScreenPaths`]: a run of its points, and whether to close it.
644	#[derive(Clone, Copy, Debug, PartialEq, Eq)]
645	pub struct PathSpan {
646	pub start: u32, // index of the first point, in points rather than floats
647	pub len: u32, // number of points
648	pub closed: bool, // join the last point back to the first
649	}
650
651	/// Screen paths as flat `f32` pairs, the shape a canvas painter or a typed array wants.
652	///
653	/// Kept by the caller and cleared between frames, so a frame allocates nothing once the
654	/// buffers have grown to fit.
655	#[derive(Clone, Debug, Default)]
656	pub struct ScreenPaths {
657	pub xy: Vec<f32>, // x0, y0, x1, y1, ... in screen pixels
658	pub spans: Vec<PathSpan>,
659	}
660
661	impl ScreenPaths {
662	pub fn new() -> Self { Self::default() }
663
664	pub fn clear(&mut self) {
665	self.xy.clear();
666	self.spans.clear();
667	}
668
669	pub fn len(&self) -> usize { self.spans.len() }
670
671	pub fn is_empty(&self) -> bool { self.spans.is_empty() }
672
673	/// The points of the `i`th path as `x, y` pairs, and whether it is closed.
674	pub fn path(&self, i: usize) -> Option<(&[f32], bool)> {
675	let s = match self.spans.get(i) {
676	Some(s) => s,
677	None => return None,
678	};
679	let a = 2 * s.start as usize;
680	let b = a + 2 * s.len as usize;
681	match self.xy.get(a..b) {
682	Some(pts) => Some((pts, s.closed)),
683	None => None,
684	}
685	}
686	}
687
688	/// Writes paths into a [`ScreenPaths`], dropping vertices within the tolerance of the last.
689	struct Pen<'a> {
690	out: &'a mut ScreenPaths,
691	eps2: f64,
692	open: bool,
693	}
694
695	impl<'a> Pen<'a> {
696	fn begin(&mut self) {
697	if self.open {
698	self.end(false);
699	}
700	let start = (self.out.xy.len() / 2) as u32;
701	self.out.spans.push(PathSpan { start, len: 0, closed: false });
702	self.open = true;
703	}
704
705	fn push(&mut self, x: f64, y: f64) {
706	if !self.open {
707	self.begin();
708	}
709	let n = self.out.xy.len();
710	if let Some(s) = self.out.spans.last() {
711	if s.len > 0 {
712	let dx = x - self.out.xy[n - 2] as f64;
713	let dy = y - self.out.xy[n - 1] as f64;
714	if dx * dx + dy * dy <= self.eps2 {
715	return;
716	}
717	}
718	}
719	self.out.xy.push(x as f32);
720	self.out.xy.push(y as f32);
721	if let Some(s) = self.out.spans.last_mut() {
722	s.len += 1;
723	}
724	}
725
726	/// Ends the open path, discarding it if it has too few points to draw.
727	fn end(&mut self, closed: bool) {
728	if !self.open {
729	return;
730	}
731	self.open = false;
732	let least = if closed { 3 } else { 2 };
733	let keep = match self.out.spans.last_mut() {
734	Some(s) => {
735	s.closed = closed;
736	s.len >= least
737	},
738	None => return,
739	};
740	if !keep {
741	if let Some(s) = self.out.spans.pop() {
742	self.out.xy.truncate(2 * s.start as usize);
743	}
744	}
745	}
746	}
747
748	// ---------------------------------------------------------------------------------------------
749	// The globe: clip to a cap and walk the rim
750	// ---------------------------------------------------------------------------------------------
751	//
752	// Lifted from Ochre's globe (`web/apps/ochre/src/web/globe.rs`, August 2026), which drew it
753	// against the horizon alone and verified it against PROJ and against the sea staying sea at
754	// every turn. Two things are new here: the clip cap may be smaller than the horizon, so that
755	// close in the rim being walked is the circle just round the screen rather than one thousands
756	// of pixels across; and the rim is walked out as points, so the result is plain polylines.
757
758	/// Which way round a ring is wound, seen from outside the sphere.
759	///
760	/// Measured by the vector area of the ring's own chords, whose component along the ring's own
761	/// centre is twice the signed area it covers. A polygon's outline and the holes in it come in
762	/// the same list of rings, and a hole winds the other way.
763	fn wound_cw(ring: &[[f64; 3]]) -> bool {
764	let mut sum = [0.0; 3];
765	let mut area = [0.0; 3];
766	let n = ring.len();
767	for i in 0..n {
768	let a = &ring[i];
769	let b = &ring[(i + 1) % n];
770	sum = [sum[0] + a[0], sum[1] + a[1], sum[2] + a[2]];
771	area = [
772	area[0] + a[1] * b[2] - a[2] * b[1],
773	area[1] + a[2] * b[0] - a[0] * b[2],
774	area[2] + a[0] * b[1] - a[1] * b[0],
775	];
776	}
777	let m = dot3(&sum, &sum).sqrt().max(f64::EPSILON);
778	dot3(&area, &sum) / m < 0.0
779	}
780
781	/// Is the view inside a ring none of whose boundary shows?
782	///
783	/// The bearing of the ring about the centre is swept right round, and a ring wound the
784	/// shapefile way sweeps a whole turn forwards when the view is inside it. A view whose
785	/// antipode is inside sweeps a whole turn backwards, so the direction is the answer and the
786	/// size alone is not.
787	fn inside(seen: &[(f64, f64, f64)], cw: bool) -> bool {
788	let mut sweep = 0.0;
789	let mut last = 0.0;
790	for i in 0..=seen.len() {
791	let p = seen[i % seen.len()];
792	let angle = p.1.atan2(p.0);
793	if i > 0 {
794	let mut step = angle - last;
795	if step > PI {
796	step -= TAU;
797	}
798	if step < -PI {
799	step += TAU;
800	}
801	sweep += step;
802	}
803	last = angle;
804	}
805	if cw { sweep > PI } else { sweep < -PI }
806	}
807
808	/// A visible stretch of a filled ring, between two crossings of the clip.
809	struct Run {
810	pts: Vec<(f64, f64)>, // on the unturned plane, the ends on the clip circle
811	entry: f64, // angle on the clip circle where it starts
812	exit: f64, // and where it ends
813	done: bool, // already drawn into a loop
814	}
815
816	/// Where the edge from `was` to `here` crosses the clip, on the clip circle.
817	fn crossing(f: &Frame, was: (f64, f64, f64), here: (f64, f64, f64)) -> (f64, f64, f64) {
818	let s = was.2 / (was.2 - here.2);
819	let (ex, ey) = (was.0 + s * (here.0 - was.0), was.1 + s * (here.1 - was.1));
820	// The crossing lies on the chord, a hair inside the sphere, and is put on the rim.
821	let m = ex.hypot(ey).max(f64::EPSILON);
822	let (ex, ey) = (ex / m * f.rho_clip, ey / m * f.rho_clip);
823	(ex, ey, ey.atan2(ex))
824	}
825
826	/// Walks the clip circle from angle `from` through `span` radians (signed), excluding both
827	/// ends.
828	fn rim(f: &Frame, from: f64, span: f64, step: f64, pen: &mut Pen) {
829	let k = (span.abs() / step).ceil() as usize;
830	for i in 1..k {
831	let a = from + span * i as f64 / k as f64;
832	let (sx, sy) = f.to_screen(f.rho_clip * a.cos(), f.rho_clip * a.sin());
833	pen.push(sx, sy);
834	}
835	}
836
837	fn ortho_fill(f: &Frame, ring: &[[f64; 3]], slack: f64, step: f64, pen: &mut Pen) {
838	let n = ring.len();
839	if n < 3 {
840	return;
841	}
842	let seen: Vec<(f64, f64, f64)> = ring.iter().map(\|p\| f.ortho(p)).collect();
843	let first = match seen.iter().position(\|p\| p.2 >= 0.0) {
844	Some(at) => at,
845	None => {
846	// Nothing of the boundary shows: the ring is out of sight, or the whole view is
847	// inside it and the clip circle is filled.
848	let cw = wound_cw(ring);
849	if inside(&seen, cw) {
850	let way = if cw { TAU } else { -TAU };
851	pen.begin();
852	let (sx, sy) = f.to_screen(f.rho_clip, 0.0);
853	pen.push(sx, sy);
854	rim(f, 0.0, way, step, pen);
855	pen.end(true);
856	}
857	return;
858	},
859	};
860	let cw = wound_cw(ring);
861
862	// Gather the visible stretches, each ending exactly where the ring meets the clip.
863	let mut runs: Vec<Run> = Vec::new();
864	let mut open = false;
865	let mut hidden = false;
866	for k in 0..=n {
867	let here = seen[(first + k) % n];
868	if k > 0 {
869	let was = seen[(first + k - 1) % n];
870	if (was.2 >= 0.0) != (here.2 >= 0.0) {
871	let (ex, ey, angle) = crossing(f, was, here);
872	if here.2 < 0.0 {
873	if let Some(run) = runs.last_mut() {
874	run.pts.push((ex, ey));
875	run.exit = angle;
876	}
877	open = false;
878	} else if k == n {
879	// The wrap: this is where the first run began.
880	runs[0].pts.insert(0, (ex, ey));
881	runs[0].entry = angle;
882	} else {
883	runs.push(Run { pts: vec![(ex, ey)], entry: angle, exit: 0.0, done: false });
884	open = true;
885	}
886	}
887	}
888	if here.2 < 0.0 {
889	hidden = true;
890	continue;
891	}
892	if k == n {
893	break;
894	}
895	if !open {
896	runs.push(Run { pts: Vec::new(), entry: 0.0, exit: 0.0, done: false });
897	open = true;
898	}
899	if let Some(run) = runs.last_mut() {
900	run.pts.push((here.0, here.1));
901	}
902	}
903	if !hidden {
904	// Wholly in view: one plain closed outline.
905	pen.begin();
906	for (x, y) in &runs[0].pts {
907	let (sx, sy) = f.to_screen(x, y);
908	pen.push(sx, sy);
909	}
910	pen.end(true);
911	return;
912	}
913	if open && runs.len() > 1 {
914	// Ended visible, so the open run wraps onto the first one.
915	let head = runs.remove(0);
916	if let Some(run) = runs.last_mut() {
917	run.pts.extend(head.pts);
918	run.exit = head.exit;
919	}
920	}
921	// A run no bigger than the rings' own tolerance is a stub from simplification, and goes
922	// with the two crossings that bound it. Entries and exits have to alternate round the rim
923	// for the nearest one to be the right one, and a ring that crosses itself breaks that.
924	runs.retain(\|run\| {
925	let (mut lox, mut hix) = (f64::MAX, f64::MIN);
926	let (mut loy, mut hiy) = (f64::MAX, f64::MIN);
927	for (x, y) in &run.pts {
928	lox = lox.min(*x);
929	hix = hix.max(*x);
930	loy = loy.min(*y);
931	hiy = hiy.max(*y);
932	}
933	hix - lox >= slack \|\| hiy - loy >= slack
934	});
935	if runs.is_empty() {
936	return;
937	}
938	// Slack as an angle on the clip circle, for deciding that two crossings are one.
939	let slack_a = slack / f.rho_clip;
940
941	// Stitch: from each run's exit the rim is walked, in the ring's own winding, to the
942	// nearest entry, and loops close where they began.
943	for start in 0..runs.len() {
944	if runs[start].done {
945	continue;
946	}
947	let mut at = start;
948	pen.begin();
949	for _ in 0..=runs.len() {
950	runs[at].done = true;
951	for (x, y) in &runs[at].pts {
952	let (sx, sy) = f.to_screen(x, y);
953	pen.push(sx, sy);
954	}
955	let exit = runs[at].exit;
956	let (mut best, mut near) = (at, f64::MAX);
957	for i in 0..runs.len() {
958	if runs[i].done && i != start {
959	continue;
960	}
961	let raw = if cw { runs[i].entry - exit } else { exit - runs[i].entry };
962	let mut d = raw.rem_euclid(TAU);
963	if d > TAU - slack_a {
964	d -= TAU;
965	}
966	if d < near {
967	near = d;
968	best = i;
969	}
970	}
971	if near > 0.0 {
972	rim(f, exit, if cw { near } else { -near }, step, pen);
973	}
974	if best == start {
975	break;
976	}
977	at = best;
978	}
979	pen.end(true);
980	}
981	}
982
983	fn ortho_line(f: &Frame, ring: &[[f64; 3]], closed: bool, pen: &mut Pen) {
984	let n = ring.len();
985	if n < 2 {
986	return;
987	}
988	let seen: Vec<(f64, f64, f64)> = ring.iter().map(\|p\| f.ortho(p)).collect();
989	let start = match seen.iter().position(\|p\| p.2 < 0.0) {
990	Some(at) => at,
991	None => {
992	// Wholly in view.
993	pen.begin();
994	for p in &seen {
995	let (sx, sy) = f.to_screen(p.0, p.1);
996	pen.push(sx, sy);
997	}
998	pen.end(closed);
999	return;
1000	},
1001	};
1002	// A closed outline is walked from a hidden vertex round to itself, so that no path is
1003	// broken where the ring happens to begin.
1004	let (from, count) = if closed { (start, n + 1) } else { (0, n) };
1005	let mut open = false;
1006	for k in 0..count {
1007	let here = seen[(from + k) % n];
1008	if k > 0 {
1009	let was = seen[(from + k - 1) % n];
1010	if (was.2 >= 0.0) != (here.2 >= 0.0) {
1011	let (ex, ey, _) = crossing(f, was, here);
1012	let (sx, sy) = f.to_screen(ex, ey);
1013	if here.2 >= 0.0 {
1014	pen.begin();
1015	pen.push(sx, sy);
1016	open = true;
1017	} else {
1018	pen.push(sx, sy);
1019	pen.end(false);
1020	open = false;
1021	}
1022	}
1023	}
1024	if here.2 >= 0.0 {
1025	if !open {
1026	pen.begin();
1027	open = true;
1028	}
1029	let (sx, sy) = f.to_screen(here.0, here.1);
1030	pen.push(sx, sy);
1031	}
1032	}
1033	if open {
1034	pen.end(false);
1035	}
1036	}
1037
1038	// ---------------------------------------------------------------------------------------------
1039	// The flat map: unwrap, repeat and clip
1040	// ---------------------------------------------------------------------------------------------
1041
1042	/// A ring on the unturned Mercator plane, in pixels from the centre with `y` north, its
1043	/// longitudes unwrapped so that it runs continuously across the antimeridian.
1044	///
1045	/// A vertex at a pole has no longitude, so it becomes two points on the map's top or bottom
1046	/// edge, under the vertices either side of it.
1047	fn merc_unwrap(f: &Frame, ring: &[[f64; 3]], closed: bool) -> Vec<(f64, f64)> {
1048	let n = ring.len();
1049	// A closed ring is started at a vertex that has a longitude.
1050	let from = if closed {
1051	match ring.iter().position(\|p\| f.merc(p).is_some()) {
1052	Some(at) => at,
1053	None => return Vec::new(),
1054	}
1055	} else {
1056	0
1057	};
1058	let mut out: Vec<(f64, f64)> = Vec::with_capacity(n + 2);
1059	let mut last: Option<(f64, f64)> = None; // unwrapped longitude, raw longitude
1060	let mut pole: Option<f64> = None; // the ordinate of a pole awaiting a longitude
1061	for k in 0..n {
1062	let p = &ring[(from + k) % n];
1063	match f.merc(p) {
1064	None => {
1065	let y = if p[2] > 0.0 { PI } else { -PI };
1066	if let Some((u, _)) = last {
1067	out.push((u, y));
1068	}
1069	pole = Some(y);
1070	},
1071	Some((lam, y)) => {
1072	let u = match last {
1073	Some((u, l)) => u + wrap_pi(lam - l),
1074	None => lam,
1075	};
1076	if let Some(py) = pole.take() {
1077	out.push((u, py));
1078	}
1079	out.push((u, y));
1080	last = Some((u, lam));
1081	},
1082	}
1083	}
1084	if closed {
1085	if let (Some(py), Some((u_last, l_last))) = (pole, last) {
1086	// A closed ring ending on a pole: the pole's second point sits under the first
1087	// vertex, reached from the last.
1088	if let Some(Some((lam0, _))) = ring.get(from).map(\|p\| f.merc(p)) {
1089	out.push((u_last + wrap_pi(lam0 - l_last), py));
1090	}
1091	}
1092	}
1093	for p in out.iter_mut() {
1094	p.0 *= f.r_px;
1095	p.1 = p.1 * f.r_px - f.y_0;
1096	}
1097	out
1098	}
1099
1100	fn merc_ring(f: &Frame, ring: &[[f64; 3]], mode: RingMode, pen: &mut Pen) {
1101	let closed = mode != RingMode::Line;
1102	let mut pts = merc_unwrap(f, ring, closed);
1103	if pts.len() < 2 {
1104	return;
1105	}
1106	let world = TAU * f.r_px;
1107	match mode {
1108	RingMode::Fill => {
1109	// A ring that winds once round the world encircles a pole, and is closed round
1110	// the side of the map where that pole is. Which pole is the one its vertices lean
1111	// towards: no land ring encloses more than a hemisphere.
1112	let first = pts[0];
1113	let last = pts[pts.len() - 1];
1114	let lam_first = first.0 / f.r_px;
1115	let lam_last = last.0 / f.r_px;
1116	let wind = lam_last + wrap_pi(lam_first - lam_last) - lam_first;
1117	if wind.abs() > PI {
1118	let lean: f64 = ring.iter().map(\|p\| p[2]).sum();
1119	let cap = if lean > 0.0 { PI * f.r_px - f.y_0 } else { -PI * f.r_px - f.y_0 };
1120	let u_end = first.0 + wind * f.r_px;
1121	pts.push((u_end, first.1));
1122	pts.push((u_end, cap));
1123	pts.push((first.0, cap));
1124	}
1125	},
1126	RingMode::Outline => pts.push(pts[0]),
1127	RingMode::Line => (),
1128	}
1129	let (mut xl, mut xh, mut yl, mut yh) = (f64::MAX, f64::MIN, f64::MAX, f64::MIN);
1130	for (x, y) in &pts {
1131	xl = xl.min(*x);
1132	xh = xh.max(*x);
1133	yl = yl.min(*y);
1134	yh = yh.max(*y);
1135	}
1136	if yh < f.y_lo \|\| yl > f.y_hi {
1137	return;
1138	}
1139	let k_lo = ((f.x_lo - xh) / world).ceil() as i64;
1140	let k_hi = ((f.x_hi - xl) / world).floor() as i64;
1141	if k_hi < k_lo \|\| k_hi - k_lo > 64 {
1142	return;
1143	}
1144	let (w, h) = (2.0 * f.cx, 2.0 * f.cy);
1145	let rect = (-CLIP_MARGIN_PX, -CLIP_MARGIN_PX, w + CLIP_MARGIN_PX, h + CLIP_MARGIN_PX);
1146	let mut screen: Vec<(f64, f64)> = Vec::with_capacity(pts.len());
1147	for k in k_lo..=k_hi {
1148	let dx = k as f64 * world;
1149	screen.clear();
1150	for (x, y) in &pts {
1151	screen.push(f.to_screen(x + dx, -y));
1152	}
1153	if mode == RingMode::Fill {
1154	let clipped = clip_polygon(&screen, rect);
1155	if clipped.len() >= 3 {
1156	pen.begin();
1157	for (sx, sy) in &clipped {
1158	pen.push(sx, sy);
1159	}
1160	pen.end(true);
1161	}
1162	} else {
1163	clip_polyline(&screen, rect, pen);
1164	}
1165	}
1166	}
1167
1168	/// Sutherland and Hodgman's polygon clip against an axis-aligned rectangle.
1169	///
1170	/// It keeps the winding, so holes stay holes under a nonzero fill, and where a polygon leaves
1171	/// and re-enters it leaves degenerate edges along the rectangle, which fill nothing.
1172	fn clip_polygon(poly: &[(f64, f64)], rect: (f64, f64, f64, f64)) -> Vec<(f64, f64)> {
1173	let (x0, y0, x1, y1) = rect;
1174	let mut cur: Vec<(f64, f64)> = poly.to_vec();
1175	for edge in 0..4 {
1176	if cur.is_empty() {
1177	break;
1178	}
1179	let keep = \|p: &(f64, f64)\| -> bool {
1180	match edge {
1181	0 => p.0 >= x0,
1182	1 => p.0 <= x1,
1183	2 => p.1 >= y0,
1184	_ => p.1 <= y1,
1185	}
1186	};
1187	let cut = \|a: &(f64, f64), b: &(f64, f64)\| -> (f64, f64) {
1188	let t = match edge {
1189	0 => (x0 - a.0) / (b.0 - a.0),
1190	1 => (x1 - a.0) / (b.0 - a.0),
1191	2 => (y0 - a.1) / (b.1 - a.1),
1192	_ => (y1 - a.1) / (b.1 - a.1),
1193	};
1194	(a.0 + t * (b.0 - a.0), a.1 + t * (b.1 - a.1))
1195	};
1196	let mut next = Vec::with_capacity(cur.len() + 4);
1197	let n = cur.len();
1198	for i in 0..n {
1199	let a = &cur[(i + n - 1) % n];
1200	let b = &cur[i];
1201	match (keep(a), keep(b)) {
1202	(true, true) => next.push(*b),
1203	(true, false) => next.push(cut(a, b)),
1204	(false, true) => {
1205	next.push(cut(a, b));
1206	next.push(*b);
1207	},
1208	(false, false) => (),
1209	}
1210	}
1211	cur = next;
1212	}
1213	cur
1214	}
1215
1216	/// Liang and Barsky's segment clip, run along a polyline, cutting it into the pieces that lie
1217	/// within the rectangle.
1218	fn clip_polyline(line: &[(f64, f64)], rect: (f64, f64, f64, f64), pen: &mut Pen) {
1219	let (x0, y0, x1, y1) = rect;
1220	let mut open = false;
1221	for pair in line.windows(2) {
1222	let (a, b) = (pair[0], pair[1]);
1223	let (dx, dy) = (b.0 - a.0, b.1 - a.1);
1224	let (mut t0, mut t1) = (0.0f64, 1.0f64);
1225	let mut gone = false;
1226	for (p, q) in [(-dx, a.0 - x0), (dx, x1 - a.0), (-dy, a.1 - y0), (dy, y1 - a.1)] {
1227	if p == 0.0 {
1228	if q < 0.0 {
1229	gone = true;
1230	break;
1231	}
1232	} else {
1233	let r = q / p;
1234	if p < 0.0 {
1235	t0 = t0.max(r);
1236	} else {
1237	t1 = t1.min(r);
1238	}
1239	}
1240	}
1241	if gone \|\| t0 > t1 {
1242	if open {
1243	pen.end(false);
1244	open = false;
1245	}
1246	continue;
1247	}
1248	if !open \|\| t0 > 0.0 {
1249	if open {
1250	pen.end(false);
1251	}
1252	pen.begin();
1253	pen.push(a.0 + t0 * dx, a.1 + t0 * dy);
1254	open = true;
1255	}
1256	pen.push(a.0 + t1 * dx, a.1 + t1 * dy);
1257	if t1 < 1.0 {
1258	pen.end(false);
1259	open = false;
1260	}
1261	}
1262	if open {
1263	pen.end(false);
1264	}
1265	}
1266
1267	#[cfg(test)]
1268	mod tests {
1269	use super::*;
1270
1271	/// How close a projected coordinate has to be to the oracle's, on a unit sphere.
1272	///
1273	/// The oracle is printed to twelve decimal places, so this is loose enough to absorb its
1274	/// own rounding and tight enough that a wrong coefficient could not pass.
1275	const TOL: f64 = 1.0e-11;
1276
1277	/// Positions and their projections on a unit sphere, Greenwich centred.
1278	///
1279	/// The expected values come from PROJ 9.7.1, an implementation outside this crate:
1280	///
1281	/// ```text
1282	/// proj -f "%.12f" +proj=eqearth +R=1 +lon_0=0
1283	/// ```
1284	///
1285	/// PROJ is the reference implementation the projection's authors' work was folded into,
1286	/// and it is not derived from anything here.
1287	const ORACLE: [(f64, f64, f64, f64); 18] = [
1288	// Perth.
1289	(-31.9535, 115.8571, 1.614621210487, -0.629372441721),
1290	// London.
1291	( 51.5074, -0.1278, -0.001565514410, 0.965105647318),
1292	// New York.
1293	( 40.7128, -74.0060, -0.981861449986, 0.787064220049),
1294	// Tokyo.
1295	( 35.6895, 139.6917, 1.909462941700, 0.697844715778),
1296	// Cape Town.
1297	(-33.9189, 18.4233, 0.254226330215, -0.665601683690),
1298	// Buenos Aires.
1299	(-34.6037, -58.3816, -0.802723755324, -0.678117876992),
1300	// Sydney.
1301	(-33.8688, 151.2093, 2.087106438222, -0.664683776716),
1302	// Honolulu.
1303	( 21.3069, -157.8580, -2.295788127853, 0.426387151159),
1304	// The origin, and the widest and tallest points of the map.
1305	( 0.0, 0.0, 0.0, 0.0),
1306	( 0.0, 180.0, 2.706629983696, 0.0),
1307	( 0.0, -180.0, -2.706629983696, 0.0),
1308	( 90.0, 0.0, 0.0, 1.317362759157),
1309	(-90.0, 0.0, 0.0, -1.317362759157),
1310	// Round numbers in all four quadrants.
1311	( 45.0, 90.0, 1.159854499103, 0.860231085522),
1312	(-45.0, -90.0, -1.159854499103, -0.860231085522),
1313	( 60.0, 30.0, 0.339843347929, 1.088300835505),
1314	(-60.0, 120.0, 1.359373391714, -1.088300835505),
1315	// A small angle, where the polynomial's high terms contribute nothing.
1316	( 1.0, 2.0, 0.030071478499, 0.020257546047),
1317	];
1318
1319	#[test]
1320	fn test_equal_earth_agrees_with_proj_00() -> Outcome<()> {
1321	for (lat, lng, x, y) in ORACLE {
1322	let p = equal_earth_unit(lat, lng);
1323	let near_x = (p.x - x).abs() < TOL;
1324	let near_y = (p.y - y).abs() < TOL;
1325	req!(near_x, true, "{}, {} came out at x = {:.12}, wanted {:.12}.", lat, lng, p.x, x);
1326	req!(near_y, true, "{}, {} came out at y = {:.12}, wanted {:.12}.", lat, lng, p.y, y);
1327	}
1328	Ok(())
1329	}
1330
1331	#[test]
1332	fn test_the_stated_extremes_are_the_extremes_01() -> Outcome<()> {
1333	let east = equal_earth_unit(0.0, 180.0);
1334	let near = (east.x - EQUAL_EARTH_HALF_WIDTH).abs() < TOL;
1335	req!(near, true, "The map is {} wide, not {}.", east.x, EQUAL_EARTH_HALF_WIDTH);
1336	let north = equal_earth_unit(90.0, 0.0);
1337	let near = (north.y - EQUAL_EARTH_HALF_HEIGHT).abs() < TOL;
1338	req!(near, true, "The map is {} tall, not {}.", north.y, EQUAL_EARTH_HALF_HEIGHT);
1339	// No point escapes the box those two describe.
1340	let mut lat = -90.0;
1341	while lat <= 90.0 {
1342	let mut lng = -180.0;
1343	while lng <= 180.0 {
1344	let p = equal_earth_unit(lat, lng);
1345	let in_x = p.x.abs() <= EQUAL_EARTH_HALF_WIDTH + TOL;
1346	let in_y = p.y.abs() <= EQUAL_EARTH_HALF_HEIGHT + TOL;
1347	req!(in_x, true, "x escaped at {}, {}.", lat, lng);
1348	req!(in_y, true, "y escaped at {}, {}.", lat, lng);
1349	lng += 3.0;
1350	}
1351	lat += 3.0;
1352	}
1353	Ok(())
1354	}
1355
1356	#[test]
1357	fn test_the_projection_is_symmetric_02() -> Outcome<()> {
1358	// A pseudocylindrical projection is symmetric about both axes, so the same
1359	// latitude north and south is the same height, and east and west the same width.
1360	for (lat, lng) in [(31.9535, 115.8571), (12.0, 5.0), (78.0, 179.0)] {
1361	let a = equal_earth_unit(lat, lng);
1362	let b = equal_earth_unit(-lat, lng);
1363	let c = equal_earth_unit(lat, -lng);
1364	let same_x = (a.x - b.x).abs() < TOL;
1365	let flip_y = (a.y + b.y).abs() < TOL;
1366	let flip_x = (a.x + c.x).abs() < TOL;
1367	let same_y = (a.y - c.y).abs() < TOL;
1368	req!(same_x, true, "x differed across the equator at {}.", lat);
1369	req!(flip_y, true, "y was not mirrored across the equator at {}.", lat);
1370	req!(flip_x, true, "x was not mirrored across Greenwich at {}.", lng);
1371	req!(same_y, true, "y differed across Greenwich at {}.", lng);
1372	}
1373	// A parallel is straight: every longitude on it has the same y.
1374	let y = equal_earth_unit(20.0, 0.0).y;
1375	for lng in [-170.0, -60.0, 45.0, 179.9] {
1376	let p = equal_earth_unit(20.0, lng);
1377	let level = (p.y - y).abs() < TOL;
1378	req!(level, true, "The 20° parallel bent at {}.", lng);
1379	}
1380	Ok(())
1381	}
1382
1383	#[test]
1384	fn test_a_central_meridian_moves_the_map_and_not_its_shape_03() -> Outcome<()> {
1385	// Recentring on 150° puts Sydney where Greenwich centring puts 1.209° east.
1386	let a = equal_earth(-33.8688, 151.2093, 150.0, 1.0);
1387	let b = equal_earth(-33.8688, 1.2093, 0.0, 1.0);
1388	let same_x = (a.x - b.x).abs() < TOL;
1389	let same_y = (a.y - b.y).abs() < TOL;
1390	req!(same_x, true, "Recentring moved x to {} rather than {}.", a.x, b.x);
1391	req!(same_y, true, "Recentring moved y at all: {} against {}.", a.y, b.y);
1392	// And a longitude that wraps past the antimeridian lands on the far side rather
1393	// than off the map.
1394	let west = equal_earth(0.0, -170.0, 20.0, 1.0);
1395	let wrapped = west.x > 0.0;
1396	req!(wrapped, true, "170° west of a map centred on 20° east came out at {}.", west.x);
1397	Ok(())
1398	}
1399
1400	#[test]
1401	fn test_a_radius_scales_the_map_uniformly_04() -> Outcome<()> {
1402	let r = equal_earth_radius_for_half_width(180.0);
1403	let east = equal_earth(0.0, 180.0, 0.0, r);
1404	let wide = (east.x - 180.0).abs() < 1.0e-9;
1405	req!(wide, true, "The map came out {} wide.", east.x);
1406	let north = equal_earth(90.0, 0.0, 0.0, r);
1407	let want = 180.0 * EQUAL_EARTH_HALF_HEIGHT / EQUAL_EARTH_HALF_WIDTH;
1408	let tall = (north.y - want).abs() < 1.0e-9;
1409	req!(tall, true, "The map came out {} tall.", north.y);
1410	// Scaling is uniform, or the projection would stop being equal-area.
1411	let one = equal_earth_unit(-31.9535, 115.8571);
1412	let big = equal_earth(-31.9535, 115.8571, 0.0, r);
1413	let scaled_x = (big.x - one.x * r).abs() < 1.0e-9;
1414	let scaled_y = (big.y - one.y * r).abs() < 1.0e-9;
1415	req!(scaled_x, true, "x did not scale by the radius.");
1416	req!(scaled_y, true, "y did not scale by the radius.");
1417	Ok(())
1418	}
1419
1420	/// Positions and their orthographic projections on a unit sphere centred on the origin.
1421	///
1422	/// From PROJ 9.7.1, an implementation outside this crate:
1423	///
1424	/// ```text
1425	/// proj -f "%.12f" +proj=ortho +R=1 +lat_0=0 +lon_0=0
1426	/// ```
1427	///
1428	/// Every position here is on the near hemisphere, because PROJ answers `*` for the far
1429	/// one rather than a coordinate; that refusal is what
1430	/// [`test_the_far_side_of_the_globe_is_known_from_the_near_07`] checks against.
1431	const ORTHO_ORACLE: [(f64, f64, f64, f64); 11] = [
1432	// The centre itself.
1433	( 0.0, 0.0, 0.0, 0.0),
1434	// London, Cape Town, Buenos Aires, New York: all within a quarter turn of Greenwich.
1435	( 51.5074, -0.1278, -0.001388311442, 0.782688550807),
1436	(-33.9189, 18.4233, 0.262254675949, -0.558018872482),
1437	(-34.6037, -58.3816, -0.700917639761, -0.567896899696),
1438	( 40.7128, -74.0060, -0.728647317901, 0.652267756393),
1439	// The rim, a tenth of a degree short of it.
1440	( 0.0, 89.9, 0.999998476913, 0.0),
1441	// Both poles, which sit at the top and bottom of the disc from the equator.
1442	( 90.0, 0.0, 0.0, 1.0),
1443	(-90.0, 0.0, 0.0, -1.0),
1444	// Round numbers, and a small angle where the foreshortening is negligible.
1445	( 1.0, 2.0, 0.034894181340, 0.017452406437),
1446	( 60.0, 30.0, 0.250000000000, 0.866025403784),
1447	(-45.0, -90.0, -0.707106781187, -0.707106781187),
1448	];
1449
1450	/// The same, seen from above Perth, so the centre is oblique in both coordinates.
1451	///
1452	/// ```text
1453	/// proj -f "%.12f" +proj=ortho +R=1 +lat_0=-31.9535 +lon_0=115.8571
1454	/// ```
1455	const ORTHO_ORACLE_PERTH: [(f64, f64, f64, f64); 8] = [
1456	// Perth, which is the centre.
1457	(-31.9535, 115.8571, 0.0, 0.0),
1458	// Sydney, Tokyo, Kuala Lumpur, Auckland, Suva.
1459	(-33.8688, 151.2093, 0.480421544470, -0.114447989448),
1460	( 35.6895, 139.6917, 0.328204336296, 0.888173527162),
1461	( 3.1390, 101.6869, -0.244435840733, 0.558819302390),
1462	(-36.8485, 174.7633, 0.685250212424, -0.290118911137),
1463	(-18.1416, 178.4419, 0.843565957319, -0.032624197369),
1464	// Fremantle, sixteen kilometres away and barely off the centre.
1465	(-32.0569, 115.7439, -0.001674457753, -0.001805544881),
1466	// The south pole, which from this latitude sits low on the disc and not on its rim.
1467	(-90.0, 0.0, 0.0, -0.848477887693),
1468	];
1469
1470	#[test]
1471	fn test_orthographic_agrees_with_proj_06() -> Outcome<()> {
1472	for (lat, lng, x, y) in ORTHO_ORACLE {
1473	let p = orthographic(lat, lng, 0.0, 0.0, 1.0);
1474	let near_x = (p.x - x).abs() < TOL;
1475	let near_y = (p.y - y).abs() < TOL;
1476	req!(near_x, true, "{}, {} came out at x = {:.12}, wanted {:.12}.", lat, lng, p.x, x);
1477	req!(near_y, true, "{}, {} came out at y = {:.12}, wanted {:.12}.", lat, lng, p.y, y);
1478	}
1479	for (lat, lng, x, y) in ORTHO_ORACLE_PERTH {
1480	let p = orthographic(lat, lng, -31.9535, 115.8571, 1.0);
1481	let near_x = (p.x - x).abs() < TOL;
1482	let near_y = (p.y - y).abs() < TOL;
1483	req!(near_x, true, "{}, {} came out at x = {:.12}, wanted {:.12}.", lat, lng, p.x, x);
1484	req!(near_y, true, "{}, {} came out at y = {:.12}, wanted {:.12}.", lat, lng, p.y, y);
1485	}
1486	Ok(())
1487	}
1488
1489	#[test]
1490	fn test_the_far_side_of_the_globe_is_known_from_the_near_07() -> Outcome<()> {
1491	// Perth is a quarter turn and more from Greenwich, so PROJ answers `*` for it on a
1492	// Greenwich-centred globe. The cosine is what says so here.
1493	let behind = orthographic_cos_c(-31.9535, 115.8571, 0.0, 0.0);
1494	let far = behind < 0.0;
1495	req!(far, true, "Perth came out in front of Greenwich, at {}.", behind);
1496	// The centre is directly under the viewer, and its antipode directly behind.
1497	let under = orthographic_cos_c(-31.9535, 115.8571, -31.9535, 115.8571);
1498	let full = (under - 1.0).abs() < TOL;
1499	req!(full, true, "The centre came out at {}.", under);
1500	let anti = orthographic_cos_c(31.9535, -64.1429, -31.9535, 115.8571);
1501	let behind_us = (anti + 1.0).abs() < TOL;
1502	req!(behind_us, true, "The antipode came out at {}.", anti);
1503	// On the horizon it is zero, and the point lands exactly on the rim.
1504	let rim = orthographic_cos_c(0.0, 90.0, 0.0, 0.0);
1505	let level = rim.abs() < TOL;
1506	req!(level, true, "A quarter turn away came out at {}.", rim);
1507	let p = orthographic(0.0, 90.0, 0.0, 0.0, 1.0);
1508	let on = (p.x.hypot(p.y) - ORTHOGRAPHIC_HALF_WIDTH).abs() < TOL;
1509	req!(on, true, "The horizon came out {} from the centre.", p.x.hypot(p.y));
1510	// And nothing escapes the disc, near side or far.
1511	let mut lat = -90.0;
1512	while lat <= 90.0 {
1513	let mut lng = -180.0;
1514	while lng <= 180.0 {
1515	let p = orthographic(lat, lng, -31.9535, 115.8571, 1.0);
1516	let inside = p.x.hypot(p.y) <= ORTHOGRAPHIC_HALF_WIDTH + TOL;
1517	req!(inside, true, "{}, {} escaped the disc at {}, {}.", lat, lng, p.x, p.y);
1518	lng += 3.0;
1519	}
1520	lat += 3.0;
1521	}
1522	Ok(())
1523	}
1524
1525	#[test]
1526	fn test_the_globe_turns_without_changing_shape_08() -> Outcome<()> {
1527	// Turning the globe by a degree of longitude and asking for a position a degree
1528	// further east is the same picture, because only the difference matters.
1529	let a = orthographic(12.0, 45.0, 0.0, 30.0, 1.0);
1530	let b = orthographic(12.0, 15.0, 0.0, 0.0, 1.0);
1531	let same_x = (a.x - b.x).abs() < TOL;
1532	let same_y = (a.y - b.y).abs() < TOL;
1533	req!(same_x, true, "Turning moved x to {} rather than {}.", a.x, b.x);
1534	req!(same_y, true, "Turning moved y to {} rather than {}.", a.y, b.y);
1535	// A radius scales the disc and nothing else.
1536	let one = orthographic(-33.8688, 151.2093, -31.9535, 115.8571, 1.0);
1537	let big = orthographic(-33.8688, 151.2093, -31.9535, 115.8571, 320.0);
1538	let scaled_x = (big.x - one.x * 320.0).abs() < 1.0e-9;
1539	let scaled_y = (big.y - one.y * 320.0).abs() < 1.0e-9;
1540	req!(scaled_x, true, "x did not scale by the radius.");
1541	req!(scaled_y, true, "y did not scale by the radius.");
1542	// A longitude past the antimeridian comes round rather than off the globe.
1543	let round = orthographic(0.0, 190.0, 0.0, 170.0, 1.0);
1544	let west = orthographic(0.0, -170.0, 0.0, 170.0, 1.0);
1545	let wrapped = (round.x - west.x).abs() < TOL;
1546	req!(wrapped, true, "190° east did not wrap to 170° west.");
1547	Ok(())
1548	}
1549
1550	#[test]
1551	fn test_a_fix_just_outside_its_range_still_projects_05() -> Outcome<()> {
1552	let over = equal_earth_unit(90.000_000_1, 0.0);
1553	let finite = over.y.is_finite();
1554	req!(finite, true, "A latitude a hair over the pole gave {}.", over.y);
1555	let clamped = (over.y - EQUAL_EARTH_HALF_HEIGHT).abs() < TOL;
1556	req!(clamped, true, "It did not clamp to the pole.");
1557	let round = equal_earth_unit(0.0, 190.0);
1558	let west = equal_earth_unit(0.0, -170.0);
1559	let wrapped = (round.x - west.x).abs() < TOL;
1560	req!(wrapped, true, "190° east did not wrap to 170° west.");
1561	Ok(())
1562	}
1563	/// Positions and their spherical Web Mercator projections on a unit sphere, at four
1564	/// central meridians, from PROJ 9.7.1:
1565	///
1566	/// ```text
1567	/// proj -f "%.12f" +proj=webmerc +R=1 +lon_0=<lon_0>
1568	/// ```
1569	const MERC_POINTS: [(f64, f64); 10] = [
1570	(-31.9535, 115.8571), // Perth
1571	( 51.5074, -0.1278), // London
1572	( 40.7128, -74.0060), // New York
1573	( 35.6895, 139.6917), // Tokyo
1574	( 21.3069, -157.8580), // Honolulu
1575	(-33.9189, 18.4233), // Cape Town
1576	( 0.0, 0.0),
1577	( 85.0, 179.9), // near the top corner
1578	(-85.0, -179.9),
1579	( 60.0, 30.0),
1580	];
1581	const MERC_ORACLE: [(f64, [(f64, f64); 10]); 4] = [
1582	(0.0, [
1583	( 2.022087856812, -0.589076140273), (-0.002230530784, 1.052273175629),
1584	(-1.291648366231, 0.779235626193), ( 2.438080102708, 0.667590039565),
1585	(-2.755141850613, 0.380755537637), ( 0.321547244083, -0.629951578195),
1586	( 0.000000000000, 0.000000000000), ( 3.139847324338, 3.131301331472),
1587	(-3.139847324338, -3.131301331472), ( 0.523598775598, 1.316957896925),
1588	]),
1589	(115.8571, [
1590	( 0.000000000000, -0.589076140273), (-2.024318387596, 1.052273175629),
1591	( 2.969449084136, 0.779235626193), ( 0.415992245896, 0.667590039565),
1592	( 1.505955599754, 0.380755537637), (-1.700540612730, -0.629951578195),
1593	(-2.022087856812, 0.000000000000), ( 1.117759467525, 3.131301331472),
1594	( 1.121250126029, -3.131301331472), (-1.498489081214, 1.316957896925),
1595	]),
1596	(-74.006, [
1597	(-2.969449084136, -0.589076140273), ( 1.289417835447, 1.052273175629),
1598	( 0.000000000000, 0.779235626193), (-2.553456838240, 0.667590039565),
1599	(-1.463493484382, 0.380755537637), ( 1.613195610314, -0.629951578195),
1600	( 1.291648366231, 0.000000000000), (-1.851689616611, 3.131301331472),
1601	(-1.848198958107, -3.131301331472), ( 1.815247141829, 1.316957896925),
1602	]),
1603	(-150.0, [
1604	(-1.643103572376, -0.589076140273), ( 2.615763347207, 1.052273175629),
1605	( 1.326345511761, 0.779235626193), (-1.227111326480, 0.667590039565),
1606	(-0.137147972622, 0.380755537637), ( 2.939541122074, -0.629951578195),
1607	( 2.617993877991, 0.000000000000), (-0.525344104850, 3.131301331472),
1608	(-0.521853446346, -3.131301331472), ( 3.141592653590, 1.316957896925),
1609	]),
1610	];
1611
1612	#[test]
1613	fn test_web_mercator_agrees_with_proj_09() -> Outcome<()> {
1614	for (lon_0, want) in MERC_ORACLE {
1615	for ((lat, lng), (x, y)) in MERC_POINTS.iter().zip(want.iter()) {
1616	let p = web_mercator(lat, lng, lon_0, 1.0);
1617	let near = (p.x - x).abs() < TOL && (p.y - y).abs() < TOL;
1618	req!(near, true, "{}, {} about {} came out at ({:.12}, {:.12}), wanted ({}, {}).",
1619	lat, lng, lon_0, p.x, p.y, x, y);
1620	}
1621	}
1622	// The map is square at its stated limit, and a pole clamps to it.
1623	let top = web_mercator(90.0, 0.0, 0.0, 1.0);
1624	let square = (top.y - PI).abs() < 1.0e-12;
1625	req!(square, true, "The map's top came out at {}, not pi.", top.y);
1626	Ok(())
1627	}
1628
1629	#[test]
1630	fn test_web_mercator_inverse_agrees_with_proj_10() -> Outcome<()> {
1631	// proj -I -f "%.12f" +proj=webmerc +R=1 +lon_0=115.8571, which answers lng, lat.
1632	for (x, y, lng, lat) in [
1633	( 0.3, -0.6, 133.045833853925, -32.483015050012),
1634	(-1.2, 1.1, 47.102164584301, 53.177781879084),
1635	( 2.5, -2.9, -100.903451217294, -83.701155006501),
1636	( 0.0, 0.0, 115.857100000000, 0.000000000000),
1637	] {
1638	let (a, b) = web_mercator_inverse(Pt::new(x, y), 115.8571, 1.0);
1639	let near = (a - lat).abs() < 1.0e-9 && (b - lng).abs() < 1.0e-9;
1640	req!(near, true, "({}, {}) came back as {}, {}; PROJ says {}, {}.", x, y, a, b, lat, lng);
1641	}
1642	// Forward and back is the identity, including across the antimeridian.
1643	for (lat, lng) in MERC_POINTS {
1644	let p = web_mercator(lat, lng, -150.0, 6.0);
1645	let (a, b) = web_mercator_inverse(p, -150.0, 6.0);
1646	let same = (a - lat).abs() < 1.0e-9 && (wrap_half_turn(b - lng)).abs() < 1.0e-9;
1647	req!(same, true, "{}, {} came back as {}, {}.", lat, lng, a, b);
1648	}
1649	Ok(())
1650	}
1651
1652	#[test]
1653	fn test_orthographic_inverse_agrees_with_proj_11() -> Outcome<()> {
1654	// proj -I -f "%.12f" +proj=ortho +R=1 +lat_0=<lat_0> +lon_0=<lon_0>, answering lng,
1655	// lat, at a centre on the equator, two oblique centres, a third far north, the north
1656	// pole, and a southern one. The last input of each is off the disc, where PROJ
1657	// answers `*`.
1658	let xy = [(0.0, 0.0), (0.3, -0.2), (0.7, 0.5), (-0.9, 0.1), (0.0, 0.99), (-0.45, -0.8),
1659	(0.001, 0.002)];
1660	let oracle: [((f64, f64), [(f64, f64); 7]); 6] = [
1661	((0.0, 0.0), [
1662	(0.000000000000, 0.000000000000), (17.829543848069, -11.536959032815),
1663	(53.929231349204, 30.000000000000), (-64.760598179321, 5.739170477267),
1664	(0.000000000000, 81.890385544006), (-48.590377890729, -53.130102354156),
1665	(0.057295903654, 0.114591635421),
1666	]),
1667	((-31.9535, 115.8571), [
1668	(115.857100000000, -31.953500000000), (139.491172018755, -41.554276226097),
1669	(160.969622305304, 8.881020396659), (50.501798095373, -8.029667432931),
1670	(115.857100000000, 49.936885544006), (14.957210228986, -62.724626645594),
1671	(115.924543725221, -31.838890517892),
1672	]),
1673	((51.5074, -0.1278), [
1674	(-0.127800000000, 51.507400000000), (22.018911603216, 37.269190760395),
1675	(95.904696130209, 45.259426785352), (-78.463343669197, 23.222821817975),
1676	(179.872200000000, 46.602214455994), (-27.392900860832, -10.795895994418),
1677	(-0.035513551106, 51.621955519567),
1678	]),
1679	((35.6895, 139.6917), [
1680	(139.691700000000, 35.689500000000), (158.631638139189, 22.439898158901),
1681	(-140.229884106383, 44.713981955783), (67.334907163601, 19.191798955953),
1682	(-40.308300000000, 62.420114455994), (109.995075616966, -24.722618420153),
1683	(139.762346389474, 35.804071028094),
1684	]),
1685	((90.0, 0.0), [
1686	(0.000000000000, 90.000000000000), (56.309932474020, 68.865707785214),
1687	(125.537677791974, 30.657298992941), (-96.340191745910, 25.104090250221),
1688	(180.000000000000, 8.109614455994), (-29.357753542791, 23.382196425509),
1689	(153.434948822922, 89.871882635420),
1690	]),
1691	((-60.0, -70.0), [
1692	(-70.000000000000, -60.000000000000), (-24.339703634894, -65.199623342309),
1693	(-24.503147524221, -11.045475113038), (-141.637585260840, -18.507178106808),
1694	(-70.000000000000, 21.890385544006), (152.308915595694, -48.046975207444),
1695	(-69.885803894134, -59.885358916100),
1696	]),
1697	];
1698	for ((lat_0, lon_0), want) in oracle {
1699	for ((x, y), (lng, lat)) in xy.iter().zip(want.iter()) {
1700	let got = match orthographic_inverse(Pt::new(x, y), lat_0, lon_0, 1.0) {
1701	Some(g) => g,
1702	None => return Err(err!("({}, {}) about {}, {} fell off the globe.",
1703	x, y, lat_0, lon_0; Test)),
1704	};
1705	// At the pole longitude is arbitrary; everywhere else it is compared too.
1706	let near_lat = (got.0 - lat).abs() < 1.0e-9;
1707	let near_lng = lat.abs() > 89.999_999 \|\| wrap_half_turn(got.1 - lng).abs() < 1.0e-9;
1708	req!(near_lat && near_lng, true,
1709	"({}, {}) about {}, {} came back as {}, {}; PROJ says {}, {}.",
1710	x, y, lat_0, lon_0, got.0, got.1, lat, lng);
1711	}
1712	let off = orthographic_inverse(Pt::new(0.8, 0.7), lat_0, lon_0, 1.0);
1713	req!(off.is_none(), true, "A point off the disc came back as {:?}.", off);
1714	}
1715	Ok(())
1716	}
1717
1718	}