第14章 - 动态渲染重采样与色彩

同一份 COG，同一张瓦片坐标，只把 rescale 从 0,10000 改成 0,3000，画面从灰蒙蒙变成层次分明；把 resampling 从 nearest 改成 average，30 cm 影像的田埂边界就糊了半米。渲染参数不是”美化选项”，它们直接改变量测结果。

这一章把渲染链拆开：数据从 GDAL 出来之后，经过哪几步、每一步在哪一层做、参数怎么传。顺序搞清楚了，调色和算法就不会互相打架。

14.1 渲染链的五个阶段

一次瓦片渲染的完整数据流，按执行顺序：

1. 读取     GDAL 按 TMS 窗口读块（受 GDAL_CACHEMAX 影响）
                ↓ 原 dtype 数组 + mask
2. 重采样   RasterIO resampling（同 CRS 缩放）
            WarpKernel reproject（跨 CRS 重投影）
                ↓ 目标分辨率的数组
3. 算法     algorithm（hillshade / normalizedIndex / contours ...）
                ↓ 可能改变波段数与 dtype
4. 拉伸     rescale（逐波段线性映射到输出 dtype 范围）
            color_formula（gamma / sigmoidal / saturation）
            colormap（查表把单波段变 RGB）
                ↓ uint8 或 uint16
5. 编码     PNG / JPEG / WEBP / TIFF

这个顺序有两个后果，值得先记住：

• rescale 作用在算法之后。写 algorithm=normalizedIndex&rescale=-1,1 时，NDVI 先算出来，再被拉伸到 0～255；写 rescale=0,10000 想先把反射率拉一下再算 NDVI，是做不到的，因为拉伸在算法后面。
• colormap 和 color_formula 都在 rescale 之后。三个都用时，先线性拉伸、再按公式调色、最后查色带。色带会覆盖 RGB，所以单波段加色带时 color_formula 基本没有意义。

render_image() 的实现把这三步写死了：先 image.rescale(rescale)，再 image.apply_color_formula(color_formula)，最后 apply_cmap。要改顺序只能换 render_func（TilerFactory 的构造参数），这是 TiTiler 留的正规扩展点。

14.2 rescale：最容易踩的四个细节

rescale 的语义是”把 [min, max] 区间的值线性映射到输出 dtype 的完整范围”。uint8 输出映射到 0～255，uint16 输出映射到 0～65535。它接受可重复参数，每个对应一个波段：

# 三个波段分别拉伸
curl -s -o out.png "http://127.0.0.1:8000/cog/tiles/WebMercatorQuad/12/3372/1552.png\
  ?url=https://example.com/s2.tif&bidx=1&bidx=2&bidx=3\
  &rescale=0,3000&rescale=0,3000&rescale=0,2000"

四个细节：

数量不匹配时全部同值。 rescale_image() 里有一句兜底：if len(in_range) != nbands: in_range = ((in_range[0]),) * nbands。给三波段只写一组 rescale=0,3000，三个波段都用 0～3000。想要”红波段 0～3000、近红外 0～5000”，必须写两组。

不给 rescale 时按 dtype 全范围拉伸。 这是官方文档明确的行为：8 位数据按 0～255、16 位按 0～65536。对于反射率乘 10000 的 Sentinel-2 数据，全范围拉伸等于把 0～10000 压到 0～255，而真实地物只占 0～4000，结果就是灰白一片。DEM 更极端，-100～8000 m 压进 0～255，平原地区几乎没有对比度。

rescale 的校验只认两段数字。 validate_rescale() 的做法是把每段里的 [ 和 ] 剥掉再按逗号切，切的段数必须正好是 2。所以 rescale=0,3000 和 rescale=[[0,3000]] 都能过（后者剥括号后还是两段），而 rescale=[[0,3000],[0,3000]] 会被切成四段、直接 422。想给两个波段设区间，只能用重复参数，不能写在一个值里。

unscale 会改变数值量级。 开了 unscale=true 之后 DN 变成物理量，0～10000 可能变成 0～1，原来的 rescale=0,3000 就完全失效。这类问题的排查顺序永远是：/info 看 scales/offsets → 决定是否 unscale → 用 /statistics 拿真实区间 → 再写 rescale。

自动拉伸的实用做法是用统计分位数，而不是手拍数字：

import requests

BASE = "http://127.0.0.1:8000/cog"
SRC = "https://example.com/s2.tif"

stats = requests.get(
    f"{BASE}/statistics",
    params={"url": SRC, "bidx": 1, "p": [2, 98]},
    timeout=60,
).json()

band = stats["1"]
lo, hi = band["percentile_2"], band["percentile_98"]
print(f"rescale={lo:.0f},{hi:.0f}")
# 直接用这个字符串去请求瓦片

/statistics 支持 p 参数多值，返回体里会带 percentile_2、percentile_98 这样的键，BandStatistics 模型的 extra 允许额外字段，所以分位数是动态加进去的。

14.3 colormap：单波段上色的两条路

色带只对单波段有意义（对 RGB 三波段用色带，会先被压成一个波段再查表，等于丢信息）。TiTiler 支持两种写法：

内置色带，用 colormap_name。可取值是 rio-tiler 的 cmap.list()，即 GET /colorMaps 返回的键集合，来源是 Matplotlib 的色带加一批遥感常用色带，全部支持加 _r 后缀反序。

# 看有哪些
curl -s "http://127.0.0.1:8000/colorMaps" | jq -r ".colormaps | keys[]" | head -20

# 用 viridis 给 DEM 上色
curl -s -o dem.png "http://127.0.0.1:8000/cog/tiles/WebMercatorQuad/11/1686/776.png\
  ?url=https://example.com/dem.tif&bidx=1&rescale=0,3000&colormap_name=viridis"

自定义色带，用 colormap 传 JSON 字符串。键是像元值（字符串形式），值是十六进制颜色或 [r,g,b,a] 数组：

import json
import urllib.parse

import requests

cmap = {
    "0": "#000004",
    "1000": "#38598c",
    "2000": "#4ec36b",
    "3000": "#f6d746",
    "4000": "#ffffff",
}
params = {
    "url": "https://example.com/dem.tif",
    "bidx": 1,
    "rescale": "0,4000",
    "colormap": json.dumps(cmap, separators=(",", ":")),
}
r = requests.get(
    "http://127.0.0.1:8000/cog/tiles/WebMercatorQuad/11/1686/776.png",
    params=params,
    timeout=30,
)
r.raise_for_status()
print(len(r.content), r.headers["content-type"])

colormap 的解析在 create_colormap_dependency 里：json.loads 之后用 object_hook 把每个值交给 parse_color，序列（list）形式也支持，会被转成区间色带。JSON 解析失败直接抛 400，信息是 Could not parse the colormap value.，不会告诉你哪个键有问题。

两个行为的边界：

• 色带是离散的。像元值 1500 落在 1000 和 2000 之间时，取哪个颜色取决于 rio-tiler 的查表实现（默认按左边界取整），想要连续渐变要用序列形式的区间色带。
• 色带的 alpha 会参与掩膜。自定义色带里把某个值写成 [0,0,0,0]，该值的像元就变透明，这是”把 nodata 藏掉”的常见手法，比调 nodata 参数更直接。

14.4 color_formula：rio-color 的三个算子

color_formula 走的是 Mapbox 的 rio-color 语法，支持三个算子：gamma、sigmoidal、saturation，可以逗号串起来依次执行。

一条实际可用的组合（官方文档给的例子）：

curl -s -o out.jpg "http://127.0.0.1:8000/cog/tiles/WebMercatorQuad/14/13491/6210.jpg\
  ?url=https://example.com/s2.tif&bidx=1&bidx=2&bidx=3\
  &rescale=0,3000&rescale=0,3000&rescale=0,3000\
  --data-urlencode "color_formula=gamma rg 1.3, sigmoidal rgb 22 0.1, saturation 1.5" -G"

注意 --data-urlencode 加 -G：公式里有空格，不编码会在 shell 层被拆成多个参数，服务端收到的 color_formula 只剩 gamma。

color_formula 有一个约束来自校验器 validate_color_formula：必须匹配它定义的正则，非法字符直接 400（InvalidColorFormat）。它不做数值范围检查，写 gamma rgb 99 也能过校验，但结果会溢出成纯白。

14.5 algorithm：改变波段数的后处理

算法是唯一能改变”波段数”和”dtype”的环节。内置 18 个，注册表在 algorithm/__init__.py 的 default_algorithms：

两个最常用的用法：

Hillshade 叠在山体阴影上，注意 buffer 默认是 3，算法会在内部裁掉这圈像素，输出仍是 256×256：

curl -s -o hs.png "http://127.0.0.1:8000/cog/tiles/WebMercatorQuad/12/3372/1552.png\
  ?url=https://example.com/dem.tif\
  --data-urlencode 'algorithm_params={"azimuth":315,"angle_altitude":45,"z_exaggeration":1.5}' \
  -G -d algorithm=hillshade -d rescale=0,255"

NDVI 用 normalizedIndex，它把前两个波段当 (b1, b2) 算 (b2-b1)/(b2+b1)：

curl -s -o ndvi.png "http://127.0.0.1:8000/cog/tiles/WebMercatorQuad/12/3372/1552.png\
  ?url=https://example.com/s2.tif&bidx=4&bidx=8\
  -d algorithm=normalizedIndex -d rescale=-1,1 -d colormap_name=rdylgn"

bidx=4&bidx=8 的顺序决定公式：NormalizedIndex.__call__ 里 b1 = array[0]、b2 = array[1]，算的是 (b2-b1)/(b2+b1)。所以要算 NDVI（NIR-Red)/(NIR+Red)，必须红在前、近红外在后。顺序写反会得到负 NDVI，画面整体反色。

algorithm_params 是 JSON 字符串，解析失败会抛 422；参数名不匹配 BaseAlgorithm 的字段会抛 400 并带 Pydantic 的字段列表。特别注意 NormalizedIndex 不接受任何参数——写了 {"b1": 3} 会报 extra fields not permitted。

14.6 重采样：nearest 与 average 的取舍

resampling 的 9 个合法值来自 rio_tiler.types.RIOResampling：nearest、bilinear、cubic、cubic_spline、lanczos、average、mode、gauss、rms。reproject 的 14 个值多出 max、min、med、q1、q3、sum。

选择规则不是”越平滑越好”，而是看数据的统计含义：

average 在降采样时是唯一保持”总量守恒”的选择。nearest 在 10 m 到 30 m 的降采样里会直接丢掉 8/9 的像元，植被指数会出现明显偏差；bilinear 会引入负值（在反射率上表现为异常黑点）。

一个实测差异：30 m DEM 降到 z=10（约 152 m/瓦片像素），nearest 与 average 的坡度图在陡崖处相差 3～5 度。如果下游拿坡度做工程判断，这个差异不能忽略。

14.7 常用渲染预设

把参数固化成几套预设，前端只传预设名，是控制缓存与效果最省事的办法。下面五套覆盖了大部分遥感出图需求，参数都可以直接抄：

同一份数据同时提供”真彩色”和”NDVI”两个预设，前端切换只改 query，不动代码；CDN 侧也能按完整 URL 缓存两份。若把参数开放给用户自由组合，缓存命中率会随组合数线性下降——第 17 章算过这笔账。

一个把预设写进服务端的做法，用 colormap_dependency 之外的自定义依赖把预设名翻译成参数：

from fastapi import HTTPException, Query

PRESETS = {
    "truecolor": {"bidx": [1, 2, 3], "rescale": ["0,3000"] * 3},
    "ndvi": {"bidx": [4, 8], "algorithm": "normalizedIndex", "rescale": ["-1,1"]},
}


def PresetParams(preset: str = Query("truecolor")) -> dict:
    """把预设名翻译成一组固定的渲染参数。"""
    if preset not in PRESETS:
        raise HTTPException(400, f"unknown preset: {preset}")
    return PRESETS[preset]

这个依赖返回的字典可以直接展开进 render_func 的调用参数，也可以在网关层改写成 query，两种做法都比让前端拼参数可控。

14.8 排错与踩坑

整幅图纯白。 三个原因按概率排序：rescale 的上界小于数据实际值（如 Sentinel-2 用了 0～3000 但数据是 0～10000 量级）；开了 unscale 又沿用旧区间；或者 color_formula 的 gamma 过大。排查顺序是先去掉 color_formula，再看 /statistics 的 p98。

整幅图纯黑。 大概率是掩膜问题：nodata_type 为 None 时没有掩膜，透明区按 0 渲染；或者 return_mask=false 让 PNG 丢了 alpha，黑色背景直接显示出来。

色带只显示一种颜色。 单波段数据加色带时，若 rescale 区间远大于实际值域，所有像元都落在色带的同一个端点。先不加色带出灰度图确认有层次，再加色带。

algorithm=normalizedIndex 报 422。 这个算法要求输入正好 2 个波段。三波段加 bidx=1&bidx=2&bidx=3 会失败，多波段会抛 400 并提示波段数不符。

hillshade 边缘出现黑边。 HillShade 默认 buffer=3，它会在内部把边缘裁掉；如果同时传了 buffer 查询参数（TileParams.buffer），两个 buffer 叠加，边缘就会出现未填充区。算法自带 buffer 时不要再加查询参数。

rescale 传二维数组报 422。 validate_rescale 要求每段切出正好两段，[[0,3000],[0,3000]] 会切成四段。多波段的写法是重复参数 rescale=0,3000&rescale=0,3000，而不是一个二维数组。

14.9 生产建议

• 把渲染参数固化成 tilejson 的固定 query。前端只拿一个 tilejson，不要在客户端拼 rescale 和 colormap，否则参数一改所有前端都要发版。
• 同一数据集准备 2～3 套”渲染预设”，比如 truecolor、falsecolor、ndvi，各自对应一条完整的 query 串，存在配置里。这比让用户自由组合参数可控得多，也能被 CDN 按 URL 直接命中。
• rescale 用 p2/p98，不要用 min/max。单景影像里一个镜面反射像元就能把 max 拉到 60000，整幅图随之变暗。
• 算法端点单独限流。hillshade 与 contours 的 CPU 开销是普通出图的 2～3 倍（numpy gradient 与模运算），和普通瓦片混在同一池里会拖慢整体 P95。
• 不要在生产暴露任意 algorithm_params。参数没有白名单，z_exaggeration 给到 1e6 会让 numpy 溢出并返回 500。如果对外开放，用固定的几组参数在网关层做映射。

14.10 本章要点

下一章处理”多景”：把几十上百个 COG 拼成一张无缝图，MosaicJSON 的 tiles 字段怎么组织，cogeo-mosaic 的六种 backend 各适合什么场景。