yooooger

Ai_tottle/_cache_tiles/tileset.json

@@ -0,0 +1 @@
+{"asset":{"gltfUpAxis":"Z","version":"1.0"},"root":{"boundingVolume":{"box":[0.019831704922580684,0.017892807236137287,-0.09323218719441684,805.0581514731019,0.0,0.0,0.0,734.5576074595929,0.0,0.0,0.0,157.5004402762163]},"children":[{"boundingVolume":{"box":[0.019831704922580684,0.017892807236137287,-0.09323218719441684,805.0581514731019,0.0,0.0,0.0,734.5576074595929,0.0,0.0,0.0,157.5004402762163]},"children":[{"boundingVolume":{"box":[0.019831704922580684,0.017892807236137287,-0.09323218719441684,805.0581514731019,0.0,0.0,0.0,734.5576074595929,0.0,0.0,0.0,157.5004402762163]},"content":{"uri":"top/Level_14/Tile_+000_+000.json"},"geometricError":60.0,"refine":"REPLACE"}],"content":{"boundingVolume":{"box":[0.019831704922580684,0.017892807236137287,-0.09323218719441684,805.0581514731019,0.0,0.0,0.0,734.5576074595929,0.0,0.0,0.0,157.5004402762163]},"uri":"top/Level_13/Tile_p0006_p0005.b3dm"},"geometricError":120.0,"refine":"REPLACE"}],"geometricError":1101.1357294254199,"refine":"REPLACE","transform":[-0.9773169501866703,-0.21178191348135814,0.0,0.0,0.10390176688103608,-0.4794788953312529,0.8713807501723461,0.0,-0.1845426826423207,0.8516151772098101,0.49060736666817356,0.0,-1178185.498415972,5437011.306290875,3111245.574712541,1.0]}}

Ai_tottle/tiles.py

@@ -1,418 +0,0 @@
-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()

Ai_tottle/tiles1.py

@@ -0,0 +1,134 @@
+import os
+import json
+import numpy as np
+import requests
+from urllib.parse import urlparse, urljoin
+from pathlib import Path
+from pyproj import Transformer
+from py3dtiles.tileset.content import B3dm
+from py3dtiles.tilers.b3dm.wkb_utils import TriangleSoup
+
+CACHE_DIR = "./_cache_tiles"
+os.makedirs(CACHE_DIR, exist_ok=True)
+
+def download_file(url, cache_dir=CACHE_DIR):
+    os.makedirs(cache_dir, exist_ok=True)
+    filename = os.path.basename(urlparse(url).path)
+    local_path = os.path.join(cache_dir, filename)
+    if not os.path.exists(local_path):
+        print(f"📥 Downloading {url}")
+        r = requests.get(url)
+        r.raise_for_status()
+        with open(local_path, "wb") as f:
+            f.write(r.content)
+    return local_path
+
+def apply_transform(points, transform):
+    n = points.shape[0]
+    homo = np.hstack([points, np.ones((n,1))])
+    transformed = (transform @ homo.T).T[:, :3]
+    return transformed
+
+def extract_positions_from_b3dm(b3dm_bytes, transform=None):
+    soup = TriangleSoup.from_glb_bytes(b3dm_bytes)
+    positions = soup.get_positions()  # Nx3 numpy
+    if transform is not None:
+        positions = apply_transform(positions, transform)
+    return positions
+
+def load_tileset_json(path_or_url):
+    if path_or_url.startswith("http"):
+        local_path = download_file(path_or_url)
+    else:
+        local_path = path_or_url
+    with open(local_path, "r", encoding="utf-8") as f:
+        js = json.load(f)
+    if "geometricError" not in js:
+        print("⚠️ tileset.json缺少geometricError，自动补 1000.0")
+        js["geometricError"] = 1000.0
+    return js, os.path.dirname(local_path), path_or_url if path_or_url.startswith("http") else None
+
+def parse_transform(tile_json):
+    # 返回 4x4 np.array
+    if "transform" in tile_json:
+        return np.array(tile_json["transform"]).reshape(4,4)
+    else:
+        return np.eye(4)
+
+def recursive_load(tile_json, base_dir, base_url, bbox, transformer, transform_accum=None):
+    if transform_accum is None:
+        transform_accum = np.eye(4)
+
+    points_list = []
+    tile_transform = parse_transform(tile_json)
+    transform_total = transform_accum @ tile_transform
+
+    print(f"\n处理 Tile，当前累计 transform 矩阵:\n{transform_total}")
+    print(f"Tile 的 bbox (输入): {bbox}")
+
+    content = tile_json.get("content")
+    if content:
+        uri = content.get("uri")
+        if uri and uri.endswith(".b3dm"):
+            if base_url is not None:
+                b3dm_url = urljoin(base_url + "/", uri)
+                b3dm_path = download_file(b3dm_url)
+            else:
+                b3dm_path = os.path.join(base_dir, uri)
+            if os.path.exists(b3dm_path):
+                with open(b3dm_path, "rb") as f:
+                    b3dm = B3dm.from_array(f.read())
+                    positions = extract_positions_from_b3dm(b3dm.glb_bytes, transform_total)
+                    print(f"加载 {b3dm_path}，原始顶点数量: {len(positions)}")
+
+                    lon, lat, elev = transformer.transform(
+                        positions[:,0], positions[:,1], positions[:,2]
+                    )
+                    print(f"转换后经度范围: {lon.min():.6f} ~ {lon.max():.6f}")
+                    print(f"转换后纬度范围: {lat.min():.6f} ~ {lat.max():.6f}")
+
+                    geo_pts = np.vstack([lon, lat, elev]).T
+                    mask = (
+                        (geo_pts[:,0] >= bbox[0]) & (geo_pts[:,0] <= bbox[2]) &
+                        (geo_pts[:,1] >= bbox[1]) & (geo_pts[:,1] <= bbox[3])
+                    )
+                    in_bbox_count = np.sum(mask)
+                    print(f"bbox内点数量: {in_bbox_count}")
+
+                    points_list.append(geo_pts[mask])
+
+        elif uri and uri.endswith(".json"):
+            if base_url is not None:
+                child_url = urljoin(base_url + "/", uri)
+                child_json, child_base_dir, child_base_url = load_tileset_json(child_url)
+            else:
+                child_path = os.path.join(base_dir, uri)
+                child_json, child_base_dir, child_base_url = load_tileset_json(child_path)
+            pts = recursive_load(child_json, child_base_dir, child_base_url, bbox, transformer, transform_total)
+            if pts.size > 0:
+                points_list.append(pts)
+
+    for child in tile_json.get("children", []):
+        pts = recursive_load(child, base_dir, base_url, bbox, transformer, transform_total)
+        if pts.size > 0:
+            points_list.append(pts)
+
+    if points_list:
+        return np.vstack(points_list)
+    else:
+        return np.empty((0,3))
+
+def load_tileset_all_points(path_or_url, bbox):
+    transformer = Transformer.from_crs("EPSG:4979", "EPSG:4326", always_xy=True)
+    tileset_json, base_dir, base_url = load_tileset_json(path_or_url)
+    points = recursive_load(tileset_json, base_dir, base_url, bbox, transformer)
+    print(f"\n总计提取 {len(points)} 个点在区域内")
+    return points
+
+if __name__ == "__main__":
+    bbox = [116.391, 39.905, 116.405, 39.915]
+    model1_url = "http://8.137.54.85:9000/300bdf2b-a150-406e-be63-d28bd29b409f/dszh/1748398014403562192_OUT/B3DM/tileset.json"
+    model2_url = "http://8.137.54.85:9000/300bdf2b-a150-406e-be63-d28bd29b409f/dszh/1748325943733189898_OUT/B3DM/tileset.json"
+
+    pts1 = load_tileset_all_points(model1_url, bbox)
+    pts2 = load_tileset_all_points(model2_url, bbox)

_cache_tiles/tileset.json

@@ -0,0 +1 @@
+{"asset":{"gltfUpAxis":"Z","version":"1.0"},"root":{"boundingVolume":{"box":[0.019831704922580684,0.017892807236137287,-0.09323218719441684,805.0581514731019,0.0,0.0,0.0,734.5576074595929,0.0,0.0,0.0,157.5004402762163]},"children":[{"boundingVolume":{"box":[0.019831704922580684,0.017892807236137287,-0.09323218719441684,805.0581514731019,0.0,0.0,0.0,734.5576074595929,0.0,0.0,0.0,157.5004402762163]},"children":[{"boundingVolume":{"box":[0.019831704922580684,0.017892807236137287,-0.09323218719441684,805.0581514731019,0.0,0.0,0.0,734.5576074595929,0.0,0.0,0.0,157.5004402762163]},"content":{"uri":"top/Level_14/Tile_+000_+000.json"},"geometricError":60.0,"refine":"REPLACE"}],"content":{"boundingVolume":{"box":[0.019831704922580684,0.017892807236137287,-0.09323218719441684,805.0581514731019,0.0,0.0,0.0,734.5576074595929,0.0,0.0,0.0,157.5004402762163]},"uri":"top/Level_13/Tile_p0006_p0005.b3dm"},"geometricError":120.0,"refine":"REPLACE"}],"geometricError":1101.1357294254199,"refine":"REPLACE","transform":[-0.9773169501866703,-0.21178191348135814,0.0,0.0,0.10390176688103608,-0.4794788953312529,0.8713807501723461,0.0,-0.1845426826423207,0.8516151772098101,0.49060736666817356,0.0,-1178185.498415972,5437011.306290875,3111245.574712541,1.0]}}