brent/geo3dp
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

improve topographical accuracy

Browse source
This commit is contained in:
Brent Schroeter 2025-05-17 16:22:08 -07:00
parent 2337591f6f
commit f57cbbfe9d
1 changed files with 32 additions and 15 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

47
src/__init__.py
View file

@ -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 = []