1. 引言
在数字孪生、三维地图可视化、智慧城市等项目中，有一个非常常见的需求：给定一个行政区划（省/市/县/区）或任意自定义区域的边界数据，生成该区域的顶面覆盖贴图——比如卫星影像纹理和地形法线贴图，以便在三维场景中将它们"披"在地形模型上，呈现出真实的地貌效果。
通常情况下，我们手头只有一份 GeoJSON 格式的区域边界数据（包含经纬度坐标的封闭多边形），而缺乏该区域的现成贴图。手动去截图、拼接、处理不仅效率低下，且难以覆盖大量区域。
为此我们搭建了一套全自动管线：给定一个 GeoJSON 边界文件，自动下载卫星瓦片并合成高清影像，同时从全球高程数据中生成法线贴图，最后根据边界裁剪输出，全程无需人工介入。 点此进入易知微孪生资产库
2. 基础知识铺垫
在深入技术细节之前，我们先快速了解几个本文会反复出现的基础概念。
2.1. GeoJSON —— 用 JSON 画地图
GeoJSON 是一种基于 JSON 的地理空间数据交换格式，可以描述点（Point）、线（LineString）、多边形（Polygon）等几何对象。你可以把它理解成"地图界的 SVG"——只不过 SVG 的坐标是屏幕像素，而 GeoJSON 的坐标是真实世界的经纬度。
例如，一个描述浙江省边界的 GeoJSON 文件结构大致如下：
{
"type": "FeatureCollection",
"features": [{
"type": "Feature",
"properties": { "name": "浙江省", "adcode": "330000" },
"geometry": {
"type": "Polygon",
"coordinates": [[[118.0, 27.0], [118.5, 28.0], ...]]
}
}]
}
本文的核心输入就是这个多边形坐标列表，我们用它来确定"需要处理哪片区域"以及"最终裁剪到什么形状"。
2.2. 瓦片地图 —— 地球被切成无数小方块
在线地图服务（谷歌、高德、OpenStreetMap 等）不是把一整张世界地图的图片发给你，而是将地球表面按照不同缩放层级切分成无数个 256×256 像素的小方块，称为瓦片（Tile） 。
瓦片的编号体系使用三个整数 (x, y, z)：
- z（Zoom Level） ：缩放层级，取值通常在 0~18 之间。z=0 时整个地球只有 1 张瓦片；z 每增加 1，瓦片数量翻 4 倍。
- x, y：该层级下瓦片的列号和行号。
这个概念对于游戏/图形开发者来说非常熟悉——类似纹理的 Mipmap 或游戏的 LOD（Level of Detail） 机制：远看用低精度，近看用高精度。
2.3. 墨卡托投影 —— 把球面"摊平"
地球是球形的，但地图是平面的。把球面坐标（经纬度）映射到平面坐标（像素）的过程叫投影。绝大多数 Web 地图（包括高德、谷歌、OpenStreetMap）使用的标准投影是 Web 墨卡托投影（Web Mercator） 。
墨卡托投影的核心特点是：
- 保持方向和形状：在地图上，北就是上，正方形还是正方形（局部保形）。
- 高纬度区域拉伸：格陵兰岛（北极圈附近）在地图上看起来比非洲还大，但实际上非洲的面积是格陵兰的 14 倍。
- 南北纬约 85° 以上的极地无法显示（数学上的奇点）。
对于我们的方案来说，必须正确处理墨卡托投影的坐标转换，才能从瓦片服务中精准下载对应区域的图片，并在法线图生成时按等面积方式采样高程数据。
2.4. DEM —— 用数字矩阵记录地形
数字高程模型（DEM, Digital Elevation Model） 是以栅格形式记录地表海拔高度的数据集——你可以把它理解成一张"灰度图"，每个像素的值代表该位置的海拔高度（单位：米）。
本文使用的是 SRTM（Shuttle Radar Topography Mission） 数据——由 NASA 的航天飞机在 2000 年通过雷达干涉测绘获得，覆盖全球北纬 60° 至南纬 56° 的陆地区域。SRTM 数据的分辨率约为 30 米（美国本土）到 90 米（全球其他区域），作为公开数据可以合法用于研究和学术目的。
合法合规提示：SRTM 数据由 NASA 和 USGS 公开发布。其他高分辨率 DEM 数据（如 ALOS AW3D30、ASTER GDEM）各有自身的许可条款，使用前请查阅相应数据提供方的使用协议。
3. 核心流程总览
整体方案包含两条独立且并行的 Pipeline：
卫星贴图生成流程：
GeoJSON 边界 → 边界框计算 → 自适应 Zoom → 瓦片下载 → 图片拼接 → 按边界裁剪 → 卫星贴图
法线贴图生成流程：
GeoJSON 边界 → 边界框计算 → DEM 高程采样 → 空洞填充 → 梯度计算 → 法线合成 → 法线贴图
两条 Pipeline 共用同一个输入（GeoJSON 区域边界数据），但后续处理链路完全独立，可以并行执行。
下面分别展开两条 Pipeline 的技术细节。
4. 卫星贴图生成详解
卫星贴图生成的本质是：根据区域边界框下载覆盖该区域的所有瓦片图片，拼成一张大图，再按边界精确裁剪。
4.1. 输入：区域边界解析
第一步是读取 GeoJSON 文件，提取其中的多边形（Polygon）或多多边形（MultiPolygon）几何体，然后遍历所有顶点坐标，计算出覆盖该区域的最小边界框（Bounding Box）：
BBox = (min_lon, min_lat, max_lon, max_lat)
注意：GeoJSON 中一个 Feature 可能是 MultiPolygon 类型（如包含岛屿的沿海省份），需要处理所有子多边形。
4.2. 自适应 Zoom 层级选择
这是整个方案中一个关键的工程决策：选多大的缩放层级？
- Zoom 太大（如 z=18）：单个瓦片覆盖范围极小，覆盖一个省可能需要上万张瓦片——下载极慢、内存爆炸
- Zoom 太小（如 z=8）：瓦片数量虽少，但分辨率太低，最终图片模糊不清
我们的策略是 "反向估算" ：对于每个候选的 Zoom 层级（从 z=18 到 z=5 依次扫描），计算如果采用该层级需要下载多少张瓦片，选择使瓦片数量落在 50~500 张 范围内的最高层级。
这一策略确保了：
- 小面积区域（如一个区县）自动匹配较高 Zoom，获得充足细节
- 大面积区域（如一个省份）自动降级 Zoom，控制下载量在可接受范围
可以通俗地理解为： "根据你的区域大小自动选择最佳分辨率，类似视频网站的自适应码率" 。
4.3. 经纬度 → 瓦片坐标转换
这是 Web 地图技术中最核心的坐标转换步骤。给定一个经纬度点 (lon, lat) 和 Zoom 层级 z，需要计算它落在哪张瓦片上。
核心公式（简化版）：
n = 2^z // 该层级下瓦片总数（每行/列 = n）
// 经纬度 → 瓦片坐标
tile_x = (lon + 180) / 360 * n // 经度线性映射到 [0, n)
lat_rad = radians(lat)
tile_y = (1 - asinh(tan(lat_rad)) / π) / 2 * n // 纬度经墨卡托投影后映射
// 取整得到瓦片行号、列号
x = floor(tile_x)
y = floor(tile_y)
反向转换（瓦片坐标 → 经纬度）同样重要，用于计算拼接后大图的实际地理跨度：
lon = tile_x / n * 360 - 180
lat = atan(sinh(π * (1 - 2 * tile_y / n)))
这里的
asinh（反双曲正弦）和atan(sinh(...))就是墨卡托投影的数学内核。tan(lat_rad)把球面坐标拉伸到投影面，asinh再把它变成适合线性分块的尺度。
4.4. 并发瓦片下载
确定瓦片范围后，生成所有需要下载的 (x, y, z) 元组列表。使用线程池并发下载，大幅提升效率。
对于省级区域（约 200~400 张瓦片），并发下载可在 30 秒内完成。
合法合规使用第三方地图服务的注意事项
在线瓦片地图服务（高德、谷歌、OSM 等）为开发者提供了便捷的地图数据获取途径，但在使用时需要特别注意：
- 遵守服务条款：各平台对数据获取方式、使用场景（商业/非商业）有明确约定，请在接入前仔细阅读。
- 控制请求频率：高频并发请求可能对服务器造成负担，建议使用合理的并发数和重试策略，并在请求间加入适当间隔。
- 商业产品需走正规渠道：如果你的产品面向商业客户，建议通过官方 API（如高德开放平台、Google Maps Platform）注册 API Key，获取稳定且有 SLA 保障的服务。
- 本文方案的技术定位：本文分享的技术方案仅用于技术研究和学习目的，实际产品化时应根据具体场景评估合规性和选择合适的服务商。
4.5. 瓦片拼接（Image Stitching）
将所有成功下载的 256×256 瓦片按行列位置"贴"到一张大图上：
合成图 = Image.new(width=cols*256, height=rows*256)
for row in range(rows):
for col in range(cols):
if (start_x+col, start_y+row) in 已下载瓦片:
大图.paste(瓦片图片, (col*256, row*256))
拼接后的像素尺寸取决于瓦片行列数，例如覆盖一个省可能需要拼接出 5000×4000 像素级别的大图。
4.6. 按 GeoJSON 边界精确裁剪
这是最后一步，也是最"GIS 味"的一步——把拼接好的矩形大图按 GeoJSON 多边形的形状裁剪：
- 多边形坐标映射：将 GeoJSON 中的每个经纬度点
(lon, lat)通过 4.3 节的反向公式映射到拼接大图中的像素坐标(px, py) - 创建遮罩（Mask） ：用 PIL 的
ImageDraw.polygon()在空画布上绘制多边形——多边形内部填充白色（255，=保留），外部黑色（0，=丢弃） - 裁剪边界框：计算多边形在图片中占据的像素边界框并裁出
- 缩放至目标尺寸：保持宽高比，将有效内容贴边缩放至 2048×2048 像素
最终输出：一张正方形 (2048×2048) 的卫星贴图，区域恰好覆盖 GeoJSON 边界内的所有内容。
5. 法线贴图生成详解
有了卫星贴图只是"画皮"——没有光影变化的地形看起来是平面的。法线贴图（Normal Map） 就是用来增强地形光影真实感的关键贴图。
5.1. 什么是法线贴图？—— 给像素"指方向"
在 3D 渲染中，每个表面都有一个"朝向"属性——叫法线（Normal） ，它是一个指向表面外侧的单位向量（长度为 1）。光照计算（如漫反射、镜面反射）必须知道法线才能算出表面该多亮多暗。
但 3D 模型的面数是有限的——一个地形模型不可能把每个小石头和微小起伏都建模出来。法线贴图 就是解决这个问题的方案：
法线贴图是一张 RGB 图片，每个像素的值不表示"颜色"，而是编码了该位置表面的"朝向"。渲染时，GPU 读取这个"假朝向"来计算光照，让平面看起来有凹凸——这就是计算机图形学中的"凹凸映射（Bump Mapping）"技术。
编码规则很简单：
| 通道 | 编码内容 | 通俗解释 |
|---|---|---|
| R（红） | 法线的 X 分量 | 表面"向左倾斜"还是"向右倾斜" |
| G（绿） | 法线的 Y 分量 | 表面"向前倾斜"还是"向后倾斜" |
| B（蓝） | 法线的 Z 分量 | 表面"直直地朝向我们"的程度 |
当一个像素的颜色是 (128, 128, 255) ——偏蓝色——表示该表面的法线是 (0, 0, 1) ，即"直直向上"，是完全平坦的地面。这也是大多数法线贴图的基色调为淡蓝色的原因。
在三维地形引擎（如 CesiumJS、Unity3D）中，给地形模型贴上法线贴图后，山脊的向阳面会变亮，背光面会变暗——就像真实光照打在真实地面上一样。
5.2. DEM 高程数据采样
法线贴图的生成流程可以概括为三步：采样高程 → 计算梯度 → 合成法线。
首先从 SRTM 数据中采样目标区域的高程。有一个重要的采样策略：按墨卡托投影像素等间距采样，而非按经纬度等间距采样。
为什么要这样做？回忆 2.3 节墨卡托投影的特点：高纬度区域被拉伸。如果按经纬度等差采样（如每隔 0.01° 取一个点），在哈尔滨（北纬 45°）的单位面积内采样点会比在海口（北纬 20°）同一单位面积内的采样点多——导致高纬度区域过度采样、法线密度不均匀。
正确的做法是：
- 将区域的经纬度范围映射到墨卡托投影像素坐标范围
- 在墨卡托投影空间中等距采样（如均分 200×200 个采样点）
- 再将每个采样点的墨卡托坐标反算回经纬度，查询 SRTM 高程值
这样保证了每个采样点在地图上覆盖的实际面积是一致的。
5.3. 数据预处理：空洞填充
SRTM 数据在某些区域可能返回无效值——水体区域（如湖泊、大型河流）在雷达测绘中无法获取有效高程，可能返回 None 或极小值。
需要识别并填充这些空洞。我们借助 SciPy 的 distance_transform_edt（欧氏距离变换）实现最近邻插值：对每个无效像素，找到离它最近的有效像素的值，以此填充。这确保了整个高程矩阵是完整且连续的。
5.4. 高度归一化策略
采样完成后，需要将所有高度值映射到 [0, 1] 区间，才能参与后续的梯度计算。但这里有个关键设计选择：线性归一化 还是 对数归一化？
| 情况 | 策略 | 理由 |
|---|---|---|
| 高差 ≤ 100m | 线性归一化 | 高度范围窄，线性映射足够保留细节 |
| 高差 > 100m | 对数归一化（log1p） | 压缩极高的山峰，为中低海拔区域留出更多"灰度空间" |
以实际地形为例：一个区域包含海平面（0m）和一座 3000m 的高峰。如果用线性归一化，海平面到 1000m 的丘陵地带只占整个灰度范围的 1/3——大量的地形成分被压缩到少数几个灰度级里，细节丢失殆尽。
对数归一化通过 log1p(height) 压缩高值：
- 0m → log1p(0) ≈ 0
- 1000m → log1p(1000) ≈ 6.9
- 3000m → log1p(3000) ≈ 8.0
可以看到低海拔区域的相对差异被拉大了，丘陵的起伏在最终的归一化结果中得到了更好的体现。
5.5. Sobel 算子计算梯度
接下来是最核心的步骤——用 Sobel 算子 计算归一化高程矩阵中每个点的梯度（即"高度在水平和垂直方向上的变化率"）。
Sobel 是一种经典的 3×3 卷积核，广泛用于图像边缘检测。这里用它来"感知"地形的陡峭程度：
Sobel X 核（检测水平方向的高度变化）：
[-1, 0, 1]
[-2, 0, 2]
[-1, 0, 1]
Sobel Y 核（检测垂直方向的高度变化）：
[-1, -2, -1]
[ 0, 0, 0]
[ 1, 2, 1]
将这两个核分别与归一化高程矩阵做二维卷积（scipy.ndimage.convolve），得到：
- dx：每个点水平方向的高度变化率——正值表示"右边比左边高"（东高西低），负值表示"左边比右边高"（西高东低）
- dy：每个点垂直方向的高度变化率——正值表示"上边比下边高"（北高南低），负值表示"下边比上边高"（南高北低）
如果你熟悉图像处理，可以把它理解为"在高度图上跑边缘检测，突出的边缘就是山脊和陡坡"。
5.6. 法线向量合成
有了梯度（dx, dy），就可以合成每个像素的法线向量了：
// 法线 = 梯度的"反向"（因为梯度指向"升高"方向，法线需要指向"表面朝向"）
normal_x = -dx × strength × scale_factor
normal_y = -dy × strength × scale_factor
normal_z = 0.3 // 保持基本向上的 Z 分量，避免法线过于"平躺"
然后将法线向量归一化（除以向量长度，使其成为单位向量），最后映射到 RGB 空间：
R = (normal_x × 0.5 + 0.5) × 255 // 将 [-1, 1] 映射到 [0, 255]
G = (normal_y × 0.5 + 0.5) × 255
B = (normal_z × 0.5 + 0.5) × 255
其中 scale_factor 是一个自适应系数：min(实际高差 / 500, 8.0)。它的含义是：
- 高差越大 → 法线起伏越剧烈 → 山脊的明暗对比越强，视觉冲击力更足
- 上限 8.0 防止在极端高差（如青藏高原）下法线过于扭曲失实
strength 是一个可调参数（默认 1.5），允许人工微调法线的表现力强度。
5.7. 输出
与卫星贴图类似，最终的法线图也需要缩放到 2048×2048。处理方式：
- 保持原始高程矩阵的宽高比贴边放置
- 空白区域填充基准色 (128, 128, 255) ——即"平坦地面"的颜色
6. 实际效果与对比
将同一区域的卫星贴图和法线贴图同时应用到三维地形模型上，与"仅贴卫星图"的效果对比：
| 对比维度 | 仅有卫星贴图 | 卫星贴图 + 法线贴图 |
|---|---|---|
| 光照感 | 平面光，无论地形高低亮度一致 | 动态光照，向阳面亮、背光面暗 |
| 地形辨识度 | 靠颜色猜地形（绿=平原，棕=山地） | 靠光影感知隆起，山谷与山脊一清二楚 |
| 视觉真实感 | 像"打印在地面上的照片" | 像"真实阳光照在真实地面上" |
特别是在山区、丘陵地带，配合法线贴图后的效果提升极为显著——原本模糊的山脊线变得清晰可辨，即使缩小到远景视角也能感受到地形起伏。
7. 真实案例
易知微基于多年在数字孪生及数据可视化领域丰富实践，沉淀了诸多经验成果，欢迎大家互相交流学习：《数字孪生+AI：行业最佳实践白皮书》下载地址：easyv.cloud/references/…
《数字孪生视觉语言白皮书》下载地址：easyv.cloud/references/…
《数字孪生世界白皮书》下载地址：easyv.cloud/references/…
《数字孪生行业方案白皮书》下载地址：easyv.cloud/references/…
《港口数智化解决方案》下载地址：easyv.cloud/references/…
想申请易知微产品免费试用的客户，欢迎点击易知微官网申请试用：易知微 - 数字孪生仿真渲染引擎与可视化应用
