324 lines
12 KiB
Python
324 lines
12 KiB
Python
"""
|
||
Lagrangian particle trajectory simulator – compact block output.
|
||
|
||
10 000 particles seeded along the Catalan coast, advected by a synthetic
|
||
Western Mediterranean current field for 72 h (3 days).
|
||
Output: 3 JSON blocks of 24 h each + a metadata file.
|
||
Each block stores per-step arrays of [lon, lat, u, v, beached] per particle.
|
||
Trails are reconstructed client-side for zero data duplication.
|
||
"""
|
||
|
||
import json
|
||
import math
|
||
import os
|
||
import numpy as np
|
||
|
||
DATA_DIR = os.path.join(os.path.dirname(__file__), "data")
|
||
|
||
NUM_STEPS = 72
|
||
DT_SECONDS = 3600
|
||
DEG_PER_M_LAT = 1.0 / 111_320
|
||
BLOCK_SIZE = 24 # steps per block
|
||
PARTICLES_PER_SEED = 333
|
||
|
||
JITTER_DEG = 0.008 # small spread so particles start in water near coast
|
||
OFFSHORE_DEG = 0.012 # push seaward if initial position is on land
|
||
|
||
SEED_POINTS = [
|
||
(3.21, 41.39), (2.82, 41.62), (2.17, 41.45), (1.82, 41.22),
|
||
(1.25, 41.12), (1.00, 40.95), (0.63, 40.62), (0.85, 40.72),
|
||
(3.12, 41.98), (3.19, 42.30), (2.82, 41.97), (2.07, 41.38),
|
||
(3.10, 41.77), (3.17, 41.70), (1.52, 41.08), (3.04, 41.87),
|
||
(0.75, 40.78), (2.93, 41.69), (1.60, 41.15), (2.40, 41.53),
|
||
(0.50, 40.55), (2.60, 41.57), (3.16, 42.12), (2.50, 41.55),
|
||
(1.40, 41.05), (3.05, 42.05), (2.92, 41.83), (1.15, 40.90),
|
||
(2.70, 41.64), (3.22, 42.20),
|
||
]
|
||
|
||
SEED_NAMES = [
|
||
"Barcelona", "Badalona", "Castelldefels", "Vilanova",
|
||
"Tarragona", "Cambrils", "Delta Ebre S", "Delta Ebre N",
|
||
"Roses", "Portbou", "L'Estartit", "Sitges",
|
||
"Blanes", "Lloret de Mar", "Salou", "Tossa de Mar",
|
||
"L'Ampolla", "St. Feliu Guíxols", "Torredembarra", "El Masnou",
|
||
"Alcanar", "Premià de Mar", "Cadaqués", "Vilassar de Mar",
|
||
"Altafulla", "L'Escala", "Palamós", "Miami Platja",
|
||
"Arenys de Mar", "Llançà",
|
||
]
|
||
|
||
# ── Land polygons ────────────────────────────────────────────────────────────
|
||
|
||
IBERIAN_PENINSULA = [
|
||
(-2.0, 36.7), (-1.5, 36.8), (-0.5, 37.4), (0.0, 38.0),
|
||
(0.2, 38.8), (-0.2, 39.5), (0.0, 39.9), (0.2, 40.4),
|
||
(0.45, 40.5), (0.5, 40.65), (0.7, 40.75), (0.8, 40.85),
|
||
(0.9, 40.9), (1.0, 41.0), (1.2, 41.1), (1.5, 41.1),
|
||
(1.8, 41.2), (2.0, 41.3), (2.1, 41.35), (2.2, 41.4),
|
||
(2.5, 41.5), (2.8, 41.6), (2.9, 41.65), (3.0, 41.7),
|
||
(3.1, 41.8), (3.15, 41.9), (3.15, 42.0), (3.2, 42.1),
|
||
(3.2, 42.3), (3.25, 42.45), (3.05, 42.5), (2.5, 42.6),
|
||
(1.8, 42.7), (1.0, 42.75), (0.0, 42.7), (-1.0, 42.8),
|
||
(-2.0, 42.5), (-2.0, 36.7),
|
||
]
|
||
FRANCE_SOUTH = [
|
||
(3.05, 42.5), (3.25, 42.45), (3.3, 42.5), (3.5, 42.6),
|
||
(3.6, 43.0), (3.5, 43.2), (3.9, 43.3), (4.1, 43.4),
|
||
(4.5, 43.3), (4.8, 43.4), (5.0, 43.3), (5.5, 43.2),
|
||
(5.8, 43.1), (6.0, 43.1), (6.3, 43.15), (6.6, 43.0),
|
||
(7.0, 43.5), (7.5, 43.8), (7.5, 44.5), (6.0, 44.5),
|
||
(4.0, 44.0), (3.0, 43.5), (3.0, 42.8), (3.05, 42.5),
|
||
]
|
||
ITALY_WEST = [
|
||
(7.5, 43.8), (7.7, 43.8), (8.0, 43.9), (8.5, 44.0),
|
||
(9.0, 44.1), (9.5, 44.3), (10.0, 44.0), (10.0, 43.5),
|
||
(10.5, 43.0), (10.5, 42.5), (11.0, 42.0), (11.5, 41.5),
|
||
(12.0, 41.2), (12.5, 41.0), (13.0, 40.5), (13.5, 40.0),
|
||
(14.0, 39.5), (14.0, 38.5), (14.0, 44.5), (7.5, 44.5),
|
||
(7.5, 43.8),
|
||
]
|
||
SARDINIA = [
|
||
(8.1, 39.0), (8.3, 38.8), (8.5, 38.85), (8.8, 38.9),
|
||
(9.2, 39.0), (9.6, 39.2), (9.7, 39.8), (9.6, 40.2),
|
||
(9.5, 40.6), (9.2, 40.9), (8.8, 41.0), (8.3, 41.0),
|
||
(8.1, 40.8), (8.1, 40.4), (8.2, 39.8), (8.1, 39.0),
|
||
]
|
||
CORSICA = [
|
||
(8.6, 41.4), (8.8, 41.4), (9.2, 41.6), (9.4, 42.0),
|
||
(9.5, 42.4), (9.4, 42.8), (9.2, 43.0), (8.8, 42.7),
|
||
(8.6, 42.3), (8.6, 41.8), (8.6, 41.4),
|
||
]
|
||
MALLORCA = [
|
||
(2.3, 39.25), (2.7, 39.3), (3.1, 39.5), (3.4, 39.7),
|
||
(3.5, 39.85), (3.4, 39.95), (3.1, 39.95), (2.7, 39.9),
|
||
(2.4, 39.75), (2.3, 39.55), (2.3, 39.25),
|
||
]
|
||
MENORCA = [
|
||
(3.8, 39.8), (3.9, 39.82), (4.1, 39.95), (4.25, 40.05),
|
||
(4.1, 40.1), (3.85, 40.0), (3.8, 39.9), (3.8, 39.8),
|
||
]
|
||
IBIZA = [
|
||
(1.2, 38.8), (1.3, 38.82), (1.45, 38.9), (1.55, 39.0),
|
||
(1.5, 39.1), (1.35, 39.1), (1.2, 39.0), (1.15, 38.9),
|
||
(1.2, 38.8),
|
||
]
|
||
NORTH_AFRICA = [
|
||
(-2.0, 36.7), (-2.0, 35.0), (10.0, 35.0), (10.0, 37.0),
|
||
(9.5, 37.2), (9.0, 37.0), (8.5, 37.0), (8.0, 36.8),
|
||
(7.5, 36.9), (7.0, 37.0), (6.0, 37.0), (5.0, 36.8),
|
||
(4.0, 36.7), (3.0, 36.7), (2.0, 36.5), (1.0, 36.3),
|
||
(0.0, 35.8), (-1.0, 35.5), (-2.0, 36.7),
|
||
]
|
||
|
||
LAND_POLYGONS = [
|
||
IBERIAN_PENINSULA, FRANCE_SOUTH, ITALY_WEST,
|
||
SARDINIA, CORSICA, MALLORCA, MENORCA, IBIZA, NORTH_AFRICA,
|
||
]
|
||
|
||
|
||
def _pip(px, py, poly):
|
||
n = len(poly)
|
||
inside = False
|
||
j = n - 1
|
||
for i in range(n):
|
||
xi, yi = poly[i]
|
||
xj, yj = poly[j]
|
||
if ((yi > py) != (yj > py)) and (px < (xj - xi) * (py - yi) / (yj - yi) + xi):
|
||
inside = not inside
|
||
j = i
|
||
return inside
|
||
|
||
|
||
def is_land(lon, lat):
|
||
for poly in LAND_POLYGONS:
|
||
if _pip(lon, lat, poly):
|
||
return True
|
||
return False
|
||
|
||
|
||
# ── Current field (NW Mediterranean / Catalan coast, realistic) ───────────────
|
||
|
||
# Typical surface currents: Northern Current ~0.1–0.25 m/s along coast (SW), weaker offshore
|
||
BASE_SPEED_MS = 0.18 # m/s, main flow magnitude
|
||
DIFFUSION_MS = 0.025 # turbulent fluctuation (m/s)
|
||
|
||
def current_field(lon, lat, rng, t_frac):
|
||
# Main flow: SW along Catalan coast (negative v = south, slight negative u = west)
|
||
angle_deg = 215 + 8 * math.sin(2 * math.pi * t_frac) # slow seasonal wobble
|
||
rad = math.radians(angle_deg)
|
||
u_base = BASE_SPEED_MS * math.cos(rad)
|
||
v_base = BASE_SPEED_MS * math.sin(rad)
|
||
|
||
# Coastal weakening: reduce speed when close to Iberian coast (lon < 2.5, lat 40–42.5)
|
||
dist_coast = 1.0
|
||
if lon < 3.2 and 40.4 < lat < 42.5:
|
||
dist_coast = max(0.15, (lon - 0.5) * 0.4 + (42.2 - lat) * 0.15)
|
||
factor = min(1.0, dist_coast)
|
||
u_main = u_base * factor
|
||
v_main = v_base * factor
|
||
|
||
# Weak cyclonic gyre east of Barcelona (3.5, 41.5), radius ~0.8 deg
|
||
cx, cy = 3.6, 41.4
|
||
dx, dy = lon - cx, lat - cy
|
||
r_deg = math.sqrt(dx * dx + dy * dy)
|
||
if r_deg < 0.9 and r_deg > 0.05:
|
||
gyre_speed = 0.06 * math.exp(-r_deg / 0.5) # m/s
|
||
# tangential: (-dy, dx) for anticlockwise
|
||
u_gyre = -gyre_speed * (dy / r_deg)
|
||
v_gyre = gyre_speed * (dx / r_deg)
|
||
else:
|
||
u_gyre, v_gyre = 0.0, 0.0
|
||
|
||
# Small-scale turbulence (Gaussian, uncorrelated)
|
||
u_noise = rng.normal(0, DIFFUSION_MS)
|
||
v_noise = rng.normal(0, DIFFUSION_MS)
|
||
|
||
return (u_main + u_gyre + u_noise, v_main + v_gyre + v_noise)
|
||
|
||
|
||
def advect_rk4(lon, lat, rng, t_frac):
|
||
cos_lat = math.cos(math.radians(lat))
|
||
dpm_lon = DEG_PER_M_LAT / max(cos_lat, 1e-6)
|
||
|
||
def vel(lo, la):
|
||
u, v = current_field(lo, la, rng, t_frac)
|
||
return u * DT_SECONDS * dpm_lon, v * DT_SECONDS * DEG_PER_M_LAT
|
||
|
||
k1 = vel(lon, lat)
|
||
k2 = vel(lon + k1[0] / 2, lat + k1[1] / 2)
|
||
k3 = vel(lon + k2[0] / 2, lat + k2[1] / 2)
|
||
k4 = vel(lon + k3[0], lat + k3[1])
|
||
|
||
new_lon = lon + (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6
|
||
new_lat = lat + (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]) / 6
|
||
u_f, v_f = current_field(new_lon, new_lat, rng, t_frac)
|
||
return new_lon, new_lat, u_f, v_f
|
||
|
||
|
||
# ── Accumulation grid ────────────────────────────────────────────────────────
|
||
|
||
ACCUM_RES = 0.1
|
||
|
||
def build_accumulation(b_lons, b_lats):
|
||
if len(b_lons) == 0:
|
||
return []
|
||
cells: dict[tuple[int, int], int] = {}
|
||
for lo, la in zip(b_lons, b_lats):
|
||
key = (round(lo / ACCUM_RES), round(la / ACCUM_RES))
|
||
cells[key] = cells.get(key, 0) + 1
|
||
return [[k[0] * ACCUM_RES, k[1] * ACCUM_RES, v] for k, v in cells.items()]
|
||
|
||
|
||
# ── Main ─────────────────────────────────────────────────────────────────────
|
||
|
||
def generate():
|
||
os.makedirs(DATA_DIR, exist_ok=True)
|
||
rng = np.random.default_rng(42)
|
||
|
||
init_lons, init_lats, origins = [], [], []
|
||
for si, (slon, slat) in enumerate(SEED_POINTS):
|
||
for _ in range(PARTICLES_PER_SEED):
|
||
lon = slon + rng.uniform(-JITTER_DEG, JITTER_DEG)
|
||
lat = slat + rng.uniform(-JITTER_DEG, JITTER_DEG)
|
||
init_lons.append(lon)
|
||
init_lats.append(lat)
|
||
origins.append(si)
|
||
|
||
n = len(init_lons)
|
||
lons = np.array(init_lons, dtype=float)
|
||
lats = np.array(init_lats, dtype=float)
|
||
|
||
# Ensure every particle starts in water: push seaward (east/northeast) if on land
|
||
for i in range(n):
|
||
for _ in range(50):
|
||
if not is_land(lons[i], lats[i]):
|
||
break
|
||
lons[i] += rng.uniform(0, OFFSHORE_DEG)
|
||
lats[i] += rng.uniform(0, OFFSHORE_DEG * 0.5)
|
||
init_lons = lons.tolist()
|
||
init_lats = lats.tolist()
|
||
|
||
# Staggered release over first 24h for more realistic dispersion
|
||
release_steps = rng.integers(0, min(24, NUM_STEPS), size=n).tolist()
|
||
init_lons_arr = np.array(init_lons, dtype=float)
|
||
init_lats_arr = np.array(init_lats, dtype=float)
|
||
|
||
beached = np.zeros(n, dtype=bool)
|
||
|
||
# Collect ALL frames: frames[step] = [[lon, lat, u, v, beached], ...]
|
||
all_frames = []
|
||
all_accum = []
|
||
|
||
for step in range(NUM_STEPS):
|
||
t_frac = step / NUM_STEPS
|
||
frame = np.zeros((n, 5))
|
||
|
||
for i in range(n):
|
||
if step < release_steps[i]:
|
||
frame[i] = [init_lons_arr[i], init_lats_arr[i], 0.0, 0.0, 0]
|
||
continue
|
||
if step == release_steps[i]:
|
||
lons[i], lats[i] = init_lons_arr[i], init_lats_arr[i]
|
||
|
||
if not beached[i]:
|
||
nl, nla, u, v = advect_rk4(float(lons[i]), float(lats[i]), rng, t_frac)
|
||
if is_land(nl, nla):
|
||
beached[i] = True
|
||
u, v = 0.0, 0.0
|
||
else:
|
||
lons[i], lats[i] = nl, nla
|
||
frame[i] = [lons[i], lats[i], u, v, 0]
|
||
continue
|
||
|
||
frame[i] = [lons[i], lats[i], 0.0, 0.0, 1 if beached[i] else 0]
|
||
|
||
all_frames.append(np.round(frame, 5).tolist())
|
||
|
||
b_mask = beached
|
||
accum = build_accumulation(lons[b_mask].tolist(), lats[b_mask].tolist())
|
||
all_accum.append(accum)
|
||
|
||
act = int(np.sum(~beached))
|
||
bch = int(np.sum(beached))
|
||
print(f" Step {step:>2d}/{NUM_STEPS - 1} | active {act:>5d} beached {bch:>5d}")
|
||
|
||
# Write blocks
|
||
num_blocks = NUM_STEPS // BLOCK_SIZE
|
||
for b in range(num_blocks):
|
||
s0 = b * BLOCK_SIZE
|
||
s1 = s0 + BLOCK_SIZE
|
||
block = {
|
||
"block": b,
|
||
"step_start": s0,
|
||
"step_end": s1,
|
||
"frames": all_frames[s0:s1],
|
||
"accumulation": all_accum[s0:s1],
|
||
}
|
||
path = os.path.join(DATA_DIR, f"block_{b}.json")
|
||
with open(path, "w") as f:
|
||
json.dump(block, f, separators=(",", ":"))
|
||
size_mb = os.path.getsize(path) / 1_048_576
|
||
print(f" Block {b} ({s0}h-{s1 - 1}h) → {path} [{size_mb:.1f} MB]")
|
||
|
||
# Metadata (release_steps: step at which each particle is released, 0..num_steps-1)
|
||
meta = {
|
||
"num_particles": n,
|
||
"num_steps": NUM_STEPS,
|
||
"num_blocks": num_blocks,
|
||
"block_size": BLOCK_SIZE,
|
||
"dt_hours": 1,
|
||
"seed_names": SEED_NAMES,
|
||
"origins": origins,
|
||
"release_steps": release_steps,
|
||
}
|
||
meta_path = os.path.join(DATA_DIR, "metadata.json")
|
||
with open(meta_path, "w") as f:
|
||
json.dump(meta, f, separators=(",", ":"))
|
||
print(f" Metadata → {meta_path}")
|
||
|
||
|
||
if __name__ == "__main__":
|
||
n = len(SEED_POINTS) * PARTICLES_PER_SEED
|
||
print(f"Simulating {n} particles for {NUM_STEPS} hours …")
|
||
generate()
|
||
print("Done.")
|