技术研究之土方计算二

标签: 无 分类: 未分类 创建时间:2025-08-06 08:50:09 更新时间:2025-09-26 04:44:39

1.前言

我想着进行土方计算,后来用了 kimi 进行分析,结果给出了方案:
(1)方案一:
OSGB(文件夹) → 点云(PotreeConverter) → DEM(GDAL) → 差值(GDAL) → 3D-Tiles(CesiumLab CLI) → Cesium(Web)。,理论上是没有什么毛病的。我一步步的按照步骤操作,结果出现各种问题,我按照问题去解决,最后还是没有解决,浪费了很长时间,后来告诉我 PotreeConverter 无法处理 osgb 格式啊,真是操蛋。

(2)方案二:
PotreeConverter 读不到点云 → OSGB 不是点云 → 用 ContextCapture 重新导出 LAS,或用 CloudCompare 采样成点云,再继续后续流程。用 CloudCompare,结果折腾了大半天,又发现 CloudCompare 也不支持 osgb 转换,又是浪费时间的一天。“CloudCompare 不支持直接读取 OSGB(倾斜摄影纹理模型);它只能处理 点云或网格,而 OSGB 属于 带纹理的三角网,因此 CLI 返回 exit 1,GUI 也弹 “unhandled file extension osgb”。CloudCompare 无法直接读取 OSGB 文件;前面提到的“File → Open OSGB → Edit → Mesh → Sample Points → Export LAS”前提是你先把 OSGB 转成 CloudCompare 能识别的格式(如 OBJ、PLY)。否则就会弹出 ”

(3)方案三:
完全开源的 OSGB → LAS 方案,使用 OpenSceneGraph 工具。我下载 OpenSceneGraph 工具,又费了很大的劲。

参考文章:
【1】.5分钟了解清楚【osgb】格式的倾斜摄影数据metadata.xml有几种规范
【2】.模方 4.0 用户手册

2.安装依赖

主要就是借助了 osgconv.exe 工具,这个工具也让我找了很久,在windows上的版本OpenSceneGraph-3.6.5-VC2022-64-2025-04.7z,反正不太好下载,下载之后,然后配置到环境变量中就可以了,linux安装如下步骤:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
# 下载cmake,https://cmake.org/download/,上传到服务器
# 设置环境变量

# 安装依赖
yum install gcc
yum install gcc-c++ libstdc++-devel

# 这步其实不用装(但是我执行了,执行了之后,其实安装了很多没用的包)
yum groupinstall "Development Tools"

# 安装依赖
yum list mesa*
yum install -y mesa*
yum install -y freeglut*
yum install -y *glew*


# 执行编译
cd OpenSceneGraph
cmake .
make
make install

# 设置软连接(这一步我没有做)
ln -s ./bin/osgconv /usr/local/bin/osgconv

# 配置环境变量
export PATH="${PATH}:/APP/osgb/openscenegraph-master/bin/"
export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:/APP/osgb/openscenegraph-master/lib/"

# 完成后执行命令使配置生效:
source /etp/profile

参考文章:
【1】.openscenegraph 官方仓库
【2】.Linux环境下OpenSceneGraph的安装和配置
【3】.Linux系统下安装Vcpkg,并使用Vcpkg安装、编译OpenSceneGraph
【4】.Linux下安装CMake的方法 这里可以设置安装
【5】.cmake 3.5.0的版本
【6】.linux 安装gcc和 g++
【7】.linux osg环境配置,5. centos下安装openscenegraph环境 这里给了依赖包,确实帮了大忙了。
【8】.【cmake 项目依赖冲突】CMake进阶:优雅解决目标依赖和安装问题
【10】.“CMAKE_CXX_COMPILER broken” while compiling with CMake 安装cmake错误
【11】.linux 安装gcc和 g++ yum install gcc-c++ libstdc++-devel (正确)
【12】.linux上安装osg_Linux环境下OpenSceneGraph的安装和配置 这里要手工安装很多的依赖,我觉得作用不大。

2.Osgb转Obj

经过不断的折腾,终于把分级的osgb合并成了 obj。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
# merge_osgb_batch.py

import os
import subprocess
import glob
import re


# ================== 配置参数 ==================
IN_ROOT = r"D:\zlc\earthwork\before"
OUT_DIR = r"D:\zlc\earthwork\single2"
OSGCONV = r"C:\Soft\OpenSceneGraph-3.6.5\bin\osgconv.exe"
OSG_PLUGIN_PATH = r"C:\Soft\OpenSceneGraph-3.6.5\bin\osgPlugins-3.6.5"

BATCH_SIZE = 10
TEMP_PREFIX = "temp_merge_"
FINAL_OUTPUT = "merged.obj"
# ==============================================

def find_leaf_osgb_files(root_dir):
"""查找所有包含几何数据的 .osgb 文件(L18~L22)"""
pattern = os.path.join(root_dir, "**", "*.osgb")
all_files = glob.glob(pattern, recursive=True)

valid_files = []
for f in all_files:
fname = os.path.basename(f)

# 排除 root.osgb
if 'root.osgb' in fname.lower():
continue

# 排除顶层结构文件:Tile_+000_+005.osgb
if re.fullmatch(r'tile_\+\d+_\+\d+\.osgb', fname, re.IGNORECASE):
continue

# 排除低层级的 _0 文件:只排除 L17 及以下的 _0
match = re.match(r'tile_\+\d+_\+\d+_l(\d+)_(0+)\.osgb', fname, re.IGNORECASE)
if match:
level = int(match.group(1))
zeros = match.group(2)
# 如果是 L17 及以下,且是 _0,排除
if level <= 17 and zeros == '0':
continue
# 如果是 L18 及以上,即使是 _0 也保留(可能是有效瓦片)
# 或者是 _00002300 这种非全零,也保留
if level >= 18:
pass # 保留
else:
continue

# ✅ 保留所有 L18 及以上层级的文件
level_match = re.search(r'_l(\d+)_', fname, re.IGNORECASE)
if level_match:
level = int(level_match.group(1))
if level >= 18:
valid_files.append(f)
continue

# 可选:如果没匹配到层级,也保留(以防漏掉)
# valid_files.append(f)

valid_files.sort(key=lambda x: x.lower())
return valid_files

def run_osgconv(args, env, cwd):
"""执行 osgconv 命令"""
try:
print(f"执行命令: {' '.join(args)}")
result = subprocess.run(
args,
env=env,
cwd=cwd,
check=True,
capture_output=True,
text=True,
encoding='utf-8'
)
print("stdout:", result.stdout.strip())
return True
except subprocess.CalledProcessError as e:
print(f"❌ 命令执行失败:{e.cmd}")
print(f"错误输出:{e.stderr}")
return False

def main():
print("� 正在检查环境...")

if not os.path.exists(OSGCONV):
print(f"❌ 错误:未找到 osgconv.exe!\n路径:{OSGCONV}")
return

if not os.path.exists(OSG_PLUGIN_PATH):
print(f"❌ 错误:未找到 OSG 插件目录!\n路径:{OSG_PLUGIN_PATH}")
return

# ✅ 创建环境变量副本
env = os.environ.copy()
env['OSG_PLUGIN_PATH'] = OSG_PLUGIN_PATH
osgconv_dir = os.path.dirname(OSGCONV)
env['PATH'] = osgconv_dir + ';' + env.get('PATH', '')

# 创建输出目录
os.makedirs(OUT_DIR, exist_ok=True)
print(f"� 已创建输出目录:{OUT_DIR}")

# 查找所有叶子节点 .osgb 文件
print("� 正在扫描 .osgb 文件...")
leaf_files = find_leaf_osgb_files(IN_ROOT)

if not leaf_files:
print("❌ 错误:未找到任何符合条件的 .osgb 文件!")
return

print(f"✅ 找到 {len(leaf_files)} 个 .osgb 文件,开始分批合并...")

intermediate_objs = []

# ========== 第一轮:分批合并为中间 OBJ 文件 ==========
for i in range(0, len(leaf_files), BATCH_SIZE):
start_idx = i + 1
end_idx = min(i + BATCH_SIZE, len(leaf_files))
temp_obj = f"{TEMP_PREFIX}{start_idx:04d}.obj"

print(f"� 合并第 {start_idx} ~ {end_idx} 个文件 → {temp_obj}")

# ✅ 构建命令:不要手动加引号!
cmd = [
OSGCONV,
"-O", "OutputTextureFiles"
]
cmd.extend(leaf_files[i:end_idx]) # ✅ 让 subprocess 自动处理路径
cmd.append(temp_obj)

success = run_osgconv(cmd, env=env, cwd=OUT_DIR)
if not success:
print(f"❌ 合并失败:{temp_obj}")
return
else:
print(f"✅ 成功生成:{temp_obj}")
intermediate_objs.append(temp_obj)

# ========== 第二轮:合并所有中间 OBJ 文件为最终 OBJ ==========
if len(intermediate_objs) == 1:
final_path = os.path.join(OUT_DIR, FINAL_OUTPUT)
os.replace(os.path.join(OUT_DIR, intermediate_objs[0]), final_path)
print(f"� 单批次完成,已重命名为:{FINAL_OUTPUT}")
else:
print(f"� 正在合并所有中间文件为 {FINAL_OUTPUT}...")
final_cmd = [
OSGCONV,
"-O", "OutputTextureFiles"
]
final_cmd.extend(intermediate_objs) # ✅ 同样不要加引号
final_cmd.append(FINAL_OUTPUT)

success = run_osgconv(final_cmd, env=env, cwd=OUT_DIR)
if not success:
print(f"❌ 最终合并失败!")
return

print(f"� === 成功生成 {FINAL_OUTPUT} === �")

if __name__ == "__main__":
main()

3.Obj转Ply

将obj转为ply,也就是提取点云,有多种方法,一种是使用 CloudCompare ,一种是使用 MeshLab。
(1)CloudCompare
将obj加载进入 CloudCompare 里面,然后找到工具栏上的 “Sample Point On a Mesh”,生成点云数据之后,然后再选择File,保存导出就好了。

参考文章:
【1】.【linux】安装meshlabserver
【2】.使用CloudCompare对obj网格模型转换为pcd/ply点云模型
【3】.使用open3d将obj格式转为pcd格式并保存(模型转点云)
【4】.Ubuntu环境编译OpenSceneGraph

4.Ply转Dem

使用 CloudCompare 进行 Ply 转为 Dem,需要用到工具栏上的 “Convert a cloud to 2D raster” 工具,工具里面还要点击 Update grid 才能导出Dem。

参考文章:
【1】.使用CloudCompare进行点云高程归一化 这里有多种方法进行处理
【2】.CloudCompare——点云转换为数字高程模型(DEM) 这里主要是使用的 Convet 转换的方法。

5.计算土方量

获取到了Dem之后,就可以计算土方量了,两期DEM进行叠加,计算一个差值。

1
pip install rasterio

参考文章:
【1】.

6.可视化

可视化这部分,后来老板提出一个三维可视化的概念:能不能实现在三维中,对挖方或者填方的地方进行选中,然后计算这部分的数值?这让我一度崩溃,我觉得这个需求不合理,但是就是没有可靠的证据能证明这个事情。所以我大量的查找资料,什么变化监测,3dtiles分区,点云聚类之类的,但是最后还是无果。我问过有些人,他们说能做,能选中变化区域,可是我始终没看到过实际的效果和例子。

(1)热力图
其实用二维好实现,就是把DEM做成一个热力图的形式,然后高的地方用红色,低的地方用蓝色等等,但是

(2)三维面片化
另外一种就是针对生成的 diff_result.tif 直接进行面片化,也就是把 dem 进行三维展示。打开 CloudCompare 选中 Dem -> Edit -> Mesh -> Delaunay 2.5D(XY Plane),可以将Dem进行三维化。

参考文章:
【1】.CloudCompare——Delaunay三角剖分【2025最新版】 通过菜单栏的’Edit > Mesh > Delaunay 2.5D (XY plane)’找到该功能。
【2】.GIS应用技巧之利用DEM制作三维立体图 这是用ArcScence进行三维可视化,还用了不同的颜色显示高度。
【3】.Blender—-利用DEM(tif)生成三维模型

7.变化检测

1
2
3
4
# 安装 pdal
conda install -c conda-forge python-pdal laspy
# 安装 numpy
conda install "blas=*=openblas" numpy scipy --force-reinstall

参考文章:
【1】.点云|CloudCompare软件使用总结
【2】.安装 pdal

8.植被区域问题

(1)变化区域
过滤掉0.2米以内的数据

(1)联通区域
将Dem差值进行联通,如果联通区域内的面积较小,那么需要忽略

(2)狭长区域
如果联通区域较为狭长,需要忽略

(3)坡度统计
如果联通区域内的坡度都较为陡峭,有多个峰谷,也需要去掉

(4)空洞统计
如果联通区域内,面积和空洞的比率较大,也需要过滤掉。

参考文章:
【1】.scipy库的label函数|标记图像连通域
【2】.scipy.ndimage.label连通算法

问题

(1)returned non-zero exit status 3221226505
使用 osgb转为 las 的时候,出现了这个问题,其实关键问题是 PotreeConverter 不支持 osgb,让我弄了很久,kimi 的大模型还真是扯淡,这么低级的错误都能犯。

1
2
3
4
5
6
7
File "D:\zlc\earthwork\scripts\osgb2las.py", line 20, in <module>
osgb_folder_to_las(osgb_root, out_las)
File "D:\zlc\earthwork\scripts\osgb2las.py", line 15, in osgb_folder_to_las
subprocess.run(cmd, check=True)
File "C:\Soft\anaconda3\envs\yolo11\Lib\subprocess.py", line 571, in run
raise CalledProcessError(retcode, process.args,
subprocess.CalledProcessError: Command '['D:\\zlc\\earthwork\\scripts\\..\\tools\\PotreeConverter.exe', '-i', 'D:\\zlc\\earthwork\\before', '-o', 'D:\\zlc\\earthwork\\output', '--output-format', 'las', '-l', '1', '--generate-page', '0']' returned non-zero exit status 3221226505.

(2) DLL load failed while importing _base: 找不到指定的程序
无法使用 conda 进行 rasterio 安装,我遇到的问题就是使用 rasterio 的时候,如果用 conda 进行安装的时候,总是会报,找不到指定程序的问题,最后我只能是用 pip 安装 rasterio。

1
2
3
4
5
File "D:\zlc\earthwork\scripts\calculare_earthwork.py", line 3, in <module>
import rasterio
File "C:\Soft\anaconda3\envs\earthwork\Lib\site-packages\rasterio\__init__.py", line 25, in <module>
from rasterio._base import DatasetBase
ImportError: DLL load failed while importing _base: 找不到指定的程序。

【尝试方案】

1
2
3
4
5
# 用pip安装
pip install rasterio

# 安装 gdal
conda install -c conda-forge gdal

参考文章:
【1】.Unable to import python Rasterio package even though it is installed 引入 osgeo

(3)ImportError: DLL load failed while importing pyexpat: 操作系统无法运行 %1。
另外还有一个问题,那就是安装了 matplotlib之后,总是报错:ImportError: DLL load failed while importing pyexpat: 操作系统无法运行 %1。

【解决方案】

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 卸载
pip uninstall matplotlib
# 重新安装
conda install -c conda-forge matplotlib

# 卸载尝试
pip uninstall matplotlib
conda install -c conda-forge matplotlib

# 重新创建环境
conda deactivate
conda env remove -n earthwork
conda create -n earthwork python=3.11
conda activate earthwork
conda config --env --add channels conda-forge
conda config --env --set channel_priority strict
# 一次性安装所有你需要的库,让 conda 来解决所有依赖关系
conda install -c conda-forge rasterio matplotlib numpy pillow

(4)Intel oneMKL FATAL ERROR: Cannot load mkl_intel_thread.2.dll.
在运行的时候,出现了一个 Intel MKL 库 问题。

【解决方案】
解决方案就是重新安装 numpy

1
2
3
4
# 卸载numpy
pip uninstall numpy
# 重新安装
conda install "blas=*=openblas" numpy --force-reinstall
小额赞助
本人提供免费与付费咨询服务,感谢您的支持!赞助请发邮件(ititchuan@gmail.com)通知,方便公布您的善意!
**光 3.01 元
Sun 3.00 元
ititchuan 3.00 元
微信公众号
广告位
诚心邀请广大金主爸爸洽谈合作
每日一省
isNaN 和 Number.isNaN 函数的区别?

1.函数 isNaN 接收参数后,会尝试将这个参数转换为数值,任何不能被转换为数值的的值都会返回 true,因此非数字值传入也会返回 true ,会影响 NaN 的判断。

2.函数 Number.isNaN 会首先判断传入参数是否为数字,如果是数字再继续判断是否为 NaN ,不会进行数据类型的转换,这种方法对于 NaN 的判断更为准确。

每日二省
为什么0.1+0.2 ! == 0.3,如何让其相等?

一个直接的解决方法就是设置一个误差范围,通常称为“机器精度”。对JavaScript来说,这个值通常为2-52,在ES6中,提供了Number.EPSILON属性,而它的值就是2-52,只要判断0.1+0.2-0.3是否小于Number.EPSILON,如果小于,就可以判断为0.1+0.2 ===0.3。

每日三省
== 操作符的强制类型转换规则?

1.首先会判断两者类型是否**相同,**相同的话就比较两者的大小。

2.类型不相同的话,就会进行类型转换。

3.会先判断是否在对比 null 和 undefined,是的话就会返回 true。

4.判断两者类型是否为 string 和 number,是的话就会将字符串转换为 number。

5.判断其中一方是否为 boolean,是的话就会把 boolean 转为 number 再进行判断。

6.判断其中一方是否为 object 且另一方为 string、number 或者 symbol,是的话就会把 object 转为原始类型再进行判断。

每日英语
Happiness is time precipitation, smile is the lonely sad.
幸福是年华的沉淀,微笑是寂寞的悲伤。