如何用Python实现2D坐标转3D空间坐标及点云目标体积计算

整体分为三个核心环节：交互框选2D目标并转地理范围、筛选对应3D点云、体积计算，全程用Python实现，依赖库均为开源常用工具。

前置依赖安装 #

先安装需要的第三方库：

pip install rasterio matplotlib open3d numpy shapely

注意：正射影像与点云必须使用同一套空间坐标系，否则需先完成坐标配准再执行后续操作

步骤1：读取正射影像，交互绘制多边形获取目标范围 #

正射影像一般为带地理变换信息的GeoTIFF格式，我们可以通过rasterio读取地理变换参数，将交互绘制的像素坐标转换为实际地理坐标：

import rasterio
import matplotlib.pyplot as plt
from shapely.geometry import Polygon, Point
import numpy as np

# 读取正射影像
with rasterio.open("你的正射影像路径.tif") as src:
    img = src.read()
    # 调整多波段影像格式为matplotlib可显示的(H,W,3)
    img = np.moveaxis(img, 0, -1)
    transform = src.transform
    crs = src.crs

# 显示影像，交互选多边形点
plt.figure(figsize=(12,8))
plt.imshow(img)
plt.title("左键点击选多边形顶点，右键点击结束选择")
# 选点，无数量上限，右键结束
points = plt.ginput(n=-1, timeout=0, show_clicks=True)
plt.close()

# 转成闭合多边形
polygon_pixel = Polygon(points)
# 像素坐标转地理坐标
polygon_geo = Polygon([
    transform * (x, y) for x, y in polygon_pixel.exterior.coords
])

步骤2：用2D地理范围筛选对应3D点云 #

用Open3D读取点云，筛选所有xy坐标落在上述多边形内的点，得到目标对应的3D点集：

import open3d as o3d

# 读取点云，支持pcd/ply/las等格式
pcd = o3d.io.read_point_cloud("你的点云路径.ply")
# 转成numpy数组，shape为(N,3)，列分别对应x,y,z
points_3d = np.asarray(pcd.points)

# 筛选xy在目标多边形内的点
mask = np.array([polygon_geo.contains(Point(x, y)) for x, y in points_3d[:, :2]])
target_points = points_3d[mask]

# 可选：可视化筛选后的点云确认范围
target_pcd = o3d.geometry.PointCloud()
target_pcd.points = o3d.utility.Vector3dVector(target_points)
o3d.visualization.draw_geometries([target_pcd])

步骤3：基于目标点云计算体积 #

根据目标类型选择对应的体积计算方案：

方案A：地面突出物（堆料、土石方等） #

需要先拟合目标区域的地面高程，计算突出部分的体积：

# 提取目标区域内的地面参考点，这里取高程最低的10%作为地面样本
z = target_points[:, 2]
ground_mask = z < np.percentile(z, 10)
ground_points = target_points[ground_mask]

# 拟合地面平面 z = ax + by + c
A = np.c_[ground_points[:, 0], ground_points[:, 1], np.ones(ground_points.shape[0])]
C, _, _, _ = np.linalg.lstsq(A, ground_points[:, 2], rcond=None)

# 计算每个点与地面的高程差，累加求体积
dists = target_points[:, 2] - (C[0]*target_points[:,0] + C[1]*target_points[:,1] + C[2])
# 过滤低于地面的无效点
dists[dists < 0] = 0
# 计算平均点距作为栅格尺寸
avg_dist = np.mean(np.sqrt(np.sum(np.diff(target_points[:, :2], axis=0)**2, axis=1)))
volume = np.sum(dists) * avg_dist **2
print(f"目标体积为: {volume:.2f} 立方米")

方案B：封闭3D物体 #

直接通过凸包或Alpha Shape重建封闭网格计算体积：

# 凸包法（适合形状规则的封闭物体）
hull, _ = target_pcd.compute_convex_hull()
hull_volume = hull.get_volume()
print(f"凸包体积为: {hull_volume:.2f} 立方米")

# Alpha Shape法（适合形状复杂的物体，alpha参数可调整精度）
alpha = 0.5 # 数值越小越贴合点云，过小将产生碎面
mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_alpha_shape(target_pcd, alpha)
alpha_volume = mesh.get_volume()
print(f"Alpha Shape体积为: {alpha_volume:.2f} 立方米")

常见优化点 #

• 点云密度较高时，可先对点云做下采样，提升筛选速度
• 地面拟合精度要求高的场景，可改用克里金插值、反距离权重插值生成地面DEM后再计算体积
• 多边形框选时可增加边界编辑功能，避免误选背景点

内容的提问来源于stack exchange，提问作者Soham Rangdal