Fix 4326 streaks in regions reaching a pole on Mercator

src/renderable-region.test.ts

@@ -31,6 +31,7 @@ const program = new Proxy(
     clipCountLoc: { name: 'u_clip_count' },
     latLookupLoc: { name: 'u_latLookup' },
     latAnchorLoc: { name: 'u_latAnchor' },
+    latShiftLoc: { name: 'u_latShift' },
   } as Record<string, unknown>,
   { get: (target, key) => (key in target ? target[key as string] : null) }
 ) as unknown as ShaderProgram
@@ -190,9 +191,10 @@ describe('latLookupAnchor', () => {
   }
   function shaderV(
     clipY: number,
-    [s0, c0, vAnchor, vPerRad]: [number, number, number, number]
+    [s0, c0, vAnchor, vPerRad]: [number, number, number, number],
+    latShift = 0
   ): number {
-    const dy = f(-2 * f(Math.PI) * f(clipY))
+    const dy = f(-2 * f(Math.PI) * f(f(clipY) - f(latShift)))
     const sh = sinhSmall(f(0.5 * dy))
     const dTan = f(f(2 * f(s0) * sh * sh) + f(f(c0) * sinhSmall(dy)))
     const dLat = f(Math.atan2(dTan, f(1 + f(f(s0) * f(f(s0) + dTan)))))
@@ -231,6 +233,33 @@ describe('latLookupAnchor', () => {
     expect(v).toBeCloseTo((72.6 - 46.7) / 50.7, 5)
   })
 
+  it('maps rows of a region whose mesh reaches past the Mercator world', () => {
+    // The top row of a global pyramid: its mesh runs nearly to the pole.
+    const edges: [number, number] = [90, 67.5]
+    const meshBounds = { x0: 0.5, y0: -2.45342, x1: 0.625, y1: 0.24298 }
+    const gl = createRecordingGl()
+    renderRegion(
+      gl,
+      program,
+      { ...region(null), meshBounds, texLatEdges: edges },
+      [0],
+      config
+    )
+    const a = gl
+      .callsTo('uniform4f')
+      .find((c) => c.args[0] === 'u_latAnchor')!
+      .args.slice(1) as [number, number, number, number]
+    const latShift = Number(
+      gl.callsTo('uniform1f').find((c) => c.args[0] === 'u_latShift')
+        ?.args[1] ?? 0
+    )
+    const meshY = (meshBounds.y0 + meshBounds.y1) / 2
+    for (let lat = 67.5; lat <= 85; lat += 0.5) {
+      const v = shaderV(latToMercatorNorm(lat) - meshY, a, latShift)
+      expect(Math.abs(v - (90 - lat) / 22.5) * 128).toBeLessThan(1e-3)
+    }
+  })
+
   it.each([
     ['crossing the equator', -30, 30, 720],
     ['reaching the Mercator limit', 60, 85, 300],

src/renderable-region.ts

@@ -115,8 +115,15 @@ export function renderRegion(
     gl.uniform1i(shaderProgram.latLookupLoc, texLatEdges ? 1 : 0)
   }
   if (texLatEdges && shaderProgram.latAnchorLoc !== null) {
-    const [s0, c0, vAnchor, vPerRad] = latLookupAnchor(shiftY, texLatEdges)
+    // A mesh reaching a pole is anchored far past the Mercator world, where
+    // the lookup's Float32 terms cancel, so the lookup anchors to the part of
+    // the region the map can show.
+    const latAnchorY = (clampUnit(meshBounds.y0) + clampUnit(meshBounds.y1)) / 2
+    const [s0, c0, vAnchor, vPerRad] = latLookupAnchor(latAnchorY, texLatEdges)
     gl.uniform4f(shaderProgram.latAnchorLoc, s0, c0, vAnchor, vPerRad)
+    if (shaderProgram.latShiftLoc !== null) {
+      gl.uniform1f(shaderProgram.latShiftLoc, latAnchorY - shiftY)
+    }
   }
 
   if (shaderProgram.clipCountLoc !== null) {
@@ -209,3 +216,5 @@ export function latLookupAnchor(
     180 / Math.PI / span,
   ]
 }
+
+const clampUnit = (y: number) => Math.min(1, Math.max(0, y))

src/shader-program.test.ts

@@ -325,6 +325,7 @@ describe('shader program uniform surface', () => {
       const program: ShaderProgram = shaderProgram
       expect(program.latLookupLoc).not.toBeNull()
       expect(program.latAnchorLoc).not.toBeNull()
+      expect(program.latShiftLoc).not.toBeNull()
     }
   })
 })

src/shader-program.ts

@@ -64,6 +64,7 @@ export interface ShaderProgram {
   // EPSG:4326 reprojection uniforms
   latLookupLoc: WebGLUniformLocation | null
   latAnchorLoc: WebGLUniformLocation | null
+  latShiftLoc: WebGLUniformLocation | null
   // Finer-data clip boxes, see FRAGMENT_CLIP
   clipBoxesLoc: WebGLUniformLocation | null
   clipCountLoc: WebGLUniformLocation | null
@@ -320,6 +321,7 @@ export function createShaderProgram(
     dataScaleLoc: gl.getUniformLocation(program, 'u_dataScale'),
     latLookupLoc: gl.getUniformLocation(program, 'u_latLookup'),
     latAnchorLoc: gl.getUniformLocation(program, 'u_latAnchor'),
+    latShiftLoc: gl.getUniformLocation(program, 'u_latShift'),
     clipBoxesLoc: gl.getUniformLocation(program, 'u_clip_boxes'),
     clipCountLoc: gl.getUniformLocation(program, 'u_clip_count'),
   }

src/shaders.test.ts

@@ -302,7 +302,7 @@ describe('mercator math parity (GPU shader <-> CPU map-utils)', () => {
     expect(vtx).toContain('float t = PI * (1.0 - 2.0 * y);')
     expect(vtx).toContain('atan(sinh(t))')
     for (const line of [
-      'float dy = -2.0 * PI * v_clipPos.y;',
+      'float dy = -2.0 * PI * (v_clipPos.y - u_latShift);',
       'float sh = zl_sinh(0.5 * dy);',
       'float dTan = 2.0 * u_latAnchor.x * sh * sh + u_latAnchor.y * zl_sinh(dy);',
       'float dLat = atan(dTan, 1.0 + u_latAnchor.x * (u_latAnchor.x + dTan));',

src/shaders.ts

@@ -404,15 +404,17 @@ const FRAG_CONST_PI = `const float PI = 3.14159265358979323846;`
  * computes v from each fragment's own latitude instead.
  *
  * v_clipPos.y is the fragment's normalized Mercator offset from the region's
- * anchor, so the latitude offset from the anchor is formed without touching
- * an absolute Float32 coordinate, which would lose rows at high zoom:
+ * anchor, and u_latShift moves that anchor to φ0 inside the Mercator world, so
+ * the latitude offset from φ0 is formed without touching an absolute Float32
+ * coordinate, which would lose rows at high zoom:
  *   tan(φ − φ0) = (tan φ − tan φ0) / (1 + tan φ tan φ0),  tan φ = sinh(y0 + dy)
  * with tan φ − tan φ0 = 2 tan φ0 sinh²(dy/2) + sec φ0 sinh(dy).
  * u_latAnchor = (tan φ0, sec φ0, v at φ0, dv/dφ), from latLookupAnchor.
  */
 const UNIFORMS_LAT_LOOKUP = `
 uniform int u_latLookup;
-uniform vec4 u_latAnchor;`
+uniform vec4 u_latAnchor;
+uniform float u_latShift;`
 
 // Builtin sinh is often exp-based and loses all relative precision near zero.
 const FUNC_SINH_SMALL = `
@@ -427,7 +429,7 @@ float zl_sinh(float x) {
 const FRAGMENT_SAMPLE_COORD = `
   vec2 sample_coord = pix_coord;
   if (u_latLookup == 1) {
-    float dy = -2.0 * PI * v_clipPos.y;
+    float dy = -2.0 * PI * (v_clipPos.y - u_latShift);
     float sh = zl_sinh(0.5 * dy);
     float dTan = 2.0 * u_latAnchor.x * sh * sh + u_latAnchor.y * zl_sinh(dy);
     float dLat = atan(dTan, 1.0 + u_latAnchor.x * (u_latAnchor.x + dTan));