从声纳到Web3D展示（TODO） - 指南

1 多波束声纳原始数据

这个数据通常来自扫测船或者其它测绘设备。因为我这边暂时没有原始数据，所以使用的公开的加州海峡群岛的声纳数据。

数据来源：XYZ format of the 2004 Multibeam Bathymetry Data in the Northeastern Channel Islands Region, Southern California [bathyxyz.zip]

​​在这里，XYZ数据​​ 是描述水下地形或物体空间位置的核心数据集，每个数据点包含 ​​经度（X）、纬度（Y）、水深（Z）​​ 三个维度信息。

-119.69927,34.31851,-89.752
-119.69912,34.31851,-89.753
-119.69897,34.31851,-89.753
-119.69882,34.31851,-89.760
-119.69942,34.31836,-90.522
-119.69927,34.31836,-90.155
-119.69912,34.31836,-89.886
-119.69897,34.31836,-89.817
-119.69882,34.31836,-89.837
-119.69868,34.31836,-89.828
-119.69853,34.31836,-89.752
-119.69838,34.31836,-89.719

2 数据建模

这一步是将xyz文件，转换成3D模型。比如DEM（Digital Elevation Model），这里有几种格式，用的最多的是GeoTIFF。

GeoTIFF​​：

是一种 ​​地理栅格数据容器格式​​（文件格式标准），可以存储 ​​任何类型的栅格数据​​（如高程、温度、卫星影像、分类图等）。 ​​内容自由度​​：可包含单波段（如纯高程）或多波段（如RGB+高程+强度）。 ​​扩展性​​：支持嵌入地理坐标、元数据、压缩算法等。

 这里可以使用OSGeo4W软件Spatial without Compromise · QGIS Web Site，这个是免费的。

安装的大小有点恐怖。。。 

使用OGIS打开的xyz图。

使用IDW进行转换，log如下：

QGIS version: 3.40.9-Bratislava
QGIS code revision: 69fa2ab303
Qt version: 5.15.13
Python version: 3.12.11
GDAL version: 3.11.3
GEOS version: 3.13.1-CAPI-1.19.2
PROJ version: Rel. 9.5.0, September 15th, 2024
PDAL version: 2.9.0 (git-version: 42b180)
Algorithm started at: 2025-07-27T03:13:25
Algorithm 'IDW interpolation' starting…
Input parameters:
{ 'DISTANCE_COEFFICIENT' : 2, 'EXTENT' : '-119.753917583,-118.385515583,33.463768291,34.374771010 [EPSG:4326]', 'INTERPOLATION_DATA' : 'file:///E:/test/bathyll.xyz?type=csv&maxFields=10000&detectTypes=yes&xField=-119.69927&yField=34.31851&zField=-89.752&crs=EPSG:4326&spatialIndex=no&subsetIndex=no&watchFile=no::~::0::~::2::~::0', 'OUTPUT' : 'TEMPORARY_OUTPUT', 'PIXEL_SIZE' : 0.1 }
Execution completed in 30.37 seconds
Results:
OUTPUT: C:/Users/Administrator/AppData/Local/Temp/processing_vPksAR/1c7f82bd4efc47c7879b5691a0a4e1ea/OUTPUT.tif
Loading resulting layers
Algorithm 'IDW interpolation' finished

最后创建的tif图。（我在这里转换的图很怪，感觉是哪里配置没对。。。）

也可以使用python来进行转换（中间要安装很多库）

import numpy as np
import pandas as pd
from scipy.interpolate import griddata
import rasterio
from rasterio.transform import from_origin
# 读取 xyz 数据
df = pd.read_csv("bathyll.xyz", header=None, names=["lon", "lat", "value"])
# 创建插值网格
grid_x, grid_y = np.mgrid[
df.lon.min():df.lon.max():500j,
df.lat.min():df.lat.max():500j
]
grid_z = griddata((df.lon, df.lat), df.value, (grid_x, grid_y), method='cubic')
# 创建 GeoTransform
res_lon = (df.lon.max() - df.lon.min()) / 500
res_lat = (df.lat.max() - df.lat.min()) / 500
transform = from_origin(df.lon.min(), df.lat.max(), res_lon, res_lat)
# 保存为 GeoTIFF
with rasterio.open("output.tif", "w",
driver="GTiff",
height=grid_z.shape[0],
width=grid_z.shape[1],
count=1,
dtype=grid_z.dtype,
crs="EPSG:4326",
transform=transform) as dst:
dst.write(grid_z, 1)

这个就是最后创建的tiff文件。

这个看着正常一些，但是和原始的xyz看着差异也不小。 

3 格栅建模的原理

 GeoTIFF 是规则栅格数据（即每个像素有一个固定大小），属于数字高程模型（DEM）。而XYZ 是散点。所以要进行转换。

规则栅格数据​​（Regular Grid Data）是一种 ​​结构化​​ 的空间数据表示形式，由 ​​均匀分布的像素（或网格单元）​​ 组成，每个像素包含一个数值（如高程、温度、浓度等）。其核心特征如下： ​​

几何规则性​​

所有像素大小相同（如 1m×1m），按行列整齐排列。 地理范围由 ​​原点坐标​​（左上角/左下角）和 ​​像素大小​​ 唯一确定。 ​

数据存储​​

以矩阵形式​​存储，每个像素值可通过行列号（row, col）快速定位。

文件格式：GeoTIFF、ASCII Grid、NetCDF、HDF5 等。

特征	点云（Point Cloud）	栅格（Raster/Grid）
数据结构	非规则、稀疏	规则、密集的网格
数据量	仅保存有值的点	保存整个区域的所有格子（无论是否有值）
空间表达	离散	连续（插值）
存储	精简，适合稀疏场景	冗余，适合密集连续场景
应用方向	三维重建、激光雷达	地图、影像、DEM、遥感

适用场景：

场景	更适合点云	更适合栅格
三维扫描/激光雷达	✅	❌
遥感影像、地形图	❌	✅
稀疏海底测深数据	✅	❌
地形分析（坡度、流向）	❌	✅
AI 图像识别训练	❌	✅

点云就是点集合的描述，栅格就是将点数据放到矩阵。两者最大的不同就是存储的数据格式不同，数据一个是不连续一个是连续。使用栅格的最大好处就是方便运算，可以直接使用很多矩阵运算，比如卷积、滤波、插值等，也可以直接用GPU来加速。

将原始点云转化成栅格需要：离散的点云 → 连续的场 → 栅格化

算法有：

IDW（反距离加权）

插值 Kriging（克里金插值）

TIN三角网再栅格化

4 LAS和LAZ

las和laz也是点云文件，但是是二进制的。laz包含压缩。

LAS的文件头（Header）：固定 375 字节

字段	字节长度	说明
文件签名（File Signature）	4 字节	固定为 “LASF”，用于识别 LAS 文件。
文件版本（Version）	2 字节	如 “1.4”（主版本号 1，次版本号 4），不同版本支持的属性字段有差异。
系统标识（System ID）	32 字节	生成数据的系统名称（如传感器型号），不足 32 字节用空字符填充。
点云总数（Point Count）	8 字节	记录文件中点的总数量（64 位整数，支持超大规模点云）。
坐标范围（Min/Max X/Y/Z）	48 字节	6 个双精度浮点数（各 8 字节），记录点云在 X、Y、Z 轴上的最小 / 最大值，用于快速定位空间范围。
坐标比例因子（Scale X/Y/Z）	24 字节	3 个双精度浮点数，用于将点数据中的整数坐标转换为实际物理坐标（见下文 “坐标存储原理”）。
坐标偏移量（Offset X/Y/Z）	24 字节	3 个双精度浮点数，与比例因子配合实现坐标转换。
点数据格式（Point Data Format）	1 字节	定义点数据记录的结构类型（如 0~6，不同类型支持的属性字段不同，例如格式 3 支持 RGB 颜色）。
点数据长度（Point Record Length）	2 字节	每个点数据记录占用的字节数（如格式 3 对应 34 字节 / 点）。

点数据记录（Point Data Records）：核心数据区

字段	字节长度	说明
X/Y/Z 坐标（整数）	4 字节 ×3	实际坐标需通过公式计算：实际坐标 = 整数坐标 × 比例因子 + 偏移量（见下文）。
反射强度（Intensity）	2 字节	0~65535 的整数，反映激光回波强度（材质越亮，强度越高）。
回波信息（Return Flags）	1 字节	二进制位标识（如第 1 位表示是否为首次回波，第 3 位表示是否为末次回波）。
分类标签（Classification）	1 字节	整数代码（如 2 = 地面、6 = 建筑物、9 = 植被），用于点云分类。
RGB 颜色	2 字节 ×3	红、绿、蓝三通道值（0~65535），记录点的颜色信息。

有两种方法转换，一种是代码一种是pdal工具。

PDAL转换

运行命令： 

pdal pipeline xyz2las.json

xyz2las.json

{
"pipeline": [
{
"type": "readers.text",
"filename": "bathyll.xyz",
"header": "X,Y,Z",
"separator": ","
},
{
"type": "writers.las",
"filename": "output.las"
}
]
}

pdal info output.las

tom@PC-20241221RKUQ:~/test$ pdal info output.las
{
"file_size": 131792591,
"filename": "output.las",
"now": "2025-07-27T22:25:29+0800",
"pdal_version": "2.6.2 (git-version: Release)",
"reader": "readers.las",
"stats":
{
"bbox":
{
"native":
{
"bbox":
{
"maxx": -119.05,
"maxy": 34.32,
"maxz": -25.66,
"minx": -119.7,
"miny": 33.89,
"minz": -1244.45
},
"boundary": { "type": "Polygon", "coordinates": [ [ [ -119.7, 33.89, -1244.45 ], [ -119.7, 34.32, -1244.45 ], [ -119.05, 34.32, -25.66 ], [ -119.05, 33.89, -25.66 ], [ -119.7, 33.89, -1244.45 ] ] ] }
}
},
"statistic":
[
{
"average": -119.3893761,
"count": 3876246,
"maximum": -119.05,
"minimum": -119.7,
"name": "X",
"position": 0,
"stddev": 0.1361889993,
"variance": 0.01854744352
},
{
"average": 34.07363823,
"count": 3876246,
"maximum": 34.32,
"minimum": 33.89,
"name": "Y",
"position": 1,
"stddev": 0.1018629814,
"variance": 0.01037606697
},
{
"average": -328.1332479,
"count": 3876246,
"maximum": -25.66,
"minimum": -1244.45,
"name": "Z",
"position": 2,
"stddev": 233.8208149,
"variance": 54672.17348
},
......

python转换

import numpy as np
import laspy
# 1. 读取 xyz 文件
xyz_file = "input.xyz"
xyz = np.loadtxt(xyz_file, usecols=(0, 1, 2))  # 只读取前三列，作为 x,y,z
# 2. 创建 LAS header
header = laspy.LasHeader(point_format=3, version="1.2")  # 可选点格式和版本
header.offsets = np.min(xyz, axis=0)
header.scales = np.array([0.001, 0.001, 0.001])  # 设置精度
# 3. 创建 LAS 文件并写入点
las = laspy.LasData(header)
las.x, las.y, las.z = xyz[:, 0], xyz[:, 1], xyz[:, 2]
# 4. 保存为 .las 文件
las.write("output.las")
print("转换完成：output.las")

LAZ则是压缩过的LAS文件。

5 在Web的显示

在web显示中，需要使用3D tile。整体转换如下：

数据格式	描述	转换路径（→ 3D Tiles）
.tif (DEM)	高程栅格	tif → terrain tiles → Cesium 可加载
.tif (影像)	遥感地图	tif → WMS/WMTS 叠加 → Cesium
.xyz	点云（文本）	.xyz → .las → 3D Tiles
.las/.laz	标准点云格式	直接 → 3D Tiles
.obj/.gltf	3D模型	直接 → 3D Tiles

生成LAS之后，可以直接在Cesium网站转换。

5.1 Cesium转换

Cesium ion

设置了位置之后，就会自动在地图上显示。（我的地图感觉设置还是有点问题。。。）

5.2 命令工具转换

另外一种方式是使用命令转换。 

流程是 LAS->EPT 点云格式（分层瓦片，结构化）->3D Tiles (.pnts + tileset.json)

LAS->EPT

Entwine Point Tile 格式5：这是一种基于八叉树的简单灵活的点云数据存储格式。EPT 数据集的组织包含 JSON 元数据部分以及二进制点数据，FME 等软件可将 ept.json 文件视为一个数据集，每个数据集包含一个点云要素。这种格式常用于处理大规模点云数据，便于数据的存储、管理和传输。

一般使用Entwine转换。 安装Entwine就可以进行转换，但是一般是使用docker转换。（目前还是报错，暂时没找到原因）

tom@PC-20241221RKUQ:~/test$ docker run --rm -v $(pwd):/data connormanning/entwine build   -i output.las   -o output
1/1: output.las
- Skipping output.las
terminate called after throwing an instance of 'nlohmann::detail::out_of_range'
what():  [json.exception.out_of_range.403] key 'bounds' not found

EPT ->3D Tiles

通常使用 entwine2tiles 或 cesium-ion 平台进一步转换

特性	EPT (Entwine Point Tile)	3D Tiles (.pnts)
是否原始点云格式	是	否（已转换）
是否支持 Cesium	是（新版）	是（广泛支持）
加载速度	快（渐进式）	快（预处理优化）
精度	高	可能略有损失
文件大小	小（支持 LAZ）	较大
可嵌入自定义属性	有限	支持更复杂属性
是否需要转换工具	否	需要使用 entwine convert