OpenStreetMap数据转3D场景【Python + PostgreSQL】

很长一段时间以来，我对 GIS 和渲染感兴趣，在分别尝试这两者之后，我决定最终尝试以 3D 方式渲染 OpenStreetMap 中的地理数据，重点关注不超过城市的小规模。

在本文中，我将介绍从建筑形状生成三角形网格、以适合 Blender 或 Godot 等游戏引擎的格式渲染和导出它的过程。我不是该领域的专家，但我确信有人面临着同样的问题，他们可能会喜欢阅读本文。

总的来说，我发现 GIS 和 3D 处理主题非常令人兴奋，因为它将计算机科学与几何和代数相结合，并且在某种意义上讲是人类如何感知和描述世界。

在这里插入图片描述

推荐：用 NSDT设计器快速搭建可编程3D场景。

1、将 OSM 数据导入 PostGIS

该过程的第一步是将 OpenStreetMap 数据放入 PostGIS 实例中。这并不是绝对必要的，因为像 osmnx 这样的库可以使用 API 在几行 Python 中动态获取这些数据，但我喜欢在关系数据库中使用它，因为它在分析和处理数据时提供了很大的灵活性，因为我可以使用 SQL。

为此，我从 Geofabrik 下载提取内容，本质上是带有单个区域数据的 protobuf 文件，并使用 pgOSM Flex 导入它们，这是我已经在 2D 可视化项目中使用的工具，并且我过去曾对此做出过一些贡献。

导入数据后，我可以使用 QGIS 直接将其可视化。
在这里插入图片描述

借助 QGIS，只需单击几下即可连接到 PostGIS 并探索数据

由于我使用的是 Python，因此我可以使用 Psycopg3 轻松提取 Shapely 表示（事实上，我添加了该功能😎）。 Asyncpg 也可以透明地执行相同的转换，但不能很好地与 Geopandas 配合使用。

Shapely 几乎涵盖了人们可能需要的有关 2D 几何（地理或非地理）的所有内容，并且非常节省内存，因此它是此类项目的宝贵资源。

2、3D 三角形网格

现在我有了一个代表建筑物等特征的 Shapely 对象，我“只”需要将其转换为 3D 形状。然后可以在 3D 图形软件或游戏引擎中对其进行处理，使用 WebGL（或 WebGPU！）在浏览器中显示并渲染为图像。

有很多技术可以表示 3D 数据，这需要另一篇文章，如果你好奇的话可以看看 Open3D 教程。现在我们只说对于这个用例有两种方法可以做到这一点：

体素看起来很可爱，可以导出到 Magica Voxels 和 Goxel，或者适应 3D 打印格式以创建小区域的触觉地图。表示的简单性使得能够以最小的努力进行大量操作。
三角形网格本质上是游戏引擎想要的，非常灵活且通常具有高性能，但处理起来可能很棘手。
我最终选择了三角形网格表示，因为使用 Open3D（或者如果我们只考虑表面，则只是一些 NumPy）很容易将它们转换为体素，而相反的则很棘手。

三角形网格只是表示网格的一种方式，但可能是 Python 中支持最好的一种方式。除了 Open3D 之外，Trimesh 库还尝试成为“3D Shapely”来处理这种格式的拓扑。

通过此设置，我开始从 Shapely 对象生成 3D 网格，使用 Trimesh 以及后来的 Open3D 对其进行渲染和体素化。

3、导出格式

我在寻找 3D 库时考虑的一个重要功能是能够处理可在其他工具中使用的格式。

情况很复杂，有许多相互竞争的格式，许多是专有的和/或记录不良的。我最终得到了两个候选者：ply 和 glTF 2。Ply 很有趣，因为它非常直接、简单，人们可以手动检查文件或在没有外部库的情况下生成它们。就这么简单，它可以在许多软件中使用；一个缺点是它不能包含动画或着色器甚至纹理等花哨的元素（或者更好的是，它可以使用尚未完全标准化的额外属性）。

glTF 2 是 Khronos 集团最新（2017 年）的标准，免版税且有文档记录，灵活且受大多数软件支持。文本形式使用易于检查的 JSON。 glTF 受到 Blender、Unreal、Unity 和 Godot 游戏引擎等的支持。当我正在研究这个 Godot 时，发布了 4.0 版本的 alpha，该版本在导入 glTF 网格时崩溃，我报告了它，并且在不到一天的时间内修复了它，真的令人印象深刻。

glTF 还受到 Three.js 和 Babylon.js 的支持，这意味着不仅可以在浏览器中可视化生成的模型（这是一个更容易跨平台部署的解决方案），而且有几个可视化工具对于像我这样的 3d n00bs 很有用调试生成的模型并对过程进行故障排除。

4、Open3D

我使用 Python 来做这件事，原因很简单，它是一种我熟悉的语言，并且有一个很好的库生态系统来处理数字和处理几何图形。鉴于 Python 的速度相对较低，3D 渲染和游戏并不是 Python 的典型用例，但该语言在数据分析和机器学习方面的成功仍然导致许多库被开发来处理 LiDAR 数据或 CT 扫描或执行分割和分割等任务。立体视觉。

我最终选择了 Open3D，因为它看起来非常活跃，支持多种表示和算法，具有包含示例和教程的不错的文档，可以以无头模式（例如从 Web 后端）以 glTF 和其他方式导出，只需最少的设置，甚至可以渲染为静态这对于自动化测试和故障排除很有用，更不用说创建视频或进行合成的可能性了。

Trimesh 因其简单性引起了我的注意，能够通过一一指定坐标来创建网格有助于原型设计，此外它还提供了许多具有合理依赖关系的开箱即用算法，并且可以使用 glTF 或直接传递 NumPy 与 Open3D 集成数组。

5、“神圣”的问题

给定一个 2D 建筑，我希望初学者能够生成 3D 棱镜。使用 Trimesh 这意味着枚举顶点的 3D 坐标，然后枚举 3 元组中的索引来指定所有三角形。

这可能是地图中最常见的对象类型之一，并且可以具有复杂的拓扑。建筑物可能是凹形的并且有孔。

所以我从这个开始：

在这里插入图片描述

有洞的建筑物的真实例子

事实上，该形状有孔（欧拉数 -6）并且是凹的，这意味着要构建它，必须忽略一些三角形，因为它们最终完全位于体积内部，并重建每个三角形的内表面（也称为法线）。三角形并不简单。

即使网格在几何上是错误的，通常也可以渲染，但如果它有效，则可以在其之上应用许多有趣的变换。

没有孔的网格称为水密网格，而修剪网格允许使用 is_watertight 属性轻松检查此属性。请注意，这与整个网格的拓扑无关，例如，环面的欧拉数为 0（Trimesh 使用 euler_number 属性显示它），但如果网格定义正确，它是无懈可击的，并且有一个定义的内部和外部体积。

6、残破的大楼

其代码在这里，我做了简化以硬编码为 wkb 的对象，通常会与 PostGIS 中的许多其他对象一起检索以重现整个区域。

第一步是将上面的形状（有 4 个孔的建筑物）变换为不重叠的三角形。这可以通过 Shapely 轻松完成：

from shapely.ops import triangulate

...

with open("check_all_triangles.svg", "w") as fw:
    fw.write(shapely.geometry.MultiPolygon(triangulate(s))._repr_svg_())

在这里插入图片描述

上面的多边形被分成三角形

然后需要使用 inside 函数分离多边形内部的三角形边：

mul_holes = shapely.geometry.MultiPolygon(
    [tri for tri in triangulate(s) if not tri.within(s)]
)

这是为了生成单个 MultiPolygon 以用于显示目的，在上面的要点中你可以看到网格创建不需要它。

在这里插入图片描述

三角形边从外部不可见，生成网格时将被忽略
在这里插入图片描述

外部三角形边，供 3D 网格使用

你会注意到内部部分缺少一侧。这稍后会产生一个问题，并且在代码的第一次迭代中，当没有生成这个有用的图像时，我没有注意到这一点。

现在，对于每个三角形，我在网格上为棱镜的两个面生成两个三角形。在这里，我需要检查三角形（而不是单边）是否是形状的一部分。

for tri_num, tri in enumerate(internal_triangles):
    # it's 4 coordinates, the last is identical to the first to form a circuit
    for x, y in tri.exterior.coords[:-1]:
        vertices.append([x, y, 0.0])
    for x, y in tri.exterior.coords[:-1]:
        vertices.append([x, y, HEIGHT])
    # now 6 points have been added
    idx = tri_num * 6
    # bottom and top faces of the prism, use as is
    faces.append([idx, idx + 1, idx + 2])
    faces.append([idx + 3, idx + 4, idx + 5])

trimesh 使用列表中顶点的索引来定义三角形面。为了简单起见，我使用列表，但内部一切都基于 Numpy，这是合并过程后使用的理想格式，它的效率要高得多，尤其是在执行向量运算时。

定义两个面（本质上是在不同高度重复的原始二维形状）后，我们需要进行横向面。

这是创建它们的逻辑

for i1, i2 in ((0, 1), (1, 2), (2, 0)):
    if LineString([tri.exterior.coords[i1], tri.exterior.coords[i2]]).within(s):
        continue
    # add the two triangles making up the exterior face
    # there are 2 equivalent ways to split it, just use one
    faces.append([idx + i1, idx + i2, idx + i1 + 3])
    faces.append([idx + i1 + 3, idx + i2 + 3, idx + i2])

为三角形的每条边创建一个 LineString 来检查它从外部是否可见，如果是，则生成两个三角形以形成连接顶边和底边的矩形。

7、（无聊的）结果

代码的第一个版本造成的混乱。旋转它就会闪烁
在这里插入图片描述

当 Pyglet 渲染的网格物体具体化时，我惊恐地看着这一团糟。老实说，我没想到它已经可以工作了，并且非常惊讶它与预期的形状有模糊的相似之处。

将网格导出到 glTF 中并在各种在线可视化工具中查看，确认它已完全损坏。我尝试通过查看单个元素的坐标并切换到更简单的形状来检查这个问题，但这并不是一件容易的事。

像上面这样的 Web 可视化工具显示线框是正确的：三角形位于它们应该在的位置，体积或孔内没有三角形（当然有一个，但这是一个单独的问题），但 trimesh 报告欧拉数接近-50，结果无法使用！

当我观察到一个关键细节时，我几乎要放弃这个项目并给自己做一个Carbonara来安慰自己：当旋转固体时，每个三角形的出现或消失取决于我观察到的一侧。

8、法线和边的顺序

在某种程度上，这是有道理的。在 3D 图形中，有一个称为法线的概念。法线只是一个向量，定义了面指向的位置、哪一侧是内侧还是外侧。当使用 Blender 或高级库时，这是自动完成的，尽管在著名的 Blender 甜甜圈教程中提到过（我强烈推荐）。我天真地认为法线对于这样一个简单的网格来说并不重要，但我错了。检查库的代码并进行一些实验，我发现法线是按顶点定义的顺序（顺时针或逆时针）定义的。从几何角度来看，三角形保持不变。

Trimesh 提供了根据相邻面重新计算法线的帮助程序，但它们基于启发式方法并产生稍微不那么混乱的结果，因此必须在首先生成网格时对其进行修复，这是内部和外部有明确定义的唯一时刻。

遗憾的是，并没有一个所有库都使用的通用标准，但 OpenGL 假设逆时针意味着面向前方。

改变顶点的顺序确实产生了更好的结果。

在这里插入图片描述