【无标题】

整合自 20260629本机操作记录.md、20260630dem处理记录.md 及 2026-06-30 编码/COG 修复结论。
案例数据：Copernicus DEM GLO-30 · 鞍山市 · 目标产物供 MapLibre raster-dem + encoding: 'terrarium' 使用。

1. 流程总览

anshan-dem.tif
  → [1] NoData 修复（dem_clean.tif）
  → [2] 重投影 + 多阶段平滑 → dem_final.tif（15 m）
  → [3] Terrarium RGB 编码（gdal_calc，INTERLEAVE=BAND）
  → [4] COG（OVERVIEW_RESAMPLING=NEAREST）
  → [5] 发布 MinIO / 可选 gdal2tiles / PMTiles
  → [6] MapLibre Demo 验证

关键约定（踩坑总结）：

Terrarium 解码公式（MapLibre / 校验脚本一致）：

height = R × 256 + G + B / 256 - 32768

编码公式（gdal_calc 实现）：

val = clip(elevation, -32768, 32767) + 32768
R = floor(val / 256)
G = floor(val) % 256
B = floor((val - floor(val)) × 256)

2. 环境准备

conda create -n gis python=3.10 -y
conda activate gis
conda install -c conda-forge gdal -y

• 需要 GDAL ≥ 3.2（gdal_calc.py 多 --calc 输出多波段）。
• 不必依赖 rio-rgbify 产出 Terrarium GeoTIFF；PyPI 0.4.0 不支持 -e terrarium，且 -b -32768 -i 1 会产出错误 RGB。

3. 预处理（Float32 高程，EPSG:3857）

输入：anshan-dem.tif（Copernicus 30 m，鞍山）
终点：dem_final.tif（Float32，EPSG:3857，像元 15 m）。多阶段 cubicspline 降采样（30→20→15 m）为生产推荐，详见 预处理数据--平滑.md §4.5。

3.1 NoData 修复

-md 10 把无数据填充为10像素内的有效数据进行填充

gdal_fillnodata.py anshan-dem.tif dem_clean.tif -md 10

3.2 重投影

gdalwarp -t_srs EPSG:3857 -r cubicspline dem_clean.tif dem_3857.tif

3.3 多阶段平滑重采样（30→20→15 m）

gdalwarp -r cubicspline -tr 30 30 dem_3857.tif dem_s1.tif
gdalwarp -r cubicspline -tr 20 20 dem_s1.tif  dem_s2.tif
gdalwarp -r cubicspline -tr 15 15 dem_s2.tif  dem_final.tif

鞍山实测：dem_final.tif 高程约 -293 – 659 m。

3.4 编码前检查

gdalinfo dem_final.tif
# 期望：Band 1 Float32；CRS EPSG:3857；像元约 15 m

4. Terrarium RGB 编码

4.1 推荐：gdal_calc.py（一条命令，三波段）

gdal_calc.py -A dem_final.tif \
  --outfile=dem_terrarium_rgb-py-band.tif \
  --type=Byte \
  --overwrite \
  --creation-option=INTERLEAVE=BAND \
  --creation-option=COMPRESS=LZW \
  --creation-option=TILED=YES \
  --calc="floor((numpy.clip(A,-32768,32767)+32768)/256)" \
  --calc="floor((numpy.clip(A,-32768,32767)+32768)%256)" \
  --calc="floor(((numpy.clip(A,-32768,32767)+32768)-floor(numpy.clip(A,-32768,32767)+32768))*256)"

务必在编码阶段指定 INTERLEAVE=BAND，避免后续 COG 金字塔 G/B 通道错位。

4.2 编码后校验（必做）

gdalinfo -stats dem_terrarium_rgb-py-band.tif

快速 Python 抽查：

python3 -c "
from osgeo import gdal
ds=gdal.Open('dem_terrarium_rgb-py-band.tif')
r=ds.GetRasterBand(1).ReadAsArray().astype(float)
g=ds.GetRasterBand(2).ReadAsArray().astype(float)
b=ds.GetRasterBand(3).ReadAsArray().astype(float)
h=r*256+g+b/256-32768
print('R', r.min(), r.max(), 'elev', round(h.min(),1), round(h.max(),1))
"

5. 转 COG（修复低 zoom 尖刺）

5.1 问题背景

INTERLEAVE=PIXEL 的 RGB GeoTIFF 在 gdal_translate -of COG 生成 overview 时，会出现 R 正确、G/B 地理错位：

• 全分辨率：相邻像元高程差 < 20 m
• 错误 overview：相邻像元可差 400+ m → MapLibre 缩小视角出现「尖刺」，放大后消失

5.2 正确命令

§4 已指定 INTERLEAVE=BAND 时，直接转 COG：

gdal_translate dem_terrarium_rgb-py-band.tif dem_terrarium_cog-py-band.tif -of COG \
  -co OVERVIEW_RESAMPLING=NEAREST

避免低 zoom 尖刺的必要条件：§4 产物为 BAND + 本步 NEAREST。仅执行本命令、而输入为 PIXEL 交错 RGB，仍可能出现尖刺。

若历史产物为 PIXEL 交错，先转 BAND 再出 COG：

gdal_translate dem_terrarium_rgb-py-band.tif dem_terrarium_band.tif \
  -co INTERLEAVE=BAND -co COMPRESS=LZW -co TILED=YES

gdal_translate dem_terrarium_band.tif dem_terrarium_cog.tif -of COG \
  -co OVERVIEW_RESAMPLING=NEAREST

5.3 COG 校验

gdalinfo dem_terrarium_cog.tif
# 期望：LAYOUT=COG；INTERLEAVE=BAND；有 Overviews

抽查 overview 与全分辨率 RGB 是否一致（同一地理位置）：

python3 -c "
from osgeo import gdal
ds=gdal.Open('dem_terrarium_cog.tif')
w,h=ds.RasterXSize,ds.RasterYSize
ov=ds.GetRasterBand(1).GetOverview(0)
ow,oh=ov.XSize,ov.YSize
ox,oy=ow//2,oh//2
fx,fy=int(ox*w/ow),int(oy*h/oh)
full=[int(ds.GetRasterBand(i).ReadAsArray(fx,fy,1,1)[0,0]) for i in (1,2,3)]
ovpx=[int(ds.GetRasterBand(i).GetOverview(0).ReadAsArray(ox,oy,1,1)[0,0]) for i in (1,2,3)]
print('full RGB',full,'overview RGB',ovpx,'match',full==ovpx)
"
# 期望：match True

6. 发布与切片（可选）

6.1 上传 MinIO

# 示例路径（与 demo 一致）
# http://yourip:9000/gis-test/terrain-test/anshan/dem_terrarium_cog.tif
mc cp dem_terrarium_cog.tif minio/gis-test/terrain-test/anshan/

6.2 XYZ 瓦片（可选）

gdal2tiles.py -z 8-15 -r near dem_terrarium_rgb-py-band.tif tiles/

Terrarium 瓦片必须 -r near，与 COG overview 同理：禁止对 RGB 做 bilinear/cubic 重采样。

6.3 PMTiles（可选）

pmtiles convert tiles/ dem_terrarium.pmtiles

需预先安装 pmtiles CLI（参见 20260629本机操作记录.md §5 故障排查）。

7. MapLibre 加载

7.1 COG（maplibre-cog-protocol）

maplibregl.addProtocol('cog', MaplibreCOGProtocol.cogProtocol);

map.addSource('dem', {
  type: 'raster-dem',
  url: 'cog://http://yourip:9000/gis-test/terrain-test/anshan/dem_terrarium_cog.tif',
  tileSize: 256,
  encoding: 'terrarium'
});

map.setTerrain({ source: 'dem', exaggeration: 1.5 });

map.addLayer({
  id: 'hillshade',
  type: 'hillshade',
  source: 'dem',
  paint: { 'hillshade-exaggeration': 0.5 }
});

• 预编码 3 波段 RGB COG：不要加 #dem 后缀。
• Demo：demo/terrarium-cog.html

7.2 PMTiles

map.addSource('terrain', {
  type: 'raster-dem',
  url: 'pmtiles://https://{host}/dem_terrarium.pmtiles',
  encoding: 'terrarium',
  tileSize: 256
});

7.3 注意事项：Demo 须使用 COG，勿直接用普通 TIF

本流程 §4 产出的 dem_terrarium_rgb-py-band.tif 仅为中间 GeoTIFF，不可作为 Web Demo 的最终发布物直接加载。

结论：面向 MapLibre Demo / 在线预览时，必须完成 §5 转 COG，并发布 dem_terrarium_cog.tif；跳过 COG 步骤、直接用编码后的普通 TIF，展示效果无法满足验收要求。参见 demo/terrarium-cog.html 中的 COG 地址配置。

8. 产物清单

dem_clean.tif              # NoData 修复（§3.1）
dem_3857.tif               # Web Mercator（§3.2）
dem_s1.tif / dem_s2.tif    # 中间尺度（§3.3）
dem_final.tif              # 预处理终点（§3.3）
dem_terrarium_rgb-py-band.tif      # Terrarium 3 波段 Byte RGB，INTERLEAVE=BAND（§4）
dem_terrarium_cog-py-band.tif      # COG 发布物，NEAREST overview（§5）
tiles/                     # 可选 XYZ（§6.2）
dem_terrarium.pmtiles      # 可选 PMTiles（§6.3）

9. 踩坑速查

对照组：同 dem_final.tif 可用 Mapbox 编码快速验证预处理是否正确：

rio rgbify -b -10000 -i 0.1 dem_final.tif terrain_rgb.tif
gdal_translate terrain_rgb.tif terrain_rgb_cog.tif -of COG
# Demo: demo/mapbox-cog.html

10. 相关文档

• 2026-06-29 本机操作记录（环境、历史踩坑、§8 详细分析）
• 2026-06-30 DEM 工程实践记录（预处理质量、PMTiles 目标）
• 预处理与平滑
• Terrain-RGB 概览
• MapLibre 加载 COG 说明
• 解码公式归档