feat(surface): Generate dense mesh surface (#997)

### Change list

- Currently, as implemented in
#976 the SurfaceLayer
only supports four vertices and two internal triangles. Over larger
areas this will run the risk of larger reprojection non-linearities.
- This PR implements a quick initial hack to this problem: offering the
user the ability to specify the height and width of the regular grid
they'd like to use.

```py
import rasterio
from sidecar import Sidecar

from lonboard import Map
from lonboard.experimental._surface import SurfaceLayer

# aws s3 cp s3://naip-visualization/ny/2022/60cm/rgb/40073/m_4007307_sw_18_060_20220803.tif ./ --request-payer
src = rasterio.open("m_4007307_sw_18_060_20220803.tif")

layer = SurfaceLayer.from_rasterio(
    src, 
    downscale=4, 
    opacity=0.8, 
    wireframe=True
)
sidecar = Sidecar(anchor="split-right")
m = Map(layer)
with sidecar:
    display(m)
```

<img width="588" height="618" alt="image"
src="https://github.com/user-attachments/assets/45943649-4669-4bef-8879-eacc4598a054"
/>
<img width="468" height="581" alt="image"
src="https://github.com/user-attachments/assets/e5dd64e4-6c59-4592-a0df-07ff90fc8f96"
/>

Author: kylebarron

lonboard/experimental/_surface.py

@@ -50,6 +50,71 @@ def load_arr_and_transform(
     return arr, overview_transform
 
 
+def apply_colormap(
+    arr: NDArray[np.uint8],
+    cmap: dict[int, tuple[int, int, int] | tuple[int, int, int, int]],
+) -> NDArray[np.uint8]:
+    """Apply rasterio colormap to single-band array."""
+    lut = np.zeros((max(cmap.keys()), 4), dtype=np.uint8)
+    for k, v in cmap.items():
+        lut[k] = v
+
+    return lut[arr]
+
+
+def generate_mesh_grid(
+    *,
+    n_rows: int = 50,
+    n_cols: int = 50,
+) -> tuple[NDArray[np.float32], NDArray[np.uint32]]:
+    """Generate a regular grid mesh with the given number of rows and columns.
+
+    Creates a mesh covering the unit square [0, 1] x [0, 1] divided into
+    n_rows x n_cols cells. Each cell is split into two triangles along the diagonal.
+
+    Args:
+        n_rows: Number of rows in the grid
+        n_cols: Number of columns in the grid
+
+    Returns:
+        positions: Array of shape ((n_rows+1) * (n_cols+1), 2) containing vertex positions
+                   in normalized image coordinates [0, 1]
+        triangles: Array of shape (n_rows * n_cols * 2, 3) containing triangle vertex indices
+
+    """
+    # Generate vertex grid
+    # We need (n_rows + 1) x (n_cols + 1) vertices to create n_rows x n_cols cells
+    x = np.linspace(0, 1, n_cols + 1, dtype=np.float32)
+    y = np.linspace(0, 1, n_rows + 1, dtype=np.float32)
+
+    # Create meshgrid and flatten to get all vertex positions
+    xx, yy = np.meshgrid(x, y)
+    positions = np.stack([xx.ravel(), yy.ravel()], axis=-1)
+
+    # Generate triangle indices
+    # For each cell (i, j), we create two triangles:
+    # Triangle 1: bottom-left, bottom-right, top-left
+    # Triangle 2: bottom-right, top-right, top-left
+    triangles = np.empty((n_rows * n_cols * 2, 3), dtype=np.uint32)
+
+    i = 0
+    for row in range(n_rows):
+        for col in range(n_cols):
+            # Vertex indices for the current cell
+            bottom_left = row * (n_cols + 1) + col
+            bottom_right = bottom_left + 1
+            top_left = (row + 1) * (n_cols + 1) + col
+            top_right = top_left + 1
+
+            triangles[i] = [bottom_left, bottom_right, top_left]
+            triangles[i + 1] = [bottom_right, top_right, top_left]
+            i += 2
+
+    triangles_array = np.array(triangles, dtype=np.uint32)
+
+    return positions, triangles_array
+
+
 def rescale_positions_to_image_crs(
     positions: NDArray[np.float32],
     *,
@@ -76,30 +141,33 @@ def from_rasterio(
         src: DatasetReader,
         *,
         downscale: int | None = None,
+        mesh_n_rows: int = 50,
+        mesh_n_cols: int = 50,
         **kwargs: Any,
     ) -> Self:
         import rasterio.plot
 
         arr, transform = load_arr_and_transform(src, downscale=downscale)
 
-        # swap axes order from (bands, rows, columns) to (rows, columns, bands)
-        image_arr: np.ndarray = rasterio.plot.reshape_as_image(arr)
+        if arr.shape[0] == 1 and src.colormap(1) is not None:
+            image_arr = apply_colormap(arr[0], src.colormap(1))
+        else:
+            # swap axes order from (bands, rows, columns) to (rows, columns, bands)
+            image_arr = rasterio.plot.reshape_as_image(arr)
 
         image_height = image_arr.shape[0]
         image_width = image_arr.shape[1]
 
-        positions = np.array(
-            [
-                [0, 0],  # bottom-left
-                [1, 0],  # bottom-right
-                [0, 1],  # top-left
-                [1, 1],  # top-right
-            ],
-            dtype=np.float32,
+        # Generate mesh grid in normalized [0, 1] coordinates
+        # These positions serve as texture coordinates
+        tex_coords, triangles = generate_mesh_grid(
+            n_rows=mesh_n_rows,
+            n_cols=mesh_n_cols,
         )
 
+        # Reproject mesh vertices from image CRS to EPSG:4326
         source_crs_coords = rescale_positions_to_image_crs(
-            positions,
+            tex_coords,
             image_height=image_height,
             image_width=image_width,
             transform=transform,
@@ -113,29 +181,13 @@ def from_rasterio(
         )
         lonlat_coords = np.stack([lons, lats], axis=-1)
 
+        # Add z-coordinate (elevation), setting to zero as we care about a flat mesh on
+        # the map surface
         final_positions = np.concatenate(
             [lonlat_coords, np.zeros((lonlat_coords.shape[0], 1), dtype=np.float32)],
             axis=-1,
         )
 
-        triangles = np.array(
-            [
-                [0, 1, 2],
-                [1, 2, 3],
-            ],
-            dtype=np.uint32,
-        )
-
-        tex_coords = np.array(
-            [
-                [0, 0],  # bottom-left
-                [1, 0],  # bottom-right
-                [0, 1],  # top-left
-                [1, 1],  # top-right
-            ],
-            dtype=np.float32,
-        )
-
         return cls(
             positions=final_positions,
             triangles=triangles,

tests/layers/surface/__init__.py

@@ -0,0 +1,87 @@
+"""Test mesh generation for SurfaceLayer."""
+
+import numpy as np
+
+from lonboard.experimental._surface import generate_mesh_grid
+
+
+def test_generate_mesh_grid_basic():
+    """Test basic mesh generation with small grid."""
+    n_rows, n_cols = 2, 2
+    positions, triangles = generate_mesh_grid(n_rows=n_rows, n_cols=n_cols)
+
+    # Check positions shape: (n_rows+1) * (n_cols+1) vertices
+    expected_num_vertices = (n_rows + 1) * (n_cols + 1)  # 3 * 3 = 9
+    assert positions.shape == (expected_num_vertices, 2)
+
+    # Check triangles shape: n_rows * n_cols * 2 triangles
+    expected_num_triangles = n_rows * n_cols * 2  # 2 * 2 * 2 = 8
+    assert triangles.shape == (expected_num_triangles, 3)
+
+    # Check all triangle indices are valid (within vertex range)
+    assert np.all(triangles >= 0)
+    assert np.all(triangles < expected_num_vertices)
+
+    # Check positions are in [0, 1] range
+    assert np.all(positions >= 0)
+    assert np.all(positions <= 1)
+
+
+def test_generate_mesh_grid_single_cell():
+    """Test with a single cell (simplest case)."""
+    positions, triangles = generate_mesh_grid(n_rows=1, n_cols=1)
+
+    # Should have 4 vertices
+    assert positions.shape == (4, 2)
+
+    # Should have 2 triangles
+    assert triangles.shape == (2, 3)
+
+    # Verify exact positions for single cell
+    expected_positions = np.array([[0, 0], [1, 0], [0, 1], [1, 1]], dtype=np.float32)
+    np.testing.assert_allclose(positions, expected_positions)
+
+    # Verify triangle indices
+    assert np.all(triangles >= 0)
+    assert np.all(triangles < 4)
+
+
+def test_generate_mesh_grid_larger():
+    """Test with larger grid to ensure scaling works."""
+    n_rows, n_cols = 50, 50
+    positions, triangles = generate_mesh_grid(n_rows=n_rows, n_cols=n_cols)
+
+    expected_num_vertices = (n_rows + 1) * (n_cols + 1)  # 51 * 51 = 2601
+    expected_num_triangles = n_rows * n_cols * 2  # 50 * 50 * 2 = 5000
+
+    assert positions.shape == (expected_num_vertices, 2)
+    assert triangles.shape == (expected_num_triangles, 3)
+
+    # Most important: all indices must be valid
+    assert np.all(triangles >= 0)
+    assert np.all(triangles < expected_num_vertices)
+
+    # Check corners are correct
+    np.testing.assert_allclose(positions[0], [0, 0])  # bottom-left
+    np.testing.assert_allclose(positions[n_cols], [1, 0])  # bottom-right
+    np.testing.assert_allclose(positions[n_cols * (n_rows + 1)], [0, 1])  # top-left
+    np.testing.assert_allclose(
+        positions[(n_rows + 1) * (n_cols + 1) - 1],
+        [1, 1],
+    )  # top-right
+
+
+def test_triangle_winding():
+    """Test that triangles have consistent winding order."""
+    positions, triangles = generate_mesh_grid(n_rows=2, n_cols=2)
+
+    # Check that all triangles have counter-clockwise winding
+    for tri in triangles:
+        v0, v1, v2 = positions[tri]
+        # Compute cross product to check winding
+        edge1 = v1 - v0
+        edge2 = v2 - v0
+        cross = edge1[0] * edge2[1] - edge1[1] * edge2[0]
+        # For counter-clockwise winding in screen space, cross product should be positive
+        # (though this depends on coordinate system)
+        assert cross > 0