From f57cbbfe9dac7b1ccaffad4a14712186eb2f8fee Mon Sep 17 00:00:00 2001
From: Brent Schroeter <brent@secondsystemtech.com>
Date: Sat, 17 May 2025 16:22:08 -0700
Subject: [PATCH] improve topographical accuracy

---
 src/__init__.py | 47 ++++++++++++++++++++++++++++++++---------------
 1 file changed, 32 insertions(+), 15 deletions(-)

diff --git a/src/__init__.py b/src/__init__.py
index 665fb2c..2d18160 100644
--- a/src/__init__.py
+++ b/src/__init__.py
@@ -32,11 +32,16 @@ def haversine_dist(lat1, lng1, lat2, lng2):
     return 2 * EARTH_RADIUS * asin(sqrt(hav_theta))
 
 
-def get_elev_grid(elev_map, map_bounds, resolution, scale_factor):
+def get_elev_grid(elev_map, map_bounds, scale_factor, resolution):
     """
     Generate a 2D array of elevation values (in model units), smoothed and
     "blocked" to avoid jagged and/or superfluous edges.
     """
+    # Are there better factoring algorithms? Yes. Do we need one here? Nope!
+    block_sizes = sorted([
+        int(resolution / x) for x in range(2, resolution) if resolution % x == 0
+    ])
+
     elev_grid = np.zeros((resolution, resolution))
     for i in range(resolution):
         for j in range(resolution):
@@ -47,7 +52,8 @@ def get_elev_grid(elev_map, map_bounds, resolution, scale_factor):
     elev_grid = (elev_grid - elev_grid.min()) * scale_factor
 
     # Convolve a uniform kernel over matrix to smooth values
-    KERNEL_SIZE = 5
+    # TODO: adjust this dynamically based on meters per grid cell
+    KERNEL_SIZE = 11
     kernel_offset = (KERNEL_SIZE - 1) / 2
     assert kernel_offset % 1 == 0
     kernel_offset = int(kernel_offset)
@@ -57,15 +63,20 @@ def get_elev_grid(elev_map, map_bounds, resolution, scale_factor):
             # Clamping to the array bounds makes the edges of the model less
             # smooth, but that's acceptable for these purposes
             windowed = elev_grid[
-                max(i - kernel_offset, 0):min(i + kernel_offset, elev_grid.shape[0]),
-                max(j - kernel_offset, 0):min(j + kernel_offset, elev_grid.shape[1]),
+                max(i - kernel_offset, 0):min(i + kernel_offset, resolution),
+                max(j - kernel_offset, 0):min(j + kernel_offset, resolution),
             ]
             smoothed[i, j] = windowed.mean()
 
-    BLOCK_TOLERANCE = 0.1
-    for block_size in range(2, resolution + 1):
-        for i in range(resolution - block_size + 1):
-            for j in range(resolution - block_size + 1):
+    BLOCK_TOLERANCE = 0.1  # Higher tolerance means more aggressive flattening
+    for block_size in block_sizes:
+        assert resolution % block_size == 0
+        block_size = int(block_size)
+        # Sweeping a window of block_size with an increment of 1 is smoother,
+        # but it may smear topographical features towards the top left of the
+        # model
+        for i in range(0, resolution - block_size + 1, block_size):
+            for j in range(0, resolution - block_size + 1, block_size):
                 block = smoothed[i:i+block_size, j:j+block_size].flatten()
                 if block.max() - block.min() < BLOCK_TOLERANCE:
                     smoothed[i:i+block_size, j:j+block_size] = block.mean()
@@ -266,18 +277,24 @@ if __name__ == "__main__":
         help="<min lon>,<min lat>,<max lon>,<max lat>",
         required=True,
     )
-    parser.add_argument(
-        "--elev-exaggeration",
-        type=float,
-        help="Additional scaling factor to apply to elevations for dramatic effect",
-        default=1.0,
-    )
     parser.add_argument(
         "--slices",
         type=int,
         help="Parallelize render and increase elevation resolution by splitting model into n tiles per side (--slices=3 splits work into 9 tiles)",
         default=ceil(sqrt(cpu_count())),
     )
+    parser.add_argument(
+        "--elev-resolution",
+        type=int,
+        help="Resolution of elevation grid, per tile",
+        default=30,
+    )
+    parser.add_argument(
+        "--elev-exaggeration",
+        type=float,
+        help="Additional scaling factor to apply to elevations for dramatic effect",
+        default=1.0,
+    )
     parser.add_argument(
         "--edge-len",
         type=float,
@@ -367,7 +384,7 @@ if __name__ == "__main__":
 
     print("Calculating elevations...")
     # Slice the model area into a rectilinear grid to approximate topography
-    elev_grid = get_elev_grid(elev_map, map_bounds, args.slices * 25, scale_factor)
+    elev_grid = get_elev_grid(elev_map, map_bounds, scale_factor, args.elev_resolution * args.slices)
     elev_grid *= args.elev_exaggeration
 
     renderers = []