1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
|
"""`mapgen.py --world DIR new-world`: a random starting world (sketch/*.png continents, config/tectonics.toml plates),
so a build runs without drawing anything. Used for the moons; edit or redraw afterwards."""
from __future__ import annotations
from pathlib import Path
import numpy as np
from PIL import Image
from scipy import ndimage
from .noise import fbm, ridged
from .sphere import latlon_to_xyz
W, H = 2000, 1000
def _spread(rng, n, min_deg, avoid=(), max_lat=70.0, tries=4000):
pts = list(avoid)
out = []
for _ in range(tries):
if len(out) == n:
break
lat = np.degrees(np.arcsin(rng.uniform(np.sin(np.radians(-max_lat)), np.sin(np.radians(max_lat)))))
lon = rng.uniform(-180, 180)
p = latlon_to_xyz(lat, lon)
if all(np.degrees(np.arccos(np.clip(p @ latlon_to_xyz(*q), -1, 1))) >= min_deg for q in pts):
pts.append((lat, lon))
out.append((float(lat), float(lon)))
return out
def _centroid(mask, LAT, LON):
w = np.cos(np.radians(LAT)) * mask
p = (latlon_to_xyz(LAT, LON) * w[..., None]).sum(axis=(0, 1))
p /= np.linalg.norm(p)
return float(np.degrees(np.arcsin(p[2]))), float(np.degrees(np.arctan2(p[1], p[0])))
def make(root: Path, seed: int, continents: int = 6, land_share: float = 0.27, toward=None, ridge: float = 0.86,
force: bool = False) -> None:
root = Path(root)
sk, cfg = root / "sketch", root / "config"
targets = [sk / "land.png", cfg / "tectonics.toml"]
if not force and any(p.exists() for p in targets):
raise SystemExit("new-world: sketch/ or config/tectonics.toml already exist (use --force to overwrite)")
rng = np.random.default_rng(seed)
w, h = W // 2, H // 2
LAT, LON = np.meshgrid(90.0 - (np.arange(h) + 0.5) * 180.0 / h, -180.0 + (np.arange(w) + 0.5) * 360.0 / w,
indexing="ij")
xyz = latlon_to_xyz(LAT, LON).reshape(-1, 3)
centres = _spread(rng, continents, 35.0, max_lat=55.0)
field = np.zeros(len(xyz))
for lat, lon in centres:
r = np.radians(rng.uniform(18.0, 34.0))
d = np.arccos(np.clip(xyz @ latlon_to_xyz(lat, lon), -1, 1))
field = np.maximum(field, np.clip(1.0 - d / r, 0.0, None) * rng.uniform(0.8, 1.2))
field += 0.55 * fbm(xyz, seed + 7, 6, 2.2)
if toward is not None:
field += 1.2 * (xyz @ latlon_to_xyz(*toward))
area = np.cos(np.radians(LAT)).ravel()
order = np.argsort(-field)
cut = field[order][np.searchsorted(np.cumsum(area[order]) / area.sum(), land_share)]
land = (field > cut).reshape(h, w)
mount = (ridged(xyz, seed + 11, 5, 3.0) > ridge).reshape(h, w) & ndimage.binary_erosion(land, iterations=6)
def save(name, m):
Image.fromarray((m * 255).astype(np.uint8), "L").resize((W, H), Image.BILINEAR).point(
lambda v: 255 if v > 127 else 0).save(sk / f"{name}.png")
sk.mkdir(parents=True, exist_ok=True)
cfg.mkdir(parents=True, exist_ok=True)
save("land", land)
save("mountains", mount)
for name in ("desert", "rainforest", "trench"):
save(name, np.zeros_like(land))
lab, n = ndimage.label(land)
for a, b in zip(lab[:, 0], lab[:, -1]):
if a and b and a != b:
lab[lab == b] = a
sizes = sorted(((int((lab == i).sum()), i) for i in np.unique(lab) if i), reverse=True)
plates, seeds = [], []
for k, (px, i) in enumerate(sizes[:continents]):
m = lab == i
if px > 0.035 * m.size: # a big landmass: two plates, a collision belt between
lat, lon = _centroid(m, LAT, LON)
ang = rng.uniform(0, np.pi)
side = (np.sin(ang) * (LAT - lat) + np.cos(ang) * ((LON - lon + 180) % 360 - 180) * np.cos(np.radians(lat))) > 0
parts = [m & side, m & ~side]
else:
parts = [m]
for j, part in enumerate(parts):
if part.sum() < 50:
continue
lat, lon = _centroid(part, LAT, LON)
seeds.append((lat, lon))
plates.append((f"continent-{k + 1}{'ab'[j] if len(parts) > 1 else ''}", lat, lon, "continental",
rng.uniform(0, 360), rng.uniform(2.0, 5.0)))
for k, (lat, lon) in enumerate(_spread(rng, 14 - len(plates) // 2, 28.0, seeds, max_lat=85.0)):
plates.append((f"ocean-{k + 1}", lat, lon, "oceanic", rng.uniform(0, 360), rng.uniform(3.0, 8.0)))
ocean = np.argwhere(~ndimage.binary_dilation(land, iterations=10))
def sea_point():
y, x = ocean[rng.integers(len(ocean))]
return float(LAT[y, x]), float(LON[y, x])
f = lambda v: f"[{v[0]:.1f}, {v[1]:.1f}]"
t = [f"# Generated by `mapgen.py new-world --seed {seed}`. Edit freely (see README.md → Tectonics).",
"# seed = [lat, lon] where the plate grows from; motion = [azimuth° clockwise from north, speed cm/yr].", ""]
for pid, lat, lon, kind, az, sp in plates:
t += ["[[plate]]", f'id = "{pid}"', f"seed = {f((lat, lon))}", f'kind = "{kind}"', f"motion = [{az:.0f}.0, {sp:.1f}]", ""]
if len(ocean):
t += ["[[hotspot]]", 'name = "island-chain"', f"center = {f(sea_point())}", "length_km = 600.0", ""]
(cfg / "tectonics.toml").write_text("\n".join(t))
print(f"new-world: {len(plates)} plates, sketch/ and config/tectonics.toml written (seed {seed})")
|
