TowardsDataScience 2023 博客中文翻译（三）

原文：TowardsDataScience

协议：CC BY-NC-SA 4.0

室内建模的 3D 点云形状检测

原文：towardsdatascience.com/3d-point-cloud-shape-detection-for-indoor-modelling-70e36e5f2511

动手教程，3D Python

10 步 Python 指南，用于自动化 3D 形状检测、分割、聚类和体素化，以实现室内点云数据集的空间占用 3D 建模。

Florent Poux, Ph.D.

·发布于 Towards Data Science ·阅读时长 28 分钟·2023 年 9 月 7 日

–

如果你有点云或数据分析的经验，你知道识别模式有多重要。识别具有相似模式的数据点或“对象”对获取有价值的见解至关重要。我们的视觉认知系统可以轻松完成这项任务，但通过计算方法复制这种人类能力是一个重大挑战。

目标是利用人类视觉系统自然的元素分组倾向。👀

3D 点云分割阶段结果示例。© F. Poux

但这有什么用呢？

首先，它通过将数据分组到不同的段中，让你轻松访问和处理数据的特定部分。其次，通过查看区域而不是单个点，它使数据处理更快。这可以节省大量时间和精力。最后，分割可以帮助你发现通过查看原始数据无法看到的模式和关系。🔍 总的来说，分割对于从点云数据中获取有用信息至关重要。如果你不确定如何做，别担心——我们会一起搞定的！🤿

策略

在我们以有效的解决方案处理项目之前，让我们框定整体方法。本教程遵循一个由十个简单步骤组成的策略，如下方的策略图所示。

本指南展示了 3D 点云室内建模的工作流程。© F. Poux

策略已列出，下面你可以找到不同步骤的快捷链接：

Step 1\. Environment Setup
Step 2\. 3D Data Preparation
Step 3\. Data Pre-Processing
Step 4\. Parameter Setting
Step 5\. RANSAC Planar Detection
Step 6\. Multi-Order RANSAC
Step 7\. Euclidean Clustering Refinement
Step 8\. Voxelization Labelling
Step 9\. Indoor Spatial Modelling
Step 10\. 3D Workflow Export

既然我们已做好准备，那就直接开始吧。

🎵读者注意：此动手指南是* UTWENTE 联合工作的一部分，由 F. Poux 和 V. Lehtola 共同作者。我们感谢来自数字双胞胎 @ITC 项目的资助，该项目由特温特大学 ITC 系提供。所有图像 © Florent Poux。

1. 设置您的 Python 环境。

在下面的上一篇文章中，我们展示了如何快速设置 Anaconda 环境并使用 IDE JupyterLab 管理代码。如果您继续以这种方式设置自己，您将成为一名完整的 Python 应用程序开发者 😆。

## 3D Python 工作流程用于 LiDAR 城市模型：逐步指南

解锁 3D 城市建模应用程序的精简工作流程的终极指南。教程涵盖了 Python…

towardsdatascience.com

🦊 Florent：我强烈推荐使用桌面设置或 IDE，避免使用 Google Colab，尤其是当您需要使用提供的库可视化 3D 点云时，因为它们在最佳情况下会不稳定，最糟糕的情况下则无法工作（不幸的是…）。

🤠 Ville：我们猜您是在 Windows 上运行？这没问题，但如果您想进入计算方法领域，Linux 是首选！

好吧，我们将采取一种“快速成果”的方法。实际上，我们将通过遵循最简化的编码方法来实现卓越的分割💻。这意味着我们对底层库非常挑剔！我们将使用三个非常强大的库，分别是numpy、matplotlib和open3d。

好的，要在一个全新的虚拟环境中安装上述库包，我们建议您从cmd终端运行以下命令：

conda install numpy
conda install matplotlib
pip install open3d==0.16.0

🍜 免责声明：我们选择了 Python，而不是 C++ 或 Julia，所以性能就是那样 😄。希望它能满足您的应用需求 😉，对于我们所谓的“离线”过程（非实时）。

在您的 Python IDE 中，请确保导入将被频繁使用的三个库：

#%% 1\. Library setup
import numpy as np
import open3d as o3d
import matplotlib.pyplot as plt

就这样！我们准备好进行室内点云建模工作流程了！

2. 点云数据准备

在之前的教程中，我们演示了如何处理点云并在使用 AHN4 LiDAR 活动获取的 3D 数据集上进行网格化。这一次，我们将使用一个通过地面激光扫描仪收集的数据集：ITC 新建筑。它被组织成三个不同的集合供您实验。紫色的是户外部分。红色是地面层，绿色是第一层，如下图所示。

使用的 ITC 数据集的 3D 点云部分。© F. Poux

你可以从这里下载数据：ITC 数据集。一旦你在本地对数据有了牢固的掌握，你可以用两行简单的代码将数据集加载到你的 Python 执行环境中：

DATANAME = "ITC_groundfloor.ply"
pcd = o3d.io.read_point_cloud("../DATA/" + DATANAME)

🦊 Florent: 这段 Python 代码片段使用 *Open3D* 库来读取位于目录“*../DATA/*”中的点云数据文件 *ITC_groundfloor.ply* ，并将其赋值给变量 *pcd*。

变量现在保存着你的点云 pcd，你将用以下代码进行操作：

一旦点云数据成功加载到 Open3D 中，下一步是应用各种预处理技术来提升其质量并提取有意义的信息。

3. 点云预处理

如果你打算在 open3d 中可视化点云，建议你将点云移动以绕过大坐标的近似，这样可以避免产生晃动的可视化效果。要对你的 pcd 点云应用这种移动，首先获取点云的中心，然后通过从原始变量中减去它来进行平移：

pcd_center = pcd.get_center()
pcd.translate(-pcd_center)

现在可以通过以下代码行进行交互式可视化：

o3d.visualization.draw_geometries([pcd])

🦚 注意: o3d.visualization.draw_geometries([pcd]) 调用了 *draw_geometries()* 函数，这个函数属于 Open3D 的 *visualization* 模块。该函数接受一个几何体列表作为参数，并在可视化窗口中显示它们。在这种情况下，列表包含一个几何体，即表示点云的 *pcd* 变量。 *draw_geometries()* 函数创建一个 3D 可视化窗口并渲染点云。你可以与可视化窗口互动以旋转、缩放，并从不同角度探索点云。

太棒了👌，我们已经准备好测试一些采样策略来统一我们的下游处理。

3.1. 点云随机采样

让我们考虑那些可以有效减少点云大小同时保持整体结构完整性和代表性的方法。如果我们将点云定义为一个矩阵 (m x n)，那么通过保留该矩阵每隔 n 行的一个行，我们可以得到一个简化的点云：

一个点云简化策略。© F. Poux

在矩阵级别，简化操作是通过根据 n 因子保留每隔 n 行的点来进行的。当然，这取决于点在文件中的存储方式。用 open3d 切片点云是相当直接的。为了简化和参数化表达式，你可以写出以下代码：

retained_ratio = 0.2
sampled_pcd = pcd.random_down_sample(retained_ratio)

🦚 注意：sampled_pcd = pcd.random_down_sample(retained_ratio) 对原始点云 *pcd* 应用随机下采样，使用的是 random_down_sample() 函数，该函数由 Open3D 提供。 retained_ratio 参数决定了下采样后保留的点的比例。例如，如果 retained_ratio 设置为 0.5，则大约 50% 的点将被随机选择并保留在采样后的点云中。

o3d.visualization.draw_geometries([sampled_pcd], window_name = "Random Sampling")

点云的下采样结果。© F. Poux

🌱 发展：在研究 3D 点云时，随机采样有其局限性，可能会导致遗漏重要信息和分析不准确。它没有考虑点之间的空间组件或关系。因此，使用其他方法以确保更全面的分析是至关重要的。

尽管这种策略快速，但随机采样可能不适合“标准化”用例。下一步是通过统计异常值去除技术来处理潜在的异常值，以确保数据的质量和可靠性，以便进行后续分析和处理。

3.2. 统计异常值去除

对 3D 点云数据使用异常值过滤器可以帮助识别和去除任何显著不同于数据集其余部分的数据点。这些异常值可能是由于测量误差或其他因素造成的，这些因素可能会影响分析。通过去除这些异常值，我们可以获得更有效的数据表示，并更好地调整算法。然而，我们需要小心不要删除有价值的点。

我们将定义一个 statistical_outlier_removal 过滤器，以去除与其邻居相比远离平均值的点。它需要两个输入参数：

• nb_neighbors，用于指定在计算给定点的平均距离时考虑多少个邻居。

• std_ratio，允许根据点云中平均距离的标准差设置阈值水平。这个数值越低，过滤器的作用就越强。

这包括以下内容：

nn = 16
std_multiplier = 10

filtered_pcd, filtered_idx = pcd.remove_statistical_outlier(nn, std_multiplier)

🦚 注意：*remove_statistical_outlier()* 函数通过指定数量的最近邻（*nn*）和一个标准差乘数（*std_multiplier*）来对点云应用统计异常值去除。该函数返回两个值：* *filtered_pcd* 和 *filtered_idx*。 *filtered_pcd* 表示经过滤的点云，其中已去除统计异常值。 *filtered_idx* 是一个索引数组，对应于在异常值去除过程中保留的原始点云 *pcd* 中的点。

为了可视化这种过滤技术的结果，我们将异常值标记为红色，并将其添加到我们希望可视化的点云对象列表中：

outliers = pcd.select_by_index(filtered_idx, invert=True)
outliers.paint_uniform_color([1, 0, 0])
o3d.visualization.draw_geometries([filtered_pcd, outliers])

异常值在点云中用红色标记和可视化。© F. Poux

在从点云中移除统计异常值之后，下一步是应用基于体素的采样技术，进一步下采样数据，以便高效处理和分析，同时保留点云的基本结构特征。

3.3. 点云体素（网格）采样

网格下采样策略基于将 3D 空间划分为规则的立方体单元，称为体素。对于这个网格的每个单元，我们只保留一个代表性点，并且可以用不同的方式选择这个点。当下采样时，我们保留该单元最接近重心的点。

体素网格采样示例。© F. Poux

具体而言，我们定义一个voxel_size，然后用它来过滤我们的点云：

voxel_size = 0.05
pcd_downsampled = filtered_pcd.voxel_down_sample(voxel_size = voxel_size)

🦚 注意：此行对过滤后的点云执行体素下采样， *filtered_pcd*，使用 *voxel_down_sample()* 函数。 *voxel_size* 参数指定体素下采样点云的每个体素的大小。较大的体素尺寸会导致点云密度显著减少。下采样操作的结果分配给变量 *pcd_downsampled*，代表下采样后的点云。

现在是仔细可视化下采样后分布的时间：

o3d.visualization.draw_geometries([pcd_downsampled])

应用在点云上的采样策略结果。© F. Poux

在此阶段，我们有一组未在进一步处理时触及的异常点和一个下采样后的点云，构成了后续过程的新主题。

3.4. 点云法线提取

点云法线指的是 3D 点云中特定点处表面的方向。它可以用于通过将点云划分为具有相似法线的区域来进行分割。例如，在我们的案例中，法线将有助于识别点云中的物体和表面，使其更容易可视化。这也是介绍一种半自动计算法线的方法的绝佳机会。我们首先定义点云中每个点与其邻居之间的平均距离：

nn_distance = 0.05

然后，我们利用这些信息在半径为radius_normals的范围内提取有限的max_nn个点，为 3D 点云中的每个点计算法线：

radius_normals=nn_distance*4
pcd_downsampled.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=radius_normals, max_nn=16), fast_normal_computation=True)

现在，pcd_downsampled点云对象很高兴地拥有法线，准备展示其最漂亮的一面😊。此时你知道该做什么：

pcd_downsampled.paint_uniform_color([0.6, 0.6, 0.6])
o3d.visualization.draw_geometries([pcd_downsampled,outliers])

用于后续实验的点云。

完成点云的体素下采样后，下一步是配置点云形状检测和聚类的参数，这在将相似点分组以及从下采样点云数据中提取有意义的结构或对象中起着关键作用。

4. 点云参数设置

在本教程中，我们为您精选了两种最有效且可靠的 3D 形状检测和聚类方法：RANSAC 和使用 DBSCAN 的欧几里得聚类。然而，在使用这些方法之前，了解参数的同时，理解基本概念是至关重要的。

RANSAC

RANSAC 算法，RANdom SAmple Consensus 的缩写，是一个强大的工具，用于处理包含异常值的数据，这在使用现实世界传感器时常常会遇到。该算法通过将数据点分为两类：内点和外点来工作。通过识别并忽略外点，您可以专注于处理可靠的内点，从而使分析更加有效。

在 3D 点云中使用 RANSAC 进行平面检测。© F. Poux

让我用一个简单的例子来说明 RANSAC 的工作原理。假设我们想通过下面的点云拟合一个平面。我们怎么做呢？

在随机点云中进行的 RANSAC 平面检测模拟。© F. Poux

首先，我们从数据中创建一个平面，为此，我们从点云中随机选择 3 个点来建立平面。然后，我们简单地检查剩余的点中有多少点落在该平面上（达到某个阈值），这将为提案打分。

RANSAC 评分系统说明。您可以看到，每次迭代都会随机抽取 3 个点，从中创建一个平面，然后选择落在平面上的点。在这里，第 159 次迭代会是最佳候选。© F. Poux

然后，我们用 3 个新的随机点重复这个过程，看看效果如何。是更好了吗？更差了吗？然后，我们反复进行这个过程，比如说 10 次、100 次、1000 次，然后选择得分最高的平面模型（即具有最佳“支持”的剩余数据点）。这将是我们的解决方案：支持点加上我们采样的三个点构成我们的内点集，其余部分是我们的外点集。够简单吧 😁？

哈哈，对于那些怀疑者，你们难道没有一个上升的问题吗？我们如何实际确定应该重复多少次过程？我们应该多频繁地尝试？好吧，这实际上是可以计算的，但我们先把它放到一边，专注于当前的问题：点云分割😉。

RANSAC 参数设置

我们希望使用 RANSAC 来检测点云中的 3D 平面形状。使用 RANSAC 时，我们需要定义三个参数：一个距离阈值distance_threshold，用于将点标记为内点或外点；一个最小点数ransac_n，用于拟合几何模型；一个迭代次数num_iterations。

距离阈值的确定：我们可能会因为需要一些领域知识来设置未来的分割和建模阈值而受到限制。因此，尝试绕过这一点以使方法对非专家开放将是令人兴奋的。我们将与你分享一个简单的想法，这可能会有用。如果我们计算数据集中点之间的平均距离，并将其作为设置阈值的基础，会怎样呢？

好吧，这是一个值得探索的想法。为了确定这样的值，我们使用 KD-Tree 来加速每个点的最近邻查询过程。从那里开始，我们可以查询点云中每个点的 k 个最近邻，然后将其打包到下面显示的 open3d 函数中：

nn_distance = np.mean(pcd.compute_nearest_neighbor_distance())

毫无意外，它接近 5 cm，因为我们将点云采样到这个值。这意味着如果我们考虑最近邻，我们会有 51 mm 的平均距离。

🌱 成长：从这里开始，我们可以将 RANSAC 参数设置为从 nn_distance 变量得出的值，这样可以根据考虑的数据集进行调整，而不依赖于领域知识。你会怎么处理这个问题？

点数的确定：这里很简单。我们想要找出平面，因此我们将取定义 3D 平面所需的最小点数：3。

迭代次数的确定：迭代次数越多，你的 3D 形状检测就越稳健，但所需时间也越长。目前，我们可以将其设置为 1000，这样可以获得良好的结果。我们将在未来的教程中探索更多巧妙的方法来找出点云的噪声比例。😉

🧙‍♂️ 专家：有一种自动获取每次迭代次数的方法。如果我们希望以概率 p（例如 99%）成功，数据中的离群点比例是 e（例如 60%），我们需要 s 个点来定义我们的模型（这里是 3）。以下公式给出了预期的迭代次数：

现在我们已经定义了 RANSAC 参数，我们可以研究第一次分割过程。

5. 使用 RANSAC 进行点云分割

首先将不同的阈值设置为非自动值进行测试：

distance_threshold = 0.1
ransac_n = 3
num_iterations = 1000

从那里开始，我们可以使用 RANSAC 对点云进行分割以检测平面，具体如下：

plane_model, inliers = pcd.segment_plane(distance_threshold=distance_threshold,ransac_n=3,num_iterations=1000)
[a, b, c, d] = plane_model
print(f”Plane equation: {a:.2f}x + {b:.2f}y + {c:.2f}z + {d:.2f} = 0")

我们将结果收集到两个变量中：plane_model，它包含平面的参数 a,b,c,d，以及 inliers 作为点索引。

这使得可以使用索引将点云分割为 inlier_cloud 点集（我们将其标记为红色）和 outlier_cloud 点集（我们将其标记为灰色）：

inlier_cloud = pcd.select_by_index(inliers)
outlier_cloud = pcd.select_by_index(inliers, invert=True)

#Paint the clouds
inlier_cloud.paint_uniform_color([1.0, 0, 0])
outlier_cloud.paint_uniform_color([0.6, 0.6, 0.6])

#Visualize the inliers and outliers
o3d.visualization.draw_geometries([inlier_cloud, outlier_cloud]) 

🦊 Florent ：*invert=True* 参数允许选择第一个参数的反义，即所有不在 *inliers* 中的索引。如果你缺少阴影，记得计算法线，如上所示。*

🌱 成长：尝试调整各种参数，并通过定性分析研究其影响。首先记住要去相关变化（一次一个变量）；否则你的分析可能会有偏差 😊。

很棒！你知道如何将点云分割为内点集和外点集 🥳！现在，让我们研究如何找到彼此接近的一些簇。让我们想象一下，一旦我们检测到大的平面部分，我们有一些“漂浮”的物体需要 delineate。怎么做呢？（是的，这是一个假问题，我们有答案给你 😀）

6. 扩展 3D 分割：多阶 RANSAC

我们的理念将非常简单。我们将首先多次运行 RANSAC（假设 n 次）以提取构成场景的不同平面区域。然后，我们将通过欧几里得聚类（DBSCAN）处理“漂浮元素”。这意味着我们必须确保有一种方法在迭代过程中存储结果。准备好了吗？

好的，让我们实例化一个空字典，以保存迭代结果（segment_models中的平面参数和segments中的点云平面区域）：

segment_models={}
segments={}

然后，我们需要确保可以影响以后迭代的次数，以检测平面。为此，让我们创建一个变量max_plane_idx来保存迭代次数：

max_plane_idx=10

🦊 Florent*：在这里，我们说我们希望迭代 10 次以找到 10 个平面，但有更聪明的方法来定义这样的参数。这实际上扩展了文章的范围，将在另一节中讨论。*

现在让我们进入一个工作循环 😁，我将首先快速说明。

执行以在 RANSAC 过程中进行分割的循环。© F. Poux

在第一次迭代（循环i=0）中，我们将内点与外点分开。我们将内点存储在segments中，然后我们只对存储在rest中的剩余点感兴趣，这将成为循环 n+1（循环i=1）的研究对象。这意味着我们希望将上一阶段的外点作为基础点云，直到达到上述的迭代阈值（与 RANSAC 迭代不同）。这转化为以下内容：

rest=pcd
for i in range(max_plane_idx):
    colors = plt.get_cmap("tab20")(i)
    segment_models[i], inliers = rest.segment_plane(
    distance_threshold=0.1,ransac_n=3,num_iterations=1000)
    segments[i]=rest.select_by_index(inliers)
    segments[i].paint_uniform_color(list(colors[:3]))
    rest = rest.select_by_index(inliers, invert=True)
    print("pass",i,"/",max_plane_idx,"done.")

就是这样！现在，为了可视化整体，我们通过循环中的第一行将每个检测到的分段涂上来自tab20的颜色（colors = plt.get_cmap("tab20")(i)），你只需写：

o3d.visualization.draw_geometries([segments[i] for i in range(max_plane_idx)]+[rest])

🦚 注意：我们传递给函数o3d.visualization.draw_geometries()的列表[segments[i] for i in range(max_plane_idx)]实际上是一个“列表推导式”🤔。它等同于编写一个for循环，将第一个元素segments[i]追加到列表中。方便的是，我们可以将[rest]添加到这个列表中，draw.geometries()方法会理解我们想要绘制一个点云。这不是很酷吗？

DBSCAN 对 3D 点云的第一次扫描结果

哈！我们以为完成了……但我们真的完成了吗？你注意到这里有什么奇怪的地方吗？如果你仔细看，会发现一些奇怪的伪影，比如实际切割一些平面元素的“红色线条/平面”。为什么？🧐

实际上，由于我们将所有点拟合到 RANSAC 平面候选者（在欧几里得空间中没有限制范围）而不考虑点的密度连续性，因此根据平面的检测顺序，我们会出现这些“线”伪影。所以下一步是防止这种行为！为此，我建议在迭代过程中包含一个基于欧几里得聚类的条件，以在连续的簇中细化内点集。为此，我们将依赖于 DBSCAN 算法。

欧几里得聚类（DBSCAN）

在点云数据集中，我们经常需要将空间上连续的点集合（即在 3D 空间中物理上接近或相邻）进行分组，如下所示。但我们如何有效地做到这一点呢？

在这张图中，很明显我们想要将彼此接近的点分组，找到 5 组点。© F. Poux

DBSCAN（基于密度的空间聚类应用与噪声）算法于 1996 年为此目的而提出。该算法被广泛使用，因此在 2014 年获得了经得起时间考验的科学贡献奖。

DBSCAN 算法涉及扫描数据集中的每个点，并基于密度构建一个可达点集合。这是通过分析每个点的邻域并在其包含足够点时将其包含在区域内来实现的。这个过程对每个邻近点重复，直到簇无法再扩展。没有足够邻居的点被标记为噪声，使得该算法对离群点具有鲁棒性。这不是很令人印象深刻吗？😆

DBSCAN 算法过程和两个参数ϵ和 min_points 对结果的影响的示意图。你可以看到，值越大，组成的簇就越少。© F. Poux

啊，我们差点忘了。参数选择（ϵ用于邻域，n_min 用于最小点数）也可能很棘手：设置参数时必须小心，以确保创建足够的内部点（如果 n_min 太大或ϵ太小，就不会发生）。特别是，这意味着 DBSCAN 在发现不同密度的簇时会遇到困难。但 DBSCAN 有一个很大的优势，就是计算效率高，不需要像 Kmeans 那样预定义簇的数量。最后，它允许发现任意形状的簇。现在我们准备深入探讨参数的黑暗面 💻

DBSCAN 用于 3D 点云聚类

让我详细说明逻辑过程（激活猛兽模式👹）。首先，我们需要定义运行 DBSCAN 的参数：

epsilon = 0.15
min_cluster_points = 5

🌱 成长：这些参数的定义是一个需要探索的问题。你必须找到一种方法来平衡过度分割和不足分割的问题。最终，你可以基于 RANSAC 的初始距离定义使用一些启发式方法。这是一个值得探索的方面 😉。

在我们之前定义的 for 循环中，我们将在分配内点（segments[i]=rest.select_by_index(inliers)）后立即运行 DBSCAN，方法是紧接着添加以下一行：

labels = np.array(segments[i].cluster_dbscan(eps=epsilon, min_points=min_cluster_points))

然后，我们将计算每个找到的簇包含多少个点，使用一种奇怪的符号表示法，该表示法利用了列表推导。结果存储在变量candidates中：

candidates=[len(np.where(labels==j)[0]) for j in np.unique(labels)]

那现在呢？我们需要找到“最佳候选者”，通常是包含最多点的簇！为此，这里是这行代码：

best_candidate=int(np.unique(labels)[np.where(candidates== np.max(candidates))[0]])

好吧，很多技巧在这里发生，但本质上，我们利用 Numpy 的熟练度来搜索并返回属于最大簇的点的索引。从这里开始，接下来的过程就简单了，我们只需确保在每次迭代中将任何剩余簇纳入后续 RANSAC 迭代中（🔥 推荐阅读 5 遍以消化！）：

rest = rest.select_by_index(inliers, invert=True) + segments[i].select_by_index(list(np.where(labels!=best_candidate)[0]))
segments[i]=segments[i].select_by_index(list(np.where(labels== best_candidate)[0]))

🦚 注意：rest 变量现在确保包含 RANSAC 和 DBSCAN 剩余的点。当然，内点现在被筛选为原始 RANSAC 内点集中最大的簇*。

当循环结束时，你会得到一组干净的段，这些段包含空间上连续的点集，符合平面形状，你可以使用以下代码行以不同颜色可视化：

colors = plt.get_cmap("tab20")(i)
segments[i].paint_uniform_color(list(colors[:3]))

结果应该类似于这个：

但这就结束了吗？不，不可能的 😄！一旦对点云应用了多阶 RANSAC 分割，下一阶段就是通过利用 DBSCAN 来细化剩余的未分割点，这有助于进一步提升点云分析的粒度，增加更多的簇。

🦚 注意：粒度是一个在数据科学中经常使用的学术词汇。数据的粒度意味着数据组织或建模的详细程度。在这里，我们可以从点云中建模的对象越小，我们的表示的粒度就越细。（很高大上，对吧？！）

7. 欧几里得聚类精炼

好了，是时候摆脱循环，处理分配给 rest 变量的剩余点，这些点尚未分配给任何段。让我们首先对我们讨论的内容有一个视觉上的了解：

o3d.visualization.draw_geometries([rest])

我们可以使用 DBSCAN 进行简单的欧几里得聚类，并将结果捕获到一个 labels 变量中。你知道怎么做的：

labels = np.array(pcd.cluster_dbscan(eps=0.1, min_points=10))

🌱 生长： 我们使用 10 cm 的半径来“生长”簇，只有在此步骤之后至少有 10 个点时才考虑一个簇。但这是正确的选择吗？ 🤔可以随意实验以找到一个好的平衡点，理想情况下，找到一种自动化的方法😀。

标签的范围在 -1 和 n 之间，其中 -1 表示“噪声”点，0 到 n 的值则是分配给相应点的簇标签。请注意，我们希望将标签作为 NumPy 数组获取。

左侧：参数定义不佳。右侧，我们对物体的轮廓有了更好的划分，以便进行后续处理。

很好。现在我们已经定义了每个点的标签组，让我们为结果着色。这是可选的，但对于迭代过程以搜索合适的参数值非常有用。为此，我们建议使用 Matplotlib 库获取特定的颜色范围，例如 tab20：

max_label = labels.max()
colors = plt.get_cmap("tab20")(labels / (max_label 
if max_label > 0 else 1))
colors[labels < 0] = 0
rest.colors = o3d.utility.Vector3dVector(colors[:, :3])
o3d.visualization.draw_geometries([rest])

🦚 注意：max_label 应该是直观的：它存储标签列表中的最大值。这允许将其用作着色方案的分母，同时处理一个特殊情况，其中聚类被扭曲，仅产生噪声+一个簇。之后，我们确保将这些噪声点的标签设置为 -1 （黑色（0）*）。然后，我们将点云 pcd 的 colors 属性设置为表示 R、G、B 的 3 列的 2D NumPy 数组。

就这样！我采用了与之前相同的方法，没有任何魔法！我只是确保使用一致的参数来进行精细的聚类，以获得你一直梦想的美丽彩虹场景 🥳！

在使用 DBSCAN 精炼点云聚类结果之后，焦点转向体素化技术，这涉及将数据组织成有意义的空间结构，从而实现对点云信息的高效建模。

8. 体素化和标记

现在我们有了一个带有分段标签的点云，查看我们是否能适应室内建模工作流将会非常有趣。一种方法是使用体素来适应 O-Space 和 R-Space。通过将点云划分为小立方体，可以更容易地理解模型的占用和空白空间。让我们深入了解一下。

9.1. 体素网格生成

创建准确且详细的空间 3D 模型意味着生成一个紧凑的体素网格。这种技术将点云划分为具有自己坐标系统的小立方体或体素。

体素网格生成以结构化带有分段信息的 3D 点云。© F. Poux

要创建这样的结构，我们首先定义我们新实体的大小：体素：

voxel_size=0.5

现在，我们想知道我们需要堆叠多少个这样的立方体才能填充由点云定义的边界框。这意味着我们首先需要计算点云的范围：

min_bound = pcd.get_min_bound()
max_bound = pcd.get_max_bound()

很好，现在我们使用o3d.geometry.VoxelGrid.create_from_point_cloud()函数对任何选择的点云进行操作，以便在其上拟合体素网格。但是等等，我们想区分哪个点云以便进一步处理？

好的，让我们举例说明你想拥有“结构性”元素的体素与不属于结构性元素的杂物体素的情况。没有标签的情况下，我们可以根据它们是否属于 RANSAC 分段或其他分段来指导我们的选择；这意味着首先将 RANSAC 处理的分段进行拼接：

pcd_ransac=o3d.geometry.PointCloud()
for i in segments:
 pcd_ransac += segments[i]

🦚 注意: 在这个阶段，你有能力用体素稍后拾取的均匀颜色为点云着色。如果这是你想要的，你可以使用： pcd_ransac.paint_uniform_color([1, 0, 0])

然后，我们可以简单地提取我们的体素网格：

voxel_grid_structural = o3d.geometry.VoxelGrid.create_from_point_cloud(pcd_ransac, voxel_size=voxel_size)

并对剩余元素做同样的处理：

rest.paint_uniform_color([0.1, 0.1, 0.8])
voxel_grid_clutter = o3d.geometry.VoxelGrid.create_from_point_cloud(rest, voxel_size=voxel_size)

最后一步是可视化我们的结果：

o3d.visualization.draw_geometries([voxel_grid_clutter,voxel_grid_structural])

使用体素表示的空间的语义表示。© F. Poux

这太棒了！它看起来像那些整个世界由立方体构成的计算机游戏！而且我们实际上可以使用分段标签来指导我们的体素建模！这开启了许多视角！但问题依然存在。使用 Open3D，提取未填充的体素是困难的；因此，完成体素化点云的结构化后，下一步涉及探索体素空间建模技术，以提供分析体素化数据空间关系和属性的替代视角，为高级点云建模开辟新途径。

🤠 Ville: 体素化是对相同点云数据的不同空间表示。如果机器人需要避开碰撞，这种表示非常有用！比如在模拟…火灾演习时也非常有用！

9. 空间建模

在室内建模应用中，基于体素的点云表示在捕捉和分析复杂环境的几何属性方面起着关键作用。随着点云数据集的规模和复杂性的增加，深入探讨体素分割技术变得至关重要，以提取有意义的结构并促进更高层次的分析。因此，让我们定义一个函数来拟合体素网格，并返回填充和空白区域：

def fit_voxel_grid(point_cloud, voxel_size, min_b=False, max_b=False):
# This is where we write what we want our function to do.
return voxel_grid, indices

在函数内部，我们将 (1) 确定点云的最小和最大坐标，(2) 计算体素网格的尺寸，(3) 创建一个空的体素网格，(4) 计算已占据体素的索引以及 (5) 将占据体素标记为 True。

从点云创建占据网格的算法工作流程。© F. Poux

这转化为我们希望在此函数中包含的以下代码：

# Determine the minimum and maximum coordinates of the point cloud
 if type(min_b) == bool or type(max_b) == bool:
   min_coords = np.min(point_cloud, axis=0)
   max_coords = np.max(point_cloud, axis=0)
 else:
   min_coords = min_b
   max_coords = max_b
 # Calculate the dimensions of the voxel grid
 grid_dims = np.ceil((max_coords - min_coords) / voxel_size).astype(int)
 # Create an empty voxel grid
 voxel_grid = np.zeros(grid_dims, dtype=bool)
 # Calculate the indices of the occupied voxels
 indices = ((point_cloud - min_coords) / voxel_size).astype(int)
 # Mark occupied voxels as True
 voxel_grid[indices[:, 0], indices[:, 1], indices[:, 2]] = True

我们已经定义了函数，现在可以使用它提取体素，按结构、杂乱和空体素进行分段：

voxel_size = 0.3

#get the bounds of the original point cloud
min_bound = pcd.get_min_bound()
max_bound = pcd.get_max_bound()

ransac_voxels, idx_ransac = fit_voxel_grid(pcd_ransac.points, voxel_size, min_bound, max_bound)
rest_voxels, idx_rest = fit_voxel_grid(rest.points, voxel_size, min_bound, max_bound)

#Gather the filled voxels from RANSAC Segmentation
filled_ransac = np.transpose(np.nonzero(ransac_voxels))

#Gather the filled remaining voxels (not belonging to any segments)
filled_rest = np.transpose(np.nonzero(rest_voxels))

#Compute and gather the remaining empty voxels
total = pcd_ransac + rest
total_voxels, idx_total = fit_voxel_grid(total.points, voxel_size, min_bound, max_bound)
empty_indices = np.transpose(np.nonzero(~total_voxels))

🦚 注意：NumPy 中的 *nonzero()* 函数查找 *ransac_voxels* 变量中非零元素的索引。*nonzero()** 函数返回一个数组的元组，每个数组对应于特定轴上非零元素的索引。然后，我们将 *np.transpose()** NumPy 函数应用于从 *np.nonzero(ransac_voxels)** 获得的结果。*transpose()** 函数最终会排列数组的轴（有效地交换行和列）。通过结合这些操作，代码行对 *ransac_voxels** 的非零元素索引进行转置，生成一个转置数组，其中每行表示原始 *ransac_voxels** 数组中非零元素的坐标或索引。😊

这种体素建模方法提供了对体素化数据空间关系和属性的宝贵见解。为了在 Python 之外以透明度可视化结果，我们需要导出数据。

占据网格匹配的结果。左侧是填充体素和空体素，中间是填充体素，右侧是空体素。© F. Poux

在实现体素化点云的结构化后，下一步涉及将点云和体素数据导出为外部格式，以便实现互操作性，并与其他软件工具和工作流程无缝集成。此导出过程确保结构化的体素数据和原始点云可以轻松共享、可视化或用于进一步分析，从而促进对点云数据利用的协作和多功能方法。

10. 导出 3D 数据集

让我们首先关注导出点云分割数据集。

10.1. 分割点云导出

要导出分割后的点云，我们必须确保可以在可读的 ASCII 文件中为每个点写入标签。为此，我们将创建一个 XYZ 片段列表，并将标签特征附加到该列表。这可以通过以下for循环完成：

xyz_segments=[]
for idx in segments:
 print(idx,segments[idx])
 a = np.asarray(segments[idx].points)
 N = len(a)
 b = idx*np.ones((N,3+1))
 b[:,:-1] = a
 xyz_segments.append(b)

从那里，我们不想忘记 DBSCAN 聚类中的剩余元素，对其应用相同的原则：

rest_w_segments=np.hstack((np.asarray(rest.points),(labels+max_plane_idx).reshape(-1, 1)))

最后，我们将其附加到片段列表中：

 xyz_segments.append(rest_w_segments)

然后我们所需要做的就是使用 numpy 导出数据集并进行外部可视化：

np.savetxt("../RESULTS/" + DATANAME.split(".")[0] + ".xyz", np.concatenate(xyz_segments), delimiter=";", fmt="%1.9f")

在成功导出分割后的点云数据集后，现在的重点是将体素数据集导出为.obj文件。

10.2. 体素模型导出

要导出体素，我们首先必须为每个体素生成一个立方体；这些立方体然后在voxel_assembly中组合在一起，堆叠生成最终文件。我们创建了voxel_modelling(filename, indices, voxel_size)函数来完成这个任务：

def voxel_modelling(filename, indices, voxel_size):
    voxel_assembly=[]
    with open(filename, "a") as f:
        cpt = 0
        for idx in indices:
            voxel = cube(idx,voxel_size,cpt)
            f.write(f"o {idx}  \n")
            np.savetxt(f, voxel,fmt='%s')
            cpt += 1
            voxel_assembly.append(voxel)
    return voxel_assembly

🦚 注意：函数 [cube()](https://drive.google.com/file/d/1kPu85YHl66gQH8Qumxlyd-Sp4PjgvBVm/view?usp=sharing) 读取提供的索引，根据索引和体素大小生成体素立方体，将体素立方体写入文件，并跟踪在 *voxel_assembly* 列表中的组装体素立方体，最终由函数返回。

然后用来导出三种不同的体素组合：

vrsac = voxel_modelling("../RESULTS/ransac_vox.obj", filled_ransac, 1)
vrest = voxel_modelling("../RESULTS/rest_vox.obj", filled_rest, 1)
vempty = voxel_modelling("../RESULTS/empty_vox.obj", empty_indices, 1)

分割和体素建模的结果。© F. Poux

• 💻 在这里获取代码访问权限：代码示例

• 🍇 在这里获取数据访问权限：3D 数据集

• 👨‍🏫3D 数据处理和 AI 课程：3D 学院

• 📖 订阅以获得 3D 教程的提前访问权限：3D AI 自动化

结论

🦊 Florent：大规模祝贺🎉！你刚刚学会了如何开发一个自动形状检测、聚类、体素化和室内建模程序，用于处理由数百万个点组成的 3D 点云，并采用不同的策略！真心地，做得很好！但道路显然不止于此，因为你刚刚解锁了一个巨大的智能过程潜力，能够在片段级别进行推理！

🤠 Ville：那么，你是否在想我们是否可以让 3D 建模成为完全无需人工干预的过程？凭借我们目前的技术，我们已经完成了一半。为什么？我们的技术可以找到，例如，平面模型的参数。这就是为什么学术界的人称之为参数建模。然而，我们仍然需要仔细选择其他一些参数，比如 RANSAC 的参数。我鼓励你尝试在不同的点云上应用你的代码！

基于地面语义的 3D 网格示例。© F. Poux

更进一步

学习之旅并未结束。我们的终身探索才刚刚开始，未来的步骤将深入到 3D Voxel 工作，探索语义和点云分析与深度学习技术。我们将解锁高级的 3D LiDAR 分析工作流程。令人兴奋的事情还有很多！

1. Lehtola, V.，Nikoohemat, S.，& Nüchter, A.（2020）。室内 3D：扫描和重建方法概述。大地空间数据手册，55–97. doi.org/10.1007/978-3-030-55462-0_3

2. Poux, F.，& Billen, R.（2019）。基于体素的 3D 点云语义分割：无监督几何和关系特征与深度学习方法。ISPRS 国际地理信息杂志。8(5)，213；doi.org/10.3390/ijgi8050213 — Jack Dangermond 奖（新闻报道链接）

3. Bassier, M., Vergauwen, M., Poux, F.，（2020）。建筑内部分类中的点云与网格特征。遥感。12，2224. https://doi:10.3390/rs12142224

《LiDAR 城市模型的 3D Python 工作流程：一步步指南》

原文：towardsdatascience.com/3d-python-workflows-for-lidar-point-clouds-100ff40e4ff0

3D Python

Florent Poux, Ph.D.

·发表于 Towards Data Science ·38 分钟阅读·2023 年 4 月 4 日

–

解锁 3D 城市建模应用的精简工作流程的终极指南。教程涵盖了结合 3D 点云、网格和体素的 Python 自动化，以进行高级分析。

《LiDAR 城市模型的 3D Python 工作流程：一步步指南》。封面来自我的另一半Marina，展示了 3D 城市建模的艺术过程。© Mimatelier。

你之前遇到过“智能城市”这个词吗？或者“智能某物”？好吧，我们会涉及这个话题！将智能城市想象成一个超能的面包师 🥐：它知道你需要什么，甚至在你提出之前就会提供最直接的建议，帮助你做出美味的选择。不，这个智能城市的比喻并不是我今天唯一要分享的内容。确实，要达到这种“智能”的水平，我们首先得从基础层面入手：3D 城市模型。

如果你曾经想要创建令人惊叹的 3D 城市模型，但发现工作流程令人生畏且复杂，那么我可以帮忙！本文探讨了如何利用 Python 和开源软件定义一个强大的 3D 工作流程，以启动你的 3D 城市建模之旅。向繁琐的手动流程说（几乎）再见，迎接高效、动态且引人注目的创作吧！

我们深入探讨了一个四步策略，描述了环境设置、3D 数据策划与准备以及 3D 几何处理，以提取关键洞察，例如使用点云数据、网格、体素以及一些灰质 🧠，了解你所在社区的建筑覆盖情况。

《LiDAR 城市模型的 3D Python 工作流程：一步步指南》。© 作者

如果你已经准备好并充满热情，现在是时候开始 3D 编程了！

🎵致读者的说明：这本实践指南是* UTWENTE 与我亲爱的同事们 Dr. Sander Oude Elberink, Dr. Mila Koeva, Dr. Ville Lehtola, Dr. Pirouz Nourian, Dr. Paulo Raposo. 和 Prof G. Vosselman* 的共同工作之一。我们感谢来自数字双胞胎* @ITC 项目的资金支持，该项目由特温特大学 ITC 学院授予。

我们将在本指南中处理的 3D Python 数据集的摘录。© 作者

引言

在进入有趣的部分之前，让我讲一个小故事，为我们将要实现的目标提供一些背景。这始于一个基本的问题：什么是 3D 城市建模，它为何有用？

在一个城市化的世界里，3D 城市建模对于高效管理我们的日常生活至关重要。通过准确地以三维方式呈现我们的城市，我们可以分析和可视化复杂的城市环境，理解提议更改的影响，并做出明智的决策，以改善居民的生活质量。这是一个基本的概念，由那句优美的话很好地表达：城市塑造生活¹。

朝向智能城市以改善我们的生活。© 作者

的确，城市是人们生活、形成社区并建立自己身份的地方。它们是如市中心和郊区这样的空间，提供了一种配置和塑造物质世界和自然环境的方式。

想象一下，如果你有一种超级能力，可以对你的城市交通网络进行建模，预测交通模式并识别拥堵区域。这将如何改变你在城市中的生活方式？这只是作为城市“用户”的一个微小例子。但从根本上讲，3D 城市建模为城市规划师、建筑师和政策制定者提供了宝贵的见解，以优化城市基础设施、减少交通拥堵并提高公共安全。

通过创建我们城市的数字化复制品，我们可以更好地规划未来，为子孙后代创造可持续、宜居的社区。²⁻³

¹ Chen, X., Orum, A. M., & Paulsen, K. E. (2018). 《城市导论：地方和空间如何塑造人类体验》。约翰·威利父子公司。

² Lehtola, V. V., Koeva, M., Elberink, S. O., Raposo, P., Virtanen, J. P., Vahdatikhaki, F., & Borsci, S. (2022). 《城市数字双胞胎：服务城市需求的技术综述》。《国际应用地球观察与地理信息杂志》，102915。

³ Nourian, P., Ohori, K. A., & Martinez-Ortiz, C. (2018). 城市计算的基本手段：基于网络的城市规划计算平台的规范，旅行者指南。城市规划，3(1)，47–57

3D 城市建模：工作流程

是时候认真起来，定义一个一致的 3D 工作流程，我们可以将其作为不同 3D 城市建模操作的灵感。我们的目标是（1）易于设置和运行，（2）为各种场景提供极大的灵活性，以及（3）足够强大以涵盖多模态应用的复杂性。

🦚 注意：不，多模态不是脏话：它只是触及到当处理 3D 城市模型时，我们遇到各种需要考虑的地理空间数据模态。在本教程中，我们将专注于一个美丽的细分领域： 点云*，* 体素*，以及* 3D 网格*。

如果我们从高层次定义工作流程，我们会遵循四个主要步骤，如下所示。

在城市模型背景下的 3D Python LiDAR 工作流程。我们从环境设置（步骤 1）和 3D 数据准备（步骤 2）开始。一旦完成这些步骤，我们进入 Python 自动化（步骤 3），其中一个特定部分处理 3D Python 挑战（步骤 4），例如地块表面或兴趣点查询。

激动了吗？仔细查看流程，你会发现我们从零开始。这允许将提出的结构适应各种场景，这些场景需要不同的环境、数据集或自动化工具。

现在让我们更深入一点。假设我们在荷兰拥有一所房子（可以靠近特文特大学），我们想更好地了解周围的区域。这是我们的起点。现在，为了更好地理解我们房子与周围环境的关系，出现了一些问题：邻里建筑区的密度有多高？房子是否可能面临洪水？我是否遵守了我所拥有地块的建筑比例？

我们应该如何解决这些问题？我们是否应该上网查找？是否应该打开地图？是否应该联系测绘服务？让我们一起探索这个过程，同时在这个背景下解锁一套新的强大技能。

在接下来的内容中，我们详细介绍了回答这些问题的四个步骤：第一步是理想的环境设置。其次，我们获取 3D 点云和作为感兴趣区域网格的城市模型。然后，我们创建一个高度关注自动化的 3D Python 笔记本。最后，我们创建一套 Python 函数以应对具有挑战性的场景。好了，让我们拿起咖啡或茶🍵，深入寻找这些问题的答案吧！

🦚 注意：我设计了下一系列操作，使其易于线性跟随，无需编码技能或可用数据。然而，如果你是经验丰富的编码人员，我提供了有用的优化技巧，以确保代码性能达到巅峰！

第一步：3D 环境设置

在开始动手编写代码和进行 3D 思考之前，一个好的做法是确保我们在适当的环境中工作。我们不会在肮脏的台面上用钝刀和过期的食材做饭（或者至少我们会避免 😁）！让我们遵循下面所示的三个子步骤。

第一步：3D 环境设置。

1. 1. 软件堆栈

为了开始，让我们首先安装一个 3D 点云和网格处理软件：CloudCompare。它是一个出色的工具，允许有效地处理点云数据的科学分析（但不限于此）。它是任何迭代实验中的一个重要部分，例如，当你想快速获得关于理论可行性的数据驱动想法时。

要获取 CloudCompare 软件，你可以前往 CloudCompare.org 的下载部分，获取适合你操作系统的最新稳定版，如下所示。

从 cloudcompare.org/ 下载 CloudCompare

下载软件后，按照线性安装步骤操作，直到获得可以打开的工作软件，启动时界面应类似于以下 GUI：

CloudCompare GUI。

然后，我们希望获得一个允许我们专注于实际代码的 Python 发行版。它被称为 Anaconda，可以在各种平台上从 Anaconda.com 获得。下载与您的操作系统匹配的软件后，您可以进行安装。提供了一个 GUI（Anaconda Navigator），您可以启动它以快速开始，如下所示。

这是 Anaconda Navigator GUI，允许你管理未来实验的独立 Python 环境。

使用 GUI Anaconda Navigator，转到左侧的“环境”标签。然后，我们通过点击“创建”来创建一个全新的隔离 Python 环境，如下所示。

在 Anaconda Navigator 中，你的左侧有四个标签。在“首页”标签中，你会发现 IDE 处于特定环境中；在“环境”标签中，你可以创建、管理和选择任何环境。

这将使我们能够更好地管理一些我们想要安装的 Python “库”（第 1.3 步），同时避免版本冲突。

🦚 注意：对于本教程，我们选择了 Python 版本 3.9.16，如上所示。创建环境后，点击它，你会看到其名称旁边有一个“播放”图标。这将开启直接在所选环境中启动 Anaconda Terminal 的可能性。这是在此环境中安装库或 IDE 的首选方式。*

如果你已经有一个工作中的 Anaconda Navigator 和一个新创建的环境，我们就可以选择一种编码方式，旨在比文本编辑器更高效。😁

1.2. Python IDE

Python IDE（集成开发环境）是一种软件应用程序，提供了一整套用于开发、测试和调试 Python 代码的工具。它们提供了一个一体化的环境，我们可以在其中编写、编辑和执行代码，管理项目文件，跟踪更改，并与其他开发人员协作。一些流行的桌面 Python IDE 包括 PyCharm、Visual Studio Code 和 Spyder，每个 IDE 都提供了适应不同编程需求的独特功能和能力。

今天，我想重点介绍一个出色的“基于网络”的 IDE：JupyterLab。JupyterLab 的主要优点之一是其笔记本界面，它提供了一个可视化的、互动的编程环境。它使得实时创建、编辑和运行代码单元变得简单，而无需在不同窗口或应用程序之间切换。此外，JupyterLab 支持广泛的数据可视化工具，是处理任何地理数据科学或 3D 机器学习项目的绝佳选择。

JupyterLab IDE 用于 3D Python。

JupyterLab 直观的界面、强大的功能集以及对多种编程语言的支持，使其成为各个技能水平的 Python 开发人员的热门选择。

🦚 注意：JupyterLab IDE 还支持多种编程语言，包括 Python、R 和 Julia，允许在一个环境中使用各种工具和库。这非常酷，因为 R 和 Julia 是很棒的语言。

在能够使用 JupyterLab 之前，我们需要在当前 Anaconda 环境中安装它。如前所述，我们必须在所选环境中打开 Anaconda Terminal（在我们的例子中是 ITC）。要通过 Anaconda Navigator 实现这一点，从你选择的环境中，（1）点击绿色箭头，（2）选择 “Open Terminal”，如下所示。*

如何在当前环境中安装依赖。

在打开的控制台（Anaconda Terminal）中，输入以下命令：conda install -c conda-forge jupyterlab，然后按“Enter”。这将直接安装 JupyterLab。

命令行在 Anaconda Terminal 中执行（我们看到我们在（ITC）环境中）。

请注意，它可能会要求你确认安装一些需要的库，你需要通过输入y并按下“Enter”键来接受。

经过 30 多秒后，安装需要我们的确认。输入‘y’并按Enter键。这将下载、解压并安装上述提到的库。

一旦过程完成，你就有一个安装在 ITC Conda 环境中的 JupyterLab IDE。不要关闭控制台；从那里，我们将通过四个简单的步骤测试一切是否顺利。

1. 启动 JupyterLab：在相同的控制台中，你可以通过输入命令：“jupyter lab”来启动 JupyterLab。这将自动在默认网页浏览器中打开 JupyterLab 界面。

2. 创建一个新的笔记本：要创建一个新的笔记本（我们在其中编写代码），请点击 JupyterLab 界面左上角的“文件”菜单，并选择“新建笔记本”。这将会在新标签页中打开一个新的笔记本。

3. 编写代码：现在你可以在笔记本中开始编写代码。要创建一个新的代码单元格，请点击工具栏中的“+”按钮或按“Ctrl + Shift + Enter”。然后你可以在单元格中编写 Python 代码（例如：‘This is working’），并通过按“Shift + Enter”运行它。

4. 保存你的工作：记得定期保存你的工作，通过点击工具栏中的“保存”按钮或使用“Ctrl + S”键盘快捷键来保存。

这些是开始使用 JupyterLab 的基本步骤。随着你对界面的熟悉，你可以探索更多高级功能，如添加 Markdown 文本、导入和导出数据以及使用 JupyterLab 扩展。

🦚 注意：对于 ITC—特温特大学的学生，我们很幸运能使用 CRIB: A Geospatial Computing Platform 进行 Jupyter Lab。 如果你的计算机在跟随本课程时显示出某些处理限制，我强烈推荐使用这个云计算服务。

1.3. 3D Python 库

在本教程中，我将介绍五个对 3D 地理空间分析至关重要的库。这些库是NumPy、Pandas、Open3D、Matplotlib和Shapely。

🦚 注意：如果你想使用上述库，我们必须确保它们在你的环境中已安装并可用。因此，在相同的环境终端中，我们使用公式 “pip install package-name==version” (其中==version是可选的，用于固定某些版本)，如下面所示：

#We install the 5 libraries with the package manager 'pip', one line at a time
pip install numpy
pip install pandas
pip install open3d==0.16.0
pip install matplotlib
pip install shapely

NumPy：这个库用于处理数组和矩阵。它提供了对大型多维数组的快速高效操作，使其成为科学计算和数据分析的强大工具。一个使用NumPy的实际示例是创建一个点云，将其作为 3D 欧几里得空间中的数据点集合。为此，你可以创建一个具有三列的 NumPy 数组，每一行表示点云中的一个点：

import numpy as np
point_cloud = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

#To print the "point_cloud" variable as [Out] from the cell
print(point_cloud) 

在这个示例中，点云包含三个坐标为(1, 2, 3)、(4, 5, 6)和(7, 8, 9)的点。简单吧，你说？😁

Pandas：这个库更侧重于数据操作和分析。它提供了强大的数据结构和工具，用于处理结构化数据，例如 CSV 文件、电子表格和数据库。虽然它并不是专门为 3D 数据处理设计的，但仍可以用来以表格格式读写点云。为此，你可以创建一个 DataFrame 对象，其中列表示每个点的 X、Y 和 Z 坐标：

import pandas as pd

# create a DataFrame with X, Y, and Z columns
points_df = pd.DataFrame({ 'X': [1, 4, 7], 'Y': [2, 5, 8], 'Z': [3, 6, 9]})

# save the DataFrame as a CSV file at the same place as your script
points_df.to_csv('point_cloud.csv', index=False)

在这个示例中，点云包含三个坐标为(1, 2, 3)、(4, 5, 6)和(7, 8, 9)的点，如使用 NumPy 获得的。然后，DataFrame 使用 to_csv 函数保存到 CSV 文件中。

🦚 注意：Pandas 可以通过另一个 Python 模块进行扩展： Geopandas。这个库使得可以直接处理存储在例如 Shapefiles 或 PostGIS 数据库中的空间数据。这扩展了当前教程的范围，但了解这些内容是有益的，因为我们在其他情况下肯定会用到它。 😉

Open3D：这个库专注于 3D 数据处理和可视化。它提供了各种工具和函数，用于处理点云、网格和其他 3D 数据格式，如体素。

🦚 注意：快速安装该库的方法是运行 Anaconda 环境终端并输入以下命令： pip install open3d==0.16.0 ，这将安装版本 0.16.0 的 Open3D，如本教程中使用的版本。

现在，在 Jupyter Lab IDE 中，让我们使用内置的 Open3D 函数创建一个点云：

import open3d as o3d

# create an empty point cloud object
point_cloud = o3d.geometry.PointCloud()

# add points to the point cloud
points = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
point_cloud.points = o3d.utility.Vector3dVector(points)

在这个示例中，点云包含三个具有坐标(1, 2, 3)、(4, 5, 6)和(7, 8, 9)的点。

🦚 注意：在上面的代码块中，我们创建了一个 Open3D PointCloud 对象，然后通过 o3d.utility.Vector3dVector 函数将点列表传递给 points 属性，该函数将点列表转换为可以添加到点云对象中的格式。

Matplotlib：这个库用于数据可视化。它提供了许多工具，用于创建高质量的图表、图形和其他可视化内容。虽然它并不是专门为 3D 数据可视化设计的，但它仍然可以创建点云的 3D 散点图。为此，你可以使用 Axes3D 类创建 3D 图，并使用 scatter 函数绘制点：

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

# create a point cloud with 1000 random points
points = np.random.rand(1000, 3)

# create a 3D plot and add the points as a scatter plot
fig = plt.figure()
ax = fig.add_subplot(111, projection=’3d’)
ax.scatter(points[:,0], points[:,1], points[:,2])

# show the plot
plt.show()

在这个示例中，使用 NumPy 生成了一个包含 1000 个随机点的点云。然后，使用 scatter 函数在 3D 图中绘制这些点。最后，使用 show 函数显示图形。

Shapely：这个库主要用于 2D 几何操作，用于创建、操作和分析几何对象。它提供了广泛的工具来处理点、线、多边形和其他几何形状。Shapely的一个常见用例是创建一个多边形，并检查一个点是否在多边形内。为此，你可以使用定义多边形的坐标列表创建一个 Polygon 对象，然后使用contains函数检查点是否在多边形内：

from shapely.geometry import Point, Polygon

# create a polygon with four vertices
polygon = Polygon([(0, 0), (0, 1), (1, 1), (1, 0)])

# create a point
point = Point(0.5, 0.5)

# check if the point is within the polygon
if polygon.contains(point):
 print(‘The point is within the polygon.’)
else:
 print(‘The point is outside the polygon.’)

这只是一个简单的例子，但 Shapely 提供了许多其他用于处理几何对象的函数，这些函数对于更复杂的应用可能会很有用。

🦚 注意：在这个例子中，使用一个坐标列表创建了一个具有四个顶点的多边形。然后，使用 Point 函数创建一个点。最后，使用 contains 函数检查点是否在多边形内，并根据检查结果将消息打印到控制台。

在这个 3D 城市建模的工作流程中，我们的第一步是巧妙地结合简单而有效且经过验证的工具和库，下面做了视觉总结。

环境设置的总结。© 作者

现在我们有了明确的软件栈、功能齐全的 Python IDE 和对 Python 库的基本理解，我们可以开始准备数据。

第 2 步：3D 数据准备

我们现在进入第二步：3D 数据准备。

3D 数据准备。我们将下载 3D 数据集，然后开发一些简单的 3D 可视化，通过下采样和导出数据为适合 Python 使用的格式来进行处理。

我们的目标是收集和准备数据集，以便我们用 Python 进行分析。因此，这一步也充当了“3D 数据可视化”阶段，我们将定性评估我们所处理的数据。如果你准备好了，就开始吧。

2.1. 下载 3D 数据集

对于 3D 城市模型分析，我们从开放数据源中收集了一些优秀的数据集。我展示了荷兰一个特定的兴趣区域，但我鼓励你研究你自己的房子或任何对你感兴趣的地点（如果你住在荷兰的话，当然，😉）。

点云数据集。首先，我们使用geotiles.nl门户网站收集点云数据集，该网站提供了一些在 CC-BY 4.0 许可下的良好数据集。为此，你可以访问该网站，缩放到感兴趣的瓦片，并获取最新版本的 AHN LiDAR 点云（在我们的案例中是 AHN4），如下图所示。文件将被下载为.laz（LasZip）文件。

从 AHN-X 项目下载 3D LiDAR 数据集，通过门户网站 geotiles.nl。

🦚 注意：AHN（Actueel Hoogtebestand Nederland）源自航拍 LiDAR 覆盖，为荷兰所有地区提供数字高程地图。它包含详细且精确的高度数据，每平方米平均有八次测量。水务局、省份和公共工程部门等组织使用 AHN 进行水务和堤坝管理。根据地面高度和高程，决定水是否能从土地上顺利排出，沟渠中的水位可以多高，河流、洪水平原和沟渠中的水是否能得到充分排放，以及堤坝是否仍然足够高且强大。AHN 还用于许多其他类型的管理，如日常管理和堤坝维护、重大维护的规格准备、3D 绘图、许可和执法。市政府、企业和研究人员也使用详细的高程数据。

如果你喜欢图形，我为你准备了 AHN LiDAR 数据的快速概述：

• AHN1: 1997–2004, 1 pt/16 m2 至 1 pt/m2

• AHN2: 2007–2012, 8 pts/m2

• AHN3: 2014–2019, 8 pts/m2

• AHN4: 2020–2022, 10–15 pts/m2，

• AHN5: 2023–2025, 10–15 pts/m2

网格数据集。

对于这一部分，我们想在与我们下载的 LiDAR 区域相同的位置找到一个 3D 网格。因此，我们使用由 TUDelft 创建的平台 3DBAG，它允许我们从 CityGML Specification 中检索荷兰的 10M 建筑物，具有 LoD1.2、LoD1.3 和 LoD2.2。在这一步，我们主要关注几何形状。虽然我们知道语义和拓扑是 CityGML / City JSON 数据模型的重要方面，但将在后续文章中探讨。😉

从门户网站 3Dbag.nl 下载数据集。

🦚 注意：3DBAG 包含多个细节层级的 3D 模型，这些模型是通过结合两个可用数据集生成的：来自 BAG 的建筑数据和来自 AHN的高度数据。3D BAG 会定期更新，保持最新的开放建筑库存和高程信息。

下载数据集后，你应该会有一个.laz 文件格式的点云和一个或多个.obj 数据集（以及其附带的.mtl 文件），这些数据集大致描述了相同的范围，但具有不同的细节层级（在我们的例子中是 LoD 1.2、1.3 和 2.2），如下面所示。

🧙‍♂️ 向导： [可选] 如果你想深入了解 3D 数据文件格式，尤其是点云网格的部分，我建议你跟随下面的教程，详细了解如何将点云网格化及其结构。

## 5 步指南：使用 Python 从点云生成 3D 网格

教程：使用 python 自动从 3D 点云生成 3D 网格（.obj, .ply, .stl, .gltf）。(附赠)…

towardsdatascience.com

为了方便起见，你可以直接从这个 Drive 文件夹 下载选择项。

2.2. 3D 数据可视化

现在我们将进入CloudCompare以确保下载的数据分析没有特殊的惊喜。😁 启动CloudCompare后，你可以从本地文件夹加载.laz点云。当导入窗口显示时，确保仅导入此应用程序的“Classification”额外字段，并接受如下面所示的全局偏移。

在 CloudCompare 中进行 3D 数据可视化。

🦚 注意：全局偏移是一个临时的偏移，以允许*CloudCompare*处理地理参考数据，该数据超出了它可以处理的可视化范围。因此，这一步是透明的，并将在保存数据时应用*。*

现在我们可以继续导入网格数据。为此，我们执行相同的导入操作，并接受建议的平移偏移，同时确保它与应用于点云的偏移相同。

从下载的资产中加载 3D 网格。

现在你可以从数据库树中看到当前项目中导入的不同对象（如果找不到 DBTree，可以查看下图）。DBTree的功能有点类似于你的操作系统资源管理器，其中包含不同的点云或网格。每个对象（例如点云或 3D 网格）可以像图层一样可视化地激活 ☑️（或停用）并选择以识别对象属性。

CloudCompare 3D 数据处理和分析界面。© 作者

🦚 注意：CloudCompare 默认不保存。为了避免崩溃，如果你担心，可以将所有数据和分析放在 DBTree 中的一个文件夹中（右键点击 > 创建新的空组），并将此文件夹保存为.bin CloudCompare 项目。

2.3. 3D 数据子采样

现在是深入 3D 数据过滤的时候了。首先，我们从DBTree中选择点云，并应用空间子采样函数以每 50 厘米保留一个点，这样可以保留足够的信息，同时不会影响计算速度。为此，我们使用下面所示的subsample函数。

3D 点云子采样

🧙‍♂️ 向导: [可选] 在我们的案例中，我们使用了空间子采样函数，每隔 50 cm 平均保留一个点。如果你想探索和深入 3D 点云采样策略，我推荐深入阅读以下文章。

## 如何使用 Python 自动化 LiDAR 点云处理

关于点云子采样的终极指南，从零开始，用 Python 编写。涵盖 LiDAR I/O、3D 体素网格处理等……

towardsdatascience.com

一旦完成子采样步骤，我们将在 DBTree 中得到一个结果子采样点云，可以用于后续处理。

在 Mesh 方面，导入后，我们可以通过点击每个网格对象旁边的小箭头图标从 DBTree 中打开它，这会显示许多不同的子网格元素。这是因为在原始的 .obj 文件中，我们有一个“对象”细化，这允许我们保留一些与几何体相关的“语义”分解。每个子网格是一个构建实体的分解。

如果你想查看这些元素，可以打开网格，通过按住 Ctrl+Shift 并在第二次鼠标点击时选择所有子元素，然后 右键点击 > 切换 以显示它们。你还可以取消选中父网格元素的 可见 属性，以确保你看到的仅是子元素，如下所示。

CloudCompare 中的 3D 数据准备。

🦚 注意: 我们不能使用子网格选择，因为它们也包含网格父项的所有顶点，这意味着，如果我们使用它，我们将不得不将网格分割为仅选中的子网格的顶点。

很棒，干得好！从那里，你可以取消选择所有元素，只保留你想要考虑的网格（在我们的例子中是 LoD 1.2），并使用例如“横截面”功能，选择你感兴趣的房屋，如下所示。

CloudCompare 中的 3D 数据选择

一旦完成，我们就准备好导出我们选择的研究站点的两种 3D 模态。

2.4. 3D 数据导出

关于点云数据，我们可以以各种文件格式导出它。由于我们希望保留分类字段以进行一些 Python 分析（例如，了解一个点是否属于地面或建筑物），我们将 3D 点云导出为 ASCII 文件，以便在 Python 中轻松处理，如下所示。

CloudCompare 中的 3D 点云数据导出。

🦚 注意：在保存文件时，我们可以选择修改文件扩展名为 *.xyz* ，并在导出对话框打开时勾选“保留列名”选项。这将允许在文件中写入列名。然而，如果你用文本编辑器打开文件，确保删除两个不必要的反斜杠，如上所示。

我们的第一个 3D 数据集已经准备好，可以作为.xyz文件用于 Python。现在，让我们继续进行 3D 网格处理。

选择你感兴趣的房屋/建筑块后，我建议你获取相关的子网格元素名称，并将其用作导出文件的名称。关于导出对话框的选项，你可以选择.obj文件扩展名，这是一个安全的选择😉。

CloudCompare 中的 3D 网格导出。

3D 网格现在已经准备好，并附带一个材料.mtl文件（在我们的案例中不太有用）。

步骤 2 已完成！干得好。我们首先收集了荷兰一个区域的 3D 点云和 3D 网格。然后我们进行了可视化，以检查它们是否符合我们的意图。接着我们过滤了点云，保留每 50 厘米左右一个点，3D 网格只保留一个代表建筑房屋的对象，并将两种数据分别导出为.xyz和.obj + .mtl文件。¹

步骤 2 的可视化总结：3D 数据准备。© F. Poux

现在让我们将这些数据集放入一个出色的 Python 设置中，以最大化自动化和 3D 分析！

¹为了方便，你可以在这个Drive 文件夹中找到这些数据集。

第 3 步. Python 自动化

现在真正有趣的部分开始了，是时候用 Python 编程了！🤓

3D Python 自动化。

如上所示，我们将遵循五个阶段的方法，包括导入库、加载数据集、设置我们的 3D Python 可视化工具、定义 3D 挑战的解决方案，然后将结果导出以供 Python 之外使用。让我们先在 Python 中建立我们自动化管道的主体，然后再处理不同的挑战。

3.1. 导入库

如步骤 1 中定义的，我们将坚持使用最少的库，以加深我们对其使用的了解。这些库包括NumPy、Pandas、Open3D、Matplotlib和Shapely。

我们将从我们的 IDE 中在 Python 笔记本（.ipynb）中编写以下代码行。

用于编写我们脚本的 IDE 视图。

我们用以下代码块导入上述提到的库：

import numpy as np
import pandas as pd
import open3d as o3d
import matplotlib.pyplot as plt
from shapely.geometry import Polygon

print(f"Open 3D Version: {o3d.__version__}")

这将返回当前的 Open 3D 版本作为字符串：Open 3D Version: 0.16.0。我们已经设置好了，可以继续加载数据集。

3.2. 加载数据集

现在我们将定义存储数据集的具体路径。我喜欢将其明确且相对于我的代码文件。这样，一切都相对表达（../ 表示转到父文件夹），使得在不同机器上编码时，浏览文件夹变得容易。

data_folder="../DATA/"
pc_dataset="30HZ1_18_sampled.xyz"
mesh_dataset="NL.IMBAG.Pand.0637100000139735.obj"
result_folder="../DATA/RESULTS/"

点云数据集。

我们可以通过首先创建一个名为 pcd_df 的 Pandas DataFrame 对象来准备点云，该对象将包含点云数据：

pcd_df= pd.read_csv(data_folder+pc_dataset, delimiter=";")
print(pcd_df.columns)

这将返回 pcd_df 数据框中的列名，它们是：[‘X’, ‘Y’, ‘Z’, ‘R’, ‘G’, ‘B’, ‘Classification’]。这对于仅选择明确的“列”而不是模糊的索引非常有用 😁。这正是我们要做的：仅选择 [‘X’, ‘Y’, ‘Z’] 坐标以创建 Open3D PointCloud 对象。这是一个很好的方法来理解当我们使用不同的库时，我们必须适应不同的机制以将数据集转换为不同的 Python 对象。在这里，我们从 Pandas DataFrame 转到 Open3D PointCloud，这是使用 Open3D 库中实现的 Open3D 函数所必需的。让我们按以下四个步骤进行：

Open3D 点云创建的概念工作流程。© 作者。

这个分解机制转化为以下内容：

创建 Open3D 点云的代码工作流程。© 作者。

不过，如果我们想要更简洁一点，这个转换可以通过一行代码完成，即：

pcd_o3d=o3d.geometry.PointCloud(o3d.utility.Vector3dVector(np.array(pcd_df[['X','Y','Z']])))

我们现在有两个重要的变量：pcd_o3d 和 pcd_df。我们还可以为 Open3D PointCloud 对象添加一些颜色：

pcd_o3d.colors = o3d.utility.Vector3dVector( np.array( pcd_df[[‘R’,’G’,’B’]]) / 255 )

我们必须注意将 R,G,B 值转换为 [0,1] 范围内的浮点值，并确保将 Vector3dVector 对象传递给 colors 属性。

网格数据集。

现在，我们可以使用 open3d 的 read_triangle_mesh() 方法将 3D 网格加载到 mesh 变量中。我们还可以使用 paint_uniform_color() 方法为网格上色：

mesh=o3d.io.read_triangle_mesh(data_folder+mesh_dataset)
mesh.paint_uniform_color([0.9,0.9,0.9])

我们现在有两个 Open3D 对象：一个包含 7 103 848 个点的 APointCloud 和一个包含 674 个点和 488 个三角形的三角网格。让我们看看这意味着什么，好吗？

3.3. Python 3D 可视化

要在 Open3D 中可视化不同的 3D 对象，我们必须传递一个包含这些 Open3D 对象的 python 列表。因此，我们的列表由一个 Open3D PointCloud 和一个 Open3D TriangleMesh 组成，即 [pcd_o3d, mesh]。让我们使用以下代码在独立窗口中可视化这个组合：

o3d.visualization.draw_geometries([pcd_o3d,mesh])

🦚 注意：上面的代码行将创建一个交互式的 Open3D 窗口，结合了 3D 点云和 3D 网格。

要调整显示的颜色，有一个实用的技巧：使用一个颜色变量，该变量将传递给 PointCloud Open3D 对象的颜色属性。此变量应包含从 0 到 1 的 R、G、B 浮点值。

🦚 注意：你是否注意到两个截图之间的细微差别？右侧我们可以更好地勾勒出边界，并且有更强的深度感（这在交互式操作中更为明显）。这是因为在第二种情况下，我们还使用了法线进行可视化。要做到这一点，你可以在运行可视化部分之前运行 pcd_o3d.estimate_normals() 和 mesh.compute_vertex_normals() 。享受吧 😉

假设我们想要根据分类属性可视化点云。我们需要遵循以下四阶段过程：

使用分类作为 Open3D 颜色的概念工作流程。

如果你仔细观察，这将转换为如下代码：

使用分类作为 Open3D 颜色的代码工作流程。

# 1=unclassified, 2=Ground, 6=building, 9=water, 26=rest
pcd_df['Classification'].unique()
colors=np.zeros((len(pcd_df), 3))
colors[pcd_df['Classification'] == 2] = [0,0,0]
colors[pcd_df['Classification'] == 6] = [1,0,0]
colors[pcd_df['Classification'] == 9] = [0,0,0]
colors[pcd_df['Classification'] == 26] = [0,0,0]
pcd_o3d.colors = o3d.utility.Vector3dVector(colors)

🦚 注意：因为我们希望对每个类别的颜色进行控制，我们可以根据唯一类别的数量调整上述代码。然后，映射是基于 LAS (1.4) 规范进行的。

LAS LiDAR 数据规格文件格式的 ASPRS 标准点类别。

我们可以可视化结果，类似于这样（根据输入的R,G,B值显示不同的颜色）：

Open3D 交互式多模式数据可视化。© 作者。

我们已经加载了变量，可以看到点云和 3D 网格，一切看起来运行顺畅！现在是定义我们想要进行的各种 3D 城市分析任务的时刻！

3.4. 3D 城市分析

3D 城市分析指的是使用城市环境的三维（3D）模型来分析和理解建筑环境的各种方面，如建筑能源性能、城市形态和行人流动。这种分析通常涉及使用专业范式和进行分析，以获得有关城市规划和设计的见解。3D 城市分析可供城市规划师、建筑师和工程师使用，以提供有关基础设施布置、公共空间设计和各种环境影响缓解的决策信息。在本教程中，我们将涉及 3D 城市分析的三个主要方面：

1. 城市形态分析：我们可以使用 Python 分析 3D 城市数据集中建筑物的形状和形式，这可以帮助提供有关城市设计和规划的决策。

2. 地理空间分析：我们可以对城市模型进行地理空间分析，例如根据可达性和环境影响等因素确定新基础设施项目的最佳位置。

3. 3D 可视化：我们可以创建交互式 3D 可视化，这可以帮助利益相关者更好地理解和参与城市规划和设计项目。

因此，这一步专门用于一个应用程序，我们在这里进行大部分分析。我们在 3D 城市分析中进行演示，并在第四部分通过各种 Python 挑战进行扩展，如下所示。

城市背景下的 3D 分析及其在第 4 步中的分解：3D Python 挑战（第四部分）。

不过，我们可以遵循一个通用的工作流程，适应于每个应用程序，通常由输入、处理管道和每个特定分析部分的输出组成。

从输入到输出的概念性工作流程，适用于点云数据。

输入可以有所不同，但在大多数情况下，它是一个包含空间信息的 NumPy 数组。输出可以是一个整数、一个列表、另一个 NumPy 数组……你需要的任何东西。😁

🦚 注意：这个特定阶段比较复杂，我们将在笔记本中留出一些空间，用于记录与第 4 步挑战相关的不同分析。现在，让我们定义一些方法来导出我们的数据。

3.5. 3D 多模态导出

在我们的 3D 管道完全功能后，我们可以将结果保存到一个或多个文件中，以便在 Python 之外处理。为此，我们将使用两种有价值的可能性。

1. Numpy（推荐用于复杂输出）： 我们可以使用以下代码行通过 Numpy 导出：np.savetxt(result_folder+pc_dataset.split(“.”)[0]+”_selection.xyz”, np.asarray(o3d_parcel_corners),delimiter=’;’, fmt=’%1.9f’)

2. Open3D（推荐仅需空间属性时使用）： 我们可以使用以下代码行通过 Open3D 导出：o3d.io.write_point_cloud(result_folder+pc_dataset.split(“.”)[0]+”_result_filtered_o3d.ply”, pcd_selection, write_ascii=False, compressed=False, print_progress=False)

🦚 注意：.split(.)允许将pc_dataset字符串对象拆分为两个字符串的列表，分别在.之前和之后，然后我们仅保留第一个元素，使用[0]。NumPy 将一个变量*o3d_parcel_corners的阶段以分隔符 ; 导出到一个.xyz* ASCII 文件中。Open3D 将从 open3d 对象*pcd_selection中写入一个.ply* 文件。*

哇，干得好！Python 自动化的批量结构已经运行起来了！恭喜！我们已经导入了库，数据集被存储在不同的显式变量中，我们确保可以处理各种可视化而无需离开 Python 的舒适环境。导出步骤已设置完成；剩下的就是解决我们在第 4 步中提出的初始问题：

我们房子周围的建筑区域有多密集？房子是否可能会受到洪水影响？我是否遵守了我拥有地块的建筑比例？

步骤 4. 3D Python 挑战

为了回答上述问题，我们将找到解决四个挑战和体素化步骤的方案，如下所示。

步骤 4：3D Python 挑战。我们研究兴趣点查询、手动边界选择、高点提取、体素化和建筑覆盖提取。

挑战 1 将允许裁剪出所需的研究区域。挑战 2 将允许提取所拥有地块的建筑比例。挑战 3 指导洪水分析。而挑战 4，经过体素化处理后，将允许获取感兴趣区域的建筑覆盖情况。让我们开始吧。

4.1. 兴趣点查询

对于这个挑战，我们从一个包含 3D PointCloud Open3D 对象和 TriangleMesh Open3D 对象的输入开始，如下所示：

一个带有 3D 网格对象的 3D 点云。© F. Poux

这个挑战的目标是只保留落在兴趣点（建筑房屋）一定距离内的数据点。我们的输入是点云和网格，输出是经过过滤的点云，符合到 POI 距离的标准，如下所示：

为此，我设置了一个六阶段的过程，如下所示：

3D 点云半径搜索的概念工作流。© 作者。

🎓 学习笔记：目标是尽量自己解决问题，如果遇到困难可以查看下面的解决方案。准备好后，可以阅读下面的解决方案。 👇

(1) 设置距离阈值时，我们将所需的半径值（例如 50）传递给新变量 dist_POI。 (2) 然后，我们使用 Mesh 对象的 get_center() Open3D 方法从 BAG 数据集中获取 POI。这允许获取网格的中心作为 POI：POI=mesh.get_center()。之后，我们可以创建 KD 树 (3)，一种用于组织空间中点的数据结构。KD 树对点云处理非常有用，因为它们允许快速进行最近邻搜索和范围查询。这是好消息，因为这正是我们要做的 😁。实际上，使用 KD 树可以快速找到离给定点最近的点，或者所有在一定距离内的点（我们的 POI），而无需遍历点云中的所有点。

# Creating a KD-Tree Data Structure
pcd_tree = o3d.geometry.KDTreeFlann(pcd_o3d)

🦚 注意：KD-Tree 通过递归地将空间划分为较小的区域，每个级别沿一个维度的中位数进行划分（X 是一个维度， Y 是另一个， Z 是另一个 😉*）。这会生成一个“树”结构，其中每个节点表示一个空间分区，每个叶子节点表示一个单独的点。*

从那里，我们可以使用 Open3D 的 search_radius_vector_3d() 方法，并从输出索引中选择点：最后可视化我们的结果（4+5+6）：

# Selecting points by distance from POI (your house) using the KD-Tree 
[k, idx, _] = pcd_tree.search_radius_vector_3d(POI, dist_POI)
pcd_selection=pcd_o3d.select_by_index(idx)

我们最终可以可视化我们的结果（6）：

o3d.visualization.draw_geometries([pcd_selection,mesh])

半径搜索下的 3D 点云和网格叠加。

总共涉及以下代码管道：

半径搜索的代码工作流。

非常好！现在你有了一个有效的解决方案，让我们从那个选择中提取地块面积。

4.2. 手动边界选择

要提取非官方边界，我们需要采用一些半自动和互动的方法。好消息是，我们可以直接在 Python 中使用 Open3D 来完成这项工作。

首先需要创建一个带有以下 draw_geometries_with_vertex_selection() 方法的互动 Open3D 窗口：

o3d.visualization.draw_geometries_with_vertex_selection([pcd_selection])

然后可以跟随下面的动画部分，这允许选择定义地块角落的点。

手动边界选择的互动过程，使用 Open3D GUI 并按住 MAJ+鼠标选择感兴趣的点。

结果将出现在你的笔记本（或 REPL）中的单元格下方，关闭窗口后即可查看。

🦚 注意：在此步骤中，使用 R、 G、 B 着色可能更容易。因此，如果您想使用正确的索引，需要在选择前返回进行更改。 😉 如果您想提取地籍边界，可以通过导入官方的 2D 矢量图形并根据此数据约束进行裁剪来实现。

从你的 REPL 中，你可以将选定点的不同索引（例如 34335、979、21544、19666、5924、21816、38008）复制并粘贴到 select_by_index() 选择方法中，以定义 o3d_parcel_corners 变量：

o3d_parcel_corners=pcd_selection.select_by_index([34335 ,979 ,21544 ,19666 ,5924 ,21816 ,38008 ])

我们仍然需要进一步准备角落，因为我们想避免考虑 Z 值。因此，我们将过滤坐标以去掉 Z 值，但要注意：这样做意味着我们认为我们在一个平坦区域（这在荷兰是事实 😉）。

o3d_parcel_corners=np.array(o3d_parcel_corners.points)[:,:2]

从那里，接下来使用 Shapely 库计算地块的面积！为此，我们可以直接使用 Polygon 函数首先从提供的角落集合中创建一个多边形：

Polygon(o3d_parcel_corners)

输出如下：

错误的几何形状

这看起来有点不对，是吧？如果你不相信我，我们可以计算面积来检查：

pgon = Polygon(o3d_parcel_corners)
print(f"This is the obtained parcel area: {pgon.area} m²")

这会输出43.583 m²。对于一个大建筑来说，这听起来很奇怪，不是吗？哈哈，一个简单的问题变得有点复杂了！实际上，这里的问题在于计算面积需要一个由多边形组成的闭合形状。这并不总是能实现，因为添加角点的顺序可能会影响数据结构。因此，问题变成了排序坐标以避免边缘交叉。

处理这个问题的一种方法是从中心点的角度来看所有坐标。然后，我们可以计算列表中每个角点与中心点之间的角度。我们这样做是为了了解每个角度的宽度，从而提供排序坐标的手段。将其翻译成代码允许定义一个排序函数，如下所示：

def sort_coordinates(XY):
    cx, cy = XY.mean(0)
    x, y = XY.T
    angles = np.arctan2(x-cx, y-cy)
    indices = np.argsort(-angles)
    return XY[indices]

🦚 注意：我使用 NumPy arctan2() 方法对每个坐标进行角度估算。这将返回以弧度表示的角度数组。剩下的工作是将角度按升序排序，以获取正确顺序的索引列表。然后，列表可以使用 argsort() 方法修正原始列表的索引。在此步骤中，它不适用于 U 形建筑。

从那里，我们可以将排序函数应用于角点，创建一个新的排序变量，计算多边形及其相关面积：

np_sorted_2D_corners=sort_coordinates(o3d_parcel_corners)
pgon = Polygon(np_sorted_2D_corners)
Polygon(np_sorted_2D_corners)
print(f"This is the parcel area: {pgon.area} m²")

这返回了一个面积为2 247.14 m²的多边形，这更为可信。😉

正确的几何形状

再次做得很好。我们现在对我们地块的面积有了一个不错的了解。现在，让我们在区域内找到高点和低点。

4.3. 查找高点和低点

要获取区域的低点和高点，一种特定的方法是使用 Open3D 的get_max_bound()和get_min_bound()方法：

pcd_selection.get_max_bound()
pcd_selection.get_min_bound()

但是，这样做不会提示哪个点是最高的，哪个点是最低的。我们需要的是点的索引，以便在之后检索坐标。为此，我建议我们这样编码：

1. 我们从 Open3D PointCloud 对象中创建一个 NumPy 数组对象np_pcd_selection，该对象仅包含X、Y和Z坐标。

2. 我们使用argmax()方法收集Z维度上的最小值和最大值的索引，并将结果存储在变量lowest_point_index和highest_point_index中。

np_pcd_selection=np.array(pcd_selection.points)
lowest_point_index=np.argmin(np_pcd_selection[:,2])
highest_point_index=np.argmax(np_pcd_selection[:,2])

让我们通过选择索引，创建TriangleMesh球体，将它们转换到点的位置，然后使用 Open3D 进行可视化来检查结果：

# We create 3D Spheres to add them to our visual scene
lp=o3d.geometry.TriangleMesh.create_sphere()
hp=o3d.geometry.TriangleMesh.create_sphere()

# We translate the 3D Spheres to the correct position
lp.translate(np_pcd_selection[lowest_point_index])
hp.translate(np_pcd_selection[highest_point_index]))

# We compute some normals and give color to each 3D Sphere
lp.compute_vertex_normals()
lp.paint_uniform_color([0.8,0.1,0.1])
hp.compute_vertex_normals()
hp.paint_uniform_color([0.1,0.1,0.8])

# We generate the scene
o3d.visualization.draw_geometries([pcd_selection,lp,hp])

如下所示，我们现在有了包含高点和低点的 3D 场景。

让我们在 350 米的扩展区域内研究建筑覆盖情况。为此，我们将重新执行代码中的某些部分，如下所述。

4.4. 点云体素化

我们想要提取已建立的覆盖范围。为此，我们将采取一种直观的方法，首先将点云的模式转换为 2D 像素的填充模拟：3D 体素。这将允许我们使用以下方法：

通过构建 3D 体素数据结构来提取高点。我们首先对点云进行体素化，然后从二进制角度对每个体素上色，最后按自上而下的索引进行筛选。© F. Poux

基本上，我们将（1）生成存在点的体素，然后（2）根据体素所持有的点的分类对体素上色，最后，我们将筛选体素，只保留和计数每个 X,Y 坐标的最高体素。这样，我们可以避免计数那些会偏向已建立覆盖范围的元素。

第一步是创建一个体素网格。这是 Open3D 的强项：通过这些简单的代码行，我们可以将一个体素网格适配到点云中，每个体素是一个 20 cm 的立方体，然后可视化结果：

voxel_grid = o3d.geometry.VoxelGrid.create_from_point_cloud(pcd_selection, voxel_size=2)
o3d.visualization.draw_geometries([voxel_grid])

3D 体素数据结构的视图。© F. Poux

🦚 注意：您看到的结果已经经过体素上色的变化。这个变化意味着，使用 Open3D 时，我们必须更改体素点的颜色。然后，VoxelGrid 方法将平均这些颜色，并将结果保留为体素的颜色。

现在我们知道如何从点云生成 3D 体素，让我们通过改变它们的颜色方案来绕过颜色平均的问题。为此，一个简单的方法是以二进制方式处理颜色。即黑色（[0,0,0]）或其他颜色（所有其他颜色）。这意味着我们可以将颜色变量初始化为 0，然后选择所有被分类为建筑的点，并赋予它们另一种颜色，例如红色（[1,0,0]）：

colors=np.zeros((len(pcd_df), 3))
colors[pcd_df['Classification'] == 6] = [1, 0, 0]
pcd_o3d.colors = o3d.utility.Vector3dVector(colors)

非常好！现在，由于我们对原始点云进行了操作，我们需要根据自己的选择重新定义 POI 和选择。为了提高效率，您会找到可以运行的代码块，以更新体素渲染到新的颜色方案：

# Defining the POI and the center of study
dist_POI=50
POI=mesh.get_center()

# Querying all the points that fall within based on a KD-Tree
pcd_tree = o3d.geometry.KDTreeFlann(pcd_o3d)
[k, idx, _] = pcd_tree.search_radius_vector_3d(POI, dist_POI)
pcd_selection=pcd_o3d.select_by_index(idx)

#Computing the voxel grid and visualizing the results
voxel_grid = o3d.geometry.VoxelGrid.create_from_point_cloud(pcd_selection, voxel_size=2)
o3d.visualization.draw_geometries([voxel_grid])

二进制着色的点云。© F. Poux

太棒了！现在，我们需要实际获取每个体素的离散整数索引，并将它们列出，还有体素的颜色和边界。这将允许我们之后遍历每个体素及其颜色，检查一个体素是否在“体素列”中是最高的。为此，一个有效的方法是使用列表推导式来向量化我们的计算，避免不必要的循环：

idx_voxels=[v.grid_index for v in voxel_grid.get_voxels()]
color_voxels=[v.color for v in voxel_grid.get_voxels()]
bounds_voxels=[np.min(idx_voxels, axis=0),np.max(idx_voxels, axis=0)]

太棒了！在 50 米的选择区域中，我们有一个 49 x 49 x 16 的体素网格，总共有 38,416 个填充的体素。现在是提取已建立的覆盖范围的时候了！

4.5. 提取已建立的覆盖范围

现在我们有了一个具有二进制颜色的体素化点云，我们专注于下图的第三阶段：

正如我们所见，选择仅顶层体素比看起来复杂一些！😁 但幸运的是，我设计了一个简单而强大的管道，可以做到这一点，见下文。

体素选择工作流以七个子步骤提取已建成的覆盖范围。© F. Poux。

让我们一步一步来。

(1) 首先，我们初始化两个字典来保存最大索引：

max_voxel={}
max_color={}

(2 到 5) 我们遍历所有填充的体素，以检查它们是否在体素列中最高或应被丢弃，这给出了以下 for 循环：

for idx, v in enumerate(idx_voxels):
    if (v[0],v[1]) in max_voxel.keys():
        if v[2]>max_voxel[(v[0],v[1])]:
            max_voxel[(v[0],v[1])]=v[2]
            max_color[(v[0],v[1])]=color_voxels[idx]
    else:
        max_voxel[(v[0],v[1])]=v[2]
        max_color[(v[0],v[1])]=color_voxels[idx] 

🦚 注意: 首先要注意的是循环定义中的 enumerate() 方法。这允许在循环 idx_voxels 变量的每个值时跟踪列表的索引。很方便！第二点是我们使用“元组”数据类型，如 (X, Y)，这给出了体素的整数位置。这使我们能够确保我们持续检查相同的 X、 Y 网格位置。最后，“if”语句允许测试表达的条件，并在条件返回 True时执行。如果条件不为 True，则执行 else 语句（如果存在）或跳过并退出条件检查。

(6) 我们初始化标记为已建成或未建成的体素的计数，并检查保留的顶部体素的颜色是否为黑色，如果是，我们更新相应类别的体素计数：

count_building_coverage,count_non_building=0,0
for col in list(max_color.values()):
    if np.all(col==0):
        count_non_building+=1
    else:
        count_building_coverage+=1

🦚 注意: np.all() 也是一个布尔检查，只会在所有值都是 True 时返回 True。在我们的例子中，黑色是* [0,0,0]，因此返回 True ，因为所有 R、 G* 和* B 都设置为 0。如果其中一个不是零，则意味着并非所有值都是* 0，np.all() 将返回 False，这会触发* else 语句。简单吧？ 😉

(7) 我们可以提取每种类型（已建成和未建成）的覆盖区域，记得将计数变量乘以实际的 2D 体素面积（在我们的例子中为 4 平方米），并计算比率：

print(f"Coverage of Buildings: {count_building_coverage*4} m²")   
print(f"Coverage of the Rest: {count_non_building*4} m²")
print(f"Built Ratio: {(count_building_coverage*4)/(count_building_coverage*2+count_non_building*2)} m²") 

这样，对于 50 米的范围，我们的覆盖率为 17.3%，相当于 1352 平方米的已建成区域和 6456 平方米的其余区域。（350 米查询的覆盖率为 19.2%，73416 平方米和 308280 平方米）。

因此，与其余部分相比，我们在选定的兴趣点上有一个很好的建筑占用平衡！

希望测试我们的 3D Python 工作流并没有让你感到太紧张！随时返回查看代码和工作流片段，掌握隐藏的代码技巧，以便（1）出色地应对挑战和（2）优化我们实现的效率。

然而，仍然有一些调整空间，例如使用分类信息结合高低 POI 来隔离水流的地面点，或使用地块表面选择作为过滤技术来提取已建成的覆盖范围。

• 💻 直接访问实操代码：Google Colab

• 🪄 直接访问点云处理：Github

• 🧙‍♂️ 直接访问点云课程：3D Academy

🔮 结论

恭喜！这在使用 LiDAR 进行城市建模的 3D Python 工作流程领域是一个真正的巨大第一步！我们下面构建的流程非常通用，可以作为你未来分析工作的参考点，这些工作可能会遵循类似的模式。

在 LiDAR 城市模型背景下的 3D Python 工作流程。我们从环境设置（步骤 1）和 3D 数据准备（步骤 2）开始。一旦完成这些步骤，我们将进入 Python 自动化（步骤 3），其中一个特定部分涉及 3D Python 挑战（步骤 4），如地块表面或兴趣点查询。© 作者

总结一下你解锁的新技能，你现在可以高效地应对以下挑战：

1. 将开源软件与 Python 结合成一个连贯的工作流程；

2. 收集开放数据集并在处理之前准备好它们；

3. 使用 3D Python 处理 3D 数据模式，特别是 3D 点云、3D 网格和 3D 体素；

4. 使用每种 3D 模式的关键方面进行高级地理空间分析，同时优化代码；

5. 作为城市规划师从分类点云中提取关键洞察，以更好地理解本地区域。

如果你感到自己能够掌握并处理这些不同的方面，那么你正走在成为伟大的 3D 地理空间专业人士的道路上！唯一需要做的就是保持努力，推动现有的边界。

🤿 更进一步

但学习旅程并未就此结束。我们的终身探索才刚刚开始，未来的步骤将深入研究 3D 体素工作，探索语义、CityGML、CityJSON，特别是如何从聪明的城市走向智慧城市。此外，我们还将利用深度学习技术分析点云，并解锁高级 3D LiDAR 分析工作流程。有很多令人兴奋的内容！

使用 Python 进行 3D 地理空间数据集成：终极指南

原文：towardsdatascience.com/3d-spatial-data-integration-with-python-7ef8ef14589a

3D Python

用多模式 Python 工作流整合地理空间数据的教程：结合 3D 点云、CityGML、体素、矢量 + 栅格数据

Florent Poux, Ph.D.

·发表于 Towards Data Science ·阅读时间 39 分钟·2023 年 11 月 7 日

–

现在的科技进步速度简直疯狂。尤其是当我们看到 3D 数据对地理空间分析和数字双胞胎的重要性时更是如此。能够以三维捕捉和分析数据意味着我们可以创建对现实世界对象和环境的精确表示。

3D 空间数据集成通过理解 3D 数据捕捉的范围来实现。© F. Poux

🦄米拉：一图胜千言。那么数字双胞胎呢？

这对城市规划、基础设施管理和灾难响应等领域尤为重要。

通过整合 3D 数据，我们可以通过依赖于精确和可靠的数据表示来提高做出明智决策的能力。此外，将这些数据整合到数字双胞胎中可以生成非常逼真的现实资产和系统的复制品，从而提高模拟和分析的效率。

但是（总是有“但是”），有效的地理空间分析和数字双胞胎创建依赖于高效整合和可视化不同的数据格式。为实现这一点，必须全面了解各种数据模式及其如何无缝集成和可视化。在数据方面，我们希望创建一个统一且全面的区域表示，以便数据能够重叠。我们真是太幸运了，因为这正是我们今天要解锁的内容！

要构建一个空间数字世界，我们必须研究 3D 数据集成。许多信息来源，如矢量数据、栅格数据、3D 点云或 3D 城市模型，可以组合形成我们星球上发生事件的统一视图。© F. Poux

在这本动手指南中，我提供了一个面向系统的 3D 数据集成工作流参考，使用 Python。因此，你不需要昂贵的软件或大规模的砖块式流水线！只需我们的 Python 朋友和精心挑选的一小组强大模块与函数。

这个倡议的终极目标是，你将拥有一个全面的指南和伴侣，陪伴你完成 3D 数据之旅！工作流被结构化为七个不同的阶段，如下所示。

使用 Python 的 3D 数据集成工作流。这是一个七步过程，用于生成统一的数据中心视图。© F. Poux

每个阶段逐步构建，以确保你可以从头开始或模块化地插入到现有系统中。由于构建得非常全面，目录将使你更容易浏览！

Chapter 1\. 3D Python Setup
1.1\. Environment Setup
1.2\. Base Libraries
1.3\. 3D Data Libraries
1.4\. Geospatial Libraries
1.5\. IDE Setup

Chapter 2\. Multi-Modal Data Curation
2.1\. 3D Data Sourcing
2.2\. Spatial Raster (GIS)
2.3\. Vector Data (GIS)
2.4\. Other Sources

Chapter 3\. Data Analysis and Profiling
3.1\. 3D Point Clouds and voxels
3.2\. 3D mesh and city models
3.3\. Spatial / Raster Imagery
3.4\. DSM, DTM, CHM

Chapter 4\. Registration / Reprojection
4.1\. Selecting a Reference System
4.2\. Data Georeferencing
4.3\. Data Reprojection
4.4\. Rigid Registration

Chapter 5\. Data Pre-Processing
5.1\. Data Cleaning
5.2\. Data Transformation
5.3\. Data Reduction
5.4\. Data Enrichment (Fusion)

Chapter 6\. Data Visualization and Validation
6.1\. 3D Data Inspection
6.2\. Point Cloud Canonical Link
6.3\. Hybrid Multi-Modal Visualization
6.4\. Projection-based Inspection

Chapter 7\. Data Sharing
7.1\. Selection Method Definition
7.2\. Data Organization
7.3\. File Format Definition
7.4 Export and External Use 

只要你准备好了，就让我们一起跳入这次精彩的 3D 数据集成探险吧，身边有杯咖啡，但不要离电脑太近 ☕（这是我亲身经历的惨痛教训）

🎵读者注意：这本动手指南是* UTWENTE 与合著者 🦊 F. Poux, 🦄* M. Koeva,* 和* 🦝 P. Nourian 的联合工作的一部分。我们感谢来自数字双胞胎 @ITC 项目的资助，由特温特大学 ITC 学院提供。所有图片均 © F. Poux

第一步：3D 数据的实施设置

第一步：3D 数据的实施设置。© F. Poux

第一个任务是快速建立一个轻量级环境，以开发我们的 3D 数据集成工作流。这是一个简单的阶段，但确保设置正确是可扩展性和复制的关键。让我们把事情做好吧！

环境设置包括基础库、3D 数据库、地理空间库和 IDE。所有这些都建立在虚拟环境管理和 Python 之上。© F. Poux

轻量级环境设置

使用 Anaconda 进行 Python 环境设置，强大的库和集成开发环境（IDE）不必痛苦。Anaconda 提供了管理 Python 包和环境的便捷方法，然后你可以使用强大的 IDE，如 Jupyter Lab 或 Spyder，让编程变得轻松。关于设置 3D Python 开发环境的详细信息，我建议你查看以下文章：

[## LiDAR 城市模型的 3D Python 工作流：逐步指南]

解锁 3D 城市建模应用程序的流畅工作流的终极指南。教程涵盖了 Python……

towardsdatascience.com

🦊 Florent：如果你不想再参加另一个会话，不用担心，我不会抛下你！作为这个终极指南的一部分，这里有一个绝佳的轻量级设置，启动时间少于 5 分钟，计时 ⌚。

首先，你可以访问 Anaconda 网站 并下载适合你操作系统（Windows、macOS 或 Linux）的 Miniconda 安装程序（一个免费的 conda 最小安装程序），并选择 Python 10 版本。然后，你可以按照 Anaconda 网站上的安装说明在你的计算机上安装 Miniconda。

就这样！你现在已经用轻量级的 miniconda 安装了最简单的 Python，这将使你非常容易隔离一个受控的虚拟环境。在继续下一步之前，我们启动 miniconda 并访问其命令行：

在 Windows 中，搜索“miniconda”应该会得到这个结果。

一旦进入 Anaconda Prompt，我们按照下面所示的简单四步骤过程进行操作。

环境创建工作流程。© F. Poux

1. 要创建一个新环境，我们写：conda create -n GEOTUTO python=3.10

2. 要切换到新创建的环境，我们写：conda activate GEOTUTO

3. 要检查 Python 版本，使用 python --version，以及已安装的包：conda list。这应该分别显示 Python 3.10 和基本库的列表。

4. 要在新的环境中安装 pip，我们写：conda install pip

就这样！我们现在准备好使用 pip 管理器安装 3D 数据集成所需的库：pip install package-name，其中你需要逐一更换包名（例如：numpy、matplotlib、laspy[lazrs,laszip]、open3d、rasterio、geopandas）。

Python 基础库

Python 基础库：Numpy 和 Matplotlib。

第一个包的安装通过提示完成：pip install numpy。不需要介绍NumPy，这是 Python 的基础数值和科学计算库。它支持大型多维矩阵，并提供了一系列数学函数，方便使用。NumPy 是许多其他科学库的基础，在数据分析、机器学习和科学研究中被广泛使用。

🦝 Nourian：NumPy 全关于线性代数。如果你需要重新学习线性代数并熟悉它，你可以从这里开始： 计算机图形学的线性代数基础 (ResearchGate)

Matplotlib 是一个流行的 Python 绘图库，支持 2D 绘图和基本的 3D 绘图功能。安装方法为：pip install matplotlib。它提供了广泛的可定制可视化选项，允许用户创建各种类型的图表，如折线图、散点图、条形图、直方图、3D 图表等。Matplotlib 在科学研究和数据科学工作流中被广泛使用。

3D Python 库

Python 3D 库：Open3D 和 Laspy。

Open3D 是一个现代的 3D 数据处理和可视化库，主要关注 3D 点云和网格。它提供了处理 3D 数据的功能，如点云配准、几何处理、网格创建和可视化。Open3D 特别适用于与 3D 计算机视觉、机器人技术和增强现实应用相关的任务。安装 Open3D 的方法是：pip install open3d。

然后，我们需要安装 Laspy：pip install laspy[lazrs,laszip]。这个 Python 库用于读取、写入和修改存储在 LAS（激光雷达数据交换格式）和 LAZ（压缩 LAS）文件格式中的 LiDAR 数据。它提供了处理来自 LiDAR 扫描仪的点云数据的工具，并在地理空间应用中（如地形建模、林业、城市规划等）被广泛使用。

地理空间库

Geopandas 是一个建立在 pandas 和 shapely 之上的库，旨在高效处理地理空间数据。安装方法为：pip install geopandas。它扩展了 pandas 的功能，支持地理空间数据类型和操作，使用户能够处理矢量数据（点、线、面）并高效地进行地理空间分析。Geopandas 在 GIS、制图和空间数据分析中被广泛使用。

Rasterio 是我们将要安装的最后一个库，安装方法为：pip install rasterio。它用于读取和写入地理空间栅格数据。支持各种标准的栅格格式，如 GeoTIFF 或 JPEG，并提供地理空间元数据、空间参考和坐标转换功能。rasterio 对于卫星影像分析、遥感和 GIS（地理信息系统）应用非常有价值。

这些库各自具有协同作用，可以处理科学计算、数据科学和地理空间应用中的各种数据类型。

设置 IDE

我们设置的最后一步是安装 IDE。我们仍在环境中的命令行界面中，输入：pip install jupyterlab，这将会在我们的环境中安装 jupyterlab。为了清晰使用，我们可以将目录更改为项目的父目录（我们称之为 INTEGRATION），其中包含 CODE 文件夹和 DATA 文件夹：cd C://COURSES/POUX/INTEGRATION。然后，我们将从这个位置启动 jupyterlab，通过在控制台输入：jupyter lab，这将会在您的网页浏览器中打开一个新的本地页面（Chrome 是首选，但 Firefox 或 Safari 也可以）。

做得好！第一阶段成功完成！🎯我们现在准备进入第二阶段：寻找可以组合以供未来使用的数据集，用于我们的 NASA 评分工作流程。🙃

步骤 2. 多模态数据策展

步骤 2. 多模态创建。© F. Poux

我们希望整合多模态数据集。但“多模态”这个词是什么意思呢？它指的是跨不同类型和背景的数据，例如图像、点云、文本和声音……

我们使用各种 2D/2.5D/3D 模态。3D 点云、3D 网格、城市模型、体素、空间栅格、360° 图像、表格数据、矢量数据。© F. Poux

🦝 Nourian：此外，这里有一个关于城市规划的基于网络计算平台的极佳指南： 必备指南 (ResearchGate)

因此，我们在此阶段的目标是识别一个感兴趣的区域，并收集尽可能多的数据，以帮助我们未来的分析。今天选择的感兴趣区域是荷兰东部 Overijssel 省（Twente 地区）Enschede 市的一部分，这里，大学的知识光芒照耀四方。🌞

确定感兴趣的区域。

🦄 Mila：通过融合 3D 地理空间和遥感数据来解锁我们城市的秘密，我们诞生了城市数字双胞胎，这些双胞胎照亮了过去，导航当前，并塑造未来。如果你想更深入地了解如何使用开源工具在网络上进行数据整合，这里有一篇研究论文： 3D 数据集成用于基于 Web 的开源 WebGL 互动可视化 (ISPRS Archives)

现在我们准备获取不同的数据集。为了清晰起见，我将这些数据源分为四类：3D 数据集、空间栅格、矢量数据集，以及其他来源，如下所示。

多模态数据策展工作流程。© F. Poux

3D 数据源

第一步是从一些开放数据仓库获取数据集，这些仓库的数据许可允许我们进行一些实验。在这方面，对于荷兰，可以从一个地方获取 LiDAR 数据、地形模型和栅格图像：geotiles.nl。你可以放大原始图块，并访问各种数据集的下载链接，如下所示。

通过 GeoTiles.nl 提取 AHN4 数据集的 34FN2 图块。© Florent Poux

🦊 Florent: 为了更好地服务于你，所有使用的数据都可以在本节末尾共享的驱动器文件夹中找到。AHN 版本是 AHN4。

你可以探索的第二个获取 CityModels 的地方是 3D BAG，即建筑和地址的 3D 注册（BAG），这是荷兰最详细、开放可用的建筑和地址数据集。它包含多个详细级别的 3D 模型，生成方式是将两个开放数据集结合起来：来自 BAG 的建筑数据和来自 AHN 的高度数据。使用 3D Bag Viewer 通过图块查询，你可以探索和访问建筑物。

BAG Viewer。开放数据的许可为 CC BY 4.0，3DBAG by tudelft3d 和 3DGI。

图块范围与我们从 geotiles.nl 获得的不同，这在我们进行集成阶段时显得有趣。在这个阶段，我们已经手头有了令人兴奋的数据集。现在我们可以继续探索空间栅格数据集。

空间栅格数据源

对于卫星和航空图像，USGS Earth Explorer 是最大的免费数据源之一。它是全球性的，具有友好的用户界面，使访问遥感数据变得简单。如果你需要下载多个数据集，它甚至有批量下载应用。如果这是你第一次使用，你需要执行一个额外的步骤：创建一个账户，但这是免费的且迅速的（我用不到 2 分钟 ⌚）。

USGS Earth Explorer。开放数据属于美国公共领域，Credits USGS。

我从 WebUI 绘制了一个多边形，然后请求获取 Landsat > Landsat Collection 2 — Level 1 组（最新的 Landsat 图像是 L8–9 OLI/TIRS 和 L7 ETM+）。这些集合之间的差异基于数据质量和处理水平。USGS 已根据质量和处理水平 将图像分类为不同级别 (来源)。一旦进入结果部分，你可以查看足迹，然后决定哪个最适合你的需求，如下所示。

提取 LandSAT 图像的多边形。© F. Poux

对于数字高程模型，我建议使用 geotiles.nl 数据服务，因为它提供了 AHN 的最新和最精确的高程模型。我们正在以令人难以置信的速度前进，下一阶段是获取一些优秀的矢量数据集！

矢量数据策展

如果你在 GIS 社区中，我希望展示 OpenStreetMap (OSM) 的力量不会冒犯你的知识基础。

OpenStreetMap 数据策展门户。数据是开放的，遵循开放数据库许可证（ODbL），OpenStreetMap 版权。

OSM 提供了不同的地图和图层，具有众包倡议，这使得它非常详尽，但精确度仍然是一个问题。实际上，OSM 对公众开放，由普通大众创建。这意味着准确性可能会根据创作者及其“制图”技能水平而有所不同。然而，OSM 是一个开放许可证街道级 GIS 数据的宝藏。

我们有几种方法可以下载 OpenStreetMap 数据，以获取一些宝贵的数据。方便的是，甚至有一个 OSM 数据维基页面 列出了所有可用的 OSM 提取数据。我的推荐是使用 Geofabrik 工具。实际上，你可以利用按语义空间范围（例如：国家、州、洲等）组织的数据。你可以快速选择一个地理位置，然后下载 OSM 数据，如下所示。

GeoFabrik 门户以收集矢量数据集。© F. Poux

🦊 Florent：我也喜欢将 OSM 数据下载为形状文件，但稍后会详细说明 😉。 Mila 让我觉得这其中浓缩了大量的知识 Geospatial Data with Python。

现在我们有一些 3D 数据集、一些栅格数据集和矢量形状文件。是时候在互联网世界中挖掘其他宝贵的资源了 💎。

其他来源

好吧，在这里，网络是您的盟友。您可以找到任何与分析相关的信息，从网页到新闻，再到声音和实时数据流。我不会夸大地说天空（或您的带宽 😁）是极限！不过，让我再实际一点，向您推荐一个平台：Mappillary。

Mappillary 是一个提供街景图像和地图数据的平台。您可以使用地图浏览器探索并下载一些 360° 图像和兴趣点。

Mapillary 数据库通过提供的门户。如果您下载图像，它们由 Mappillary 按 CC-BY-SA 许可。

如果您想仅选择一些感兴趣的元素或与 OSM 数据进行交叉验证，您还可以按类别过滤元素。

现在，好消息来了！为了跟随即将到来的代码行，我为您简化了获取这些数据的过程，您可以直接在这个 数据驱动仓库 中找到它们。一旦您在您的 DATA 文件夹中拥有所需的数据，我们可以开始一次愉快的探索性数据分析。

第 3 步：探索性数据分析

第 3 步：探索性数据分析。

在这一阶段，我们已经设置好了 Python 代码和数据，因此我们准备好进入深度集中模式，定时器设置为 15 分钟 ⌚。实际上，我们将使用 Python 探索我们收集的各种数据集。

加载和读取 3D 数据。我们使用 Open3D（点云、网格、体素）、RasterIO 和 Geopandas。© F. Poux

这意味着，对于每种模式，我们将进行读取、分析、协调和分类其内容。通常，这意味着处理库的双向通信。为了保持简洁，我擅自将一些模式归为一组，如下所示。

涵盖的多种模式。© F. Poux

在加载任何内容之前，让我们通过在您的笔记本/脚本中输入以下九行来导入所有已安装的库：

#base libraries
import numpy as np
import matplotlib.pyplot as plt

#3D libraries
import open3d as o3d
import laspy as lp

#geospatial libraries
import rasterio as ra
import geopandas as gpd

3D 点云

让我们从我的门生开始：3D 点云。它们本质上是一组在 3D 空间中表示物理对象或环境的点。这些点来自各种来源，包括 LiDAR、摄影测量、人工智能（是的，你没看错）、扫描设备等。我们手中的第一个数据集来自 Aerial LiDAR AHN 活动，格式为 LAZ 文件。它从 LAS LiDAR 数据交换格式（LAS）扩展而来，是存储 LiDAR 点云数据的压缩对等格式。它支持 2D 和 3D 点数据及其属性，如强度和分类。为了用 Python 读取该文件，我们使用laspy库及其扩展，如下所示：

#Load a las file
las = lp.read("../DATA/34FN2_13.laz")

太好了！我们现在已经将文件加载到las变量中，可以使用laspy函数快速探索，获取可能的属性、red通道的max value或投影信息：

print([dimension.name for dimension in las.point_format.dimensions])
print(np.max(las.red))
print(las.header.vlrs[2].string)

这将产生以下结果：

['X', 'Y', 'Z', 'intensity', ‘return_number’, ‘number_of_returns’, ‘synthetic’, ‘key_point’, ‘withheld’, ‘overlap’, ‘scanner_channel’, ‘scan_direction_flag’, ‘edge_of_flight_line’, ‘classification’, ‘user_data’, ‘scan_angle’, ‘point_source_id’, ‘gps_time’, ‘red’, ‘green’, ‘blue’, ‘nir’, ‘Amplitude’, ‘Reflectance’, ‘Deviation’]
255
COMPD_CS["Amersfoort / RD New + NAP height",PROJCS[“Amersfoort / RD New”,GEOGCS[“Amersfoort”,DATUM[“Amersfoort”,SPHEROID[“Bessel 1841”,6377397.155,299.1528128,AUTHORITY[“EPSG”,”7004"]],AUTHORITY[“EPSG”,”6289"]],PRIMEM[“Greenwich”,0,AUTHORITY[“EPSG”,”8901"]],UNIT[“degree”,0.0174532925199433,AUTHORITY[“EPSG”,”9122"]],AUTHORITY[“EPSG”,”4289"]],PROJECTION[“Oblique_Stereographic”],PARAMETER[“latitude_of_origin”,52.1561605555556],PARAMETER[“central_meridian”,5.38763888888889],PARAMETER[“scale_factor”,0.9999079],PARAMETER[“false_easting”,155000],PARAMETER[“false_northing”,463000],UNIT[“metre”,1,AUTHORITY[“EPSG”,”9001"]],AXIS[“Easting”,EAST],AXIS[“Northing”,NORTH],AUTHORITY[“EPSG”,”28992"]],VERT_CS[“NAP height”,VERT_DATUM[“Normaal Amsterdams Peil”,2005,AUTHORITY[“EPSG”,”5109"]],UNIT[“metre”,1,AUTHORITY[“EPSG”,”9001"]],AXIS[“Gravity-related height”,UP],AUTHORITY[“EPSG”,”5709"]],AUTHORITY[“EPSG”,”7415"]]

第三行给我们一个有趣的概况：数据以Amersfoort / RD New + NAP height参考系统表示。我们记下来以备后用。

🦊 Florent：如你所见，第一行给了我们几个后续使用的属性。很高兴默认情况下我们可以以半结构化的方式存储这么多信息。然而，我们需要整理这些信息的相关性；通常我们默认使用*X*、*Y*、*Z*、*intensity*、*red*、*green*、*blue*。我会把其他的称为额外的奖励😁。

在这个阶段，我们没有自然的方法来可视化数据是否正确；我们必须将其转换为我们喜爱的numpy对象，然后再转换为open3d点云对象，以完成库之间的切换。首先，我们将一些laspy对象属性转换为coords和colors numpy 对象，以保存我们的点云坐标X、Y、Z以及点云颜色：

coords = np.vstack((las.x, las.y, las.z))
colors = np.vstack((las.red, las.green, las.blue))

然后我们将 numpy 对象转换为open3d对象以可视化点云。不幸的是，从laspy到open3d没有直接的方法，尽管有😁。

pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(coords.transpose())
pcd.colors = o3d.utility.Vector3dVector(colors.transpose()/255)
o3d.visualization.draw_geometries([pcd])

🦊 Florent：这是一个很好的练习，可以看到切换库几乎总是会导致在转换之间丢失一些处理时间，并且也暴露出可能的内存错误。因此，当处理复杂的数据集和工作流程时，我建议总是使用有限数量的库。*

在 Open3D 中可视化的 3D 点云数据集。© F. Poux

现在你可以在 Python 中很漂亮地可视化完整的 LiDAR 点云，媲美专业软件！让我们加载另一个点云，以探索另一种广泛使用的文件格式：.ply。PLY 格式用于存储 3D 网格数据及其属性，如颜色和法线。它也适用于带有附加信息的点云。

使用open3d，这非常简单；你可以执行以下代码，将点云填充到pcd_itc变量中，并使用o3d.visualization函数绘制点云：

#Load a point cloud and visualize
pcd_itc = o3d.io.read_point_cloud("../DATA/ITC_outdoor.ply")
o3d.visualization.draw_geometries([pcd_itc])

在 Open3D 中可视化的 3D 室内点云。© F. Poux

这个点云没有显示任何元数据，这意味着我们将认为它以本地参考系统表示，需要以某种方式与参考数据集进行配准。

🦊Florent: 剧透警告！点云数据集就像一个特殊的中心岩石，尤其在处理多模态数据时。实际上，正如我稍后将展示的，它可以作为标准参考，将所有其他数据集与之关联。我已经提醒过你有剧透了！

是时候转向 3D 体素🧊

3D 体素

体素是另一种 3D 数据格式。体素本质上是表示空间体积的 3D 像素。它们非常有趣，你可以在以下案例研究中查看它们的实际用途：

## 3D Python 工作流用于 LiDAR 城市模型：一步一步的指南

解锁 3D 城市建模应用的最终指南。该教程涵盖了 Python……

towardsdatascience.com ## 如何使用 Python 自动化 3D 点云的体素建模

实用教程，将大型点云转换为 3D 体素🧊，使用 Python 和 open3d。解锁自动化工作流……

towardsdatascience.com

🦝 Pirouz: 体素非常像乐高模型。查看这篇关于空间数据体素化并深入研究数字乐高的论文： 地理空间应用的体素化算法 (ResearchGate)。如果你对体素情有独钟？那么在 Google 上搜索“体素艺术”，并查看库 topoGenesis (GitHub)。

现在，我们将使用以下代码行读取和可视化体素数据集：

vox_read = o3d.io.read_voxel_grid("../DATA/34FN2_13_vox.ply", format='auto')
o3d.visualization.draw_geometries([vox_read])

在 Open3D 中可视化的 3D 体素数据集。© F. Poux

在这个阶段，我们阅读、分析，并确保有一个单一的库来处理所有这些数据（numpy）和一个用于 3D 可视化的库（open3d）。

是时候在我们的 Python 实验中加入城市模型了😁

城市模型

城市模型通常遵循一种表达标准：CityGML 模型。CityGML 是一种基于 XML 的数据格式，用于表示城市环境的三维几何。CityGML 模型可以包含有关建筑物、道路、桥梁和其他基础设施的信息。城市规划师、建筑师和工程师常常使用这些模型来模拟和分析城市环境。一个优秀的获取地点，我之前故意没有提到的，是由我的同事 Olaf 精心策划的 GitHub 仓库（任何冰雪奇缘的引用禁止 ⛄）：CityGML Models。

🦊 Florent: 基本上，在 Python 中，我们必须确保将这些模型转换为三维网格，同时保留所需的语义/拓扑或其他信息。一个友好的 Python 工具，同样由 TUM* 提供，可以在这里找到：* citygml2obj。但如果你使用数据，我已经为你完成了繁重的工作，因此继续阅读下一个阅读方式：三维网格。

3D 网格

接下来，让我们谈谈网格。三角网格由顶点、连接顶点的边缘和完成对象“包络”的三角面组成。它们使我们能够描绘三维物体的形状。网格可以使用三维建模软件创建，也可以通过点云、城市模型和专用算法生成，甚至是人工智能！数据集采用 OBJ 格式，通常用于导出表示几何体和纹理坐标的三维网格模型。

🦝 Pirouz: 如果你想知道为什么这些花哨的词汇被用来代替点、线和多边形，你需要多了解一点另一个花哨的词汇：拓扑。*

🦊 Florent: 你可以阅读 Pirouz 写的关于拓扑的故事，这可能是最简单的拓扑故事之一： 几何和拓扑基础 (ResearchGate)

要加载网格并获取其范围，我们将编写以下几行：

mesh = o3d.io.read_triangle_mesh("../DATA/10–976–660-LoD22–3D.obj")
mesh.get_axis_aligned_bounding_box()

从输出中我们看到坐标接近于 Amersfoort / RD New 所提供的。这令人欣慰；我们现在在这个参考系统中有了另一个数据集。

从那里开始，为了可视化网格，我们首先计算一些顶点法线，并使用 open3d 的相同可视化函数来进行可视化：

mesh.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh])

在三维方面，我们可以开始了！让我们用 Python 和 rasterio 探索我们的空间栅格。

空间影像

我们的影像采用 GeoTIFF 格式，这是 TIFF 格式的扩展，包含地理空间元数据，适用于栅格高程数据（DSM/DTM）和影像。要打开和分析我们拥有的空间影像，我们只需使用 rasterio 的以下代码行：

sat_image = ra.open("../DATA/RGB_34FN2.tiff")
sat_image_array = sat_image.read(1)
print(sat_image.meta)

这将产生以下输出：

{‘driver’: ‘GTiff’, ‘dtype’: ‘uint8’, ‘nodata’: None, ‘width’: 20002, ‘height’: 25002, ‘count’: 3, ‘crs’: CRS.from_epsg(28992), ‘transform’: Affine(0.25, 0.0, 254999.75,
 0.0, -0.25, 475000.25)}

有趣的是，我们可以清楚地看到这是一个20002x25002的图像，具有三个通道（R、G和B），其表达在CRS 28992中，这再次对应于 Amersfoort！这太棒了！

现在，为了使用 numpy 和 matplotlib 绘图，下面是获取图像的最小代码行，以避免内存错误，结果如图所示：

plt.imshow(sat_image_array)
plt.axis("off")
plt.show()

与此同时，rasterio 还提供了[show()](https://rasterio.readthedocs.io/en/latest/api/rasterio.plot.html#rasterio.plot.show)函数，以执行日常任务，例如将多波段图像显示为 RGB，并用适当的地理参考范围标记坐标轴。这使得绘图变得更加简单，如下所示：

from rasterio.plot import show
cir_image = ra.open("../DATA/CIR_34FN2.tiff")
show(cir_image)
show(sat_image)

使用 matplotlib 和 rasterio 的绘图函数可视化的各种光栅。© F. Poux

这真是太棒了！我们有多个 3D 数据集（3D 点云、体素、3D 网格）和卫星图像（包括红外和 R、G、B），我们可以用 numpy 处理这些数据，并且大多数可以用一个通用 CRS 表示。让我们转向其他光栅，重点关注高程数据（2.5D）。

高程光栅（DSM，DTM）

高程模型通常以光栅文件的形式出现，其中每个像素都有一个可以转换为其高程的值。这就是为什么我们称之为 2.5D 数据，因为它是一个纯粹的俯视图或单点视角，当它不在任何 GIS 投影系统中时也称为深度图像。我们可以用相同的代码行导入高程模型：

dtm = ra.open("../DATA/DTM5_34FN2.TIF")

现在我们有了一个保存了高程数据的dtm变量，我们可以通过dtm.meta提取其整体元数据，或者更专业地使用dtm.shape、dtm.crs、dtm.bounds、dtm.overviews(1)：

print(dtm.meta)
print(dtm.shape, dtm.crs, dtm.bounds, dtm.overviews(1))

这给我们带来了以下结果：

{‘driver’: ‘GTiff’, ‘dtype’: ‘float32’, ‘nodata’: 3.4028234663852886e+38, ‘width’: 1000, ‘height’: 1250, ‘count’: 1, ‘crs’: CRS.from_epsg(28992), ‘transform’: Affine(5.0, 0.0, 255000.0,
 0.0, -5.0, 475000.0)}
(1250, 1000) EPSG:28992 BoundingBox(left=255000.0, bottom=468750.0, right=260000.0, top=475000.0)

再次令人感兴趣的是，我们看到这个数据集使用的是相同的 CRS Amersfoort / RD New。

为了绘制图形并了解我们正在处理的内容，我们可以将数据集转换为 numpy 数组，并使用索引选择其中的一部分到选择变量中：

dtm_array = dtm.read(1).astype(‘float64’)
selection = dtm_array[300:800,300:800]

最后，我们使用两种方式绘图：rasterio和matplotlib，分别绘制全尺度数据集和放大的选择区域：

#For the dtm plot with rasterio
show(dtm)

#For the dtm zoomed-in plot with matplotlib
fig, ax = plt.subplots(1, figsize=(15, 15))
ra.plot.show(selection, cmap=’Greys_r’, ax=ax)
plt.axis(‘off’)
plt.show()

使用 RasterIO 和 matplotlib 绘制的 DTM 图。© F. Poux

最后，我们不能在不深入探讨矢量空间数据的情况下讨论光栅数据集。

空间矢量数据

在这个阶段，必须一次性明确，如果你看到的是像素，那么你并不是在谈论矢量数据。相反，矢量数据集由顶点和路径组成，这些顶点和路径以点（X，Y 坐标）、线（连接这些标记的顶点）和多边形（连接顶点以闭合由线组成的路径）的形式存在。

GIS 系统中的矢量数据。© F. Poux

🌱 成长：考虑到这一点，你的直觉把 3D 点云放在哪里？体素模型？3D 网格？这些都是有深远意义的有趣问题，但有助于更好地定位 2D 与 3D 的协同作用。

有趣的是，制图师（不仅仅是制图师）将这些用作“符号”来在地图上描绘现实世界的组件。这意味着他们必须不断确定实体所代表的“详细程度” (LoD)，这赋予这些实体如此多的符号意义（一个点可以是一个国家、一个城市、一个公民或一个特定地点，……）。

为了处理这些矢量数据集，我们主要使用一种开放规范——shapefile 格式。这允许我们空间描述几何形状并附加一些额外的信息。在我们的计算机上，shapefile 文件格式实际上是多个文件的集合，用于表示地理数据的不同方面（参考）：

• .shp — 形状文件，包含特征几何本身。

• .shx — 形状索引格式，包含特征几何的定位索引。这对于大型文件非常有用，因为它利用巧妙的“索引”实现快速搜索。

• .dbf — 属性格式，以表格方式（dBase IV 格式）存储每个形状的各种属性。

在顶部，我们可能会有与 shapefile 格式配套的可选文件。最值得注意的是“.prj”文件。这个额外的文件实际上定义了坐标系统和任何必要的投影信息。因此，拥有这些文件使我们能够使用 geopandas 加载 Overijssel 区域的街道数据：

vector_data = gpd.read_file("../DATA/gis_osm_roads_free_1.shp")

矢量数据输出。© F. Poux

现在让我们通过探索其坐标参考系统（CRS）来进一步了解其投影系统：

print(vector_data.crs)

这导致了以下结果：

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

如我们所见，该数据集以 ESPG 4326 表达，这代表 WGS84. 因此，与我们之前的数据集有所不同。

Florent：geopandas 是 pandas 和 shapely 的完美结合。因此，如果你对这两者都很熟悉，使用 geopandas 将会非常轻松！一些有用的命令是 vector_data.columns 或 vector_data.describe() ，可以让你快速了解其内容。

但现在让我们通过可视化我们的数据集来查看其内容：

vector_data.plot()

这导致了以下图示：

来自 OpenStreetMap 的矢量数据。© F. Poux

这非常有趣！我们还可以使用 Matplotlib 来显示矢量数据。该区域人口密集，道路网络广泛且非常密集！现在让我们从另一个数据源中提取一些兴趣点。

矢量：兴趣点

数据集包中包含一个 GeoJSON 文件，这是一个轻量级且流行的矢量数据格式，具有几何体和属性，适用于基于 Web 的应用程序和数据交换。要加载该文件，我们也使用 geopandas，代码如下：

trashcan_data = gpd.read_file("../DATA/mapillary-trashcan_points.json")

我们分析其 CRS：

trashcan_data.crs

这给我们的答案是它仍然在 WGS84 CRS 中。最后，我们绘制以检查任何异常：

trashcan_data.plot()

关注点。© F. Poux

没有更多上下文，这可能有点平淡；因此，首先选择叠加栅格影像。但在此之前，让我们探索另一个可以尝试的数据集：360° 影像！

360° 影像

从 mappillary 平台，我们可以提取一些附有特定许可的 360° 图像。你仍然可以使用 geopandas 来读取这些影像：

image = ra.open("../DATA/An8XJdYJkSXycbyDWiVyw8dfVi2IzZ-jr9z7IxneeEXnOJPt0K1C89cMdXNkl5FZT0x6aVuPFRpda6eWV9fcnpgJkPfWjrd7k9harUP48csXGFf2azE1qF7FhkYK4h__j9t6Vd3aUfKcEdlqF_rsVO4.jpg")
print(image.meta)

这将输出：

Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
{'driver': 'JPEG', 'dtype': 'uint8', 'nodata': None, 'width': 2048, 'height': 1024, 'count': 3, 'crs': None, 'transform': Affine(1.0, 0.0, 0.0,
       0.0, 1.0, 0.0)}

在这里，我们可以看到（当然）这张图像没有地理变换数据。这是正常的：这不是空间栅格！然而，如果我们愿意，我们可以添加一个位置，以便数据点表示拍摄照片的位置，这暗示了通过数据整合来叠加数据模态。但在去到那里之前，让我们绘制图像：

ra.plot.show(image)

结果如下：

来自 Mapillary 的 360° 影像。© F. Poux

我们现在已经成功加载了 3D 点云、作为 3D 网格的 3D 城市模型、3D 体素数据集、空间栅格（卫星影像、航拍影像和高程栅格）、矢量数据集（线条和点）以及 360° 影像。第一个目标现在显然已经成功！但我们需要检查如何将这些数据整合在一起，首先需要检查的是它们是否都表示在相同的坐标参考系统（CRS）中。为此，我将我们的分析步骤的结果汇总在下面的表格中。

| Data | Lidar Point Cloud                | TLS Point Cloud | 3D City Mesh                     | Satellite Imagery   | Digital Terrain Model | Streets (Vector)         | Trashcans                      |
|------|----------------------------------|-----------------|----------------------------------|---------------------|-----------------------|--------------------------|--------------------------------|
| Name | 34FN2_13.laz                     | ITC_outdoor.ply | 10-976-660-LoD22-3D.obj          | RGB_34FN2.tiff      | DTM5_34FN2.TIF        | gis_osm_roads_free_1.shp | mapillary-trashcan_points.json |
| CRS  | Amersfoort / RD New + NAP height | None            | Amersfoort / RD New + NAP height | Amersfoort / RD New | Amersfoort / RD New   | WGS 84                   | WGS 84                         |
| CRSC | 7415                             |                 | 7415                             | 28992               | 28992                 | 4326                     | 4326                           |

我们看到混合了 EPSG:7415、28992、4326 或缺失 EPSG。下一阶段是澄清并为所有数据集使用统一的系统。

数据注册和重投影

步骤 4：注册和重投影

在整合之前，我们必须确保所有数据集都处于兼容的格式和坐标参考系统（CRS）中。如果不是这样，我们会执行数据重投影，将它们带入一个共同的 CRS。这通过四个主要阶段实现：（1）选择参考系统，（2）地理配准主要数据集，（3）重投影其他数据集，以及（4）严格对齐本地数据集，如下所示。

已覆盖 4 个步骤工作流。© F. Poux

让我们逐步了解这些阶段。

选择参考系统

好的，让我们深入探索参考系统的激动人心的世界吧！想象一下参考系统就像是引导你的数据穿越地球广阔景观的 GPS 坐标。参考系统定义了一组规则和参数，以一种对我们这些凡人有意义的方式来表示地球表面。这就像选择一种数据流利使用的秘密语言，让它找到自己在世界上的位置。那么，我们如何选择正确的参考系统呢？让我们务实一些，首先尝试勾画项目的范围和关注区域。

投影经典案例：一个标准线，一个标准圆锥和一个标准点。© F. Poux

不同的参考系统在全球不同区域表现更佳。例如，Amersfoort / RD New 对于荷兰的本地项目可能是一个绝佳的选择。值得信赖的 WGS84 或全球横轴墨卡托（UTM）可能是你全球项目的最佳伙伴。因此，我们使用 Amersfoort / RD New + NAP 高程（用于高程信息）作为我们的 CRS。

记住，在选择参考系统时，让你的项目性质和地区来指导你的决策。选择一个能与数据语言相通的系统，并为你的地理空间工作带来和谐。🗺️

数据地理参考

好的，我们定义了 CRS；这一步确保数据标记为 EPSG:7415 的 3D 数据和 EPSG:28992 的 2D 数据。在分析时，我们可以检查这些数据集是否已经一致：

• LiDAR 3D 点云

• 栅格影像（包括空间栅格和高程栅格）

然而，体素数据集和城市模型数据集中的网格看起来似乎在相同的 CRS 中，但没有元数据。因此，我们迅速叠加这些数据集以检查任何可能的不一致：

o3d.visualization.draw_geometries([pcd,vox_read,mesh])

结果如下：

点云、建筑物和体素数据集被合并在一起并可视化。© F. Poux

🦊弗洛朗：一切都对齐得很好，如图所示。各种数据集被故意裁剪以显示它们之间的重叠。你还可以看到点云在表示建筑物的 3D 模型上有些许粗糙和“分辨率”差异。

剩余的数据集需要处理如下：

• 表示为 WGS84 的矢量数据集（包括 OSM 和 Mappillary）

• ITC 点云数据集没有地理参考。

让我们首先处理数据投影。

数据重投影

空间数据重投影是一个基本过程，涉及将空间数据从一个坐标参考系统转换到另一个系统，通常是为了匹配特定空间分析的投影和坐标单位。重投影确保具有不同投影的数据集可以准确地叠加、集成或分析。

该过程涉及将地理坐标（纬度和经度）从一个基准转换到另一个基准的数学计算，考虑参数如比例、旋转和畸变。空间数据重新投影对于实现数据一致性、促进不同数据集之间的互操作性，以及在各种制图系统和应用中进行准确的地理空间分析至关重要。

在参考系统之间重新投影空间数据时，必须仔细考虑几个重要因素，以确保结果的准确性和意义。以下是一些关键考虑因素：

1. 坐标系统和投影：了解源系统和目标系统的坐标系统和投影。确保它们兼容；如果不兼容，选择合适的转换方法。

2. 基准和椭球体：检查源系统和目标系统中使用的基准和椭球体。基准的差异可能会导致坐标的显著偏移。如有需要，应用基准转换以正确对齐数据。如果你感到有些困惑，这里有一个不错的基准和椭球体讲座 (tamia.edu)

3. 准确性/精度：虽然这些指向两个不同的特征（这扩展了我们文章的范围），了解你特定分析所需的“水平”很重要。重新投影可能会引入一些错误，特别是在大规模转换中。确实，任何数据“错误”都有可能影响我们分析的结果。

4. 畸变：不同的地图投影可能会引入形状、面积、距离或角度的畸变。注意这些畸变及其对数据解释的影响。

5. 覆盖区域：某些投影适合特定区域，但可能不适合全球数据集。选择一种在整个覆盖区域内保持数据属性的投影。

6. 元数据和文档：记录重新投影过程，并记录应用于数据的转换。妥善记录源坐标系统和目标投影以备将来参考。

通过仔细考虑这些因素，我们确保重新投影准确执行，得到的数据适用于我们的预期应用。必须注意在重新投影过程中可能出现的错误和伪影，并验证结果以保持数据的质量和可靠性。在我们的案例中，我们需要重新投影两个矢量数据集。为此，我们使用以下步骤：

trashcan_data_georeferenced = trashcan_data.to_crs(‘epsg:28992’)
trashcan_data_georeferenced.plot(color=’red’)
vector_data_georeferenced = vector_data.to_crs(‘epsg:28992’)
vector_data_georeferenced.plot(edgecolor=’green’)

这导致了以下结果：

生成的图表。© F. Poux

🦊 Florent: 3D 重投影也是可能的，但目前超出了本指南的范围。然而，您可以浏览课程库，其中会找到一些很好的代码和示例，用于在“Segment Anything”语境下进行 3D 重投影。

注意我们图表坐标轴上的变化！在这一阶段，我们似乎只有一个数据需要注册：ITC 的 3D 点云。

数据刚性注册

3D 数据注册经典工作流程。© F. Poux

在 3D 集成中，数据注册变得至关重要。点云注册方法通常由两个阶段组成：粗略对齐，将两个点云相对快速地定位得较近。然后是精细注册，如迭代最近点（ICP）或基于特征的注册，以更高的精度对齐多个点云。为了给您一个关于全局注册的提示，我们可以通过选择三对公共点的列表来对齐 ITC 点云，以获得初步估计，然后用 ICP 进行精细化。这允许我们获得叠加的点云。这些内容大大扩展了当前指南的范围，将在另一个会议中讨论。😉

我们现在拥有一个在一致的坐标参考系统（CRS）中对齐的数据集。下一阶段是看看我们是否能“提炼数据”，以便从中建立一个智能的分析工作流程。

第 5 步。数据处理、转换和融合

第 5 步。数据预处理

这一阶段也可以在注册和重投影部分之前完成。实际上，由于我们主要按数据集工作，这是一种解决方案，并可能影响（无论是好是坏）前面步骤的结果。

空间数据集成将不同类型的空间数据（2D、3D 或 2.5D）结合起来，以创建一个统一且全面的区域数据重叠表示。这一过程使我们能够利用每种数据类型的优点，从集成的数据集中获得更有价值的见解。以下是空间数据集成的一些处理策略的快速视图：

3D 数据预处理的主要领域：数据清理、数据转换、数据减少和数据丰富。© F. Poux

现在，让我们探讨一下这四个主要阶段的实际代码实现。

数据清理（栅格数据集）

第一阶段是仔细处理缺失或错误的数据点。可以使用插值或外推等数据清理技术来填补空白或替换错误值。在我们的案例中，我们可以从数字地形模型（DTM）中填补空白的栅格值。

我通过对每个空白区域的相邻像素进行插值，试图用以下代码片段为我们的数据集提供更好的结构：

#We import a specific function
from rasterio.fill import fillnodata

#We transform to a numpy array
dtm_array = dtm.read(1)

#We interpolate
interpolated_dtm = fillnodata(dtm_array, mask=dtm.read_masks(1), max_search_distance=100, smoothing_iterations=0)

#We plot the results
fig, ax = plt.subplots(1, figsize=(15, 15))
ra.plot.show(interpolated_dtm, cmap='Greys_r', ax=ax)
plt.axis('off')
plt.show()

🌱 Growing: 你会如何使用这些来提取你区域的等高线？

这会导致以下结果：

DTM 可视化结果。© F. Poux

数据转换（3D 点云）

除了文件格式之间的转换，数据转换还指的是确保结构和内容完整性的各种过程，以便于一个或多个处理步骤。例如，在许多特征工程任务中，这对于提升 3D 机器学习模型的性能至关重要。

一种关键的数据转换方法是缩放，确保数据集中的所有值都在特定范围内，例如 0 到 1。要对点云执行此操作，我们按以下步骤进行：

#We compute a MinMaxScaler bounds
coords_itc = np.array(pcd_itc.points)
min_itc = np.min(coords_itc, axis=0)
max_itc = np.max(coords_itc, axis=0)

#MinMaxScaling
coords_itc_mmscaling = (coords_itc - min_itc) / (max_itc - min_itc)

这将我们原始的点云转换为所有点的坐标在 0 和 1 之间。因此，如果我们要绘制它们沿三个轴的分布，我们会使用以下内容：

n_bins = 20
fig, ax = plt.subplots(1, 3, sharey=True, tight_layout=True)
ax[0].hist(coords_itc_mmscaling[:,0], bins=n_bins, color = "skyblue")
ax[1].hist(coords_itc_mmscaling[:,1], bins=n_bins, color = "salmon")
ax[2].hist(coords_itc_mmscaling[:,2], bins=n_bins,color = "purple")

获得：

3D 点云分布分析。© F. Poux

如你所见，我们的点云现在拥有 X、Y 和 Z 坐标，这些坐标的范围在 0 和 1 之间。

🍇 注意：这在另一种口味上很有趣。实际上，我们会根据 Z 分布的模式来寻找线索，这样的 Z 驱动算法是有意义的，例如，用于区分屋顶和地面点。

数据减少：裁剪掩码，点云采样

数据减少包括各种技术，通过这些技术数据被简化到最简单的形式，以便进行最佳的分析任务。

其中之一是数据裁剪：减少数据集的空间范围。一个明确的例子是在栅格数据集上执行，以限制内存占用，当我们想要关注更小的区域时加载一个巨大的区域。为此，使用 rasterio，我们首先需要创建一个 geopandas 边界框，作为过滤掩码：

from shapely.geometry import Polygon

BBlidar=[np.min(coords, axis=1),np.max(coords, axis=1)]

minx=BBlidar[0][0]
miny=BBlidar[0][1]
maxx=BBlidar[1][0]
maxy=BBlidar[1][1]
areabbox = gpd.GeoDataFrame({'geometry':Polygon([(minx,maxy),(maxx,maxy),(maxx,miny),(minx,miny),(minx,maxy)])},index=[0],crs="EPSG:28992")

然后，我们得到一个漂亮的区域框，将其用作原始文件（或任何栅格文件）的掩码，同时确保将原始文件的元数据复制到裁剪后的版本中：

from rasterio.mask import mask

in_raster = ra.open("../DATA/CIR_34FN2.tiff")

# Do the clip operation
out_raster, out_transform = mask(in_raster, areabbox.geometry, filled=False, crop=True)

# Copy the metadata from the source and update the new clipped layer 
out_meta=in_raster.meta.copy() 
out_meta.update({
    "driver":"GTiff",
    "height":out_raster.shape[1], # height starts with shape[1]
    "width":out_raster.shape[2], # width starts with shape[2]
    "transform":out_transform})

# Plot the CIR raster
ra.plot.show(out_raster)

这会导致以下结果：

裁剪 CIR 图像。© F. Poux

另一种明确的数据减少技术是数据采样。在我们的案例中，就是试图减少点云中的点数。这可以通过几种方法来实现（这些方法也超出了本文的范围），其中一种方法是使用预定义的体素数据结构：我们希望每个体素保留一个最佳候选点：

voxel_size = 1

pcd_downsampled = pcd.voxel_down_sample(voxel_size = voxel_size)
o3d.visualization.draw_geometries([pcd_downsampled])

这会导致以下下采样点云：

3D 点云下采样结果。© F. Poux

最后，我们可以通过各种方式丰富我们的数据集。

数据丰富：将 POI 近似值注入点云

根据应用需求，可以采用各种融合技术来合并数据集。对于 2D 和 2.5D 集成，可以使用基于栅格的方法，如加权平均或多数投票，来结合多个数据层。在 3D 集成中，使用点云融合技术，如合并、平均或体素化，来创建一个代表整个区域的单一合成 3D 点云。

🦊 Florent: 使用数据增强（融合）技术，我们将来自多个来源的信息结合起来（这些来源可以是基于网络的、来自数据采集任务的，或任何相关的数据收集方式）。这使得我们能够创建更全面、更准确的基础现象表示。这些技术包括栅格加权平均（涉及对来自多个栅格层的像素值进行加权平均，其中权重代表每个栅格层的重要性或可靠性）、多数投票（通过将多数类别分配给每个像素位置来结合分类栅格数据，对土地覆盖分类非常有帮助）。但这些内容留到以后再讲。 😉

扩展数据集的快速方法是计算额外的特征。例如，这可以通过从点云中提取从邻域计算出的法线来完成，这样可以提取最佳拟合平面，然后计算法线。我们使用以下代码自动完成此操作：

nn_distance = np.mean(pcd.compute_nearest_neighbor_distance())
print(nn_distance)

#setting the radius search to compute normals
radius_normals=nn_distance*4
pcd_downsampled.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=radius_normals, max_nn=16), fast_normal_computation=True)

# Visualizing the point cloud
pcd_downsampled.paint_uniform_color([0.6, 0.6, 0.6])
o3d.visualization.draw_geometries([pcd_downsampled])

如下所示，我们可以在互动窗口查看器中按下“n”键来可视化我们新的点云，无论是否显示法线：

从 3D 点云中计算和可视化法线。© F. Poux

步骤 6. 多模态数据可视化

步骤 6. 多模态数据可视化。

多模态数据可视化是一种强大的方法，可以无缝集成我们多样的地理空间数据，包括矢量数据、空间栅格影像、数字表面模型（DSM）、数字地形模型（DTM）、点云和城市模型。

用于 3D 数据可视化的库：Matplotlib、Open3D 和 RasterIO（2D/2.5D）。© F. Poux

数据类型的全面融合使我们能够对环境有更全面的理解，为更明智的决策和复杂的分析奠定基础。现在，让我们深入探讨一个 Python 解决方案，展示如何实现这一卓越的可视化效果！

诀窍在于了解我们每个库的优点和限制！事实上，我们将专注于 Open3D，以检查数据模态（网格、3D 点云和体素）之间的一致性，并使用点云数据集作为标准参考，然后将所有其他数据集链接到它。

点云作为标准参考框架。© F. Poux

因此，第一步是分析 3D 模态：

o3d.visualization.draw_geometries([pcd,vox_read,mesh])

组合的 3D 数据集。© F. Poux

🦊 Florent：虽然图像可能令人困惑，但请注意，当你在一张图中可视化这三者时，什么也没有出错。事实上，数据集的范围不同，使我们能够更好地在视觉上区分模态。

太好了！接下来，我们使用 Numpy 对点云数据集进行自上而下的视图，以检查 X-Y 一致性，例如，与矢量数据一起：

# Plot the Georeferenced Shapefile and Point Cloud together using matplotlib
vector_data_clipped_pc = vector_data_georeferenced.cx[257000:258000, 471200:472400]
fig, ax = plt.subplots(figsize=(15, 15))
vector_data_clipped_pc.plot(ax=ax, color='yellow', edgecolor='yellow')
# Plot the point cloud
ax.scatter(coords[0][::100], coords[1][::100], c=coords[2][::100])
# Customize the plot (add labels, titles, legends, etc. if needed)
ax.set_xlabel('X (meters)')
ax.set_ylabel('Y (meters)')
ax.set_title('Shapefile and Point Cloud Visualization')
# Show the plot
plt.show()

这明显地将矢量数据集叠加到我们的点云上，为使这两个数据集相互交流开辟了可能性：

集成矢量和栅格数据集。© F. Poux

然后，我们通过将矢量数据集（重投影）与空间影像连接，继续我们的链接：

#Clipping vector data
vector_data_clipped = vector_data_georeferenced.cx[255000:260000, 469000:475000]

#Plot the Raster Imagery with Vector Dataset
fig, ax = plt.subplots(figsize=(15, 15))
ra.plot.show(sat_image, ax=ax)
vector_data_clipped.plot(ax=ax, facecolor='none', edgecolor='green')

#Zooming-in
fig, ax = plt.subplots(figsize=(15, 15))
plt.axis([258000, 258500, 471000, 471500])
#The line above actually define the extent of the plot using the coordinates in the CRS.
ra.plot.show(sat_image, ax=ax)
vector_data_clipped.plot(ax=ax, facecolor='none', edgecolor='cyan')

这导致了可视化结果，需要稍微放大，如代码中所示，以更好地定位两个数据集的适配情况。

栅格和矢量数据集的集成结果，带有 LiDAR 重投影。© F. Poux

对 Mappillary 的位置重复相同的过程，使用如下代码。

#Overlaying Raster with Vector data
fig, ax = plt.subplots(figsize=(15, 15))
plt.axis([257600, 260000, 471000, 471500])
ra.plot.show(sat_image, ax=ax)
trashcan_data_georeferenced.plot(ax=ax, facecolor='none', edgecolor='green')

Mappillary 的兴趣点视图。© F. Poux

最终，我们希望使用 Numpy 将栅格、矢量和点云模态叠加在一个图中：

# Plot the Shapefile and Point Cloud together using matplotlib
vector_data_clipped_pc = vector_data_georeferenced.cx[257000:258000, 471200:472400]
fig, ax = plt.subplots(figsize=(15, 15))

# plot the sattelite imagery
ra.plot.show(sat_image, ax=ax)

# Plot the shapefile
vector_data_clipped_pc.plot(ax=ax, color='blue', edgecolor='blue')

# Plot the point cloud
ax.scatter(coords[0][::100], coords[1][::100], c=coords[2][::100])

# Customize the plot (add labels, titles, legends, etc. if needed)
ax.set_xlabel('X (meters)')
ax.set_ylabel('Y (meters)')
ax.set_title('Shapefile and Point Cloud Visualization')

# Show the plot
plt.show()

这段小代码（注意其中使用的快捷方式和聪明的技巧以缩短代码）使我们能够对我们各种模态的数据集成做出最终的说明。

使用 DTM 的 Shapefile 和点云可视化。© F. Poux

我们接近我们系统方法的结束！

通过仔细应用这些处理技术，空间数据集成使我们能够创建环境的整体表示，促进更好的决策和洞察力，适用于广泛的应用。不过，我们还必须考虑最后一步：分享新处理的数据。

第 7 步：数据共享

第 7 步：数据共享。© F. Poux

空间数据导出在地理空间数据工作流中至关重要，使我们能够无缝地共享、可视化和分析 2D 和 3D 数据模态。在这方面，我们可以从更高的角度，几乎是空中视角，以更好地定位数据共享的影响范围。我选择与学术界和发布新科学发现的证据相关联。确实，我下面所示的这个过程是一个道德且强健的研发周期的关键。正如你所看到的，标准格式的数据导出，我们在这里介绍，是在推动自动化解决方案之前数据浏览的核心环节。

3D 数据共享工作流。我们从实验开始，以填充一个数据库。然后使用该数据库来探索或导出结果，以进行筛选和发布。结果随后用于重新填充和更新 3D 数据库。© F. Poux

因此，导出过程必须创建各种适用于应用程序和软件环境的文件格式。每种文件格式都有其特定性，满足不同的使用案例和先前表达的需求。对于深入的 3D 数据导出，我建议参考以下文章：

## LiDAR 城市模型的 3D Python 工作流：逐步指南

解锁 3D 城市建模应用程序的精简工作流的终极指南。教程涵盖了 Python…

towardsdatascience.com

了解如何在矢量和栅格数据集上进行相同的操作是至关重要的。像往常一样，越简单越好：让我们专注于最后的导出步骤。要导出栅格图像，你可以使用 rasterio 执行以下操作：

with ra.open("../RESULTS/CIR_34FN2_cropped.tiff",'w',**out_meta) as wf:
    wf.write(out_raster)

关于矢量数据集，geopandas 使得从我们的处理阶段导出地理参考和裁剪结果变得轻而易举：

vector_data_clipped_pc.to_file("../RESULTS/roads.shp")

最后，你可以将这些数据导入到像 QGIS 或 CloudCompare 这样的软件中，结果如下：

3D 点云、3D 体素、3D 城市模型、3D DTM、栅格和矢量数据集的叠加。它们与提出的方法论进行了整合。© F. Poux

这是 3D、栅格和矢量数据集的完美层叠！

结论

我们开始了一段激动人心的旅程，探索使用 Python 进行 3D 数据集成的世界。在这份全面的指南中，我们揭示了将各种 3D 数据模态（如点云、网格、城市模型、DSM、DTM 和体素）结合起来，以创建空间环境的统一和全面表示的巨大潜力。

本 3D 数据集成指南中介绍的工作流。© F. Poux

总结一下，我在下面列出了我们所涵盖的七个步骤的关键要点：

• 在 Python 中设置多模态编码环境需要配备尽可能少的强大库。

• 在获取数据集时，最好仔细考虑通过网络接口提供的具有明确许可选项的开放数据集，这样您可以轻松收集相关样本。

• 独立地对每种模态进行分析并评估其主要特性（坐标参考系统、精度、分辨率等）是一项优秀的实践，建议通过快速分析方案来实现。

• 掌握坐标参考系统，并了解如何从一个系统转换到另一个系统，以及特定的转换意味着什么，对于具有空间数据集的可扩展工作流是强制性的。

• 预处理算法和技术通常是优化和组织的数据集与具有混乱管理系统的数据仓库之间的关键所在。

• nD 数据可视化对成功的数据集成工作流至关重要，并要求您仔细考虑哪个数据集/数据模态充当您的标准锚点。

• 数据共享作为最终阶段，允许确保一个从 A 到 Z 的系统，考虑实际和操作性的因素，即使在学术界或研发活动中，生产阶段也可能被忽视。

在结束本指南时，我无法不对 3D 数据处理的未来感到兴奋。目前，凭借 Python 的多功能性和不断增长的地理空间技术，我们有了更多的机会来更好地理解复杂的空间现象（并对此采取行动）。

因此，您可以期待更深入的高效算法、创新可视化技术和与前沿技术无缝集成的指南。通过探索我们能够利用集成的 3D 数据工作流完成的任务，我们为新应用铺平了道路，其中一些应用很快就会到来。这使我说，您正走在塑造空间分析和可视化未来的正确轨道上，改变我们感知和互动的方式。

参考文献

1. Cárdenas, Ivan L., Morales, Luis Rodrigoandrés, Koeva, Mila, Atun, Funda, & Pfeffer, Karin. (2023 年 8 月 31 日). 《用于生理等效温度计算的数字双胞胎指南》。Zenodo. doi.org/10.5281/zenodo.83064562.

2. Kumalasari, D.; Koeva, M.; Vahdatikhaki, F.; Petrova Antonova, D.; Kuffer, M. 《规划可步行城市：面向数字双胞胎实施的生成设计方法》。遥感. 2023, 15, 1088. doi.org/10.3390/rs150410883.

3. Rajan, V.; Koeva, M.; Kuffer, M.; Da Silva Mano, A.; Mishra, S. 《通过多时相遥感数据对过去和现在的沙贾汉纳巴德进行三维建模》。遥感. 2023, 15, 2924. doi.org/10.3390/rs151129244.

4. Ying, Y.; 科娃，M.; Kuffer, M.; Zevenbergen, J. 关于三维属性估值——城市三维建模方法在数字孪生创建中的综述。《国际地理信息科学期刊》，2023，12，2. doi.org/10.3390/ijgi120100025.

5. La Guardia, M.; 科娃，M. 朝向网络上的数字孪生：基于开源结构的异构三维数据融合。《遥感》，2023，15，721. doi.org/10.3390/rs150307216.

6. Khawte, S. S., 科娃，M. N., Gevaert, C. M., Oude Elberink, S., 和 Pedro, A. A.：基于无人机数据的巴西贫民窟数字孪生创建，《国际摄影测量与遥感档案》，XLVIII-4/W4–2022，75–81，doi.org/10.5194/isprs-archives-XLVIII-4-W4-2022-75-2022, 2022