diff --git a/examples/explore_perlin.py b/examples/explore_perlin.py
new file mode 100644
index 0000000..37dc83c
--- /dev/null
+++ b/examples/explore_perlin.py
@@ -0,0 +1,345 @@
+"""Explore procedural Perlin terrain with hydro flow — no data files needed.
+
+Generates a synthetic island archipelago using xarray-spatial's ridged
+multi-octave Perlin noise with domain warping, all on GPU via dask + cupy.
+Hydraulic erosion carves realistic drainage channels and hydrological flow
+direction, accumulation, stream order, and stream links are recomputed for
+every window — including dynamically loaded ones.
+
+The terrain is infinite: fly to the edge and a new Perlin window is
+generated on the fly.  The fixed normalization (clip → scale) and fixed
+sea level ensure adjacent windows tile seamlessly.
+
+Controls:
+    WASD / arrows  Move camera
+    Shift+Y        Toggle hydro flow particle animation
+    G              Cycle data layers (elevation, slope, stream_link, ...)
+    R / Shift+R    Increase / decrease resolution
+    Y              Cycle color stretch
+
+Usage:
+    python explore_perlin.py
+    python explore_perlin.py --size 2048 --chunks 512 --seed 42
+    python explore_perlin.py --noise fbm --warp 0.5
+    python explore_perlin.py --erosion-iters 500000
+
+Requirements:
+    pip install rtxpy[all] xarray-spatial dask cupy scipy
+"""
+
+import argparse
+import time
+
+import cupy as cp
+import dask.array as da
+import numpy as np
+import xarray as xr
+from xrspatial import (
+    generate_terrain,
+    slope,
+    aspect,
+    fill as xrs_fill,
+    flow_direction,
+    flow_accumulation,
+    stream_order,
+    stream_link,
+)
+from xrspatial.erosion import erode
+from scipy.ndimage import uniform_filter
+
+import rtxpy  # noqa: F401 — registers .rtx accessor
+
+WINDOW_HALF = 10_000  # half-width of each window in metres (20 km total)
+ZFACTOR = 4000  # raw elevation multiplier
+
+# full_extent == window size → each window covers 1.0 noise-space units,
+# giving proper multi-octave detail.  Windows beyond this range still work
+# because Perlin noise is defined for all coordinates.
+FULL_EXTENT = (-WINDOW_HALF, -WINDOW_HALF, WINDOW_HALF, WINDOW_HALF)
+
+# Sea level as a fraction of ZFACTOR.  Ridged noise after fixed
+# normalisation sits in roughly [0.5·zf, zf].  Subtracting 0.55·zf
+# puts the coastline where the noise dips, creating an archipelago.
+SEA_LEVEL_FRAC = 0.55
+
+# Noise parameters (set from CLI args at startup, shared with loader).
+# NOTE: worley_blend is 0 because the dask+cupy path normalises Worley
+# per-chunk (min/max), creating visible grid seams.  Domain warping is
+# seamless (uses deterministic Perlin displacement).
+NOISE_PARAMS: dict = {}
+
+
+def _apply_sea_level(terrain):
+    """Subtract a fixed sea level and zero out ocean pixels.
+
+    Because the sea level is an absolute constant (not per-window), the
+    transform is identical for every window, preserving tile coherence.
+    """
+    sea = ZFACTOR * SEA_LEVEL_FRAC
+    terrain.data[:] = cp.maximum(terrain.data - sea, 0)
+    return terrain
+
+
+def generate_window(size, seed, x_range, y_range, chunks=None):
+    """Generate a cupy-backed elevation DataArray from Perlin noise.
+
+    When *chunks* is given, the template is dask+cupy so noise is generated
+    in parallel across GPU chunks before being materialised.
+    """
+    if chunks is not None:
+        data = da.zeros((size, size), dtype=np.float32,
+                        chunks=(chunks, chunks))
+        data = data.map_blocks(cp.asarray, dtype=np.float32,
+                               meta=cp.array((), dtype=np.float32))
+        template = xr.DataArray(data, dims=['y', 'x'])
+        print(f"Dask template: {size}x{size}, "
+              f"chunks={chunks}x{chunks} "
+              f"({(size // chunks) ** 2} tasks)")
+    else:
+        template = xr.DataArray(
+            cp.zeros((size, size), dtype=cp.float32), dims=['y', 'x'],
+        )
+
+    t0 = time.time()
+    terrain = generate_terrain(
+        template,
+        x_range=x_range,
+        y_range=y_range,
+        seed=seed,
+        zfactor=ZFACTOR,
+        full_extent=FULL_EXTENT,
+        **NOISE_PARAMS,
+    )
+
+    # Materialise dask graph → cupy
+    if isinstance(terrain.data, da.Array):
+        print("Computing dask graph on GPU...")
+        terrain = terrain.compute()
+    if not hasattr(terrain.data, 'device'):
+        terrain = terrain.copy(data=cp.asarray(terrain.data))
+
+    # Fixed sea-level cutoff (same for every window → seamless tiling)
+    terrain = _apply_sea_level(terrain)
+
+    dt = time.time() - t0
+    elev_min = float(cp.nanmin(terrain.data))
+    elev_max = float(cp.nanmax(terrain.data))
+    water_pct = float((terrain.data == 0).sum()) / terrain.data.size * 100
+    print(f"Generated {terrain.shape[0]}x{terrain.shape[1]} terrain "
+          f"in {dt:.2f}s  (elev {elev_min:.0f}–{elev_max:.0f} m, "
+          f"{water_pct:.0f}% water)")
+    return terrain
+
+
+def erode_terrain(terrain, iterations, seed):
+    """Apply hydraulic erosion to carve drainage channels."""
+    print(f"Eroding terrain ({iterations:,} droplets)...")
+    t0 = time.time()
+    terrain = erode(terrain, iterations=iterations, seed=seed)
+    if not hasattr(terrain.data, 'device'):
+        terrain = terrain.copy(data=cp.asarray(terrain.data))
+    # Re-clamp: erosion can push values slightly below 0
+    terrain.data[:] = cp.maximum(terrain.data, 0)
+    dt = time.time() - t0
+    print(f"  Erosion done in {dt:.2f}s")
+    return terrain
+
+
+def compute_hydro(terrain):
+    """Condition DEM and compute D8 hydro flow layers.
+
+    Returns a hydro dict ready for explore(), plus a stream_link DataArray
+    suitable for adding to a Dataset overlay.
+    """
+    print("Conditioning DEM for hydrological flow...")
+    t0 = time.time()
+
+    elev = cp.asnumpy(terrain.data).astype(np.float32)
+
+    # Mark water (sea-level pixels are 0) as ocean sentinel
+    ocean = (elev == 0.0) | np.isnan(elev)
+    elev[ocean] = -100.0
+
+    # Smooth to fill noise pits
+    smoothed = uniform_filter(elev, size=15, mode='nearest')
+    smoothed[ocean] = -100.0
+
+    sm_gpu = cp.asarray(smoothed)
+    filled = xrs_fill(terrain.copy(data=sm_gpu))
+
+    # Resolve flats: small perturbation toward drainage
+    delta = filled.data - sm_gpu
+    resolved = filled.data + delta * 0.01
+    cp.random.seed(0)
+    resolved += cp.random.uniform(0, 0.001, resolved.shape, dtype=cp.float32)
+    resolved[cp.asarray(ocean)] = -100.0
+
+    conditioned = terrain.copy(data=resolved)
+    fd = flow_direction(conditioned)
+    fa = flow_accumulation(fd)
+    so = stream_order(fd, fa, threshold=50)
+    sl = stream_link(fd, fa, threshold=50)
+
+    # Mask ocean back to NaN
+    ocean_gpu = cp.asarray(ocean)
+    fd.data[ocean_gpu] = cp.nan
+    fa.data[ocean_gpu] = cp.nan
+    so.data[ocean_gpu] = cp.nan
+    sl.data[ocean_gpu] = cp.nan
+
+    dt = time.time() - t0
+    n_streams = int(cp.nanmax(sl.data))
+    print(f"  Hydro computed in {dt:.2f}s — "
+          f"{n_streams} stream segments found")
+
+    # Clean stream_link for Dataset overlay (NaN → 0)
+    sl_clean = cp.nan_to_num(sl.data, nan=0.0).astype(cp.float32)
+
+    hydro = {
+        'flow_dir': fd.data,
+        'flow_accum': fa.data,
+        'stream_order': so.data,
+        'stream_link': sl.data,
+        'accum_threshold': 50,
+        'enabled': False,
+    }
+    return hydro, terrain.copy(data=sl_clean)
+
+
+def make_terrain_loader(size, seed, chunks, erosion_iters=0, do_hydro=True):
+    """Create a callback that generates new Perlin terrain at the camera.
+
+    No clamping — Perlin noise is defined everywhere, so exploration is
+    unlimited.  The full_extent mapping just sets the scale factor so each
+    window covers 1.0 noise-space units.
+
+    Each new window gets its own erosion and hydro computation so flow
+    patterns are consistent with the local terrain.  Returns a
+    ``(terrain_da, hydro_dict)`` tuple that the engine's terrain reload
+    handler picks up to reinitialise hydro particles.
+    """
+
+    def loader(cam_x, cam_y):
+        x_range = (cam_x - WINDOW_HALF, cam_x + WINDOW_HALF)
+        y_range = (cam_y - WINDOW_HALF, cam_y + WINDOW_HALF)
+        try:
+            terrain = generate_window(size, seed, x_range, y_range,
+                                      chunks=chunks)
+            if erosion_iters > 0:
+                terrain = erode_terrain(terrain, erosion_iters, seed)
+
+            hydro = None
+            if do_hydro:
+                try:
+                    hydro, _ = compute_hydro(terrain)
+                except Exception as e:
+                    print(f"Loader hydro error: {e}")
+
+            return (terrain, hydro)
+        except Exception as e:
+            print(f"Terrain loader error: {e}")
+            return None
+
+    return loader
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(
+        description="Explore procedural Perlin terrain with hydro flow.",
+    )
+    parser.add_argument("--size", type=int, default=2048,
+                        help="Grid size in pixels (default: 2048)")
+    parser.add_argument("--chunks", type=int, default=512,
+                        help="Dask chunk size (default: 512, 0=no dask)")
+    parser.add_argument("--seed", type=int, default=42,
+                        help="Perlin noise seed (default: 42)")
+    parser.add_argument("--noise", type=str, default='ridged',
+                        choices=['fbm', 'ridged'],
+                        help="Noise mode (default: ridged)")
+    parser.add_argument("--warp", type=float, default=0.4,
+                        help="Domain warp strength (default: 0.4)")
+    parser.add_argument("--worley", type=float, default=0.0,
+                        help="Worley cellular noise blend (default: 0.0, "
+                             "causes chunk seams in dask+cupy path)")
+    parser.add_argument("--octaves", type=int, default=6,
+                        help="Noise octaves (default: 6)")
+    parser.add_argument("--erosion-iters", type=int, default=200_000,
+                        help="Hydraulic erosion droplets (default: 200000)")
+    parser.add_argument("--no-erode", action="store_true",
+                        help="Skip hydraulic erosion")
+    parser.add_argument("--no-hydro", action="store_true",
+                        help="Skip hydro flow computation")
+    args = parser.parse_args()
+
+    chunks = args.chunks if args.chunks > 0 else None
+
+    # Noise params shared between initial window and terrain loader.
+    # Erosion is handled separately (not via generate_terrain's erode param)
+    # so it runs after sea-level subtraction.
+    NOISE_PARAMS = {
+        'noise_mode': args.noise,
+        'warp_strength': args.warp,
+        'worley_blend': args.worley,
+        'octaves': args.octaves,
+    }
+    desc = [args.noise]
+    if args.warp > 0:
+        desc.append(f"warp={args.warp}")
+    if args.worley > 0:
+        desc.append(f"worley={args.worley}")
+    if not args.no_erode:
+        desc.append(f"erode({args.erosion_iters:,})")
+    print(f"Noise: {', '.join(desc)}")
+
+    # ---- Generate initial terrain from Perlin noise ----------------------
+    x_range = (-WINDOW_HALF, WINDOW_HALF)
+    y_range = (-WINDOW_HALF, WINDOW_HALF)
+    terrain = generate_window(args.size, args.seed, x_range, y_range,
+                              chunks=chunks)
+
+    # ---- Hydraulic erosion ------------------------------------------------
+    if not args.no_erode:
+        terrain = erode_terrain(terrain, args.erosion_iters, args.seed)
+
+    # ---- Build Dataset with analysis layers ------------------------------
+    print("Computing terrain analysis layers...")
+    ds = xr.Dataset({
+        'elevation': terrain.rename(None),
+        'slope': slope(terrain),
+        'aspect': aspect(terrain),
+    })
+
+    # ---- Hydro flow ------------------------------------------------------
+    hydro = None
+    if not args.no_hydro:
+        try:
+            hydro, sl_layer = compute_hydro(terrain)
+            ds['stream_link'] = sl_layer.rename(None)
+        except Exception as e:
+            print(f"Skipping hydro: {e}")
+
+    # ---- Terrain loader for infinite Perlin exploration ------------------
+    # The loader runs in a background thread (engine does not block the
+    # render loop).  Cap erosion at 50k for reloaded windows so the load
+    # completes in a few seconds instead of 30–60s.
+    loader_erosion = 0 if args.no_erode else min(args.erosion_iters, 50_000)
+    loader = make_terrain_loader(
+        args.size, args.seed, chunks,
+        erosion_iters=loader_erosion,
+        do_hydro=not args.no_hydro,
+    )
+
+    # ---- Launch ----------------------------------------------------------
+    print(f"\nLaunching explore "
+          f"(G = cycle layers, Shift+Y = hydro flow)...\n")
+    ds.rtx.explore(
+        z='elevation',
+        width=2048,
+        height=1600,
+        render_scale=0.5,
+        color_stretch='cbrt',
+        terrain_loader=loader,
+        hydro_data=hydro,
+        repl=True,
+    )
+    print("Done")
diff --git a/examples/port_jefferson.py b/examples/port_jefferson.py
new file mode 100644
index 0000000..89190c1
--- /dev/null
+++ b/examples/port_jefferson.py
@@ -0,0 +1,28 @@
+"""Port Jefferson, NY — 1 m USGS 3DEP LiDAR DEM exploration.
+
+Downloads USGS 1-meter LiDAR-derived DEM tiles for Port Jefferson
+harbor and the surrounding village on Long Island's north shore,
+then launches the interactive rtxpy viewer with buildings, roads,
+and satellite imagery.
+
+First run downloads the DEM tiles from USGS (~30 MB per 10 km tile)
+and caches them as a zarr store. Subsequent runs load from cache.
+
+Usage:
+    python examples/port_jefferson.py
+"""
+from rtxpy import quickstart
+
+# Port Jefferson harbor + village — tight bbox to keep 1m DEM manageable
+# ~3.5 km E-W x ~3 km N-S
+quickstart(
+    name='port_jefferson',
+    bounds=(-73.10, 40.92, -73.04, 40.97),
+    crs='EPSG:32618',
+    source='usgs_1m',
+    features=['buildings', 'roads', 'water'],
+    tiles='satellite',
+    hydro=True,
+    wind=True,
+    weather=True,
+)
diff --git a/examples/trinidad.py b/examples/trinidad.py
index 1fd05a8..c3dd543 100644
--- a/examples/trinidad.py
+++ b/examples/trinidad.py
@@ -7,6 +7,7 @@
     crs='EPSG:32620',
     features=['buildings', 'roads', 'water', 'fire', 'places',
               'infrastructure', 'land_use'],
+    weather=True,
     hydro=True,
     coast_distance=True,
     ao_samples=1,
diff --git a/rtxpy/__init__.py b/rtxpy/__init__.py
index f8f52e4..ca3c137 100644
--- a/rtxpy/__init__.py
+++ b/rtxpy/__init__.py
@@ -10,6 +10,8 @@
 from .mesh import (
     triangulate_terrain,
     voxelate_terrain,
+    add_terrain_skirt,
+    build_terrain_skirt,
     write_stl,
     load_glb,
     load_mesh,
@@ -28,7 +30,7 @@
 
 # Optional convenience — network helpers with lazy dependency checks
 try:
-    from .remote_data import fetch_dem, fetch_lidar, fetch_osm, fetch_buildings, fetch_roads, fetch_water, fetch_wind, fetch_firms, fetch_places, fetch_infrastructure, fetch_land_use, fetch_restaurant_grades, fetch_gtfs
+    from .remote_data import fetch_dem, fetch_lidar, fetch_osm, fetch_buildings, fetch_roads, fetch_water, fetch_wind, fetch_weather, fetch_firms, fetch_places, fetch_infrastructure, fetch_land_use, fetch_restaurant_grades, fetch_gtfs
 except ImportError:
     pass
 
diff --git a/rtxpy/accessor.py b/rtxpy/accessor.py
index fdab96e..f39719a 100644
--- a/rtxpy/accessor.py
+++ b/rtxpy/accessor.py
@@ -2147,7 +2147,8 @@ def _safe_id(s):
         return results
 
     def triangulate(self, geometry_id='terrain', scale=1.0,
-                     pixel_spacing_x=None, pixel_spacing_y=None):
+                     pixel_spacing_x=None, pixel_spacing_y=None,
+                     skirt=True, skirt_depth=None):
         """Triangulate the terrain and add it to the scene.
 
         Creates a triangle mesh from the raster elevation data and adds it
@@ -2165,6 +2166,13 @@ def triangulate(self, geometry_id='terrain', scale=1.0,
         pixel_spacing_y : float, optional
             Y spacing between pixels in world units. If None, uses the
             auto-computed value from the DataArray coordinates.
+        skirt : bool, optional
+            If True, add a vertical skirt around the terrain edges to
+            give it a solid "slice of land" appearance. Default is True.
+        skirt_depth : float, optional
+            Depth of the skirt below the minimum elevation (in raw
+            elevation units before scaling). If None, defaults to 15%
+            of the elevation range.
 
         Returns
         -------
@@ -2180,7 +2188,7 @@ def triangulate(self, geometry_id='terrain', scale=1.0,
         >>> # Triangulate with real-world spacing (e.g., 25m per pixel)
         >>> verts, indices = dem.rtx.triangulate(pixel_spacing_x=25.0, pixel_spacing_y=25.0)
         """
-        from .mesh import triangulate_terrain
+        from .mesh import triangulate_terrain, add_terrain_skirt
         import numpy as np
 
         # Fall back to auto-computed spacing from __init__
@@ -2200,6 +2208,11 @@ def triangulate(self, geometry_id='terrain', scale=1.0,
         # Triangulate the terrain (creates vertices in pixel coordinates)
         triangulate_terrain(vertices, indices, self._obj, scale=scale)
 
+        # Add skirt (vertical walls + bottom face around edges)
+        if skirt:
+            vertices, indices = add_terrain_skirt(
+                vertices, indices, H, W, skirt_depth=skirt_depth)
+
         # Scale x,y coordinates to world units if pixel spacing != 1.0
         if pixel_spacing_x != 1.0 or pixel_spacing_y != 1.0:
             # Vertices are stored as [x0,y0,z0, x1,y1,z1, ...]
@@ -2211,9 +2224,11 @@ def triangulate(self, geometry_id='terrain', scale=1.0,
         self._pixel_spacing_y = pixel_spacing_y
         self._terrain_mesh_type = 'tin'
 
-        # Add to scene (pass grid dims for cluster-accelerated BVH)
+        # Add to scene — skip cluster grid_dims when skirt is present
+        # (skirt breaks the regular grid topology clusters require)
+        gd = (H, W) if not skirt else None
         self._rtx.add_geometry(geometry_id, vertices, indices,
-                               grid_dims=(H, W))
+                               grid_dims=gd)
 
         return vertices, indices
 
@@ -2295,7 +2310,7 @@ def voxelate(self, geometry_id='terrain', scale=1.0, base_elevation=None,
 
     def heightfield(self, geometry_id='terrain', scale=1.0,
                     pixel_spacing_x=None, pixel_spacing_y=None,
-                    tile_size=32):
+                    tile_size=32, skirt=True, skirt_depth=None):
         """Add the terrain as a heightfield using custom intersection.
 
         Instead of materializing an explicit triangle mesh, this uploads
@@ -2319,6 +2334,12 @@ def heightfield(self, geometry_id='terrain', scale=1.0,
         tile_size : int, optional
             Number of cells per tile dimension for AABB grouping.
             Smaller tiles = tighter BVH but more AABBs. Default is 32.
+        skirt : bool, optional
+            If True, add a vertical skirt around the terrain edges to
+            give it a solid "slice of land" appearance. Default is True.
+        skirt_depth : float, optional
+            Depth of the skirt below the minimum elevation. If None,
+            defaults to 15% of the elevation range.
 
         Returns
         -------
@@ -2326,6 +2347,7 @@ def heightfield(self, geometry_id='terrain', scale=1.0,
             The 2-D elevation array that was uploaded to the GPU.
         """
         import numpy as np
+        from .mesh import build_terrain_skirt
 
         if pixel_spacing_x is None:
             pixel_spacing_x = self._pixel_spacing_x
@@ -2356,6 +2378,14 @@ def heightfield(self, geometry_id='terrain', scale=1.0,
             tile_size=tile_size,
         )
 
+        # Add skirt as separate triangle geometry alongside the heightfield
+        if skirt:
+            sv, si = build_terrain_skirt(
+                elev, H, W, scale=1.0, skirt_depth=skirt_depth,
+                pixel_spacing_x=pixel_spacing_x,
+                pixel_spacing_y=pixel_spacing_y)
+            self._rtx.add_geometry('terrain_skirt', sv, si)
+
         self._pixel_spacing_x = pixel_spacing_x
         self._pixel_spacing_y = pixel_spacing_y
         self._terrain_mesh_type = 'heightfield'
@@ -2563,7 +2593,8 @@ def explore(self, width=800, height=600, render_scale=0.5,
                 start_position=None, look_at=None, key_repeat_interval=0.05,
                 pixel_spacing_x=None, pixel_spacing_y=None,
                 mesh_type='heightfield', color_stretch='linear', title=None,
-                subsample=1, wind_data=None, hydro_data=None, gtfs_data=None,
+                subsample=1, wind_data=None, weather_data=None,
+                hydro_data=None, gtfs_data=None,
                 terrain_loader=None,
                 scene_zarr=None, ao_samples=0, gi_bounces=1, denoise=False,
                 fog_density=0.0, fog_color=(0.7, 0.8, 0.9),
@@ -2705,6 +2736,7 @@ def explore(self, width=800, height=600, render_scale=0.5,
             baked_meshes=self._baked_meshes if self._baked_meshes else None,
             subsample=subsample,
             wind_data=wind_data,
+            weather_data=weather_data,
             hydro_data=hydro_data,
             gtfs_data=gtfs_data,
             accessor=self,
@@ -3034,8 +3066,10 @@ def explore(self, z, width=800, height=600, render_scale=0.5,
                 pixel_spacing_x=None, pixel_spacing_y=None,
                 mesh_type='heightfield', color_stretch='linear', title=None,
                 subtitle=None, legend=None,
-                subsample=1, wind_data=None, hydro_data=None,
+                subsample=1, wind_data=None, weather_data=None,
+                hydro_data=None,
                 gtfs_data=None,
+                terrain_loader=None,
                 scene_zarr=None,
                 ao_samples=0, gi_bounces=1, denoise=False,
                 fog_density=0.0, fog_color=(0.7, 0.8, 0.9),
@@ -3144,9 +3178,11 @@ def explore(self, z, width=800, height=600, render_scale=0.5,
             baked_meshes=terrain_da.rtx._baked_meshes if terrain_da.rtx._baked_meshes else None,
             subsample=subsample,
             wind_data=wind_data,
+            weather_data=weather_data,
             hydro_data=hydro_data,
             gtfs_data=gtfs_data,
             accessor=terrain_da.rtx,
+            terrain_loader=terrain_loader,
             scene_zarr=scene_zarr,
             ao_samples=ao_samples,
             gi_bounces=gi_bounces,
diff --git a/rtxpy/engine.py b/rtxpy/engine.py
index 1a4559d..d17bc58 100644
--- a/rtxpy/engine.py
+++ b/rtxpy/engine.py
@@ -173,6 +173,7 @@ def _hydro_splat_kernel(
         min_vis_age,  # int32 scalar — minimum visible age
         ref_depth,    # float32 scalar — depth-scaling reference distance
         terrain,      # (tH, tW) float32 — terrain elevation
+        depth_t,      # (sh, sw) float32 — ray-trace t-values for occlusion
         output,       # (sh, sw, 3) float32 — frame buffer (atomic add)
         # Camera basis — scalar args to avoid tiny GPU allocations
         cam_x, cam_y, cam_z,
@@ -278,6 +279,20 @@ def _hydro_splat_kernel(
         if sx < 0 or sx >= sw or sy < 0 or sy >= sh:
             return
 
+        # Depth test: cull particles occluded by terrain.
+        # Convert ray t-value at this pixel to forward depth, then compare
+        # to the particle's forward depth (already computed as `depth`).
+        if depth_t.shape[0] > 0:
+            t_val = depth_t[sy, sx]
+            if t_val > 0.0 and t_val < 1.0e20:
+                # Forward depth = t / sqrt(1 + u_cam^2 + v_cam^2)
+                u_px = (2.0 * float(sx) / float(sw) - 1.0) * fov_scale * aspect_ratio
+                v_px = (1.0 - 2.0 * float(sy) / float(sh)) * fov_scale
+                inv_cos = math.sqrt(1.0 + u_px * u_px + v_px * v_px)
+                terrain_fwd = t_val / inv_cos
+                if depth > terrain_fwd:
+                    return
+
         # Per-particle color and radius
         color_r = colors[pidx, 0]
         color_g = colors[pidx, 1]
@@ -1303,7 +1318,8 @@ def __init__(self, raster, width: int = 800, height: int = 600,
                  title: str = None,
                  subtitle: str = None,
                  legend: dict = None,
-                 subsample: int = 1):
+                 subsample: int = 1,
+                 skirt: bool = True):
         """
         Initialize the interactive viewer.
 
@@ -1342,7 +1358,7 @@ def __init__(self, raster, width: int = 800, height: int = 600,
         self.terrain = TerrainState(
             raster, pixel_spacing_x=pixel_spacing_x,
             pixel_spacing_y=pixel_spacing_y, mesh_type=mesh_type,
-            subsample=subsample,
+            subsample=subsample, skirt=skirt,
         )
 
         self.rtx = rtx
@@ -1585,6 +1601,12 @@ def __init__(self, raster, width: int = 800, height: int = 600,
                     spacing_y=self.pixel_spacing_y,
                     ve=1.0,
                 )
+                if self.terrain_skirt:
+                    sv, si = mesh_mod.build_terrain_skirt(
+                        terrain_np, H, W, scale=1.0,
+                        pixel_spacing_x=self.pixel_spacing_x,
+                        pixel_spacing_y=self.pixel_spacing_y)
+                    rtx.add_geometry('terrain_skirt', sv, si)
                 cache_key = (self.subsample_factor, mesh_type)
                 self._terrain_mesh_cache[cache_key] = (
                     None, None, terrain_np.copy(),
@@ -1605,6 +1627,11 @@ def __init__(self, raster, width: int = 800, height: int = 600,
                     idxs = np.zeros(nt * 3, dtype=np.int32)
                     mesh_mod.triangulate_terrain(verts, idxs, raster, scale=1.0)
 
+                    # Add skirt for TIN meshes
+                    if self.terrain_skirt:
+                        verts, idxs = mesh_mod.add_terrain_skirt(
+                            verts, idxs, H, W)
+
                 if self.pixel_spacing_x != 1.0 or self.pixel_spacing_y != 1.0:
                     verts[0::3] *= self.pixel_spacing_x
                     verts[1::3] *= self.pixel_spacing_y
@@ -1614,9 +1641,9 @@ def __init__(self, raster, width: int = 800, height: int = 600,
                     verts.copy(), idxs.copy(), terrain_np.copy(),
                 )
 
-                # Only pass grid_dims for TIN meshes — cluster GAS
-                # partitioning assumes regular grid triangle layout.
-                gd = (H, W) if mesh_type != 'voxel' else None
+                # Only pass grid_dims for TIN meshes without skirt —
+                # cluster GAS requires regular grid triangle layout.
+                gd = (H, W) if mesh_type != 'voxel' and not self.terrain_skirt else None
                 rtx.add_geometry('terrain', verts, idxs,
                                  grid_dims=gd)
 
@@ -2908,6 +2935,14 @@ def mesh_type(self):
     def mesh_type(self, value):
         self.terrain.mesh_type = value
 
+    @property
+    def terrain_skirt(self):
+        return self.terrain.terrain_skirt
+
+    @terrain_skirt.setter
+    def terrain_skirt(self, value):
+        self.terrain.terrain_skirt = value
+
     @property
     def _water_mask(self):
         return self.terrain._water_mask
@@ -3427,6 +3462,14 @@ def _rebuild_at_resolution(self, factor):
                     spacing_y=self.pixel_spacing_y,
                     ve=ve,
                 )
+                if self.terrain_skirt:
+                    sv, si = mesh_mod.build_terrain_skirt(
+                        terrain_np, H, W, scale=ve,
+                        pixel_spacing_x=self.pixel_spacing_x,
+                        pixel_spacing_y=self.pixel_spacing_y)
+                    self.rtx.add_geometry('terrain_skirt', sv, si)
+                elif self.rtx.has_geometry('terrain_skirt'):
+                    self.rtx.remove_geometry('terrain_skirt')
         else:
             if cache_key in self._terrain_mesh_cache:
                 # Cache hit — reuse pre-built mesh (stored at scale=1.0)
@@ -3457,6 +3500,10 @@ def _rebuild_at_resolution(self, factor):
                     indices = np.zeros(num_tris * 3, dtype=np.int32)
                     mesh_mod.triangulate_terrain(vertices, indices, sub, scale=1.0)
 
+                    if self.terrain_skirt:
+                        vertices, indices = mesh_mod.add_terrain_skirt(
+                            vertices, indices, H, W)
+
                 # Scale x,y to world units
                 if self.pixel_spacing_x != 1.0 or self.pixel_spacing_y != 1.0:
                     vertices[0::3] *= self.pixel_spacing_x
@@ -3474,7 +3521,7 @@ def _rebuild_at_resolution(self, factor):
             # 4. Replace terrain geometry (add_geometry overwrites existing key
             #    in-place, preserving dict insertion order and instance IDs)
             if self.rtx is not None:
-                gd = (H, W) if self.mesh_type != 'voxel' else None
+                gd = (H, W) if self.mesh_type != 'voxel' and not self.terrain_skirt else None
                 self.rtx.add_geometry('terrain', vertices, indices,
                                       grid_dims=gd)
 
@@ -3683,6 +3730,14 @@ def _rebuild_vertical_exaggeration(self, ve):
                     spacing_y=self.pixel_spacing_y,
                     ve=ve,
                 )
+                if self.terrain_skirt:
+                    sv, si = mesh_mod.build_terrain_skirt(
+                        terrain_np, H, W, scale=ve,
+                        pixel_spacing_x=self.pixel_spacing_x,
+                        pixel_spacing_y=self.pixel_spacing_y)
+                    self.rtx.add_geometry('terrain_skirt', sv, si)
+                elif self.rtx.has_geometry('terrain_skirt'):
+                    self.rtx.remove_geometry('terrain_skirt')
         else:
             if cache_key in self._terrain_mesh_cache:
                 verts_base, indices, terrain_np = self._terrain_mesh_cache[cache_key]
@@ -3712,6 +3767,10 @@ def _rebuild_vertical_exaggeration(self, ve):
                     mesh_mod.triangulate_terrain(vertices, indices, self.raster,
                                                  scale=1.0)
 
+                    if self.terrain_skirt:
+                        vertices, indices = mesh_mod.add_terrain_skirt(
+                            vertices, indices, H, W)
+
                 if self.pixel_spacing_x != 1.0 or self.pixel_spacing_y != 1.0:
                     vertices[0::3] *= self.pixel_spacing_x
                     vertices[1::3] *= self.pixel_spacing_y
@@ -3725,7 +3784,7 @@ def _rebuild_vertical_exaggeration(self, ve):
 
             # Replace terrain geometry (preserves dict insertion order)
             if self.rtx is not None:
-                gd = (H, W) if self.mesh_type != 'voxel' else None
+                gd = (H, W) if self.mesh_type != 'voxel' and not self.terrain_skirt else None
                 self.rtx.add_geometry('terrain', vertices, indices,
                                       grid_dims=gd)
 
@@ -4218,6 +4277,45 @@ def _toggle_firms(self):
         print(f"FIRMS fire: {'ON' if self._firms_visible else 'OFF'}")
         self._update_frame()
 
+    def _init_weather(self, weather_data):
+        """Interpolate weather variables from lat/lon grid onto terrain pixels.
+
+        Registers each variable (temperature, precipitation, cloud_cover,
+        pressure) as an overlay layer accessible via G-key cycling.
+        """
+        from .tiles import _build_latlon_grids
+        from scipy.interpolate import RegularGridInterpolator
+
+        raster = self._base_raster
+        H, W = raster.shape
+
+        lats_grid, lons_grid = _build_latlon_grids(raster)
+        points = np.stack([lats_grid.ravel(), lons_grid.ravel()], axis=-1)
+
+        w_lats = weather_data['lats']   # (ny,)
+        w_lons = weather_data['lons']   # (nx,)
+
+        variables = ['temperature', 'precipitation', 'cloud_cover', 'pressure']
+        units = {'temperature': '°C', 'precipitation': 'mm',
+                 'cloud_cover': '%', 'pressure': 'hPa'}
+        added = []
+
+        for var in variables:
+            if var not in weather_data:
+                continue
+            grid = weather_data[var]  # (ny, nx)
+            interp = RegularGridInterpolator(
+                (w_lats, w_lons), grid,
+                method='linear', bounds_error=False, fill_value=np.nan,
+            )
+            layer = interp(points).reshape(H, W).astype(np.float32)
+            _add_overlay(self, var, layer)
+            mean_val = np.nanmean(layer)
+            added.append(f"{var} {mean_val:.1f}{units.get(var, '')}")
+
+        if added:
+            print(f"  Weather: {len(added)} layers added ({', '.join(added)})")
+
     def _init_wind(self, wind_data):
         """Interpolate wind U/V from lat/lon grid onto the terrain pixel grid.
 
@@ -5266,6 +5364,21 @@ def _draw_hydro_on_frame(self, img):
 
         on_screen = valid & (sx_all >= 0) & (sx_all < sw) & (sy_all >= 0) & (sy_all < sh)
 
+        # Depth test against ray-traced terrain
+        depth_t_np = None
+        if hasattr(self, '_d_depth_t') and self._d_depth_t is not None and self._d_depth_t.shape[0] > 0:
+            depth_t_np = cp.asnumpy(self._d_depth_t) if hasattr(self._d_depth_t, 'get') else np.asarray(self._d_depth_t)
+        if depth_t_np is not None:
+            sx_c = np.clip(sx_all, 0, sw - 1)
+            sy_c = np.clip(sy_all, 0, sh - 1)
+            t_vals = depth_t_np[sy_c, sx_c]
+            u_px = (2.0 * sx_c / sw - 1.0) * fov_scale * aspect_ratio
+            v_px = (1.0 - 2.0 * sy_c / sh) * fov_scale
+            inv_cos = np.sqrt(1.0 + u_px**2 + v_px**2)
+            terrain_fwd = t_vals / inv_cos
+            occluded = (t_vals > 0) & (t_vals < 1e20) & (depth > terrain_fwd)
+            on_screen = on_screen & ~occluded
+
         ages = self._hydro_ages
         lifetimes = self._hydro_lifetimes
 
@@ -5390,6 +5503,11 @@ def _splat_hydro_gpu(self, d_frame):
         if not isinstance(terrain_data, cp.ndarray):
             terrain_data = cp.asarray(terrain_data)
 
+        # Depth buffer for occlusion culling (populated by _update_frame)
+        depth_t = getattr(self, '_d_depth_t', None)
+        if depth_t is None:
+            depth_t = cp.empty((0, 0), dtype=cp.float32)
+
         # Single kernel launch — alpha computed on GPU
         threadsperblock = 256
         blockspergrid = (total + threadsperblock - 1) // threadsperblock
@@ -5405,6 +5523,7 @@ def _splat_hydro_gpu(self, d_frame):
             int(self._hydro_min_visible_age),
             float(self._hydro_ref_depth),
             terrain_data,
+            depth_t,
             d_frame,
             float(cam_pos[0]), float(cam_pos[1]), float(cam_pos[2]),
             float(forward[0]), float(forward[1]), float(forward[2]),
@@ -5939,10 +6058,36 @@ def _sync_drone_from_pos_for(self, obs, pos):
                 self._calculate_viewshed(quiet=True)
 
     def _check_terrain_reload(self):
-        """Check if camera is near terrain edge and reload a new window if needed."""
+        """Check if camera is near terrain edge and reload a new window.
+
+        The terrain loader runs in a background thread so it doesn't block
+        the render loop (erosion/hydro can take many seconds).  Each tick
+        we either (a) submit a new loader job if near-edge, or (b) poll
+        for a completed result and swap in the new terrain.
+        """
         if self._terrain_loader is None:
             return
 
+        # --- Phase 2: check for completed background load ---
+        future = self.terrain._terrain_reload_future
+        if future is not None:
+            if not future.done():
+                return  # still computing — keep rendering
+            # Harvest result
+            self.terrain._terrain_reload_future = None
+            try:
+                result, cam_lon, cam_lat = future.result()
+            except Exception as e:
+                print(f"Terrain loader error: {e}")
+                self._last_reload_time = time.time()
+                return
+            if result is None:
+                self._last_reload_time = time.time()
+                return
+            self._apply_terrain_reload(result, cam_lon, cam_lat)
+            return
+
+        # --- Phase 1: decide whether to submit a new load ---
         now = time.time()
         if now - self._last_reload_time < self._reload_cooldown:
             return
@@ -5966,11 +6111,29 @@ def _check_terrain_reload(self):
         cam_lon = self._coord_origin_x + cam_col * self._coord_step_x
         cam_lat = self._coord_origin_y + cam_row * self._coord_step_y
 
-        # Call the terrain loader
-        new_raster = self._terrain_loader(cam_lon, cam_lat)
-        if new_raster is None:
-            self._last_reload_time = now
-            return
+        # Submit loader to background thread
+        from concurrent.futures import ThreadPoolExecutor
+        pool = self.terrain._terrain_reload_pool
+        if pool is None:
+            pool = ThreadPoolExecutor(max_workers=1)
+            self.terrain._terrain_reload_pool = pool
+
+        loader = self._terrain_loader
+
+        def _bg_load(lon, lat):
+            return loader(lon, lat), lon, lat
+
+        self.terrain._terrain_reload_future = pool.submit(_bg_load, cam_lon, cam_lat)
+        # Prevent re-submission while the load is in progress
+        self._last_reload_time = now + 999999
+
+    def _apply_terrain_reload(self, result, cam_lon, cam_lat):
+        """Apply a completed terrain reload result (runs on main thread)."""
+        new_hydro = None
+        if isinstance(result, tuple):
+            new_raster, new_hydro = result
+        else:
+            new_raster = result
 
         cam_z = self.position[2]
 
@@ -5991,6 +6154,7 @@ def _check_terrain_reload(self):
         self._hydro_terrain_np = None
         self._d_base_frame = None     # invalidate GPU wind/hydro buffers
         self._d_wind_scratch = None
+        self._d_depth_t = None        # invalidate depth buffer
 
         # Update coordinate tracking
         self._coord_origin_x = new_origin_x
@@ -6055,6 +6219,14 @@ def _check_terrain_reload(self):
                     spacing_y=self.pixel_spacing_y,
                     ve=ve,
                 )
+                if self.terrain_skirt:
+                    sv, si = mesh_mod.build_terrain_skirt(
+                        terrain_np, H, W, scale=ve,
+                        pixel_spacing_x=self.pixel_spacing_x,
+                        pixel_spacing_y=self.pixel_spacing_y)
+                    self.rtx.add_geometry('terrain_skirt', sv, si)
+                elif self.rtx.has_geometry('terrain_skirt'):
+                    self.rtx.remove_geometry('terrain_skirt')
             self._terrain_mesh_cache[cache_key] = (None, None, terrain_np.copy())
         else:
             if self.mesh_type == 'voxel':
@@ -6072,6 +6244,10 @@ def _check_terrain_reload(self):
                 indices = np.zeros(num_tris * 3, dtype=np.int32)
                 mesh_mod.triangulate_terrain(vertices, indices, new_raster, scale=1.0)
 
+                if self.terrain_skirt:
+                    vertices, indices = mesh_mod.add_terrain_skirt(
+                        vertices, indices, H, W)
+
             # Scale x,y to world units
             if self.pixel_spacing_x != 1.0 or self.pixel_spacing_y != 1.0:
                 vertices[0::3] *= self.pixel_spacing_x
@@ -6089,10 +6265,20 @@ def _check_terrain_reload(self):
 
             # Replace terrain geometry
             if self.rtx is not None:
-                gd = (H, W) if self.mesh_type != 'voxel' else None
+                gd = (H, W) if self.mesh_type != 'voxel' and not self.terrain_skirt else None
                 self.rtx.add_geometry('terrain', vertices, indices,
                                       grid_dims=gd)
 
+        # Reinitialize hydro if the loader provided new flow data
+        if new_hydro is not None and self._hydro_data is not None:
+            was_enabled = self._hydro_enabled
+            flow_dir = new_hydro['flow_dir']
+            flow_accum = new_hydro['flow_accum']
+            hydro_opts = {k: v for k, v in new_hydro.items()
+                          if k not in ('flow_dir', 'flow_accum', 'enabled')}
+            self._init_hydro(flow_dir, flow_accum, **hydro_opts)
+            self._hydro_enabled = was_enabled
+
         # Reposition camera in new window
         self.position = np.array([
             new_col * self.pixel_spacing_x,
@@ -7409,6 +7595,18 @@ def _update_frame(self):
         d_output = self._render_frame()
         self.frame_count += 1
 
+        # Extract depth buffer from primary hits for hydro occlusion culling.
+        # primary_hits[:, 0] holds ray t-values; persists until next render.
+        if self._hydro_enabled:
+            from .analysis.render import _render_buffers as _depth_bufs
+            if _depth_bufs.primary_hits is not None:
+                h, w = self.render_height, self.render_width
+                _ddt = getattr(self, '_d_depth_t', None)
+                if _ddt is None or _ddt.shape != (h, w):
+                    self._d_depth_t = cp.empty((h, w), dtype=cp.float32)
+                t_src = _depth_bufs.primary_hits.reshape(h * w, 4)[:, 0]
+                self._d_depth_t[:] = t_src.reshape(h, w)
+
         # Progressive accumulation (needed for AO convergence and/or DOF)
         needs_accum = self.ao_enabled or self.dof_enabled
         if needs_accum:
@@ -8478,6 +8676,12 @@ def _run_repl():
             sys.stdout.write('\033[?25h')  # show cursor
             sys.stdout.flush()
 
+        # Clean up terrain reload thread pool
+        pool = self.terrain._terrain_reload_pool
+        if pool is not None:
+            pool.shutdown(wait=False)
+            self.terrain._terrain_reload_pool = None
+
         # Clean up tile service
         if self._tile_service is not None:
             self._tile_service.shutdown()
@@ -8506,6 +8710,7 @@ def explore(raster, width: int = 800, height: int = 600,
             baked_meshes=None,
             subsample: int = 1,
             wind_data=None,
+            weather_data=None,
             hydro_data=None,
             gtfs_data=None,
             accessor=None,
@@ -8525,6 +8730,7 @@ def explore(raster, width: int = 800, height: int = 600,
             minimap_layer: str = None,
             minimap_colors: dict = None,
             info_text: str = None,
+            skirt: bool = True,
             repl: bool = False,
             tour=None):
     """
@@ -8670,6 +8876,7 @@ def explore(raster, width: int = 800, height: int = 600,
         subtitle=subtitle,
         legend=legend,
         subsample=subsample,
+        skirt=skirt,
     )
     viewer._geometry_colors_builder = geometry_colors_builder
     viewer._baked_meshes = baked_meshes or {}
@@ -8698,6 +8905,10 @@ def explore(raster, width: int = 800, height: int = 600,
     if wind_data is not None:
         viewer._init_wind(wind_data)
 
+    # Weather overlay initialization
+    if weather_data is not None:
+        viewer._init_weather(weather_data)
+
     # Hydro flow initialization
     if hydro_data is not None:
         hydro_start_enabled = hydro_data.get('enabled', True)
diff --git a/rtxpy/mesh.py b/rtxpy/mesh.py
index 4e9a260..8102e93 100644
--- a/rtxpy/mesh.py
+++ b/rtxpy/mesh.py
@@ -133,6 +133,193 @@ def triangulate_terrain(verts, triangles, terrain, scale=1.0):
     return 0
 
 
+def add_terrain_skirt(vertices, indices, H, W, skirt_depth=None):
+    """Add a vertical skirt around the edges of a TIN terrain mesh.
+
+    Creates vertical wall faces dropping from the terrain perimeter down
+    to a base level, plus a flat bottom face, giving the terrain a solid
+    "slice of land" appearance.
+
+    Parameters
+    ----------
+    vertices : np.ndarray
+        Flat float32 vertex buffer from triangulate_terrain, shape (H*W*3,).
+    indices : np.ndarray
+        Flat int32 index buffer from triangulate_terrain.
+    H : int
+        Number of rows in the terrain grid.
+    W : int
+        Number of columns in the terrain grid.
+    skirt_depth : float, optional
+        Depth of the skirt below the minimum terrain elevation.
+        If None, defaults to 15% of the elevation range.
+
+    Returns
+    -------
+    new_vertices : np.ndarray
+        Extended vertex buffer with skirt vertices appended.
+    new_indices : np.ndarray
+        Extended index buffer with skirt triangles appended.
+    """
+    z_vals = vertices[2::3]
+    z_min = float(np.nanmin(z_vals))
+    z_max = float(np.nanmax(z_vals))
+
+    if skirt_depth is None:
+        z_range = z_max - z_min
+        skirt_depth = max(z_range * 0.15, 1.0)
+
+    skirt_z = z_min - skirt_depth
+
+    # Build perimeter vertex indices in clockwise order (viewed from above).
+    # Top row L→R, right col T→B, bottom row R→L, left col B→T.
+    top = np.arange(W, dtype=np.int32)
+    right = (np.arange(1, H, dtype=np.int32)) * W + (W - 1)
+    bottom = (H - 1) * W + np.arange(W - 2, -1, -1, dtype=np.int32)
+    left = np.arange(H - 2, 0, -1, dtype=np.int32) * W
+
+    perim = np.concatenate([top, right, bottom, left])
+    n_perim = len(perim)
+    n_orig_verts = H * W
+
+    # Skirt vertices: same x,y as perimeter, z = skirt_z
+    skirt_verts = np.empty(n_perim * 3, dtype=np.float32)
+    skirt_verts[0::3] = vertices[perim * 3]
+    skirt_verts[1::3] = vertices[perim * 3 + 1]
+    skirt_verts[2::3] = skirt_z
+
+    # Wall triangles — two per perimeter edge segment.
+    # CW perimeter traversal with this winding gives outward-facing normals:
+    #   tri1: top_a, bot_b, top_b
+    #   tri2: top_a, bot_a, bot_b
+    idx = np.arange(n_perim, dtype=np.int32)
+    idx_next = np.roll(idx, -1)
+
+    top_a = perim
+    top_b = perim[idx_next]
+    bot_a = (n_orig_verts + idx).astype(np.int32)
+    bot_b = (n_orig_verts + idx_next).astype(np.int32)
+
+    wall_tris = np.empty(n_perim * 6, dtype=np.int32)
+    wall_tris[0::6] = top_a
+    wall_tris[1::6] = bot_b
+    wall_tris[2::6] = top_b
+    wall_tris[3::6] = top_a
+    wall_tris[4::6] = bot_a
+    wall_tris[5::6] = bot_b
+
+    # Bottom face — two triangles using the four corner skirt vertices.
+    c00 = np.int32(n_orig_verts)
+    cW1 = np.int32(n_orig_verts + W - 1)
+    cHW = np.int32(n_orig_verts + W + H - 2)
+    cH0 = np.int32(n_orig_verts + 2 * W + H - 3)
+
+    bottom_tris = np.array([
+        c00, cHW, cW1,
+        c00, cH0, cHW,
+    ], dtype=np.int32)
+
+    new_vertices = np.concatenate([vertices, skirt_verts])
+    new_indices = np.concatenate([indices, wall_tris, bottom_tris])
+
+    return new_vertices, new_indices
+
+
+def build_terrain_skirt(terrain_data, H, W, scale=1.0, skirt_depth=None,
+                        pixel_spacing_x=1.0, pixel_spacing_y=1.0):
+    """Build standalone skirt geometry for a terrain grid.
+
+    Creates a triangle mesh of just the vertical walls and bottom face
+    around the terrain perimeter. Useful for adding a skirt alongside
+    heightfield terrain (which uses custom intersection, not triangles).
+
+    Parameters
+    ----------
+    terrain_data : np.ndarray
+        2D elevation array of shape (H, W).
+    H : int
+        Number of rows.
+    W : int
+        Number of columns.
+    scale : float, optional
+        Scale factor for elevation values. Default is 1.0.
+    skirt_depth : float, optional
+        Depth below minimum elevation. If None, defaults to 15% of range.
+    pixel_spacing_x : float, optional
+        X spacing in world units. Default is 1.0.
+    pixel_spacing_y : float, optional
+        Y spacing in world units. Default is 1.0.
+
+    Returns
+    -------
+    vertices : np.ndarray
+        Flat float32 vertex buffer.
+    indices : np.ndarray
+        Flat int32 index buffer.
+    """
+    # Build perimeter indices (clockwise from above)
+    top = np.arange(W, dtype=np.int32)
+    right = (np.arange(1, H, dtype=np.int32)) * W + (W - 1)
+    bottom = (H - 1) * W + np.arange(W - 2, -1, -1, dtype=np.int32)
+    left = np.arange(H - 2, 0, -1, dtype=np.int32) * W
+    perim = np.concatenate([top, right, bottom, left])
+    n_perim = len(perim)
+
+    perim_h = perim // W
+    perim_w = perim % W
+
+    # Elevation at perimeter vertices
+    perim_z = np.array(
+        [terrain_data[h, w] * scale for h, w in zip(perim_h, perim_w)],
+        dtype=np.float32)
+
+    z_min = float(np.nanmin(terrain_data * scale))
+    z_max = float(np.nanmax(terrain_data * scale))
+    if skirt_depth is None:
+        skirt_depth = max((z_max - z_min) * 0.15, 1.0)
+    skirt_z = z_min - skirt_depth
+
+    # Replace NaN edge elevations with skirt_z so walls still render
+    nan_mask = np.isnan(perim_z)
+    perim_z[nan_mask] = skirt_z
+
+    # Vertices: first n_perim = terrain edge, next n_perim = skirt base
+    vertices = np.empty(n_perim * 2 * 3, dtype=np.float32)
+    vertices[0:n_perim * 3:3] = perim_w.astype(np.float32) * pixel_spacing_x
+    vertices[1:n_perim * 3:3] = perim_h.astype(np.float32) * pixel_spacing_y
+    vertices[2:n_perim * 3:3] = perim_z
+    off = n_perim * 3
+    vertices[off + 0::3] = perim_w.astype(np.float32) * pixel_spacing_x
+    vertices[off + 1::3] = perim_h.astype(np.float32) * pixel_spacing_y
+    vertices[off + 2::3] = skirt_z
+
+    # Wall triangles
+    idx = np.arange(n_perim, dtype=np.int32)
+    idx_next = np.roll(idx, -1)
+    top_a = idx
+    top_b = idx_next
+    bot_a = (n_perim + idx).astype(np.int32)
+    bot_b = (n_perim + idx_next).astype(np.int32)
+
+    wall_tris = np.empty(n_perim * 6, dtype=np.int32)
+    wall_tris[0::6] = top_a
+    wall_tris[1::6] = bot_b
+    wall_tris[2::6] = top_b
+    wall_tris[3::6] = top_a
+    wall_tris[4::6] = bot_a
+    wall_tris[5::6] = bot_b
+
+    # Bottom face
+    c00 = np.int32(n_perim)
+    cW1 = np.int32(n_perim + W - 1)
+    cHW = np.int32(n_perim + W + H - 2)
+    cH0 = np.int32(n_perim + 2 * W + H - 3)
+    bottom_tris = np.array([c00, cHW, cW1, c00, cH0, cHW], dtype=np.int32)
+
+    indices = np.concatenate([wall_tris, bottom_tris])
+    return vertices, indices
+
+
 @cuda.jit
 def _voxelate_terrain_gpu(verts, triangles, data, H, W, scale, base_z, stride):
     """GPU kernel for terrain voxelation — one box-column per cell."""
diff --git a/rtxpy/quickstart.py b/rtxpy/quickstart.py
index d7a6301..58525dc 100644
--- a/rtxpy/quickstart.py
+++ b/rtxpy/quickstart.py
@@ -400,6 +400,7 @@ def quickstart(
     tiles='satellite',
     tile_zoom=None,
     wind=True,
+    weather=False,
     hydro=False,
     coast_distance=False,
     cache_dir=None,
@@ -574,6 +575,15 @@ def quickstart(
         except Exception as e:
             print(f"Skipping GTFS realtime: {e}")
 
+    # -- weather --------------------------------------------------------------
+    weather_data = None
+    if weather:
+        try:
+            from .remote_data import fetch_weather as _fetch_weather
+            weather_data = _fetch_weather(bounds, grid_size=15)
+        except Exception as e:
+            print(f"Skipping weather: {e}")
+
     # -- wind -----------------------------------------------------------------
     wind_data = None
     if wind:
@@ -923,7 +933,8 @@ def quickstart(
 
     print("\nLaunching explore...\n")
     ds.rtx.explore(z='elevation', scene_zarr=zarr_path,
-                   wind_data=wind_data, hydro_data=hydro_data,
+                   wind_data=wind_data, weather_data=weather_data,
+                   hydro_data=hydro_data,
                    gtfs_data=gtfs_data, **defaults)
 
 
diff --git a/rtxpy/remote_data.py b/rtxpy/remote_data.py
index 0a23cc4..7ef2781 100644
--- a/rtxpy/remote_data.py
+++ b/rtxpy/remote_data.py
@@ -1941,6 +1941,97 @@ def fetch_wind(bounds, grid_size=20):
     }
 
 
+def fetch_weather(bounds, grid_size=20):
+    """Fetch current weather conditions from Open-Meteo for a bounding box.
+
+    Queries the Open-Meteo forecast API for temperature, precipitation,
+    cloud cover, and surface pressure on a regular lat/lon grid.
+
+    Parameters
+    ----------
+    bounds : tuple of float
+        (west, south, east, north) in WGS84 degrees.
+    grid_size : int
+        Number of grid points along each axis (default 20).
+        Total API points = grid_size².  Open-Meteo allows up to
+        ~1 000 points per request.
+
+    Returns
+    -------
+    dict
+        ``'temperature'`` : ndarray (ny, nx) — 2 m temperature (°C).
+        ``'precipitation'`` : ndarray (ny, nx) — precipitation (mm).
+        ``'cloud_cover'`` : ndarray (ny, nx) — cloud cover (%).
+        ``'pressure'`` : ndarray (ny, nx) — surface pressure (hPa).
+        ``'lats'`` : ndarray (ny,) — latitude values.
+        ``'lons'`` : ndarray (nx,) — longitude values.
+
+    Examples
+    --------
+    >>> from rtxpy import fetch_weather
+    >>> wx = fetch_weather((-73.10, 40.92, -73.04, 40.97), grid_size=5)
+    >>> wx['temperature'].shape
+    (5, 5)
+    """
+    try:
+        import requests
+    except ImportError:
+        raise ImportError(
+            "requests is required for fetch_weather(). "
+            "Install it with: pip install requests"
+        )
+    import numpy as np
+
+    west, south, east, north = bounds
+
+    lons = np.linspace(west, east, grid_size)
+    lats = np.linspace(south, north, grid_size)
+    grid_lons, grid_lats = np.meshgrid(lons, lats)
+
+    lat_str = ",".join(f"{v:.4f}" for v in grid_lats.ravel())
+    lon_str = ",".join(f"{v:.4f}" for v in grid_lons.ravel())
+
+    print(f"Fetching weather data ({grid_size}x{grid_size} grid)...")
+    resp = requests.get(
+        _OPEN_METEO_URL,
+        params={
+            "latitude": lat_str,
+            "longitude": lon_str,
+            "current": "temperature_2m,precipitation,cloud_cover,surface_pressure",
+        },
+        timeout=30,
+    )
+    resp.raise_for_status()
+    data = resp.json()
+
+    ny, nx = grid_size, grid_size
+    temperature = np.empty((ny, nx), dtype=np.float32)
+    precipitation = np.empty((ny, nx), dtype=np.float32)
+    cloud_cover = np.empty((ny, nx), dtype=np.float32)
+    pressure = np.empty((ny, nx), dtype=np.float32)
+
+    for i, point in enumerate(data):
+        row = i // nx
+        col = i % nx
+        current = point.get("current", point)
+        temperature[row, col] = current["temperature_2m"]
+        precipitation[row, col] = current["precipitation"]
+        cloud_cover[row, col] = current["cloud_cover"]
+        pressure[row, col] = current["surface_pressure"]
+
+    mean_temp = float(np.mean(temperature))
+    print(f"  Mean temperature: {mean_temp:.1f}°C")
+
+    return {
+        "temperature": temperature,
+        "precipitation": precipitation,
+        "cloud_cover": cloud_cover,
+        "pressure": pressure,
+        "lats": lats,
+        "lons": lons,
+    }
+
+
 # ---------------------------------------------------------------------------
 # NASA FIRMS fire detection footprints
 # ---------------------------------------------------------------------------
diff --git a/rtxpy/viewer/terrain.py b/rtxpy/viewer/terrain.py
index 400b442..c57e1b3 100644
--- a/rtxpy/viewer/terrain.py
+++ b/rtxpy/viewer/terrain.py
@@ -18,14 +18,16 @@ class TerrainState:
         '_gpu_terrain', '_gpu_base_terrain',
         'mesh_type', '_water_mask',
         'vertical_exaggeration', '_land_color_range',
+        'terrain_skirt',
         '_terrain_loader',
         '_coord_origin_x', '_coord_origin_y',
         '_coord_step_x', '_coord_step_y',
         '_reload_cooldown', '_last_reload_time',
+        '_terrain_reload_future', '_terrain_reload_pool',
     )
 
     def __init__(self, raster, pixel_spacing_x=1.0, pixel_spacing_y=1.0,
-                 mesh_type='heightfield', subsample=1):
+                 mesh_type='heightfield', subsample=1, skirt=True):
         self.raster = raster
         self._base_raster = raster
         self.pixel_spacing_x = pixel_spacing_x
@@ -35,6 +37,7 @@ def __init__(self, raster, pixel_spacing_x=1.0, pixel_spacing_y=1.0,
         self.mesh_type = mesh_type
         self.subsample_factor = max(1, int(subsample))
         self.vertical_exaggeration = 1.0
+        self.terrain_skirt = skirt
 
         # Elevation stats (set by viewer __init__ after ocean-fill)
         self.terrain_shape = (0, 0)
@@ -58,3 +61,5 @@ def __init__(self, raster, pixel_spacing_x=1.0, pixel_spacing_y=1.0,
         self._coord_step_y = -1.0
         self._reload_cooldown = 2.0
         self._last_reload_time = 0.0
+        self._terrain_reload_future = None
+        self._terrain_reload_pool = None