Address terrain globe Mercator review comments

Author: charlieforward9

modules/geo-layers/src/terrain-layer/terrain-layer.ts

@@ -18,6 +18,11 @@ import {SimpleMeshLayer} from '@deck.gl/mesh-layers';
 import {COORDINATE_SYSTEM} from '@deck.gl/core';
 import type {Mesh} from '@loaders.gl/schema';
 import {TerrainWorkerLoader} from '@loaders.gl/terrain';
+import {
+  MAX_LATITUDE as MAX_WEB_MERCATOR_LATITUDE,
+  lngLatToWorld,
+  worldToLngLat
+} from '@math.gl/web-mercator';
 import TileLayer, {TileLayerProps} from '../tile-layer/tile-layer';
 import type {
   Bounds,
@@ -33,8 +38,6 @@ const TILE_OVERLAP_PIXELS = 1;
 const MIN_TERRAIN_MESH_MAX_ERROR = 1;
 const MAX_LATITUDE = 90;
 const MAX_LONGITUDE = 180;
-const DEGREES_TO_RADIANS = Math.PI / 180;
-const RADIANS_TO_DEGREES = 180 / Math.PI;
 
 const defaultProps: DefaultProps<TerrainLayerProps> = {
   ...TileLayer.defaultProps,
@@ -110,6 +113,7 @@ type TerrainLoadProps = {
   elevationData: string | null;
   elevationDecoder: ElevationDecoder;
   meshMaxError: number;
+  remapToWebMercatorTile?: boolean;
   signal?: AbortSignal;
 };
 
@@ -205,6 +209,7 @@ export default class TerrainLayer<ExtraPropsT extends {} = {}> extends Composite
     bounds,
     elevationDecoder,
     meshMaxError,
+    remapToWebMercatorTile,
     signal
   }: TerrainLoadProps): Promise<Mesh> | null {
     if (!elevationData) {
@@ -223,7 +228,16 @@ export default class TerrainLayer<ExtraPropsT extends {} = {}> extends Composite
       }
     };
     const {fetch} = this.props;
-    return fetch(elevationData, {propName: 'elevationData', layer: this, loadOptions, signal});
+    const terrain = fetch(elevationData, {
+      propName: 'elevationData',
+      layer: this,
+      loadOptions,
+      signal
+    });
+
+    return remapToWebMercatorTile
+      ? terrain.then(mesh => (mesh ? remapMeshToWebMercatorTile(mesh, bounds) : mesh))
+      : terrain;
   }
 
   getTiledTerrainData(tile: TileLoadProps): Promise<MeshAndTexture> {
@@ -245,22 +259,18 @@ export default class TerrainLayer<ExtraPropsT extends {} = {}> extends Composite
       topRight = [bbox.right, bbox.top];
     }
     const bounds: Bounds = [bottomLeft[0], bottomLeft[1], topRight[0], topRight[1]];
-    const overlappedBounds = getOverlappedBounds(
-      bounds,
-      this.props.tileSize,
-      Boolean(viewport.resolution && viewport.resolution > 0)
-    );
+    const isGlobe = Boolean(viewport.resolution && viewport.resolution > 0);
+    const overlappedBounds = getOverlappedBounds(bounds, this.props.tileSize, isGlobe);
 
     const terrain =
       this.loadTerrain({
         elevationData: dataUrl,
         bounds: overlappedBounds,
         elevationDecoder,
         meshMaxError,
+        remapToWebMercatorTile: isGlobe,
         signal
-      })?.then(mesh =>
-        viewport.resolution && mesh ? remapMeshToWebMercatorTile(mesh, overlappedBounds) : mesh
-      ) ?? Promise.resolve(null);
+      }) ?? Promise.resolve(null);
     const surface = textureUrl
       ? // If surface image fails to load, the tile should still be displayed
         fetch(textureUrl, {propName: 'texture', layer: this, loaders: [], signal}).catch(_ => null)
@@ -431,14 +441,14 @@ function remapMeshToWebMercatorTile(mesh: Mesh, bounds: Bounds): Mesh {
   }
 
   const [, south, , north] = bounds;
-  const northY = lngLatToMercatorY(north);
-  const southY = lngLatToMercatorY(south);
+  const northY = lngLatToMercatorWorldY(north);
+  const southY = lngLatToMercatorWorldY(south);
   const remappedPositions = new Float32Array(positions);
 
   for (let i = 0; i < texCoords.length / 2; i++) {
     const v = texCoords[i * 2 + 1];
     const mercatorY = northY + (southY - northY) * v;
-    remappedPositions[i * 3 + 1] = mercatorYToLat(mercatorY);
+    remappedPositions[i * 3 + 1] = worldToLngLat([0, mercatorY])[1];
   }
 
   return {
@@ -453,12 +463,10 @@ function remapMeshToWebMercatorTile(mesh: Mesh, bounds: Bounds): Mesh {
   };
 }
 
-function lngLatToMercatorY(latitude: number): number {
-  const clampedLatitude = Math.max(-85.051129, Math.min(85.051129, latitude));
-  const sin = Math.sin(clampedLatitude * DEGREES_TO_RADIANS);
-  return 0.5 - Math.log((1 + sin) / (1 - sin)) / (4 * Math.PI);
-}
-
-function mercatorYToLat(y: number): number {
-  return Math.atan(Math.sinh(Math.PI * (1 - 2 * y))) * RADIANS_TO_DEGREES;
+function lngLatToMercatorWorldY(latitude: number): number {
+  const clampedLatitude = Math.max(
+    -MAX_WEB_MERCATOR_LATITUDE,
+    Math.min(MAX_WEB_MERCATOR_LATITUDE, latitude)
+  );
+  return lngLatToWorld([0, clampedLatitude])[1];
 }

test/modules/geo-layers/terrain-layer.spec.ts

@@ -8,6 +8,11 @@ import {TerrainLayer, TileLayer} from '@deck.gl/geo-layers';
 import {_GlobeView as GlobeView} from '@deck.gl/core';
 import {SimpleMeshLayer} from '@deck.gl/mesh-layers';
 import {TerrainLoader} from '@loaders.gl/terrain';
+import {
+  MAX_LATITUDE as MAX_WEB_MERCATOR_LATITUDE,
+  lngLatToWorld,
+  worldToLngLat
+} from '@math.gl/web-mercator';
 
 test('TerrainLayer', async () => {
   const testCases = generateLayerTests({
@@ -50,15 +55,27 @@ test('TerrainLayer', async () => {
 });
 
 test('TerrainLayer#globe remaps WebMercator tile rows to lng/lat mesh positions', async () => {
+  const sourcePositions = new Float32Array([0, 80, 0, 0.5, 40, 0, 1, 0, 0]);
+  const tileSize = 512;
+  const bbox = {west: 0, south: 0, east: 1, north: 80};
+  const yPad = ((bbox.north - bbox.south) / tileSize) * 1;
+  const overlappedSouth = bbox.south - yPad;
+  const overlappedNorth = bbox.north + yPad;
+  const expectedMiddleLatitude = worldToLngLat([
+    0,
+    (lngLatToMercatorWorldY(overlappedNorth) + lngLatToMercatorWorldY(overlappedSouth)) / 2
+  ])[1];
+
   const sourceMesh = {
     attributes: {
-      POSITION: {value: new Float32Array([0, 80, 0, 0.5, 40, 0, 1, 0, 0]), size: 3},
+      POSITION: {value: sourcePositions, size: 3},
       TEXCOORD_0: {value: new Float32Array([0, 0, 0.5, 0.5, 1, 1]), size: 2}
     }
   };
   const layer = new TerrainLayer({
     id: 'terrain-globe-mercator',
     elevationData: 'terrain/{z}/{x}/{y}.png',
+    tileSize,
     fetch: () => Promise.resolve(sourceMesh)
   });
   layer.context = {
@@ -77,15 +94,34 @@ test('TerrainLayer#globe remaps WebMercator tile rows to lng/lat mesh positions'
   const [mesh] = await layer.getTiledTerrainData({
     index: {x: 0, y: 0, z: 1},
     id: '0-0-1',
-    bbox: {west: 0, south: 0, east: 1, north: 80},
+    bbox,
     zoom: 1
   });
-  const positions = mesh.attributes.POSITION.value;
+  const positions = mesh!.attributes.POSITION.value;
+
+  expect(positions, 'remap copies the loader positions').not.toBe(sourcePositions);
+  expect(sourcePositions[1], 'source top row is unchanged').toBe(80);
+  expect(sourcePositions[4], 'source middle row is unchanged').toBe(40);
+  expect(sourcePositions[7], 'source bottom row is unchanged').toBe(0);
 
-  expect(positions[1], 'top row latitude is preserved').toBeGreaterThan(80);
-  expect(
-    positions[4],
-    'middle row uses Mercator latitude instead of linear latitude'
-  ).toBeGreaterThan(40);
-  expect(positions[7], 'bottom row latitude is preserved').toBeLessThan(0);
+  expect(positions[1], 'top row latitude follows the overlapped tile north').toBeCloseTo(
+    overlappedNorth,
+    5
+  );
+  expect(positions[4], 'middle row uses Mercator latitude instead of linear latitude').toBeCloseTo(
+    expectedMiddleLatitude,
+    5
+  );
+  expect(positions[7], 'bottom row latitude follows the overlapped tile south').toBeCloseTo(
+    overlappedSouth,
+    5
+  );
 });
+
+function lngLatToMercatorWorldY(latitude: number): number {
+  const clampedLatitude = Math.max(
+    -MAX_WEB_MERCATOR_LATITUDE,
+    Math.min(MAX_WEB_MERCATOR_LATITUDE, latitude)
+  );
+  return lngLatToWorld([0, clampedLatitude])[1];
+}