| … | … | @@ -340,7 +340,57 @@ def to_png(rgb): |
| 340 | 340 | return buf.getvalue() |
| 341 | 341 | |
| 342 | 342 | |
| 343 | | def colourise_array(h, mpp, view, vmin, vmax): |
| 343 | # HYPSOMETRIC / BATHYMETRIC TINTING. |
| 344 | # |
| 345 | # A single ramp across the whole range is right for consistency between tiles |
| 346 | # and wrong for contrast within one. Seafloor at -2200..-1000 m occupies about |
| 347 | # 15% of a +/-4000 m ramp, so it arrives on screen almost flat -- which read as |
| 348 | # "the model has no detail down there" when the measurement says it has 5x less |
| 349 | # detail, not none (land 8.41 vs seafloor 1.65, and the drop lands exactly on |
| 350 | # the 0 m contour, not on any tile boundary). |
| 351 | # |
| 352 | # Real maps solve this by tinting land and sea on SEPARATE scales. Both then get |
| 353 | # their full range, and sea level is a hard, readable edge instead of the middle |
| 354 | # of a gradient. The scales are absolute, so tiles still agree with each other. |
| 355 | LAND_STOPS = [ |
| 356 | (0, (0.29, 0.48, 0.25)), # coastal green |
| 357 | (300, (0.49, 0.60, 0.29)), |
| 358 | (800, (0.72, 0.63, 0.35)), # tan |
| 359 | (1600, (0.66, 0.46, 0.25)), # brown |
| 360 | (2600, (0.55, 0.44, 0.38)), # bare rock |
| 361 | (3600, (0.91, 0.91, 0.94)), # snow |
| 362 | (6500, (1.00, 1.00, 1.00)), |
| 363 | ] |
| 364 | SEA_STOPS = [ |
| 365 | (0, (0.47, 0.68, 0.74)), # shore |
| 366 | (-200, (0.31, 0.53, 0.63)), # shelf |
| 367 | (-1000, (0.19, 0.37, 0.52)), |
| 368 | (-2500, (0.11, 0.24, 0.38)), |
| 369 | (-4500, (0.05, 0.13, 0.24)), |
| 370 | (-9000, (0.02, 0.06, 0.14)), # trench |
| 371 | ] |
| 372 | |
| 373 | |
| 374 | def _ramp(h, stops): |
| 375 | xs = np.array([s[0] for s in stops], dtype=np.float64) |
| 376 | cs = np.array([s[1] for s in stops], dtype=np.float64) |
| 377 | if xs[0] > xs[-1]: |
| 378 | xs = xs[::-1]; cs = cs[::-1] |
| 379 | out = np.empty(h.shape + (3,)) |
| 380 | for k in range(3): |
| 381 | out[..., k] = np.interp(h, xs, cs[:, k]) |
| 382 | return out |
| 383 | |
| 384 | |
| 385 | def hypsometric(h): |
| 386 | """Land and sea tinted on separate scales, so each gets its full range.""" |
| 387 | rgb = _ramp(h, LAND_STOPS) |
| 388 | sea = _ramp(h, SEA_STOPS) |
| 389 | m = (h <= 0)[..., None] |
| 390 | return np.where(m, sea, rgb) |
| 391 | |
| 392 | |
| 393 | def colourise_array(h, mpp, view, vmin, vmax, exag=1.0): |
| 344 | 394 | """The colour computation without the PNG encode, so a caller can crop an |
| 345 | 395 | apron off before encoding.""" |
| 346 | 396 | finite = h[np.isfinite(h)] |
| … | … | @@ -354,8 +404,14 @@ def colourise_array(h, mpp, view, vmin, vmax): |
| 354 | 404 | elif view == "height": |
| 355 | 405 | t = np.nan_to_num(np.clip((h - lo) / (hi - lo), 0, 1)) |
| 356 | 406 | rgb = np.stack([0.15 + 0.85 * t, 0.25 + 0.6 * t, 0.35 + 0.3 * (1 - t)], axis=-1) |
| 407 | elif view == "hypsometric": |
| 408 | # Exaggeration affects the SHADING only, never the colour -- the tint |
| 409 | # keeps telling the truth about elevation while the relief becomes |
| 410 | # legible. Submarine relief needs 5-10x before it reads at all. |
| 411 | sh = hillshade(h, mpp, exaggeration=exag) |
| 412 | rgb = hypsometric(h) * (0.30 + 0.70 * sh[..., None]) |
| 357 | 413 | else: |
| 358 | | sh = hillshade(h, mpp) |
| 414 | sh = hillshade(h, mpp, exaggeration=exag) |
| 359 | 415 | if view == "shaded-height": |
| 360 | 416 | t = np.nan_to_num(np.clip((h - lo) / (hi - lo), 0, 1)) |
| 361 | 417 | base = np.stack([0.20 + 0.75 * t, 0.30 + 0.55 * t, |
| … | … | @@ -574,8 +630,9 @@ def api_tile(): |
| 574 | 630 | # Shade WITH the apron, then crop it. Shading the bare tile is what put |
| 575 | 631 | # a seam on every join even though the elevation is continuous. |
| 576 | 632 | full = a.astype(np.float64) |
| 577 | | shaded = colourise_array(full, 30.0, request.args.get("view", "shaded-height"), |
| 578 | | vmin, vmax) |
| 633 | shaded = colourise_array(full, 30.0, request.args.get("view", "hypsometric"), |
| 634 | vmin, vmax, |
| 635 | float(request.args.get("exag", 1.0))) |
| 579 | 636 | ap = tiles().APRON |
| 580 | 637 | png = to_png(shaded[ap:-ap, ap:-ap] if ap else shaded) |
| 581 | 638 | r = Response(png, mimetype="image/png") |
| … | … | @@ -609,11 +666,12 @@ def api_coarse(): |
| 609 | 666 | elev = coarse_to_metres(a[..., 0].astype(np.float64)) |
| 610 | 667 | vmin = float(request.args.get("vmin", -4000)) |
| 611 | 668 | vmax = float(request.args.get("vmax", 4000)) |
| 612 | | view = request.args.get("view", "shaded-height") |
| 669 | view = request.args.get("view", "hypsometric") |
| 613 | 670 | if view == "terrain-class": |
| 614 | 671 | rgb = colourise_terrain_class(elev) |
| 615 | 672 | else: |
| 616 | | rgb = colourise_array(elev, 7680.0, view, vmin, vmax) |
| 673 | rgb = colourise_array(elev, 7680.0, view, vmin, vmax, |
| 674 | float(request.args.get("exag", 1.0))) |
| 617 | 675 | resp = Response(to_png(rgb), mimetype="image/png") |
| 618 | 676 | resp.headers["Cache-Control"] = "public, max-age=3600" |
| 619 | 677 | resp.headers["X-Height"] = str(h); resp.headers["X-Width"] = str(w) |