419 lines
14 KiB
Python
419 lines
14 KiB
Python
import os
|
||
import requests
|
||
import numpy as np
|
||
from py3dtiles.tileset.content import B3dm
|
||
from pyproj import Proj, Transformer
|
||
from shapely.geometry import Polygon, Point
|
||
from shapely.ops import unary_union
|
||
import math
|
||
from urllib.parse import urljoin
|
||
from pygltflib import GLTF2
|
||
|
||
# 目标区域经纬度坐标（转换为多边形）
|
||
region_coords = [
|
||
[102.22321717600258, 29.384100779345513],
|
||
[102.22612442019208, 29.384506810595088],
|
||
[102.22638603372953, 29.382061071072794],
|
||
[102.22311237980807, 29.38186133280733],
|
||
[102.22321717600258, 29.384100779345513] # 闭合多边形
|
||
]
|
||
|
||
# 创建多边形对象
|
||
region_polygon = Polygon(region_coords)
|
||
|
||
# 两个3D Tiles模型的URL
|
||
tileset_urls = [
|
||
"http://8.137.54.85:9000/300bdf2b-a150-406e-be63-d28bd29b409f/dszh/1748398014403562192_OUT/B3DM/tileset.json",
|
||
"http://8.137.54.85:9000/300bdf2b-a150-406e-be63-d28bd29b409f/dszh/1748325943733189898_OUT/B3DM/tileset.json"
|
||
]
|
||
|
||
# 坐标系转换
|
||
wgs84 = Proj(init='epsg:4326') # WGS84经纬度
|
||
web_mercator = Proj(init='epsg:3857') # Web墨卡托投影
|
||
|
||
def adjust_z_in_transform(tileset_path, output_path=None, delta_z=0):
|
||
import json
|
||
import numpy as np
|
||
|
||
if not os.path.exists(tileset_path):
|
||
print(f"❌ tileset.json 文件不存在: {tileset_path}")
|
||
return
|
||
|
||
with open(tileset_path, 'r', encoding='utf-8') as f:
|
||
data = json.load(f)
|
||
|
||
root = data.get('root', {})
|
||
|
||
# 插入默认 transform
|
||
if 'transform' not in root:
|
||
print("⚠️ 未找到 transform 字段，使用单位矩阵")
|
||
root['transform'] = [
|
||
1, 0, 0, 0,
|
||
0, 1, 0, 0,
|
||
0, 0, 1, 0,
|
||
0, 0, 0, 1
|
||
]
|
||
|
||
transform = np.array(root['transform']).reshape(4, 4)
|
||
print(f"原始 Z 平移: {transform[3, 2]}")
|
||
transform[3, 2] += delta_z
|
||
print(f"修正后 Z 平移: {transform[3, 2]}")
|
||
|
||
root['transform'] = transform.flatten().tolist()
|
||
data['root'] = root
|
||
|
||
if output_path is None:
|
||
output_path = tileset_path.replace(".json", f"_adjusted_z{int(delta_z)}.json")
|
||
|
||
with open(output_path, 'w', encoding='utf-8') as f:
|
||
json.dump(data, f, indent=2)
|
||
|
||
print(f"✅ 高度调整完成，输出文件: {output_path}")
|
||
|
||
def download_tileset(tileset_url):
|
||
"""下载tileset.json数据"""
|
||
try:
|
||
response = requests.get(tileset_url)
|
||
response.raise_for_status()
|
||
return response.json()
|
||
except Exception as e:
|
||
print(f"下载tileset失败: {e}")
|
||
return None
|
||
|
||
def extract_vertices(gltf_bytes):
|
||
try:
|
||
with open("temp.glb", "wb") as f:
|
||
f.write(gltf_bytes)
|
||
|
||
gltf = GLTF2().load("temp.glb")
|
||
|
||
for mesh in gltf.meshes:
|
||
for primitive in mesh.primitives:
|
||
if not hasattr(primitive.attributes, "POSITION"):
|
||
continue
|
||
|
||
accessor_idx = primitive.attributes.POSITION
|
||
accessor = gltf.accessors[accessor_idx]
|
||
buffer_view = gltf.bufferViews[accessor.bufferView]
|
||
buffer = gltf.buffers[buffer_view.buffer]
|
||
|
||
byte_offset = (buffer_view.byteOffset or 0) + (accessor.byteOffset or 0)
|
||
byte_length = accessor.count * 3 * 4 # 3 floats per vertex
|
||
|
||
data_bytes = gltf.binary_blob()[byte_offset: byte_offset + byte_length]
|
||
vertices = np.frombuffer(data_bytes, dtype=np.float32).reshape((accessor.count, 3))
|
||
|
||
return vertices
|
||
|
||
except Exception as e:
|
||
print(f"提取顶点数据失败: {e}")
|
||
|
||
return np.array([])
|
||
|
||
def find_closest_vertex(vertices, lon, lat):
|
||
"""找到离目标点最近的顶点"""
|
||
if not vertices:
|
||
return None
|
||
|
||
# 计算距离并找到最近的顶点
|
||
min_distance = float('inf')
|
||
closest_vertex = None
|
||
|
||
for vertex in vertices:
|
||
v_lon, v_lat, v_z = vertex
|
||
# 计算经纬度距离（简化为平面距离）
|
||
distance = math.hypot(v_lon - lon, v_lat - lat)
|
||
if distance < min_distance:
|
||
min_distance = distance
|
||
closest_vertex = vertex
|
||
|
||
return closest_vertex
|
||
|
||
def compare_heights(heights1, heights2, tolerance=0.5):
|
||
"""比较两个高度数据集，找出差异"""
|
||
# 找到所有点的并集
|
||
all_points = set(heights1.keys()).union(set(heights2.keys()))
|
||
differences = []
|
||
|
||
for point in all_points:
|
||
h1 = heights1.get(point, None)
|
||
h2 = heights2.get(point, None)
|
||
|
||
# 检查是否有一个模型在该点没有数据
|
||
if h1 is None or h2 is None:
|
||
differences.append({
|
||
'point': point,
|
||
'height1': h1,
|
||
'height2': h2,
|
||
'difference': None,
|
||
'type': 'missing_data'
|
||
})
|
||
else:
|
||
# 检查高度差异是否超过容忍度
|
||
diff = abs(h1 - h2)
|
||
if diff > tolerance:
|
||
differences.append({
|
||
'point': point,
|
||
'height1': h1,
|
||
'height2': h2,
|
||
'difference': diff,
|
||
'type': 'height_difference'
|
||
})
|
||
|
||
return differences
|
||
|
||
def get_b3dm_from_tile_json(json_url):
|
||
try:
|
||
response = requests.get(json_url)
|
||
response.raise_for_status()
|
||
data = response.json()
|
||
|
||
# 递归查找 b3dm uri
|
||
def find_b3dm_uri(node):
|
||
if 'content' in node and 'uri' in node['content']:
|
||
uri = node['content']['uri']
|
||
if uri.endswith('.b3dm'):
|
||
return uri
|
||
if 'children' in node:
|
||
for child in node['children']:
|
||
result = find_b3dm_uri(child)
|
||
if result:
|
||
return result
|
||
return None
|
||
|
||
root = data.get('root', {})
|
||
b3dm_uri = find_b3dm_uri(root)
|
||
if not b3dm_uri:
|
||
print(f"{json_url} 中找不到 content.uri")
|
||
return None
|
||
|
||
base_url = os.path.dirname(json_url)
|
||
full_b3dm_url = urljoin(base_url + '/', b3dm_uri)
|
||
return full_b3dm_url
|
||
|
||
except Exception as e:
|
||
print(f"解析 JSON {json_url} 时出错: {e}")
|
||
return None
|
||
|
||
|
||
def get_heights_in_region(tileset_url, sample_density=10):
|
||
"""获取区域内的高度数据"""
|
||
tileset_json = download_tileset(tileset_url)
|
||
if not tileset_json:
|
||
return {}
|
||
|
||
tiles_in_region = get_tiles_in_region(tileset_json, tileset_url)
|
||
if not tiles_in_region:
|
||
print(f"在{tileset_url}中未找到区域内的瓦片")
|
||
return {}
|
||
|
||
min_lon, min_lat = min(p[0] for p in region_coords), min(p[1] for p in region_coords)
|
||
max_lon, max_lat = max(p[0] for p in region_coords), max(p[1] for p in region_coords)
|
||
lon_steps = np.linspace(min_lon, max_lon, sample_density)
|
||
lat_steps = np.linspace(min_lat, max_lat, sample_density)
|
||
|
||
heights = {}
|
||
|
||
for tile_info in tiles_in_region:
|
||
try:
|
||
response = requests.get(tile_info['url'])
|
||
response.raise_for_status()
|
||
b3dm_data = response.content
|
||
|
||
# ✅ 尝试解析为 b3dm
|
||
try:
|
||
gltf_bytes = parse_b3dm(b3dm_data)
|
||
except Exception:
|
||
# 可能 tile_info['url'] 是 JSON，不是真 b3dm
|
||
print(f"尝试从 {tile_info['url']} 获取真实 b3dm 地址...")
|
||
actual_b3dm_url = get_b3dm_from_tile_json(tile_info['url'])
|
||
if not actual_b3dm_url:
|
||
print(f"跳过：无法从 {tile_info['url']} 获取有效 b3dm")
|
||
continue
|
||
response = requests.get(actual_b3dm_url)
|
||
response.raise_for_status()
|
||
b3dm_data = response.content
|
||
gltf_bytes = parse_b3dm(b3dm_data)
|
||
|
||
# ✅ 模拟解析 glb
|
||
vertices = extract_vertices(gltf_bytes)
|
||
if not vertices.size:
|
||
continue
|
||
|
||
# ✅ 应用变换
|
||
transformed_vertices = [transform_point(v, tile_info['transform']) for v in vertices]
|
||
transformer = Transformer.from_crs("EPSG:3857", "EPSG:4326", always_xy=True)
|
||
wgs84_vertices = []
|
||
for x, y, z in transformed_vertices:
|
||
lon, lat = transformer.transform(x, y)
|
||
wgs84_vertices.append((lon, lat, z))
|
||
|
||
for lon in lon_steps:
|
||
for lat in lat_steps:
|
||
if point_in_region(lon, lat):
|
||
closest_vertex = find_closest_vertex(wgs84_vertices, lon, lat)
|
||
if closest_vertex:
|
||
key = (round(lon, 6), round(lat, 6))
|
||
heights[key] = closest_vertex[2]
|
||
|
||
except Exception as e:
|
||
print(f"处理瓦片 {tile_info['url']} 时出错: {e}")
|
||
continue
|
||
|
||
return heights
|
||
|
||
def get_tiles_in_region(tileset_json, tileset_base_url):
|
||
"""获取区域内的所有瓦片"""
|
||
tiles_in_region = []
|
||
|
||
# 去除 tileset.json 得到根路径
|
||
tileset_root_url = tileset_base_url.rsplit('/', 1)[0]
|
||
|
||
def recursive_search(tile, parent_transform=None):
|
||
tile_transform = tile.get('transform', [1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1])
|
||
combined_transform = multiply_matrices(parent_transform, tile_transform) if parent_transform else tile_transform
|
||
|
||
if 'boundingVolume' in tile and is_bounding_volume_intersects_region(tile['boundingVolume']):
|
||
if 'content' in tile and 'uri' in tile['content']:
|
||
# 修复URL拼接
|
||
tile_url = urljoin(tileset_root_url + '/', tile['content']['uri'])
|
||
tiles_in_region.append({
|
||
'url': tile_url,
|
||
'transform': combined_transform
|
||
})
|
||
|
||
if 'children' in tile:
|
||
for child in tile['children']:
|
||
recursive_search(child, combined_transform)
|
||
|
||
if 'root' in tileset_json:
|
||
recursive_search(tileset_json['root'])
|
||
|
||
return tiles_in_region
|
||
|
||
def is_bounding_volume_intersects_region(bounding_volume):
|
||
"""检查边界体是否与区域相交"""
|
||
# 简化实现，实际需要根据不同边界体类型实现
|
||
if 'region' in bounding_volume:
|
||
# region格式: [west, south, east, north, minHeight, maxHeight]
|
||
region = bounding_volume['region']
|
||
bv_polygon = Polygon([
|
||
[region[0], region[1]],
|
||
[region[2], region[1]],
|
||
[region[2], region[3]],
|
||
[region[0], region[3]],
|
||
[region[0], region[1]]
|
||
])
|
||
return region_polygon.intersects(bv_polygon)
|
||
elif 'box' in bounding_volume:
|
||
# 对于box类型，需要转换到经纬度后再判断
|
||
# 这里简化处理，返回True让更细致的检查在后续进行
|
||
return True
|
||
elif 'sphere' in bounding_volume:
|
||
# 对于sphere类型，简化处理
|
||
return True
|
||
return False
|
||
|
||
def multiply_matrices(a, b):
|
||
"""计算两个4x4矩阵的乘积"""
|
||
result = [0.0] * 16
|
||
for i in range(4):
|
||
for j in range(4):
|
||
result[i*4 + j] = a[i*4 + 0] * b[0*4 + j] + \
|
||
a[i*4 + 1] * b[1*4 + j] + \
|
||
a[i*4 + 2] * b[2*4 + j] + \
|
||
a[i*4 + 3] * b[3*4 + j]
|
||
return result
|
||
|
||
def parse_b3dm(b3dm_data: bytes):
|
||
"""
|
||
解析 b3dm 文件，返回 glb 二进制数据
|
||
"""
|
||
import struct
|
||
|
||
if b3dm_data[:4] != b'b3dm':
|
||
raise ValueError("不是有效的 b3dm 文件")
|
||
|
||
# 读取 header（28 字节）
|
||
header = struct.unpack('<4sIIIIII', b3dm_data[:28])
|
||
_, version, byte_length, ft_json_len, ft_bin_len, bt_json_len, bt_bin_len = header
|
||
|
||
glb_start = 28 + ft_json_len + ft_bin_len + bt_json_len + bt_bin_len
|
||
glb_bytes = b3dm_data[glb_start:]
|
||
|
||
return glb_bytes
|
||
|
||
def point_in_region(lon, lat):
|
||
"""判断点是否在目标区域内"""
|
||
return region_polygon.contains(Point(lon, lat))
|
||
|
||
def transform_point(point, matrix):
|
||
"""应用变换矩阵到点"""
|
||
x, y, z = point
|
||
x_out = x * matrix[0] + y * matrix[4] + z * matrix[8] + matrix[12]
|
||
y_out = x * matrix[1] + y * matrix[5] + z * matrix[9] + matrix[13]
|
||
z_out = x * matrix[2] + y * matrix[6] + z * matrix[10] + matrix[14]
|
||
return (x_out, y_out, z_out)
|
||
|
||
def main():
|
||
sample_density = 20
|
||
|
||
print("正在从第一个3D Tiles模型提取区域高度数据...")
|
||
heights1 = get_heights_in_region(tileset_urls[0], sample_density)
|
||
|
||
print("正在从第二个3D Tiles模型提取区域高度数据...")
|
||
heights2 = get_heights_in_region(tileset_urls[1], sample_density)
|
||
|
||
if not heights1 or not heights2:
|
||
print("无法获取足够的高度数据进行比较")
|
||
return
|
||
|
||
# 计算平均高度
|
||
avg1 = np.mean(list(heights1.values()))
|
||
avg2 = np.mean(list(heights2.values()))
|
||
|
||
print(f"\n模型1 平均高度: {avg1:.2f} 米")
|
||
print(f"模型2 平均高度: {avg2:.2f} 米")
|
||
|
||
delta = avg1 - avg2
|
||
print(f"高度差: {delta:.2f} 米")
|
||
|
||
# 🔧 自动统一高度
|
||
if abs(delta) > 0.5:
|
||
print("\n⚙️ 正在统一高度基准（修改模型2的 transform）...")
|
||
# tileset_urls[1] 是远程 URL，下载后调整
|
||
try:
|
||
ts2_url = tileset_urls[1]
|
||
response = requests.get(ts2_url)
|
||
response.raise_for_status()
|
||
with open("tileset_model2.json", "w", encoding="utf-8") as f:
|
||
f.write(response.text)
|
||
|
||
adjust_z_in_transform("tileset_model2.json", "tileset_model2_adjusted.json", delta_z=delta)
|
||
except Exception as e:
|
||
print(f"❌ 调整高度失败: {e}")
|
||
else:
|
||
print("\n✅ 高度差异在容忍范围内，无需调整")
|
||
|
||
# 🔍 差异分析
|
||
print("\n正在分析详细差异点...")
|
||
differences = compare_heights(heights1, heights2, 0.5)
|
||
|
||
if differences:
|
||
print(f"共发现 {len(differences)} 处显著高度差异:")
|
||
for i, diff in enumerate(differences[:10], 1): # 仅显示前10条
|
||
lon, lat = diff['point']
|
||
print(f"\n位置 {i}: 经度 {lon}, 纬度 {lat}")
|
||
print(f"模型1高度: {diff['height1']:.2f}米")
|
||
print(f"模型2高度: {diff['height2']:.2f}米")
|
||
if diff['difference'] is not None:
|
||
print(f"差异: {diff['difference']:.2f}米")
|
||
else:
|
||
print("差异: 一个模型在该位置没有数据")
|
||
else:
|
||
print("两个模型在指定区域高度基本一致 ✅")
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|