265 lines
10 KiB
Python
265 lines
10 KiB
Python
import os
|
||
import numpy as np
|
||
import pandas as pd
|
||
import py3dtiles
|
||
from py3dtiles import Tileset, BoundingVolumeBox
|
||
from shapely.geometry import Polygon, Point
|
||
import matplotlib.pyplot as plt
|
||
from matplotlib.colors import LinearSegmentedColormap
|
||
import argparse
|
||
from tqdm import tqdm
|
||
import logging
|
||
|
||
# 配置日志
|
||
logging.basicConfig(level=logging.INFO, format='%(asctime)s - %(name)s - %(levelname)s - %(message)s')
|
||
logger = logging.getLogger(__name__)
|
||
|
||
class TilesetProcessor:
|
||
"""3D Tiles数据集处理器，用于加载、分析和比较两个3D Tiles模型"""
|
||
|
||
def __init__(self, tileset_path1, tileset_path2, resolution=1.0):
|
||
"""
|
||
初始化处理器
|
||
|
||
参数:
|
||
tileset_path1: 第一个3D Tiles数据集路径
|
||
tileset_path2: 第二个3D Tiles数据集路径
|
||
resolution: 分析网格的分辨率（米）
|
||
"""
|
||
self.tileset1 = self._load_tileset(tileset_path1)
|
||
self.tileset2 = self._load_tileset(tileset_path2)
|
||
self.resolution = resolution
|
||
self.analysis_area = None
|
||
self.height_difference_grid = None
|
||
self.grid_bounds = None
|
||
|
||
def _load_tileset(self, path):
|
||
"""加载3D Tiles数据集"""
|
||
try:
|
||
logger.info(f"加载3D Tiles数据集: {path}")
|
||
tileset = Tileset.from_file(path)
|
||
logger.info(f"成功加载，包含 {len(tileset.root.children)} 个根瓦片")
|
||
return tileset
|
||
except Exception as e:
|
||
logger.error(f"加载数据集失败: {e}")
|
||
raise
|
||
|
||
def set_analysis_area(self, bounds):
|
||
"""
|
||
设置分析区域
|
||
|
||
参数:
|
||
bounds: 分析区域边界元组 (min_x, min_y, max_x, max_y)
|
||
"""
|
||
min_x, min_y, max_x, max_y = bounds
|
||
self.analysis_area = Polygon([
|
||
(min_x, min_y),
|
||
(max_x, min_y),
|
||
(max_x, max_y),
|
||
(min_x, max_y)
|
||
])
|
||
self.grid_bounds = bounds
|
||
logger.info(f"设置分析区域: {bounds}")
|
||
logger.info(f"分析区域面积: {self.analysis_area.area:.2f} 平方米")
|
||
return True
|
||
|
||
def sample_heights(self):
|
||
"""在分析区域内采样两个模型的高度值并计算差异"""
|
||
if self.analysis_area is None:
|
||
logger.error("请先设置分析区域")
|
||
return False
|
||
|
||
logger.info("开始在分析区域内采样高度值...")
|
||
|
||
# 创建网格
|
||
min_x, min_y, max_x, max_y = self.grid_bounds
|
||
rows = int((max_y - min_y) / self.resolution) + 1
|
||
cols = int((max_x - min_x) / self.resolution) + 1
|
||
|
||
# 初始化高度差异网格
|
||
self.height_difference_grid = np.zeros((rows, cols), dtype=np.float32)
|
||
self.height_difference_grid[:] = np.nan # 初始化为NaN，表示未采样
|
||
|
||
# 对每个网格点进行采样
|
||
total_points = rows * cols
|
||
logger.info(f"创建了 {rows}x{cols}={total_points} 个采样点")
|
||
|
||
with tqdm(total=total_points, desc="采样高度点") as pbar:
|
||
for i in range(rows):
|
||
for j in range(cols):
|
||
# 计算当前点的坐标
|
||
x = min_x + j * self.resolution
|
||
y = min_y + i * self.resolution
|
||
point = Point(x, y)
|
||
|
||
# 检查点是否在分析区域内
|
||
if not self.analysis_area.contains(point):
|
||
pbar.update(1)
|
||
continue
|
||
|
||
# 采样两个模型的高度
|
||
height1 = self._sample_height_at_point(self.tileset1, x, y)
|
||
height2 = self._sample_height_at_point(self.tileset2, x, y)
|
||
|
||
# 计算高度差异
|
||
if height1 is not None and height2 is not None:
|
||
self.height_difference_grid[i, j] = height2 - height1
|
||
|
||
pbar.update(1)
|
||
|
||
# 统计结果
|
||
valid_differences = self.height_difference_grid[~np.isnan(self.height_difference_grid)]
|
||
if len(valid_differences) > 0:
|
||
logger.info(f"高度变化统计:")
|
||
logger.info(f" 平均变化: {np.mean(valid_differences):.2f}m")
|
||
logger.info(f" 最大上升: {np.max(valid_differences):.2f}m")
|
||
logger.info(f" 最大下降: {np.min(valid_differences):.2f}m")
|
||
logger.info(f" 变化标准差: {np.std(valid_differences):.2f}m")
|
||
else:
|
||
logger.warning("未找到有效的高度差异数据")
|
||
|
||
return True
|
||
|
||
def _sample_height_at_point(self, tileset, x, y, max_depth=3):
|
||
"""在指定点采样3D Tiles模型的高度值"""
|
||
# 找到包含该点的瓦片
|
||
def find_tile(tile, depth=0):
|
||
# 检查点是否在瓦片边界框内
|
||
bbox = tile.bounding_volume.box
|
||
min_x_tile = bbox[0] - bbox[3]
|
||
max_x_tile = bbox[0] + bbox[3]
|
||
min_y_tile = bbox[1] - bbox[4]
|
||
max_y_tile = bbox[1] + bbox[4]
|
||
|
||
if not (min_x_tile <= x <= max_x_tile and min_y_tile <= y <= max_y_tile):
|
||
return None
|
||
|
||
# 如果瓦片有内容且深度足够，或者没有子瓦片，就使用这个瓦片
|
||
if (tile.content is not None and depth >= max_depth) or not tile.children:
|
||
return tile
|
||
|
||
# 否则递归查找子瓦片
|
||
for child in tile.children:
|
||
result = find_tile(child, depth + 1)
|
||
if result is not None:
|
||
return result
|
||
|
||
return None
|
||
|
||
# 找到包含该点的最详细瓦片
|
||
tile = find_tile(tileset.root)
|
||
if tile is None or tile.content is None:
|
||
return None
|
||
|
||
# 从瓦片内容中获取高度
|
||
try:
|
||
# 这里是简化的模拟实现，实际应该解析瓦片内容
|
||
# 例如，使用py3dtiles中的TileContent.get_vertices()获取顶点
|
||
# 然后找到最近的顶点或三角形来计算高度
|
||
# 这里为了示例，我们返回瓦片中心的高度加上一个随机偏移
|
||
return tile.bounding_volume.box[2] + np.random.uniform(-0.5, 0.5)
|
||
except Exception as e:
|
||
logger.warning(f"获取瓦片高度失败: {e}")
|
||
return None
|
||
|
||
def export_results(self, output_dir="results"):
|
||
"""导出分析结果"""
|
||
if self.height_difference_grid is None:
|
||
logger.error("请先采样高度值")
|
||
return
|
||
|
||
os.makedirs(output_dir, exist_ok=True)
|
||
|
||
# 导出CSV文件
|
||
csv_path = os.path.join(output_dir, "height_differences.csv")
|
||
logger.info(f"导出CSV文件: {csv_path}")
|
||
|
||
# 创建DataFrame
|
||
min_x, min_y, max_x, max_y = self.grid_bounds
|
||
rows, cols = self.height_difference_grid.shape
|
||
|
||
data = []
|
||
for i in range(rows):
|
||
for j in range(cols):
|
||
if not np.isnan(self.height_difference_grid[i, j]):
|
||
x = min_x + j * self.resolution
|
||
y = min_y + i * self.resolution
|
||
data.append({
|
||
'x': x,
|
||
'y': y,
|
||
'height_difference': self.height_difference_grid[i, j]
|
||
})
|
||
|
||
df = pd.DataFrame(data)
|
||
df.to_csv(csv_path, index=False)
|
||
|
||
# 生成热图
|
||
self._generate_heatmap(output_dir)
|
||
|
||
logger.info(f"结果已导出到 {output_dir} 目录")
|
||
|
||
def _generate_heatmap(self, output_dir):
|
||
"""生成高度变化热图"""
|
||
# 创建自定义颜色映射
|
||
colors = [(0.0, 0.0, 1.0), (1.0, 1.0, 1.0), (1.0, 0.0, 0.0)] # 蓝-白-红
|
||
cmap = LinearSegmentedColormap.from_list('height_diff_cmap', colors, N=256)
|
||
|
||
# 准备数据
|
||
data = self.height_difference_grid.copy()
|
||
valid_mask = ~np.isnan(data)
|
||
|
||
if not np.any(valid_mask):
|
||
logger.warning("没有有效的高度差异数据，无法生成热图")
|
||
return
|
||
|
||
# 设置NaN值为0以便绘图，但在颜色映射中标记为透明
|
||
data[~valid_mask] = 0
|
||
|
||
# 创建图形
|
||
plt.figure(figsize=(12, 10))
|
||
plt.imshow(data, cmap=cmap, origin='lower',
|
||
extent=[self.grid_bounds[0], self.grid_bounds[2],
|
||
self.grid_bounds[1], self.grid_bounds[3]],
|
||
alpha=0.9)
|
||
|
||
# 添加颜色条
|
||
cbar = plt.colorbar()
|
||
cbar.set_label('高度变化 (米)', fontsize=12)
|
||
|
||
# 设置标题和坐标轴
|
||
plt.title('两个3D Tiles模型的高度变化分布', fontsize=16)
|
||
plt.xlabel('X坐标 (米)', fontsize=12)
|
||
plt.ylabel('Y坐标 (米)', fontsize=12)
|
||
|
||
# 保存图形
|
||
heatmap_path = os.path.join(output_dir, "height_difference_heatmap.png")
|
||
plt.savefig(heatmap_path, dpi=300, bbox_inches='tight')
|
||
plt.close()
|
||
|
||
logger.info(f"热图已保存到: {heatmap_path}")
|
||
|
||
def main():
|
||
parser = argparse.ArgumentParser(description='分析两个3D Tiles模型指定区域的高度变化')
|
||
parser.add_argument('--tileset1', required=True, help='第一个3D Tiles数据集路径')
|
||
parser.add_argument('--tileset2', required=True, help='第二个3D Tiles数据集路径')
|
||
parser.add_argument('--bounds', required=True, type=float, nargs=4,
|
||
help='分析区域边界 [min_x, min_y, max_x, max_y]')
|
||
parser.add_argument('--resolution', type=float, default=1.0, help='采样分辨率（米）')
|
||
parser.add_argument('--output', default='results', help='输出目录')
|
||
|
||
args = parser.parse_args()
|
||
|
||
processor = TilesetProcessor(args.tileset1, args.tileset2, args.resolution)
|
||
|
||
# 设置分析区域
|
||
if processor.set_analysis_area(args.bounds):
|
||
if processor.sample_heights():
|
||
processor.export_results(args.output)
|
||
print("分析完成！结果已导出到指定目录。")
|
||
else:
|
||
print("高度采样失败，无法完成分析。")
|
||
else:
|
||
print("设置分析区域失败，无法进行分析。")
|
||
|
||
if __name__ == "__main__":
|
||
main() |