C# 结合 GDAL 把 .tif 高程图转换为 .obj 地形模型

把一张 .tif 高程图转成 .obj，核心并不是“把图片另存为模型”，而是先把像素里的高度值读出来，再按网格关系重建顶点和三角面。下面这篇会按实际流程拆开讲：先用 GDAL 读取高程数据，再在 C# 中生成顶部、侧壁和底面，最后写出标准 .obj 文件。读完后，你可以判断这套方案适合直接做地形模型，还是只适合作为后续 OpenTK / OpenGL 渲染前的数据预处理。

读取 .tif 高程图：先拿到可用的高度矩阵

这一步使用的是 GDAL。目标很明确：从 .tif 中读取单波段像素值，转成后续建模可用的二维高程数组，并同时标记哪些点不能参与建模。

GDAL 读取流程

原始代码里先完成 GDAL 初始化，再打开数据集、读取第一波段，并把像素值灌入一维缓冲区，最后再整理成二维数组：

Gdal.AllRegister(); //初始化GDAL

//获取高程图像.tif的图像信息dem
using Dataset ds = Gdal.Open(".tif图片路径", Access.GA_ReadOnly);
width = ds.RasterXSize;
height = ds.RasterYSize;
Band band = ds.GetRasterBand(1);
band.GetNoDataValue(out double noData, out int hasNoData);

float[] buffer = new float[width * height];
band.ReadRaster(0, 0, width, height, buffer, width, height, 0, 0);

float[,] dem = new float[height, width];
isNoData = new bool[height, width];

for (int y = 0; y < height; y++)
{
    for (int x = 0; x < width; x++)
    {
        float v = buffer[y * width + x];
        bool nd = (hasNoData == 1 && Math.Abs(v - (float)noData) < 0.001f) || v == 0f;
        dem[y, x] = nd ? 0 : v;
        isNoData[y, x] = nd;
    }
}

为什么要单独处理 NoData

高程图不是每个像素都一定有效。代码里把两类值视为不可用点：

• GDAL 明确标出的 NoData 值；
• 值为 0f 的像素。

随后把这些点在 dem 中置为 0，同时在 isNoData 里记录下来。这样做的实际意义是：后面生成网格时，可以直接跳过无效区域，避免把空洞也拼成三角面。

这里也要注意一个判断前提：把 0f 直接当成无效值，只适合“0 本身不代表真实海拔”的数据。如果你的数据里 0 就是合法高程，这条规则需要调整。

把高程矩阵变成 Mesh：顶点和三角面怎么组织

有了二维高程数组后，下一步就是把每个有效像素映射成三维顶点，再按网格相邻关系生成三角面。原文使用了一个简化的 ObjMesh 类来存储顶点和面：

internal class ObjMesh //自定义的模型类
{
    public List<Vector3> Vertices { get; } = new();
    public List<(int, int, int)> Faces { get; } = new();
}

先确定坐标系和中心偏移

如果直接按像素坐标建模，模型原点通常会落在角上，不利于后续在 OpenTK / OpenGL 里旋转或缩放。这里先计算中心偏移量，把几何中心移动到坐标原点附近：

int rows = dem.GetLength(0);
int cols = dem.GetLength(1);

// 计算中心点偏移量
float centerX = (cols - 1) * 0.5f * xyScale;
float centerZ = -(rows - 1) * 0.5f * xyScale;  // 注意：Z轴为负方向

ObjMesh mesh = new ObjMesh();
int[,] topIndexMap = new int[rows, cols];
int[,] bottomIndexMap = new int[rows, cols];

这里有两个关键参数：

• xyScale：控制平面上相邻采样点的间距；
• zScale：控制高程值映射到 Y 轴后的拉伸倍数。

另外，代码把 Z 轴写成负方向，这样二维图像从上到下扫描时，三维空间里的前后方向也更统一。

先生成顶部顶点，再铺顶部三角面

顶部表面是整个地形模型的主体。处理方法是：按 step 间隔扫描高程矩阵，为每个有效点创建一个三维顶点；无效点则在映射表中记为 0。

// ---------- 顶部顶点 ----------
for (int y = 0; y < rows; y += step)
{
    for (int x = 0; x < cols; x += step)
    {
        if (outerNoData[y, x])
        {
            topIndexMap[y, x] = 0;
            continue;
        }

        topIndexMap[y, x] = mesh.Vertices.Count + 1;
        // 将坐标原点移动到中心点
        float xPos = x * xyScale - centerX;
        float zPos = -y * xyScale - centerZ;  // 注意：Z轴为负方向
        mesh.Vertices.Add(new Vector3(
            xPos,
            dem[y, x] * zScale,
            zPos
        ));
    }
}

之所以要维护 topIndexMap，是因为 OBJ 顶点索引从 1 开始，后面拼三角面时不能再反查坐标，只能依赖索引表快速找到每个采样点对应的顶点编号。

随后再把相邻四个点组成两个三角面：

// ---------- 顶部三角面 ----------
for (int y = 0; y < rows - step; y += step)
{
    for (int x = 0; x < cols - step; x += step)
    {
        int v00 = topIndexMap[y, x];
        int v10 = topIndexMap[y, x + step];
        int v01 = topIndexMap[y + step, x];
        int v11 = topIndexMap[y + step, x + step];

        if (v00 == 0 || v10 == 0 || v01 == 0 || v11 == 0)
            continue;

        mesh.Faces.Add((v00, v01, v10));
        mesh.Faces.Add((v10, v01, v11));
    }
}

这段逻辑的重点是：只要四个角里有一个点无效，当前网格就直接跳过。这样能避免空洞边缘出现错误连面。

为什么还要补底面和侧壁：让地形从曲面变成实体

如果只输出顶部三角面，得到的其实更接近“一张起伏曲面”，而不是一个适合展示、打印或继续加工的封闭实体。原文后半段的重点，就是补齐底部和四周侧壁。

先给每个有效顶部点生成对应的底部点

做法很直接：只要顶部有有效点，就在固定的 bottomY 高度再生成一个底部顶点，并记录到 bottomIndexMap。

// ---------- 创建底部顶点映射 ----------
for (int y = 0; y < rows; y += step)
{
    for (int x = 0; x < cols; x += step)
    {
        // 只有顶部有效的点才有底部顶点
        if (topIndexMap[y, x] != 0)
        {
            Vector3 topV = mesh.Vertices[topIndexMap[y, x] - 1];
            bottomIndexMap[y, x] = mesh.Vertices.Count + 1;
            mesh.Vertices.Add(new Vector3(
                topV.X,
                bottomY,
                topV.Z
            ));
        }
    }
}

这样顶部和底部在 X、Z 平面上一一对应，只是 Y 值不同，后面补墙和补底都更容易处理。

侧壁如何处理边界和缺口

侧壁生成是这篇代码里最复杂的部分。它分别沿 X 方向和 Y 方向扫描，如果相邻两个顶部点都存在，并且至少有一个点是边缘点，就创建两片三角形，把顶部边和底部边连起来。

当一边有点、另一边没点时，代码不会直接放弃，而是临时补出一个新的底部点，再创建单个三角形去封住缺口。这一段正是为了处理 NoData 区域造成的不规则边界。

原文保留了完整实现，核心片段如下：

// ---------- 创建侧边墙 ----------
// 垂直边（X方向）
for (int y = 0; y < rows; y += step)
{
    for (int x = 0; x < cols - step; x += step)
    {
        int idx00 = topIndexMap[y, x];
        int idx10 = topIndexMap[y, x + step];

        // 如果两个顶点都存在，且至少有一个是边缘点
        if (idx00 != 0 && idx10 != 0)
        {
            bool isEdge = IsEdgePoint(topIndexMap, x, y, rows, cols, step) ||
                         IsEdgePoint(topIndexMap, x + step, y, rows, cols, step);

            if (isEdge)
            {
                int b00 = bottomIndexMap[y, x];
                int b10 = bottomIndexMap[y, x + step];
                mesh.Faces.Add((idx00, b00, b10));
                mesh.Faces.Add((idx00, b10, idx10));
            }
        }
        // 如果一个存在一个不存在，创建三角形连接
        else if (idx00 != 0 && idx10 == 0)
        {
            int b00 = bottomIndexMap[y, x];
            Vector3 topV = mesh.Vertices[idx00 - 1];
            int b10 = mesh.Vertices.Count + 1;
            // 计算新的底部顶点坐标（使用中心偏移）
            float xPos = (x + step) * xyScale - centerX;
            float zPos = -y * xyScale - centerZ;
            mesh.Vertices.Add(new Vector3(
                xPos,
                bottomY,
                zPos
            ));
            mesh.Faces.Add((idx00, b00, b10));
            // 更新映射
            bottomIndexMap[y, x + step] = b10;
        }
        else if (idx00 == 0 && idx10 != 0)
        {
            int b10 = bottomIndexMap[y, x + step];
            Vector3 topV = mesh.Vertices[idx10 - 1];
            int b00 = mesh.Vertices.Count + 1;
            // 计算新的底部顶点坐标（使用中心偏移）
            float xPos = x * xyScale - centerX;
            float zPos = -y * xyScale - centerZ;
            mesh.Vertices.Add(new Vector3(
                xPos,
                bottomY,
                zPos
            ));
            mesh.Faces.Add((idx10, b10, b00));
            // 更新映射
            bottomIndexMap[y, x] = b00;
        }
    }
}

Y 方向的侧壁处理逻辑与此相同，只是扫描方向换成了纵向相邻点。

底面怎么补，才能尽量形成封闭体

底面部分分成两层处理：

• 如果一个网格四个角都有有效顶部点，就生成对应的两片底部三角面；
• 如果某些点是孤立点，附近又凑不出完整网格，就额外生成一个小三角形底面，避免底部破洞。

这种写法的目的很直接：尽可能把模型补成一个封闭壳体，而不只是保留上表面。

其中常规底面生成如下：

// ---------- 创建底面（只覆盖有顶部顶点的地方） ----------
// 只在顶部面存在的区域创建底面三角形
for (int y = 0; y < rows - step; y += step)
{
    for (int x = 0; x < cols - step; x += step)
    {
        int t00 = topIndexMap[y, x];
        int t10 = topIndexMap[y, x + step];
        int t01 = topIndexMap[y + step, x];
        int t11 = topIndexMap[y + step, x + step];

        // 如果这个网格的顶部存在（四个角都有顶部顶点），则创建对应的底面
        if (t00 != 0 && t10 != 0 && t01 != 0 && t11 != 0)
        {
            int b00 = bottomIndexMap[y, x];
            int b10 = bottomIndexMap[y, x + step];
            int b01 = bottomIndexMap[y + step, x];
            int b11 = bottomIndexMap[y + step, x + step];

            // 底面三角面（注意法线方向向下，顶点顺序要逆时针）
            mesh.Faces.Add((b00, b10, b01));
            mesh.Faces.Add((b10, b11, b01));
        }
    }
}

这里最容易出错的是顶点顺序。原文特别强调底面法线要朝下，所以三角面索引顺序必须与顶部不同。

最后一步：把 Mesh 写成 .obj 文件

当前面的顶点和面都准备好之后，导出 .obj 就是标准文本写入过程。OBJ 文件里最关键的两类记录分别是：

• v x y z：顶点坐标；
• f i j k：三角面引用的顶点索引。

原始实现如下：

private static void WriteObj(string path, ObjMesh mesh)
{
    using StreamWriter sw = new StreamWriter(path);
    sw.WriteLine("# DEM Terrain with Walls + Bottom");
    foreach (var v in mesh.Vertices)
        sw.WriteLine($"v {v.X:F6} {v.Y:F6} {v.Z:F6}");
    foreach (var f in mesh.Faces)
        sw.WriteLine($"f {f.Item1} {f.Item2} {f.Item3}");

}

这里保留了 F6 精度格式，说明作者希望导出的坐标保持 6 位小数。对于大多数地形网格导出，这已经够用了。

这套 .tif 转 .obj 方案适合什么场景

从代码结构看，这套实现更像是“高程数据转封闭地形网格”的工具逻辑，适合接在 OpenTK / OpenGL 渲染流程之前使用，也可以单独当作地形模型预处理步骤。

它的优点是路径清晰：读取高程、过滤无效值、生成顶部、补侧壁、补底面、导出 OBJ，整个过程完整闭环。需要额外留意的地方主要有三个：

• 0f 是否应当被视为无效高程；
• step 会直接影响模型精度与顶点数量；
• xyScale、zScale、bottomY 会决定最终模型比例和厚度。

如果你的目标只是可视化地形表面，顶部网格可能已经够用；但如果要导入建模软件、做封闭展示模型，或者继续用于打印、布尔运算，补齐侧壁和底面就是必要步骤。