From 5c865a4418c67bf89aca7be44ca7291479766a2b Mon Sep 17 00:00:00 2001 From: liyubo Date: Wed, 14 Jan 2026 11:37:35 +0800 Subject: [PATCH] =?UTF-8?q?=E5=9D=A1=E5=BA=A6=E5=9D=A1=E5=90=91=E5=9C=9F?= =?UTF-8?q?=E6=96=B9=E9=87=8F=E8=AE=A1=E7=AE=97?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- b3dm/b3dm_api.py | 60 ++ b3dm/data_3dtiles_manager.py | 566 +++++++++++++++++ b3dm/data_3dtiles_to_dem.py | 757 +++++++++++++++++++++++ b3dm/earthwork_api.py | 377 +++++++++++ b3dm/earthwork_calculator_3d_tiles.py | 526 ++++++++++++++++ b3dm/earthwork_calculator_point_cloud.py | 691 +++++++++++++++++++++ b3dm/glb_with_draco.py | 469 ++++++++++++++ b3dm/terrain3d_analyzer_color.py | 380 ++++++++++++ b3dm/terrain_api.py | 342 ++++++++++ b3dm/terrain_calculator.py | 361 +++++++++++ b3dm/tileset_data_source.py | 167 +++++ b3dm/tileset_to_ply.py | 395 ++++++++++++ b3dm/volume_calculator.py | 684 ++++++++++++++++++++ 13 files changed, 5775 insertions(+) create mode 100644 b3dm/b3dm_api.py create mode 100644 b3dm/data_3dtiles_manager.py create mode 100644 b3dm/data_3dtiles_to_dem.py create mode 100644 b3dm/earthwork_api.py create mode 100644 b3dm/earthwork_calculator_3d_tiles.py create mode 100644 b3dm/earthwork_calculator_point_cloud.py create mode 100644 b3dm/glb_with_draco.py create mode 100644 b3dm/terrain3d_analyzer_color.py create mode 100644 b3dm/terrain_api.py create mode 100644 b3dm/terrain_calculator.py create mode 100644 b3dm/tileset_data_source.py create mode 100644 b3dm/tileset_to_ply.py create mode 100644 b3dm/volume_calculator.py diff --git a/b3dm/b3dm_api.py b/b3dm/b3dm_api.py new file mode 100644 index 0000000..4147fcf --- /dev/null +++ b/b3dm/b3dm_api.py @@ -0,0 +1,60 @@ +from sanic import Sanic, Request, json +from sanic_cors import CORS +import logging +import time +from earthwork_api import earthwork_bp +from terrain_api import terrain_bp + +# 配置日志 +logging.basicConfig(level=logging.INFO) +logger = logging.getLogger(__name__) + +# 创建Sanic应用 +app = Sanic("TerrainAnalysisAPI") +# 显式注册蓝图 +app.blueprint(earthwork_bp) +app.blueprint(terrain_bp) + +CORS(app, automatic_options=True) + +# 中间件：请求计时 +@app.middleware("request") +async def add_start_time(request: Request): + request.ctx.start_time = time.time() + +@app.middleware("response") +async def add_response_time(request: Request, response): + if hasattr(request.ctx, "start_time"): + process_time = (time.time() - request.ctx.start_time) * 1000 + response.headers["X-Process-Time"] = f"{process_time:.2f}ms" + +@app.get("/api/v1/health") +async def health_check(request: Request): + """健康检查""" + return json({ + "status": "healthy", + "timestamp": time.time(), + "service": "terrain-analysis-api", + "version": "1.0.0" + }) + +# 错误处理 +@app.exception(Exception) +async def handle_exception(request: Request, exception): + """全局异常处理""" + logger.error(f"未处理的异常: {exception}") + return json({ + "error": "服务器内部错误", + "message": str(exception) if app.debug else "请稍后重试", + "success": False + }, status=500) + +if __name__ == "__main__": + # 启动服务器 + app.run( + host="0.0.0.0", + port=8000, + debug=True, # 生产环境设为False + access_log=True, + auto_reload=True + ) \ No newline at end of file diff --git a/b3dm/data_3dtiles_manager.py b/b3dm/data_3dtiles_manager.py new file mode 100644 index 0000000..18ae940 --- /dev/null +++ b/b3dm/data_3dtiles_manager.py @@ -0,0 +1,566 @@ +from minio import Minio +from minio.error import S3Error +import json +import os +import numpy as np +from urllib.parse import urlparse +import hashlib +import time +import re +import pickle +from datetime import datetime + +class MinIO3DTilesManager: + def __init__(self, endpoint_url, access_key, secret_key, secure=False, + mapping_file="minio_path_mapping.pkl"): + """ + 初始化MinIO客户端 + + Args: + endpoint_url: MinIO服务地址 (如: 222.212.85.86:9001) + access_key: 访问密钥 + secret_key: 秘密密钥 + secure: 是否使用HTTPS + mapping_file: 路径映射文件名 + """ + if endpoint_url.startswith('http://'): + endpoint_url = endpoint_url.replace('http://', '') + elif endpoint_url.startswith('https://'): + endpoint_url = endpoint_url.replace('https://', '') + secure = True + + self.endpoint_url = endpoint_url + self.access_key = access_key + self.secret_key = secret_key + + self.minio_client = Minio( + endpoint_url, + access_key=access_key, + secret_key=secret_key, + secure=secure + ) + + # 获取脚本所在目录 + self.script_dir = os.path.dirname(os.path.abspath(__file__)) + + # 映射文件路径 + self.mapping_file = os.path.join(self.script_dir, mapping_file) + + # 加载现有的路径映射 + self.path_mapping = self.load_path_mapping() + + def load_path_mapping(self): + """加载路径映射数据""" + if os.path.exists(self.mapping_file): + try: + with open(self.mapping_file, 'rb') as f: + mapping = pickle.load(f) + return mapping + except Exception as e: + return {} + else: + return {} + + def save_path_mapping(self): + """保存路径映射数据""" + try: + with open(self.mapping_file, 'wb') as f: + pickle.dump(self.path_mapping, f) + return True + except Exception as e: + return False + + def get_cache_key(self, tileset_url, save_dir=None): + """生成缓存键""" + # 基于URL和保存目录生成缓存键 + cache_data = f"{tileset_url}|{save_dir}" + return hashlib.md5(cache_data.encode()).hexdigest() + + def get_cached_tileset_info(self, tileset_url, save_dir=None): + """获取缓存的tileset信息""" + cache_key = self.get_cache_key(tileset_url, save_dir) + + # 检查缓存映射中是否有这个tileset + for file_id, info in self.path_mapping.items(): + if info.get('cache_key') == cache_key and info.get('is_tileset_root'): + # 检查入口文件是否存在 + local_path = info.get('local_path') + if local_path and os.path.exists(local_path): + return local_path + return None + + def update_tileset_cache(self, tileset_url, save_dir, local_path): + """更新tileset缓存信息""" + cache_key = self.get_cache_key(tileset_url, save_dir) + + # 将tileset根文件标记为缓存 + entry_bucket, entry_path = self.parse_minio_url(tileset_url) + file_id = f"{entry_bucket}/{entry_path}" + + if file_id in self.path_mapping: + self.path_mapping[file_id]['cache_key'] = cache_key + self.path_mapping[file_id]['is_tileset_root'] = True + self.path_mapping[file_id]['tileset_url'] = tileset_url + self.path_mapping[file_id]['save_dir'] = save_dir + self.path_mapping[file_id]['cache_time'] = datetime.now().isoformat() + + def download_full_tileset(self, tileset_url, save_dir=None, region_filter=None, use_cache=True): + """ + 下载完整的3D Tiles数据集，支持缓存功能 + + Args: + tileset_url: MinIO上的tileset.json URL + save_dir: 本地保存目录 + region_filter: 区域过滤器 + use_cache: 是否使用缓存 + + Returns: + tuple: (success, result) + - success: True/False + - result: 如果success=True且use_cache=True，返回本地路径；否则返回True/False + """ + if save_dir is None: + save_dir = os.path.join(self.script_dir, "data_3dtiles") + + # 清理保存目录名称 + save_dir = self.clean_file_path(save_dir) + + # 检查缓存：只需检查入口文件是否存在 + if use_cache: + cached_path = self.get_cached_tileset_info(tileset_url, save_dir) + if cached_path: + # 入口文件存在，默认缓存完备 + return True, cached_path + + # 解析URL + entry_bucket, entry_path = self.parse_minio_url(tileset_url) + if not entry_bucket or not entry_path: + return False, "无法解析URL" + + entry_dir = os.path.dirname(entry_path) + + # 创建保存目录 + os.makedirs(save_dir, exist_ok=True) + + visited = set() + + # 下载入口文件 + entry_local_path = self.get_local_path( + entry_bucket, entry_path, + entry_bucket, entry_dir, + save_dir + ) + + success, result = self.download_file(entry_bucket, entry_path, entry_local_path) + if not success: + return False, f"入口文件下载失败: {result}" + + entry_id = f"{entry_bucket}/{entry_path}" + visited.add(entry_id) + + # 加载tileset数据 + tileset_data = self.load_json_from_minio(entry_bucket, entry_path) + if not tileset_data or "root" not in tileset_data: + return False, "无效的tileset.json文件" + + # 遍历下载所有文件 + self.traverse_and_download_tileset( + tileset_data["root"], + entry_bucket, + entry_dir, + entry_bucket, + entry_dir, + save_dir, + region_filter, + None, + visited + ) + + # 更新缓存信息 + self.update_tileset_cache(tileset_url, save_dir, entry_local_path) + + # 保存路径映射 + self.save_path_mapping() + + if use_cache: + return True, entry_local_path + else: + return True, True + + def get_tileset_local_path(self, tileset_url, save_dir=None): + """ + 获取已缓存的tileset本地路径 + + Args: + tileset_url: tileset的URL + save_dir: 保存目录 + + Returns: + str: 本地路径，如果未缓存则返回None + """ + if save_dir is None: + save_dir = os.path.join(self.script_dir, "data_3dtiles") + + return self.get_cached_tileset_info(tileset_url, save_dir) + + def clear_tileset_cache(self, tileset_url=None, save_dir=None): + """ + 清除tileset缓存 + + Args: + tileset_url: 指定要清除的tileset URL，如果为None则清除所有 + save_dir: 保存目录 + + Returns: + bool: 成功/失败 + """ + try: + if tileset_url: + # 清除指定tileset的缓存 + cache_key = self.get_cache_key(tileset_url, save_dir) + + # 找出所有相关的缓存条目 + to_remove = [] + for file_id, info in self.path_mapping.items(): + if info.get('cache_key') == cache_key: + to_remove.append(file_id) + + # 删除这些条目 + for file_id in to_remove: + del self.path_mapping[file_id] + + print(f"已清除tileset缓存: {tileset_url}") + else: + # 清除所有缓存 + self.path_mapping = {} + if os.path.exists(self.mapping_file): + os.remove(self.mapping_file) + print("已清除所有缓存") + + return True + except Exception as e: + return False + + # 以下是原有的辅助方法 + def clean_filename(self, filename): + """清理文件名中的特殊字符""" + if not filename: + return "" + cleaned = re.sub(r'[<>:"/\\|?*\x00-\x1F]', '_', filename) + cleaned = re.sub(r'_+', '_', cleaned) + cleaned = cleaned.strip(' _') + return cleaned + + def parse_minio_url(self, url): + """解析MinIO URL""" + if url.startswith('http://') or url.startswith('https://'): + parsed = urlparse(url) + path = parsed.path.lstrip('/') + parts = path.split('/', 1) + if len(parts) == 2: + bucket, key = parts + else: + bucket = parts[0] + key = "" + return bucket, key + else: + parts = url.split('/', 1) + if len(parts) == 2: + bucket, key = parts + else: + bucket = parts[0] + key = "" + return bucket, key + + def download_file(self, bucket_name, object_name, file_path): + """从MinIO下载文件""" + try: + os.makedirs(os.path.dirname(file_path), exist_ok=True) + + # 清理文件名 + clean_file_path = self.clean_file_path(file_path) + + # 检查是否已下载 + file_id = f"{bucket_name}/{object_name}" + if file_id in self.path_mapping: + mapped_path = self.path_mapping[file_id]['local_path'] + if os.path.exists(mapped_path): + return True, mapped_path + + # 下载文件 + self.minio_client.fget_object( + bucket_name, + object_name, + clean_file_path + ) + + # 更新路径映射 + self.path_mapping[file_id] = { + 'local_path': clean_file_path, + 'bucket': bucket_name, + 'object': object_name, + 'download_time': datetime.now().isoformat(), + 'size': os.path.getsize(clean_file_path) + } + + return True, clean_file_path + + except S3Error as e: + return False, str(e) + except Exception as e: + return False, str(e) + + def clean_file_path(self, file_path): + """清理文件路径中的所有特殊字符""" + dir_name = os.path.dirname(file_path) + file_name = os.path.basename(file_path) + + if dir_name: + dir_parts = dir_name.split(os.sep) + cleaned_parts = [] + for part in dir_parts: + cleaned_part = self.clean_filename(part) + if cleaned_part: + cleaned_parts.append(cleaned_part) + cleaned_dir = os.sep.join(cleaned_parts) + else: + cleaned_dir = "" + + cleaned_file = self.clean_filename(file_name) + + if cleaned_dir: + cleaned_path = os.path.join(cleaned_dir, cleaned_file) + else: + cleaned_path = cleaned_file + + return cleaned_path + + def load_json_from_minio(self, bucket_name, object_name): + """从MinIO加载JSON文件""" + try: + self.minio_client.stat_object(bucket_name, object_name) + + response = self.minio_client.get_object(bucket_name, object_name) + content = response.read().decode('utf-8') + response.close() + response.release_conn() + + return json.loads(content) + + except S3Error as e: + return None + except Exception as e: + return None + + def get_local_path(self, bucket_name, object_name, base_bucket, base_object, save_dir): + """生成保持目录结构的本地路径""" + clean_bucket = self.clean_filename(bucket_name) + bucket_dir = clean_bucket + + if bucket_name == base_bucket and base_object: + base_dir = os.path.dirname(base_object) + + if base_dir: + if object_name.startswith(base_dir): + relative_path = object_name[len(base_dir):].lstrip('/\\') + else: + relative_path = object_name + else: + relative_path = object_name + else: + relative_path = object_name + + if relative_path: + path_parts = relative_path.split('/') + cleaned_parts = [] + for part in path_parts: + cleaned_part = self.clean_filename(part) + if cleaned_part: + cleaned_parts.append(cleaned_part) + + if cleaned_parts: + cleaned_relative = '/'.join(cleaned_parts) + local_path = os.path.join(save_dir, bucket_dir, cleaned_relative) + else: + local_path = os.path.join(save_dir, bucket_dir) + else: + local_path = os.path.join(save_dir, bucket_dir) + + return os.path.normpath(local_path) + + def traverse_and_download_tileset(self, tile_obj, current_bucket, current_dir, + base_bucket, base_dir, save_dir, + region_filter=None, parent_transform=None, + visited=None): + """递归遍历并下载3D Tiles文件""" + if visited is None: + visited = set() + + current_transform = parent_transform + if "transform" in tile_obj: + tile_mat = tile_obj["transform"] + if current_transform is None: + current_transform = tile_mat + else: + mat1 = np.array(current_transform).reshape(4, 4) + mat2 = np.array(tile_mat).reshape(4, 4) + combined_mat = np.dot(mat1, mat2).flatten().tolist() + current_transform = combined_mat + + skip_current_tile = False + if region_filter and "boundingVolume" in tile_obj: + if not region_filter.check_tile_bounding_volume(tile_obj["boundingVolume"]): + skip_current_tile = True + + if not skip_current_tile and "content" in tile_obj and "uri" in tile_obj["content"]: + tile_uri = tile_obj["content"]["uri"] + + file_bucket = current_bucket + file_path = "" + + if tile_uri.startswith('http://') or tile_uri.startswith('https://'): + parsed_bucket, parsed_path = self.parse_minio_url(tile_uri) + if parsed_bucket: + file_bucket = parsed_bucket + file_path = parsed_path + else: + if current_dir: + file_path = os.path.join(current_dir, tile_uri).replace('\\', '/') + else: + file_path = tile_uri + + file_path = file_path.lstrip('/') + + file_id = f"{file_bucket}/{file_path}" + + if file_id not in visited: + print(f"下载文件:{file_id}") + visited.add(file_id) + + local_path = self.get_local_path( + file_bucket, file_path, + base_bucket, base_dir, + save_dir + ) + + self.download_file(file_bucket, file_path, local_path) + + if file_path.lower().endswith('.json'): + sub_tileset = self.load_json_from_minio(file_bucket, file_path) + if sub_tileset and "root" in sub_tileset: + sub_dir = os.path.dirname(file_path) if file_path else "" + self.traverse_and_download_tileset( + sub_tileset["root"], + file_bucket, + sub_dir, + base_bucket, + base_dir, + save_dir, + region_filter, + current_transform, + visited + ) + + if "children" in tile_obj: + for child_tile in tile_obj["children"]: + self.traverse_and_download_tileset( + child_tile, + current_bucket, + current_dir, + base_bucket, + base_dir, + save_dir, + region_filter, + current_transform, + visited + ) + + def upload_file(self, bucket_name, object_name, file_path): + """上传文件到MinIO""" + try: + if not os.path.exists(file_path): + return False, f"文件不存在: {file_path}" + + file_size = os.path.getsize(file_path) + self.minio_client.fput_object(bucket_name, object_name, file_path) + + return True, f"{bucket_name}/{object_name}" + + except S3Error as e: + return False, f"MinIO上传错误: {e}" + except Exception as e: + return False, f"上传失败: {str(e)}" + + def upload_directory(self, bucket_name, local_dir, remote_prefix=""): + """上传目录到MinIO""" + if not os.path.exists(local_dir): + return [], [f"目录不存在: {local_dir}"] + + uploaded_files = [] + failed_files = [] + + for root, dirs, files in os.walk(local_dir): + for file in files: + local_path = os.path.join(root, file) + rel_path = os.path.relpath(local_path, local_dir) + if remote_prefix: + remote_path = os.path.join(remote_prefix, rel_path).replace('\\', '/') + else: + remote_path = rel_path.replace('\\', '/') + + success, message = self.upload_file(bucket_name, remote_path, local_path) + if success: + uploaded_files.append(remote_path) + else: + failed_files.append((remote_path, message)) + + return uploaded_files, failed_files + + def check_and_create_bucket(self, bucket_name): + """检查并创建bucket""" + try: + if not self.minio_client.bucket_exists(bucket_name): + self.minio_client.make_bucket(bucket_name) + return True, f"创建bucket: {bucket_name}" + return True, f"bucket已存在: {bucket_name}" + except S3Error as e: + return False, f"创建bucket失败: {e}" + + +# 使用示例 +if __name__ == "__main__": + # 配置参数 + ENDPOINT_URL = "222.212.85.86:9000" + ACCESS_KEY = "WuRenJi" + SECRET_KEY = "WRJ@2024" + + # 初始化管理器 + manager = MinIO3DTilesManager( + endpoint_url=ENDPOINT_URL, + access_key=ACCESS_KEY, + secret_key=SECRET_KEY, + secure=False + ) + + # 使用缓存下载tileset + tileset_url = "http://222.212.85.86:9000/300bdf2b-a150-406e-be63-d28bd29b409f/model/石棉0908/terra_b3dms/tileset.json" + + # 第一次下载（会下载到本地） + success, result = manager.download_full_tileset(tileset_url, use_cache=True) + if success: + print(f"下载成功，本地路径: {result}") + + # 第二次下载相同URL（直接从缓存返回） + success, result = manager.download_full_tileset(tileset_url, use_cache=True) + if success: + print(f"从缓存获取，本地路径: {result}") + + # 强制重新下载（忽略缓存） + success, result = manager.download_full_tileset(tileset_url, use_cache=False) + if success: + print("强制重新下载成功") + + # 获取缓存的本地路径 + local_path = manager.get_tileset_local_path(tileset_url) + if local_path: + print(f"缓存的本地路径: {local_path}") \ No newline at end of file diff --git a/b3dm/data_3dtiles_to_dem.py b/b3dm/data_3dtiles_to_dem.py new file mode 100644 index 0000000..f4b387d --- /dev/null +++ b/b3dm/data_3dtiles_to_dem.py @@ -0,0 +1,757 @@ +# data_3dtiles_to_dem52 +import os +import json +import numpy as np +import pandas as pd +import pyproj +import struct +from osgeo import gdal, osr +import uuid +from glb_with_draco import DracoGLBParser + +# 解决GDAL中文路径/警告问题（生产必加） +gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES") +gdal.SetConfigOption("CPL_ZIP_ENCODING", "UTF-8") +gdal.PushErrorHandler('CPLQuietErrorHandler') + +class RegionFilter: + """完整的区域过滤器，支持所有格式""" + + def __init__(self, region_coords=None, enable_tile_filter=True, debug=False): + """ + 初始化区域过滤器 + :param debug: 是否输出调试信息 + """ + self.region_coords = region_coords + self.enable_tile_filter = enable_tile_filter + self.debug = debug + + if region_coords: + # 提取区域边界 + lons = [coord[0] for coord in region_coords] + lats = [coord[1] for coord in region_coords] + self.min_lon = min(lons) + self.max_lon = max(lons) + self.min_lat = min(lats) + self.max_lat = max(lats) + + # 扩展边界（避免边缘误差） + self.expand_factor = 0.1 + lon_expand = (self.max_lon - self.min_lon) * self.expand_factor + lat_expand = (self.max_lat - self.min_lat) * self.expand_factor + + self.filter_min_lon = self.min_lon - lon_expand + self.filter_max_lon = self.max_lon + lon_expand + self.filter_min_lat = self.min_lat - lat_expand + self.filter_max_lat = self.max_lat + lat_expand + + if self.debug: + print(f"[DEBUG] 区域过滤器初始化:") + print(f" 目标区域: Lon[{self.min_lon:.6f}, {self.max_lon:.6f}], " + f"Lat[{self.min_lat:.6f}, {self.max_lat:.6f}]") + print(f" 过滤区域: Lon[{self.filter_min_lon:.6f}, {self.filter_max_lon:.6f}], " + f"Lat[{self.filter_min_lat:.6f}, {self.filter_max_lat:.6f}]") + else: + self.min_lon = self.min_lat = -180 + self.max_lon = self.max_lat = 180 + self.filter_min_lon = self.filter_min_lat = -180 + self.filter_max_lon = self.filter_max_lat = 180 + if self.debug: + print("[DEBUG] 区域过滤器: 未指定区域，将处理所有数据") + + def check_tile_bounding_volume(self, bounding_volume): + """检查瓦片的包围体是否与指定区域相交""" + if not self.region_coords or not self.enable_tile_filter: + return True + + try: + if 'region' in bounding_volume: + return self._check_region(bounding_volume['region']) + + elif 'box' in bounding_volume: + box = bounding_volume['box'] + if self.debug: + print(f"[DEBUG] 检查box包围体，长度={len(box)}") + + if len(box) == 12: + result = self._check_box_12(box) + elif len(box) == 15: + result = self._check_box_15(box) + else: + if self.debug: + print(f"[DEBUG] 异常box长度 {len(box)}，默认通过") + return True + + if self.debug and not result: + print(f"[DEBUG] Box被过滤") + return result + + elif 'sphere' in bounding_volume: + return self._check_sphere(bounding_volume['sphere']) + + return True + + except Exception as e: + if self.debug: + print(f"[DEBUG] 包围体检查出错: {e}") + return True + + def _check_region(self, region): + """检查region格式 [west, south, east, north, minHeight, maxHeight]""" + if len(region) != 6: + return True + + west, south, east, north, min_h, max_h = region + + # 检查是否完全在过滤区域外 + if (east < self.filter_min_lon or west > self.filter_max_lon or + north < self.filter_min_lat or south > self.filter_max_lat): + if self.debug: + print(f"[DEBUG] Region过滤: [{west:.3f},{south:.3f},{east:.3f},{north:.3f}]") + return False + + return True + + def _check_box_12(self, box): + """检查12值box格式""" + # 提取参数 + cx, cy, cz = box[0], box[1], box[2] + halfX = box[3] + halfY = box[7] + halfZ = box[11] + + if self.debug: + print(f"[DEBUG] 12值box: center=({cx:.1f},{cy:.1f},{cz:.1f}), " + f"halfs=({halfX},{halfY},{halfZ})") + + try: + # 获取转换器 + transformer = self._get_transformer() + + # 转换中心点 + center_lon, center_lat, _ = transformer.transform( + cx, cy, cz, radians=False + ) + + # 计算最大偏移（简化方法，避免计算所有角点） + max_half = max(halfX, halfY, halfZ) + + # 转换为经纬度偏移 + earth_radius = 6378137.0 + lon_offset = np.degrees(max_half / (earth_radius * np.cos(np.radians(center_lat)))) + lat_offset = np.degrees(max_half / earth_radius) + + # 计算box范围 + box_min_lon = center_lon - lon_offset + box_max_lon = center_lon + lon_offset + box_min_lat = center_lat - lat_offset + box_max_lat = center_lat + lat_offset + + if self.debug: + print(f"[DEBUG] 中心: ({center_lon:.6f}, {center_lat:.6f})") + print(f"[DEBUG] 范围: Lon[{box_min_lon:.6f}, {box_max_lon:.6f}], " + f"Lat[{box_min_lat:.6f}, {box_max_lat:.6f}]") + + # 检查相交 + if (box_max_lon < self.filter_min_lon or box_min_lon > self.filter_max_lon or + box_max_lat < self.filter_min_lat or box_min_lat > self.filter_max_lat): + return False + + return True + + except Exception as e: + if self.debug: + print(f"[DEBUG] Box12检查失败: {e}") + return True # 失败时默认通过 + + def _check_box_15(self, box): + """检查标准15值box格式""" + if len(box) < 15: + return True + + cx, cy, cz = box[0], box[1], box[2] + halfX, halfY, halfZ = box[12], box[13], box[14] + + # 简化处理：只检查中心点 + transformer = self._get_transformer() + center_lon, center_lat, _ = transformer.transform(cx, cy, cz, radians=False) + + # 计算最大偏移 + max_half = max(halfX, halfY, halfZ) + earth_radius = 6378137.0 + lon_offset = np.degrees(max_half / (earth_radius * np.cos(np.radians(center_lat)))) + lat_offset = np.degrees(max_half / earth_radius) + + box_min_lon = center_lon - lon_offset + box_max_lon = center_lon + lon_offset + box_min_lat = center_lat - lat_offset + box_max_lat = center_lat + lat_offset + + if (box_max_lon < self.filter_min_lon or box_min_lon > self.filter_max_lon or + box_max_lat < self.filter_min_lat or box_min_lat > self.filter_max_lat): + return False + + return True + + def _check_sphere(self, sphere): + """检查sphere格式 [centerX, centerY, centerZ, radius]""" + if len(sphere) < 4: + return True + + cx, cy, cz, radius = sphere[0], sphere[1], sphere[2], sphere[3] + + transformer = self._get_transformer() + center_lon, center_lat, _ = transformer.transform(cx, cy, cz, radians=False) + + # 计算半径对应的经纬度偏移 + earth_radius = 6378137.0 + lon_offset = np.degrees(radius / (earth_radius * np.cos(np.radians(center_lat)))) + lat_offset = np.degrees(radius / earth_radius) + + sphere_min_lon = center_lon - lon_offset + sphere_max_lon = center_lon + lon_offset + sphere_min_lat = center_lat - lat_offset + sphere_max_lat = center_lat + lat_offset + + if (sphere_max_lon < self.filter_min_lon or sphere_min_lon > self.filter_max_lon or + sphere_max_lat < self.filter_min_lat or sphere_min_lat > self.filter_max_lat): + return False + + return True + + def _get_transformer(self): + """获取或创建坐标转换器""" + if not hasattr(self, '_transformer'): + ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84') + lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84') + self._transformer = pyproj.Transformer.from_proj(ecef, lla) + return self._transformer + + def filter_points(self, points): + """ + 过滤点集，只保留区域内的点 + :param points: 点列表或numpy数组，每行 [lon, lat, height] + :return: 过滤后的点列表 + """ + if not self.region_coords or len(points) == 0: + return points + + # 转换为numpy数组处理 + if isinstance(points, list): + points_array = np.array(points) + else: + points_array = points + + if len(points_array) == 0 or points_array.shape[1] < 2: + return points_array.tolist() if isinstance(points, list) else points_array + + # 检查每个点是否在区域内 + in_region_mask = ( + (points_array[:, 0] >= self.min_lon) & + (points_array[:, 0] <= self.max_lon) & + (points_array[:, 1] >= self.min_lat) & + (points_array[:, 1] <= self.max_lat) + ) + + filtered_points = points_array[in_region_mask] + + print(f"区域过滤: {len(points_array)} 个点 -> {len(filtered_points)} 个点 " + f"(过滤掉 {len(points_array) - len(filtered_points)} 个点)") + + return filtered_points.tolist() if isinstance(points, list) else filtered_points + +# ========== 核心工具函数：矩阵变换 ========== +def apply_transform_matrix(vertices, transform_matrix): + """ + 将模型的局部相对顶点坐标，通过transform矩阵转换为绝对ECEF坐标 + :param vertices: 原始局部顶点 (n,3) numpy数组 + :param transform_matrix: 瓦片的transform矩阵 一维列表/数组，长度16（4x4矩阵） + :return: 绝对ECEF坐标 (n,3) numpy数组 + """ + if transform_matrix is None or len(transform_matrix) != 16: + return vertices + + # reshape为4x4，然后转置（因为glTF是列主序） + mat = np.array(transform_matrix).reshape(4, 4).astype(np.float64).T + + # 顶点齐次坐标化 (n,3) -> (n,4) 最后一列补1 + ones = np.ones((vertices.shape[0], 1), dtype=np.float64) + vertices_hom = np.hstack([vertices, ones]) + + # 矩阵乘法：顶点坐标 * 变换矩阵 = 绝对坐标 + vertices_ecef_hom = np.dot(vertices_hom, mat.T) + + # 还原为三维坐标 (n,4) -> (n,3) + vertices_ecef = vertices_ecef_hom[:, :3] + + return vertices_ecef + +def parse_b3dm_to_points(b3dm_path, region_filter=None, transform_matrix=None): + """ + 解析B3DM文件,提取顶点的经纬度+高程 + 【关键修改】区域过滤移到读取顶点后进行 + """ + # 获取脚本所在目录 + script_dir = os.path.dirname(os.path.abspath(__file__)) + temp_dir = os.path.join(script_dir, "temp_glb") + + # 创建临时目录（如果不存在） + os.makedirs(temp_dir, exist_ok=True) + + # 1. 读取B3DM二进制文件 + with open(b3dm_path, "rb") as f: + b3dm_data = f.read() + + # 跳过头部 + header = struct.unpack('<4sIIIIII', b3dm_data[:28]) + ft_json_len, ft_bin_len, bt_json_len, bt_bin_len = header[3:7] + + offset = 28 + offset += ft_json_len # 跳过Feature Table JSON + offset += ft_bin_len # 跳过Feature Table Binary + offset += bt_json_len # 跳过Batch Table JSON + offset += bt_bin_len # 跳过Batch Table Binary + + # 提取glb数据 + glb_data = b3dm_data[offset:] + if len(glb_data) < 12: + return [] + + # 2. 使用脚本目录下的临时文件 + temp_file_path = None + try: + # 生成唯一临时文件名 + temp_filename = f"temp_{uuid.uuid4().hex[:8]}.glb" + temp_file_path = os.path.join(temp_dir, temp_filename) + + # 将GLB数据写入临时文件 + with open(temp_file_path, "wb") as tmp_glb: + tmp_glb.write(glb_data) + + # 使用DracoGLBParser解析 + parser = DracoGLBParser(temp_file_path) + # 解析 GLB 结构 + parser.parse_glb_structure() + # 分析结构 + parser.analyze_structure() + # 解码 Draco 网格 + mesh = parser.decode_draco_meshes() + + except Exception as e: + print(f"读取GLB数据失败 {b3dm_path}: {e}") + return [] + finally: + # 清理临时文件 + if temp_file_path and os.path.exists(temp_file_path): + try: + os.unlink(temp_file_path) + except Exception as e: + pass + + if mesh is None: + print(f"无法加载模型: {b3dm_path}") + return [] + + # 获取顶点数据 + vertices = parser.get_all_vertices() + + if vertices.size == 0 or len(vertices.shape) < 2 or vertices.shape[1] != 3: + print(f"顶点数据格式无效: {b3dm_path}") + return [] + + # 3. 应用transform矩阵，局部坐标 → 绝对ECEF坐标 + vertices = apply_transform_matrix(vertices, transform_matrix) + + # 4. ECEF坐标转WGS84经纬度+高程 + try: + ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84') + lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84') + transformer = pyproj.Transformer.from_proj(ecef, lla, always_xy=True) + + lons, lats, heights = transformer.transform( + vertices[:, 0], vertices[:, 1], vertices[:, 2], radians=False + ) + + # 组合成点集 + points = np.column_stack([lons, lats, heights]) + + # 5. 基本数据清洗（过滤异常值） + # 过滤nan/inf + valid_mask = np.isfinite(points).all(axis=1) + points = points[valid_mask] + + if len(points) == 0: + print(f"B3DM文件 {os.path.basename(b3dm_path)} 转换后无有效点") + return [] + + # 过滤经纬度超限值 + geo_mask = ( + (points[:, 0] >= -180) & (points[:, 0] <= 180) & + (points[:, 1] >= -90) & (points[:, 1] <= 90) + ) + points = points[geo_mask] + + if len(points) == 0: + print(f"B3DM文件 {os.path.basename(b3dm_path)} 经纬度超限") + return [] + + print(f"从 {os.path.basename(b3dm_path)} 提取到 {len(points)} 个原始顶点") + + # 6. 【关键修改】应用区域过滤器（如果提供） + # 在读取顶点后进行区域过滤，而不是在瓦片级别过滤 + if region_filter: + points = region_filter.filter_points(points) + if len(points) == 0: + print(f"B3DM文件 {os.path.basename(b3dm_path)} 的所有点都在区域外，跳过") + return [] + + print(f"最终保留 {len(points)} 个在区域内的顶点") + return points.tolist() + + except Exception as e: + print(f"坐标转换失败 {b3dm_path}: {e}") + return [] + +def traverse_nested_tiles(tile_obj, base_dir, b3dm_paths, tile_transforms, region_filter=None, parent_transform=None): + """ + 深度递归遍历瓦片,自动识别「子JSON」和「B3DM」 + 【修改】移除瓦片级别的区域过滤，完全依赖顶点级别的过滤 + """ + # 1. 计算当前瓦片的最终变换矩阵 + current_transform = parent_transform + if "transform" in tile_obj: + tile_mat = tile_obj["transform"] + if current_transform is None: + current_transform = tile_mat + else: + # 合并变换矩阵 + mat1 = np.array(current_transform).reshape(4, 4) + mat2 = np.array(tile_mat).reshape(4, 4) + combined_mat = np.dot(mat1, mat2).flatten().tolist() + current_transform = combined_mat + + # 2. 【修改】不再检查瓦片包围体，直接处理内容 + # 这样可以避免因粗略的包围体判断而漏掉部分在区域内的顶点 + + # 3. 处理当前瓦片的内容 + if "content" in tile_obj and "uri" in tile_obj["content"]: + tile_uri = tile_obj["content"]["uri"] + tile_abs_path = os.path.join(base_dir, tile_uri) + + if tile_uri.lower().endswith(".json"): + # 情况1:uri是子JSON文件 → 递归解析这个子JSON + if os.path.exists(tile_abs_path): + print(f"解析嵌套子JSON文件: {tile_abs_path}") + with open(tile_abs_path, "r", encoding="utf-8") as f: + sub_tileset = json.load(f) + sub_base_dir = os.path.dirname(tile_abs_path) + traverse_nested_tiles(sub_tileset["root"], sub_base_dir, b3dm_paths, tile_transforms, region_filter, current_transform) + else: + print(f"嵌套子JSON文件不存在,跳过: {tile_abs_path}") + + elif tile_uri.lower().endswith(".b3dm"): + # 情况2:uri是B3DM文件 → 收集路径+对应transform矩阵 + if os.path.exists(tile_abs_path): + b3dm_paths.append(tile_abs_path) + tile_transforms.append(current_transform) + print(f"收集到B3DM文件: {tile_abs_path}") + else: + print(f"B3DM文件不存在,跳过: {tile_abs_path}") + + # 4. 递归遍历当前tile的子节点children + if "children" in tile_obj and len(tile_obj["children"]) > 0: + for child_tile in tile_obj["children"]: + traverse_nested_tiles(child_tile, base_dir, b3dm_paths, tile_transforms, region_filter, current_transform) + +def parse_tileset(tileset_path, region_coords=None): + """重构主解析函数,支持区域过滤+矩阵变换""" + if not os.path.exists(tileset_path): + raise FileNotFoundError(f"根tileset.json文件不存在: {tileset_path}") + + # 初始化区域过滤器（将在顶点级别使用） + region_filter = RegionFilter(region_coords) + + # 读取根tileset.json + with open(tileset_path, "r", encoding="utf-8") as f: + tileset_json = json.load(f) + + root_dir = os.path.dirname(tileset_path) + b3dm_paths = [] + tile_transforms = [] + + print(f"开始遍历tileset结构...") + # 调用深度递归函数 + traverse_nested_tiles(tileset_json["root"], root_dir, b3dm_paths, tile_transforms, region_filter, None) + + print(f"\n遍历完成，共发现 {len(b3dm_paths)} 个B3DM文件:") + for i, b3dm_path in enumerate(b3dm_paths): + print(f" {i+1}. {os.path.basename(b3dm_path)}") + + # 批量解析所有B3DM文件,合并点云 + all_points = [] + if len(b3dm_paths) == 0: + print("未提取到任何有效的B3DM文件！") + return all_points + + print(f"\n===== 开始解析B3DM文件 =====") + for i, (b3dm_path, transform_mat) in enumerate(zip(b3dm_paths, tile_transforms), 1): + print(f"解析文件 {i}/{len(b3dm_paths)}: {os.path.basename(b3dm_path)}") + points = parse_b3dm_to_points(b3dm_path, region_filter, transform_mat) + if points: + all_points.extend(points) + print(f" 提取到 {len(points)} 个点") + else: + print(f" 未提取到有效点") + + # 点云去重+优化 + if all_points: + all_points = np.array(all_points) + original_count = len(all_points) + all_points = np.unique(all_points.round(decimals=6), axis=0) + print(f"\n最终提取点云数量: {len(all_points)} 个 (已去重, 去除了 {original_count - len(all_points)} 个重复点)") + + # ========== 新增：输出整个地图文件的经纬度高程范围 ========== + print("\n" + "=" * 60) + print("地图文件总体范围统计:") + print("-" * 60) + + # 计算经纬度范围 + min_lon = np.min(all_points[:, 0]) + max_lon = np.max(all_points[:, 0]) + min_lat = np.min(all_points[:, 1]) + max_lat = np.max(all_points[:, 1]) + min_height = np.min(all_points[:, 2]) + max_height = np.max(all_points[:, 2]) + avg_height = np.mean(all_points[:, 2]) + std_height = np.std(all_points[:, 2]) + + # 计算中心点 + center_lon = (min_lon + max_lon) / 2 + center_lat = (min_lat + max_lat) / 2 + center_height = (min_height + max_height) / 2 + + print(f"经度范围: {min_lon:.6f}° ~ {max_lon:.6f}° (跨度: {max_lon - min_lon:.6f}°)") + print(f"纬度范围: {min_lat:.6f}° ~ {max_lat:.6f}° (跨度: {max_lat - min_lat:.6f}°)") + print(f"高程范围: {min_height:.2f}m ~ {max_height:.2f}m (总高差: {max_height - min_height:.2f}m)") + print(f"平均高程: {avg_height:.2f}m (±{std_height:.2f}m)") + print(f"\n中心点坐标:") + print(f" 位置: ({center_lon:.6f}°, {center_lat:.6f}°)") + print(f" 高程: {center_height:.2f}m") + + # 输出边界坐标（用于复制使用） + print(f"\n边界坐标:") + print(f" 西北角: ({min_lon:.6f}, {max_lat:.6f})") + print(f" 东北角: ({max_lon:.6f}, {max_lat:.6f})") + print(f" 西南角: ({min_lon:.6f}, {min_lat:.6f})") + print(f" 东南角: ({max_lon:.6f}, {min_lat:.6f})") + + # 如果有区域过滤，显示过滤效果 + if region_coords: + region_min_lon = min(coord[0] for coord in region_coords) + region_max_lon = max(coord[0] for coord in region_coords) + region_min_lat = min(coord[1] for coord in region_coords) + region_max_lat = max(coord[1] for coord in region_coords) + + print(f"\n区域过滤效果:") + print(f" 原始地图范围: Lon[{region_min_lon:.6f}, {region_max_lon:.6f}], Lat[{region_min_lat:.6f}, {region_max_lat:.6f}]") + print(f" 提取数据范围: Lon[{min_lon:.6f}, {max_lon:.6f}], Lat[{min_lat:.6f}, {max_lat:.6f}]") + + # 计算覆盖率 + region_width = region_max_lon - region_min_lon + region_height = region_max_lat - region_min_lat + extracted_width = max_lon - min_lon + extracted_height = max_lat - min_lat + + width_coverage = extracted_width / region_width * 100 if region_width > 0 else 0 + height_coverage = extracted_height / region_height * 100 if region_height > 0 else 0 + + print(f" 经度方向覆盖率: {width_coverage:.1f}%") + print(f" 纬度方向覆盖率: {height_coverage:.1f}%") + + print("=" * 60) + + return all_points + +def points_to_dem(points, output_dem_path, pixel_size=0.0001): + """将离散点云插值为DEM（GeoTIFF格式）- 使用Scipy插值优化版本""" + if len(points) == 0: + raise ValueError("无有效点云数据,无法生成DEM") + + # 转换为numpy数组 + points_array = np.array(points) + lons = points_array[:, 0] + lats = points_array[:, 1] + heights = points_array[:, 2] + + min_lon, max_lon = lons.min(), lons.max() + min_lat, max_lat = lats.min(), lats.max() + + print(f"DEM范围: Lon[{min_lon:.6f}, {max_lon:.6f}], Lat[{min_lat:.6f}, {max_lat:.6f}]") + print(f"点云数量: {len(points)} 个") + print(f"高程范围: {heights.min():.2f} ~ {heights.max():.2f} 米") + + # 计算网格尺寸 + width = int((max_lon - min_lon) / pixel_size) + 1 + height = int((max_lat - min_lat) / pixel_size) + 1 + + # 限制网格大小，避免过大 + max_grid_size = 5000 # 最大网格尺寸 + if width > max_grid_size or height > max_grid_size: + print(f"警告: 网格尺寸过大 ({width}x{height})，自动调整像素大小...") + # 重新计算像素大小 + larger_dim = max(width, height) + pixel_size = pixel_size * (larger_dim / max_grid_size) + width = int((max_lon - min_lon) / pixel_size) + 1 + height = int((max_lat - min_lat) / pixel_size) + 1 + print(f"调整后像素大小: {pixel_size:.6f}°") + + print(f"DEM网格: {width}x{height} (像素大小: {pixel_size:.6f}°)") + + # 创建网格坐标 + x_grid = np.linspace(min_lon, max_lon, width) + y_grid = np.linspace(max_lat, min_lat, height) # 纬度从上到下递减 + xi, yi = np.meshgrid(x_grid, y_grid) + + # 使用scipy进行插值 + from scipy.interpolate import griddata + + print("开始插值计算...") + + # 方法1: 先尝试线性插值 + try: + zi = griddata((lons, lats), heights, (xi, yi), method='linear') + nan_count = np.isnan(zi).sum() + nan_percent = nan_count / (width * height) * 100 + + print(f"线性插值完成，空白区域: {nan_count} 像素 ({nan_percent:.1f}%)") + + # 如果有空白区域，使用最近邻方法填充 + if nan_count > 0: + print("使用最近邻方法填充空白区域...") + zi_nn = griddata((lons, lats), heights, (xi, yi), method='nearest') + mask = np.isnan(zi) + zi[mask] = zi_nn[mask] + + # 再次检查 + nan_count = np.isnan(zi).sum() + if nan_count > 0: + print(f"仍有 {nan_count} 个空白像素，填充为最低高程值") + min_height = heights.min() + zi[np.isnan(zi)] = min_height + + except Exception as e: + print(f"线性插值失败: {e}") + print("尝试使用最近邻插值...") + zi = griddata((lons, lats), heights, (xi, yi), method='nearest') + + # 创建GeoTIFF + print("创建GeoTIFF文件...") + driver = gdal.GetDriverByName("GTiff") + dem_ds = driver.Create( + output_dem_path, width, height, 1, gdal.GDT_Float32, + options=["COMPRESS=LZW", "TILED=YES", "PREDICTOR=2", "ZLEVEL=9"] + ) + + if dem_ds is None: + raise RuntimeError(f"无法创建DEM文件: {output_dem_path}") + + # 设置投影和地理变换 + srs = osr.SpatialReference() + srs.ImportFromEPSG(4326) # WGS84 + dem_ds.SetProjection(srs.ExportToWkt()) + + geotransform = [ + min_lon, pixel_size, 0, + max_lat, 0, -pixel_size + ] + dem_ds.SetGeoTransform(geotransform) + + # 写入数据 + band = dem_ds.GetRasterBand(1) + band.WriteArray(zi) + band.SetNoDataValue(-9999.0) + band.SetDescription("Elevation") + band.SetUnitType("meters") + + # 计算统计信息 + print("计算统计信息...") + band.FlushCache() + band.ComputeStatistics(False) + + # 设置颜色表（可选） + try: + import matplotlib.pyplot as plt + cmap = plt.cm.terrain + colors = cmap(np.linspace(0, 1, 256)) + colors = (colors[:, :3] * 255).astype(np.uint8) + + color_table = gdal.ColorTable() + for i in range(256): + color_table.SetColorEntry(i, (colors[i, 0], colors[i, 1], colors[i, 2], 255)) + + band.SetColorTable(color_table) + band.SetColorInterpretation(gdal.GCI_PaletteIndex) + except: + pass # 如果设置颜色表失败，继续 + + dem_ds = None # 关闭文件 + + print(f"DEM生成成功: {output_dem_path}") + print(f" 文件大小: {os.path.getsize(output_dem_path) / (1024*1024):.2f} MB") + + # 验证文件 + if os.path.exists(output_dem_path): + ds = gdal.Open(output_dem_path, gdal.GA_ReadOnly) + if ds: + band = ds.GetRasterBand(1) + stats = band.GetStatistics(True, True) + print(f" DEM统计: 最小值={stats[0]:.2f}m, 最大值={stats[1]:.2f}m") + print(f" 平均值={stats[2]:.2f}m, 标准差={stats[3]:.2f}m") + ds = None + else: + print("警告: DEM文件可能未成功生成") + +def generate_dem(REGION_COORDS=None, tileset_path=None, dem_path=None): + # 配置参数 + script_dir = os.path.dirname(os.path.abspath(__file__)) + if tileset_path: + TILESET_PATH = tileset_path + else: + TILESET_PATH = os.path.dirname(script_dir) + "/data/3dtiles/tileset.json" + if dem_path : + OUTPUT_DEM_PATH = dem_path + else : + OUTPUT_DEM_PATH = os.path.join(script_dir, f"o_dem_{uuid.uuid4().hex[:8]}.tif") + PIXEL_SIZE = 0.0001 + + # 执行流程 + print("=" * 60) + print("开始解析3D Tiles...") + if REGION_COORDS: + print(f"启用区域过滤: {REGION_COORDS}") + else: + print("未启用区域过滤，将处理所有数据") + + points = parse_tileset(TILESET_PATH, REGION_COORDS) + print(f"解析完成,共提取点云: {len(points)}个") + + if len(points) > 0: + points_to_dem(points, OUTPUT_DEM_PATH, pixel_size=PIXEL_SIZE) + return OUTPUT_DEM_PATH + else: + print("无点云数据,无法生成DEM") + return None + +if __name__ == "__main__": + # 测试示例 + # 石棉县核心区域 + # SHIMIAN_CORE = [(100.22476304, 29.38340151), (110.32476304, 31.28340151)] + SHIMIAN_CORE = [(100.22476304, 29.18340151), (110.32476304, 31.28340151)] + + # 可以根据需要启用或禁用区域过滤 + REGION_COORDS = SHIMIAN_CORE # 启用区域过滤 + # REGION_COORDS = None # 禁用区域过滤 + + dem_path = generate_dem(REGION_COORDS) + if dem_path: + print(f"\nDEM文件已生成: {dem_path}") \ No newline at end of file diff --git a/b3dm/earthwork_api.py b/b3dm/earthwork_api.py new file mode 100644 index 0000000..246fb97 --- /dev/null +++ b/b3dm/earthwork_api.py @@ -0,0 +1,377 @@ +# pip install fastapi uvicorn pdal pyvista numpy +from sanic import Blueprint, Request, json +from sanic.response import text +from typing import Dict, Any +import logging +from earthwork_calculator_point_cloud import EarthworkCalculatorPointCloud +# 导入计算模块 +from earthwork_calculator_3d_tiles import EarthworkCalculator3dTiles, AlgorithmType, EarthworkResult3dTiles +from tileset_data_source import TilesetDataSource + +earthwork_bp = Blueprint("earthwork", url_prefix="") + +# 配置日志 +logging.basicConfig(level=logging.INFO) +logger = logging.getLogger(__name__) + +# 全局变量 +_calculator_point_cloud = None +_calculator_3d_tiles = None +_data_source_3d_tiles = None + +# 初始化函数 +async def init_app(): + """初始化应用""" + global _data_source_3d_tiles, _calculator_3d_tiles, _calculator_point_cloud + + try: + # 配置数据源 + tileset_path = "./data/3dtiles/tileset.json" + + # 初始化数据源 + _data_source_3d_tiles = TilesetDataSource(tileset_path) + await _data_source_3d_tiles.initialize() + + # 初始化计算器-3dTiles + _calculator_3d_tiles = EarthworkCalculator3dTiles(_data_source_3d_tiles) + + # 初始化计算器-点云 + point_cloud_path = "./data/pointCloud/simulated_points.laz" + _calculator_point_cloud = EarthworkCalculatorPointCloud(point_cloud_path) + + logger.info("土方量计算器初始化完成") + + except ImportError as e: + logger.error(f"依赖库缺失: {str(e)}") + raise + except Exception as e: + logger.error(f"初始化失败: {str(e)}") + raise + +# 土方量计算接口-3dTiles +@earthwork_bp.post("/api/v1/calc/earthwork3dTiles") +async def calc_earthwork(request: Request): + """ + 土方量计算接口 + + 请求参数示例: + { + "polygonCoords": [[120.1, 30.1], [120.2, 30.1], [120.2, 30.2], [120.1, 30.2]], + "designElevation": 50.0, + "algorithm": "tin", + "resolution": 1.0, + "crs": "EPSG:4326", + "interpolationMethod": "linear" + } + """ + try: + # 1. 接收前端传参 + data = request.json + if not data: + return _error_response("请求参数不能为空", 400) + + # 2. 提取参数 + polygon_coords = data.get("polygonCoords") + design_elevation = data.get("designElevation") + + if not polygon_coords: + return _error_response("多边形坐标不能为空", 400) + if design_elevation is None: + return _error_response("设计高程不能为空", 400) + + # 3. 可选参数 + algorithm = data.get("algorithm", "tin") + resolution = data.get("resolution", 1.0) + crs = data.get("crs", "EPSG:4326") + interpolation_method = data.get("interpolationMethod", "linear") + + # 4. 参数验证 + if len(polygon_coords) < 3: + return _error_response("多边形至少需要3个点", 400) + + # 检查多边形是否闭合，如不闭合则自动闭合 + if polygon_coords[0] != polygon_coords[-1]: + polygon_coords.append(polygon_coords[0]) + + # 算法验证 + if algorithm not in ["grid", "tin", "prism"]: + return _error_response("算法必须是 grid, tin 或 prism", 400) + + # 分辨率验证 + if resolution <= 0 or resolution > 100: + return _error_response("分辨率必须在0-100米之间", 400) + + # 5. 确保计算器已初始化 + if _calculator_3d_tiles is None: + await init_app() + + # 6. 执行计算 + algorithm_type = AlgorithmType(algorithm) + + result = await _calculator_3d_tiles.calculate( + polygon_coords=polygon_coords, + design_elevation=design_elevation, + algorithm=algorithm_type, + resolution=resolution, + target_crs=crs, + interpolation_method=interpolation_method + ) + + # 7. 返回成功响应 + return _success_response(result) + + except ValueError as e: + logger.warning(f"参数验证失败: {str(e)}") + return _error_response(f"参数错误: {str(e)}", 400) + except Exception as e: + logger.error(f"计算失败: {str(e)}") + return _error_response(f"服务器内部错误: {str(e)}", 500) + + +# 验证接口 +@earthwork_bp.post("/api/v1/calc/earthwork3dTiles/validate") +async def validate_earthwork(request: Request): + """验证计算参数接口""" + try: + # 1. 接收前端传参 + data = request.json + if not data: + return _error_response("请求参数不能为空", 400) + + # 2. 提取参数 + polygon_coords = data.get("polygonCoords") + + if not polygon_coords: + return _error_response("多边形坐标不能为空", 400) + + # 3. 参数验证 + if len(polygon_coords) < 3: + return _error_response("多边形至少需要3个点", 400) + + # 检查多边形是否闭合，如不闭合则自动闭合 + if polygon_coords[0] != polygon_coords[-1]: + polygon_coords.append(polygon_coords[0]) + + # 4. 确保计算器已初始化 + if _calculator_3d_tiles is None: + await init_app() + + # 5. 执行验证 + validation_result = await _calculator_3d_tiles.validate(polygon_coords) + + # 6. 返回结果 + return _success_response(validation_result) + + except Exception as e: + logger.error(f"验证失败: {str(e)}") + return _error_response(f"验证失败: {str(e)}", 400) + +# 获取算法列表接口 +@earthwork_bp.get("/api/v1/calc/earthwork3dTiles/algorithms") +async def get_algorithms(request: Request): + """获取支持的算法列表接口""" + try: + algorithms = [ + { + "id": "grid", + "name": "格网法", + "description": "将计算区域划分为规则格网，通过插值计算每个格网的高程变化，适合平坦或规则地形", + "accuracy": "中等", + "performance": "快速", + "parameters": { + "resolution": { + "name": "格网分辨率", + "description": "格网大小（米），影响计算精度和性能", + "default": 1.0, + "range": [0.1, 10.0] + } + } + }, + { + "id": "tin", + "name": "三角网法", + "description": "基于不规则三角网（TIN）构建地形表面，计算每个三角形的体积变化，适合复杂地形", + "accuracy": "高", + "performance": "中等", + "parameters": { + "resolution": { + "name": "不适用", + "description": "三角网法不使用固定的分辨率参数", + "default": None + } + } + }, + { + "id": "prism", + "name": "三棱柱法", + "description": "结合三角网和垂直棱柱的高精度算法，计算每个三棱柱的体积，精度最高", + "accuracy": "最高", + "performance": "较慢", + "parameters": { + "resolution": { + "name": "棱柱宽度", + "description": "棱柱宽度（米），影响计算精度", + "default": 1.0, + "range": [0.1, 5.0] + } + } + } + ] + + return _success_response(algorithms) + + except Exception as e: + logger.error(f"获取算法列表失败: {str(e)}") + return _error_response(f"获取算法列表失败: {str(e)}", 500) + +# 批量计算接口 +@earthwork_bp.post("/api/v1/calc/earthwork3dTiles/batch") +async def batch_calc_earthwork(request: Request): + """批量土方量计算接口""" + try: + # 1. 接收前端传参 + data = request.json + if not data: + return _error_response("请求参数不能为空", 400) + + calculations = data.get("calculations", []) + + if not calculations: + return _error_response("计算任务列表不能为空", 400) + + if len(calculations) > 100: + return _error_response("批量计算数量超过限制（最多100个）", 400) + + # 2. 确保计算器已初始化 + if _calculator_3d_tiles is None: + await init_app() + + # 3. 执行批量计算 + results = [] + errors = [] + + for i, calc_data in enumerate(calculations): + try: + # 提取参数 + polygon_coords = calc_data.get("polygonCoords") + design_elevation = calc_data.get("designElevation") + + if not polygon_coords or design_elevation is None: + errors.append({ + "index": i, + "error": "缺少必要参数" + }) + continue + + # 参数验证 + if len(polygon_coords) < 3: + errors.append({ + "index": i, + "error": "多边形至少需要3个点" + }) + continue + + # 检查多边形是否闭合 + if polygon_coords[0] != polygon_coords[-1]: + polygon_coords.append(polygon_coords[0]) + + # 可选参数 + algorithm = calc_data.get("algorithm", "tin") + resolution = calc_data.get("resolution", 1.0) + crs = calc_data.get("crs", "EPSG:4326") + interpolation_method = calc_data.get("interpolationMethod", "linear") + + # 执行计算 + algorithm_type = AlgorithmType(algorithm) + + result = await _calculator_3d_tiles.calculate( + polygon_coords=polygon_coords, + design_elevation=design_elevation, + algorithm=algorithm_type, + resolution=resolution, + target_crs=crs, + interpolation_method=interpolation_method + ) + + results.append(result) + + except Exception as e: + errors.append({ + "index": i, + "error": str(e), + "polygon": polygon_coords if 'polygon_coords' in locals() else None + }) + continue + + # 4. 返回结果 + batch_result = { + "results": results, + "errors": errors, + "summary": { + "total": len(calculations), + "success": len(results), + "failed": len(errors), + "successRate": f"{(len(results)/len(calculations)*100):.1f}%" if calculations else "0%" + } + } + + message = f"批量计算完成，成功 {len(results)} 个，失败 {len(errors)} 个" + return _success_response(batch_result) + + except Exception as e: + logger.error(f"批量计算失败: {str(e)}") + return _error_response(f"批量计算失败: {str(e)}", 500) + +# 核心接口：土方量计算-点云 +@earthwork_bp.post("/api/v1/calc/earthworkPointCloud") +async def calc_earthwork_point_cloud(request: Request): + try: + # 1. 接收前端传参 + data = request.json + if not data: + return json({ + "code": 400, + "msg": "请求参数不能为空", + "data": None + }, status=400) + + polygon_coords = data.get("polygonCoords") # 计算区域多边形坐标 + design_elev = data.get("designElevation") # 设计高程 + crs = data.get("crs", "EPSG:4326") # 坐标系，默认WGS84 + + # 2. 确保计算器已初始化 + if _calculator_point_cloud is None: + await init_app() + + result = _calculator_point_cloud.calculate_earthwork(polygon_coords=polygon_coords, design_elev=design_elev, crs=crs) + + # 3. 处理结果 + if not result["success"]: + return _error_response(result["error"], 400) + + # 4. 格式化结果 + formatted_result = _calculator_point_cloud.format_result(result) + + # 5. 返回成功响应 + return _success_response(formatted_result) + + except Exception as e: + logger.error(f"服务器错误: {str(e)}") + return _error_response(f"服务器内部错误: {str(e)}", 500) + +def _success_response(data: Dict[str, Any]) -> json: + """成功响应""" + + return json({ + "code": 200, + "msg": "计算成功", + "data": data + }) + +def _error_response(message: str, status_code: int = 400) -> json: + """错误响应""" + return json({ + "code": status_code, + "msg": message, + "data": None + }, status=status_code) \ No newline at end of file diff --git a/b3dm/earthwork_calculator_3d_tiles.py b/b3dm/earthwork_calculator_3d_tiles.py new file mode 100644 index 0000000..090e0d7 --- /dev/null +++ b/b3dm/earthwork_calculator_3d_tiles.py @@ -0,0 +1,526 @@ +# earthwork_calculator.py +import numpy as np +from pyproj import Transformer, CRS +from scipy.spatial import Delaunay +from dataclasses import dataclass +from typing import List, Dict, Optional, Tuple, Union +import logging +from enum import Enum +from abc import ABC, abstractmethod +import math + +logger = logging.getLogger(__name__) + +class AlgorithmType(Enum): + """计算算法类型""" + GRID = "grid" + TIN = "tin" + PRISM = "prism" + +@dataclass +class EarthworkResult3dTiles: + """土方量计算结果""" + cut_volume: float # 挖方量 (m³) + fill_volume: float # 填方量 (m³) + net_volume: float # 净方量 (m³) + area: float # 计算区域面积 (m²) + avg_elevation: float # 平均高程 + min_elevation: float # 最低高程 + max_elevation: float # 最高高程 + points_count: int # 使用的点数 + bounding_box: Dict[str, List[float]] # 边界框 + volume_accuracy: float # 计算精度 + algorithm: str # 使用的算法 + resolution: float # 计算分辨率 + + def to_dict(self) -> Dict: + """转换为字典""" + return { + "volume": { + "cut": round(self.cut_volume, 3), + "fill": round(self.fill_volume, 3), + "net": round(self.net_volume, 3), + "unit": "m³" + }, + "area": { + "value": round(self.area, 3), + "unit": "m²" + }, + "elevation": { + "average": round(self.avg_elevation, 3), + "min": round(self.min_elevation, 3), + "max": round(self.max_elevation, 3), + "unit": "m" + }, + "statistics": { + "points_count": self.points_count, + "accuracy": round(self.volume_accuracy, 3), + "algorithm": self.algorithm + }, + "bounding_box": self.bounding_box, + "calculation_params": { + "resolution": self.resolution, + "accuracy": self.volume_accuracy + } + } + +class TerrainDataSource(ABC): + """地形数据源抽象类""" + + @abstractmethod + async def get_points_in_polygon(self, + polygon_coords: List[List[float]], + z_range: Optional[Tuple[float, float]] = None) -> np.ndarray: + """ + 获取多边形区域内的点云数据 + + Args: + polygon_coords: 多边形坐标 [[x1,y1], [x2,y2], ...] + z_range: 高程范围 (min_z, max_z) + + Returns: + Nx3的numpy数组 [x, y, z] + """ + pass + + @abstractmethod + async def get_data_bounds(self) -> Dict[str, List[float]]: + """获取数据范围""" + pass + + @abstractmethod + def get_crs(self) -> str: + """获取数据坐标系""" + pass + +class GeometryUtils: + """几何计算工具类""" + + @staticmethod + def calculate_polygon_area(polygon_coords: List[List[float]]) -> float: + """计算多边形面积（平面面积）""" + polygon_np = np.array(polygon_coords) + x = polygon_np[:, 0] + y = polygon_np[:, 1] + return 0.5 * np.abs(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1))) + + @staticmethod + def is_point_in_polygon(point: np.ndarray, polygon: np.ndarray) -> bool: + """判断点是否在多边形内""" + from matplotlib.path import Path + path = Path(polygon) + return path.contains_point(point) + + @staticmethod + def calculate_triangle_area(points: np.ndarray) -> float: + """计算三角形面积""" + a = np.linalg.norm(points[0] - points[1]) + b = np.linalg.norm(points[1] - points[2]) + c = np.linalg.norm(points[2] - points[0]) + s = (a + b + c) / 2 + return np.sqrt(s * (s - a) * (s - b) * (s - c)) + + @staticmethod + def create_grid(polygon: np.ndarray, resolution: float) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """创建规则格网""" + x_min, y_min = polygon.min(axis=0) + x_max, y_max = polygon.max(axis=0) + + # 扩展一个格网单元 + x_min -= resolution + x_max += resolution + y_min -= resolution + y_max += resolution + + x_grid = np.arange(x_min, x_max + resolution, resolution) + y_grid = np.arange(y_min, y_max + resolution, resolution) + xx, yy = np.meshgrid(x_grid, y_grid) + + return xx, yy, x_grid, y_grid + + @staticmethod + def interpolate_grid(xx: np.ndarray, yy: np.ndarray, + points: np.ndarray, method: str = 'linear') -> np.ndarray: + """格网插值""" + from scipy.interpolate import LinearNDInterpolator, CloughTocher2DInterpolator + + grid_points = np.column_stack([xx.ravel(), yy.ravel()]) + + if method == 'linear': + interpolator = LinearNDInterpolator(points[:, :2], points[:, 2]) + elif method == 'cubic': + interpolator = CloughTocher2DInterpolator(points[:, :2], points[:, 2]) + else: + raise ValueError(f"不支持的插值方法: {method}") + + elevations = interpolator(grid_points) + return elevations.reshape(xx.shape) + +class EarthworkCalculator3dTiles: + """土方量计算器""" + + def __init__(self, data_source: TerrainDataSource): + """ + 初始化计算器 + + Args: + data_source: 地形数据源 + """ + self.data_source = data_source + self._transformer_cache = {} + + async def calculate(self, + polygon_coords: List[List[float]], + design_elevation: float, + algorithm: AlgorithmType = AlgorithmType.TIN, + resolution: float = 1.0, + target_crs: str = "EPSG:4326", + interpolation_method: str = 'linear') -> EarthworkResult3dTiles: + """ + 计算土方量 + + Args: + polygon_coords: 多边形坐标 + design_elevation: 设计高程 + algorithm: 计算算法 + resolution: 格网分辨率（米） + target_crs: 目标坐标系 + interpolation_method: 插值方法 + + Returns: + EarthworkResult: 计算结果 + """ + try: + # 1. 获取数据 + points = await self.data_source.get_points_in_polygon(polygon_coords) + + if points.size == 0: + raise ValueError("区域内没有找到高程数据") + + # 2. 坐标转换 + points = await self._transform_coordinates(points, target_crs) + + # 3. 执行计算 + if algorithm == AlgorithmType.GRID: + result = await self._calculate_grid(points, polygon_coords, + design_elevation, resolution, + interpolation_method) + elif algorithm == AlgorithmType.TIN: + result = await self._calculate_tin(points, polygon_coords, + design_elevation) + elif algorithm == AlgorithmType.PRISM: + result = await self._calculate_prism(points, polygon_coords, + design_elevation, resolution) + else: + raise ValueError(f"不支持的算法: {algorithm}") + + return result + + except Exception as e: + logger.error(f"土方量计算失败: {str(e)}") + raise + + async def _transform_coordinates(self, points: np.ndarray, + target_crs: str) -> np.ndarray: + """坐标转换""" + source_crs = self.data_source.get_crs() + + if source_crs == target_crs: + return points + + cache_key = f"{source_crs}->{target_crs}" + if cache_key not in self._transformer_cache: + self._transformer_cache[cache_key] = Transformer.from_crs( + CRS.from_string(source_crs), + CRS.from_string(target_crs), + always_xy=True + ) + + transformer = self._transformer_cache[cache_key] + points_2d = transformer.transform(points[:, 0], points[:, 1]) + + return np.column_stack([points_2d[0], points_2d[1], points[:, 2]]) + + async def _calculate_grid(self, points: np.ndarray, + polygon_coords: List[List[float]], + design_elevation: float, + resolution: float, + interpolation_method: str) -> EarthworkResult3dTiles: + """格网法计算""" + polygon_np = np.array(polygon_coords) + + # 创建格网 + xx, yy, x_grid, y_grid = GeometryUtils.create_grid(polygon_np, resolution) + + # 插值 + natural_elevations = GeometryUtils.interpolate_grid(xx, yy, points, interpolation_method) + + # 初始化挖填量 + cut_volume = 0.0 + fill_volume = 0.0 + total_area = 0.0 + + # 遍历每个格网单元 + for i in range(len(x_grid) - 1): + for j in range(len(y_grid) - 1): + # 格网四个角点 + cell_corners = np.array([ + [x_grid[i], y_grid[j]], + [x_grid[i+1], y_grid[j]], + [x_grid[i+1], y_grid[j+1]], + [x_grid[i], y_grid[j+1]] + ]) + + # 检查格网中心点是否在多边形内 + cell_center = cell_corners.mean(axis=0) + if not GeometryUtils.is_point_in_polygon(cell_center, polygon_np): + continue + + # 获取格网四个角点的高程 + cell_elevations = [ + natural_elevations[j, i], + natural_elevations[j, i+1], + natural_elevations[j+1, i+1], + natural_elevations[j+1, i] + ] + + # 检查是否有无效数据 + if any(np.isnan(elev) for elev in cell_elevations): + continue + + # 计算格网平均高程 + avg_elevation = np.mean(cell_elevations) + cell_area = resolution * resolution + total_area += cell_area + + # 计算挖填量 + height_diff = design_elevation - avg_elevation + if height_diff > 0: + fill_volume += height_diff * cell_area + else: + cut_volume += abs(height_diff) * cell_area + + # 计算统计信息 + area = GeometryUtils.calculate_polygon_area(polygon_coords) + mask = ~np.isnan(natural_elevations) + valid_elevations = natural_elevations[mask] + + return EarthworkResult3dTiles( + cut_volume=cut_volume, + fill_volume=fill_volume, + net_volume=cut_volume - fill_volume, + area=area, + avg_elevation=np.mean(valid_elevations) if valid_elevations.size > 0 else 0, + min_elevation=np.min(valid_elevations) if valid_elevations.size > 0 else 0, + max_elevation=np.max(valid_elevations) if valid_elevations.size > 0 else 0, + points_count=points.shape[0], + bounding_box={ + "min": [x_grid[0], y_grid[0]], + "max": [x_grid[-1], y_grid[-1]] + }, + volume_accuracy=self._calculate_accuracy(points, resolution), + algorithm=AlgorithmType.GRID.value, + resolution=resolution + ).to_dict() + + async def _calculate_tin(self, points: np.ndarray, + polygon_coords: List[List[float]], + design_elevation: float) -> EarthworkResult3dTiles: + """三角网法计算""" + polygon_np = np.array(polygon_coords) + + # 创建Delaunay三角网 + triangulation = Delaunay(points[:, :2]) + + # 筛选多边形内的三角形 + cut_volume = 0.0 + fill_volume = 0.0 + total_area = 0.0 + + for simplex in triangulation.simplices: + triangle_points = points[simplex] + triangle_center = triangle_points.mean(axis=0)[:2] + + # 检查三角形中心是否在多边形内 + if not GeometryUtils.is_point_in_polygon(triangle_center, polygon_np): + continue + + # 计算三角形面积 + area = GeometryUtils.calculate_triangle_area(triangle_points[:, :2]) + total_area += area + + # 计算平均高程（使用三个顶点的高程） + avg_elevation = triangle_points[:, 2].mean() + + # 计算挖填量 + height_diff = design_elevation - avg_elevation + if height_diff > 0: + fill_volume += height_diff * area + else: + cut_volume += abs(height_diff) * area + + # 计算统计信息 + area = GeometryUtils.calculate_polygon_area(polygon_coords) + + return EarthworkResult3dTiles( + cut_volume=cut_volume, + fill_volume=fill_volume, + net_volume=cut_volume - fill_volume, + area=area, + avg_elevation=points[:, 2].mean(), + min_elevation=points[:, 2].min(), + max_elevation=points[:, 2].max(), + points_count=points.shape[0], + bounding_box={ + "min": points.min(axis=0)[:2].tolist(), + "max": points.max(axis=0)[:2].tolist() + }, + volume_accuracy=self._calculate_accuracy(points, 0), + algorithm=AlgorithmType.TIN.value, + resolution=0 + ).to_dict() + + async def _calculate_prism(self, points: np.ndarray, + polygon_coords: List[List[float]], + design_elevation: float, + resolution: float) -> EarthworkResult3dTiles: + """三棱柱法计算""" + # 先创建TIN + polygon_np = np.array(polygon_coords) + triangulation = Delaunay(points[:, :2]) + + cut_volume = 0.0 + fill_volume = 0.0 + total_area = 0.0 + + for simplex in triangulation.simplices: + triangle_points = points[simplex] + triangle_center = triangle_points.mean(axis=0)[:2] + + if not GeometryUtils.is_point_in_polygon(triangle_center, polygon_np): + continue + + # 计算三角形面积 + area = GeometryUtils.calculate_triangle_area(triangle_points[:, :2]) + total_area += area + + # 对于每个三角形，计算三棱柱体积 + # 这里简化处理，使用梯形公式 + for i in range(3): + j = (i + 1) % 3 + # 计算边的长度 + edge_length = np.linalg.norm(triangle_points[i, :2] - triangle_points[j, :2]) + + # 计算边两端点的挖填高度 + height_i = design_elevation - triangle_points[i, 2] + height_j = design_elevation - triangle_points[j, 2] + + # 计算边的平均挖填高度 + avg_height = (abs(height_i) + abs(height_j)) / 2 + + # 计算边的面积（假设边宽度为resolution） + edge_area = edge_length * resolution + + if height_i > 0 or height_j > 0: + fill_volume += avg_height * edge_area + else: + cut_volume += avg_height * edge_area + + area = GeometryUtils.calculate_polygon_area(polygon_coords) + + return EarthworkResult3dTiles( + cut_volume=cut_volume, + fill_volume=fill_volume, + net_volume=cut_volume - fill_volume, + area=area, + avg_elevation=points[:, 2].mean(), + min_elevation=points[:, 2].min(), + max_elevation=points[:, 2].max(), + points_count=points.shape[0], + bounding_box={ + "min": points.min(axis=0)[:2].tolist(), + "max": points.max(axis=0)[:2].tolist() + }, + volume_accuracy=self._calculate_accuracy(points, resolution), + algorithm=AlgorithmType.PRISM.value, + resolution=resolution + ).to_dict() + + def _calculate_accuracy(self, points: np.ndarray, resolution: float) -> float: + """计算精度评估""" + if points.shape[0] < 10: + return 0.5 + + # 基于点密度和分辨率估算精度 + if points.shape[0] > 0: + # 计算点密度 + bounds = points[:, :2] + area = (bounds[:, 0].max() - bounds[:, 0].min()) * \ + (bounds[:, 1].max() - bounds[:, 1].min()) + + if area > 0: + point_density = points.shape[0] / area + + if point_density > 100 and resolution <= 0.5: # 高密度，高分辨率 + return 0.05 + elif point_density > 10 and resolution <= 1.0: + return 0.1 + elif point_density > 1: + return 0.2 + + return 0.3 + + async def validate(self, polygon_coords: List[List[float]]) -> Dict: + """验证计算参数""" + try: + points = await self.data_source.get_points_in_polygon(polygon_coords) + + validation_result = { + "polygon_valid": len(polygon_coords) >= 3, + "area": GeometryUtils.calculate_polygon_area(polygon_coords), + "points_available": points.size > 0, + "points_count": points.shape[0] if points.size > 0 else 0, + "data_quality": "good" if points.shape[0] > 100 else "poor", + "suggested_algorithm": self._suggest_algorithm(points), + "estimated_accuracy": self._calculate_accuracy(points, 1.0) if points.size > 0 else None + } + + if points.size > 0: + validation_result.update({ + "elevation_range": { + "min": float(points[:, 2].min()), + "max": float(points[:, 2].max()), + "avg": float(points[:, 2].mean()) + }, + "bounding_box": { + "min": points.min(axis=0).tolist(), + "max": points.max(axis=0).tolist() + } + }) + + return validation_result + + except Exception as e: + logger.error(f"验证失败: {str(e)}") + return { + "polygon_valid": False, + "error": str(e) + } + + def _suggest_algorithm(self, points: np.ndarray) -> str: + """建议计算算法""" + if points.shape[0] < 100: + return AlgorithmType.TIN.value + + # 计算点的规则性 + bounds = points[:, :2] + x_range = bounds[:, 0].max() - bounds[:, 0].min() + y_range = bounds[:, 1].max() - bounds[:, 1].min() + + # 如果点分布相对规则，建议使用格网法 + if x_range > 0 and y_range > 0: + aspect_ratio = max(x_range, y_range) / min(x_range, y_range) + if aspect_ratio < 2: # 相对规则 + return AlgorithmType.GRID.value + + return AlgorithmType.TIN.value \ No newline at end of file diff --git a/b3dm/earthwork_calculator_point_cloud.py b/b3dm/earthwork_calculator_point_cloud.py new file mode 100644 index 0000000..3d4793e --- /dev/null +++ b/b3dm/earthwork_calculator_point_cloud.py @@ -0,0 +1,691 @@ +# earthwork_calculator.py +import pdal +import pyvista as pv +import numpy as np +import json +from typing import List, Tuple, Dict, Any, Optional +from dataclasses import dataclass, field +from enum import Enum +import traceback + + +class EarthworkAlgorithm(Enum): + """土方量计算算法枚举""" + GRID = "grid" + TIN = "tin" + PRISM = "prism" + + +@dataclass +class EarthworkResultPointCloud: + """土方量计算结果""" + cut_volume: float # 挖方量 (m³) + fill_volume: float # 填方量 (m³) + net_volume: float # 净方量 (m³) + area: float # 计算区域面积 (m²) + avg_elevation: float # 平均高程 + min_elevation: float # 最低高程 + max_elevation: float # 最高高程 + points_count: int # 使用的点数 + triangle_count: int = 0 # 三角形数量 + grid_count: int = 0 # 网格数量（仅GRID算法使用） + prism_count: int = 0 # 棱柱体数量（仅PRISM算法使用） + bounding_box: Dict[str, List[float]] = field(default_factory=dict) # 边界框 + volume_accuracy: float = 0.95 # 计算精度 + algorithm: str = "TIN三角网法" # 使用的算法 + resolution: float = 1.0 # 计算分辨率 + algorithm_params: Dict[str, Any] = field(default_factory=dict) # 算法参数 + + def to_dict(self) -> Dict: + """转换为字典""" + result = { + "volume": { + "cut": round(self.cut_volume, 3), + "fill": round(self.fill_volume, 3), + "net": round(self.net_volume, 3), + "unit": "m³" + }, + "area": { + "value": round(self.area, 3), + "unit": "m²" + }, + "elevation": { + "average": round(self.avg_elevation, 3), + "min": round(self.min_elevation, 3), + "max": round(self.max_elevation, 3), + "unit": "m" + }, + "statistics": { + "points_count": self.points_count, + "accuracy": round(self.volume_accuracy, 3), + "algorithm": self.algorithm + }, + "bounding_box": self.bounding_box, + "calculation_params": { + "resolution": self.resolution, + "accuracy": self.volume_accuracy, + **self.algorithm_params + } + } + + # 根据算法类型添加特定的统计信息 + if self.algorithm.startswith("GRID"): + result["statistics"]["grid_count"] = self.grid_count + elif self.algorithm.startswith("TIN"): + result["statistics"]["triangle_count"] = self.triangle_count + elif self.algorithm.startswith("PRISM"): + result["statistics"]["prism_count"] = self.prism_count + + return result + + +class EarthworkCalculatorPointCloud: + """土方量计算核心类（支持多种算法）""" + + def __init__(self, point_cloud_path: str = "./data/your_point_cloud.laz"): + """ + 初始化土方量计算器 + + Args: + point_cloud_path: 点云数据文件路径 + """ + self.point_cloud_path = point_cloud_path + + def validate_inputs(self, polygon_coords: List[List[float]], design_elev: float) -> Tuple[bool, str]: + """验证输入参数""" + if not polygon_coords or len(polygon_coords) < 3: + return False, "多边形坐标至少需要3个点" + + try: + design_elev = float(design_elev) + except (TypeError, ValueError): + return False, "设计高程必须是有效数字" + + return True, "" + + def create_polygon_string(self, polygon_coords: List[List[float]]) -> str: + """创建PDAL多边形字符串""" + coords_list = [] + for coord in polygon_coords: + if len(coord) >= 2: + coords_list.append(f"{coord[0]} {coord[1]}") + + # 确保多边形闭合 + if coords_list and coords_list[0] != coords_list[-1]: + coords_list.append(coords_list[0]) + + return "POLYGON((" + ", ".join(coords_list) + "))" + + def calculate_bounding_box(self, points: np.ndarray) -> Dict[str, List[float]]: + """ + 计算边界框 + + Args: + points: 点云坐标数组 + + Returns: + Dict: 边界框信息 + """ + if len(points) == 0: + return {"min": [0, 0, 0], "max": [0, 0, 0]} + + min_vals = np.min(points, axis=0) + max_vals = np.max(points, axis=0) + + return { + "min": [float(min_vals[0]), float(min_vals[1]), float(min_vals[2])], + "max": [float(max_vals[0]), float(max_vals[1]), float(max_vals[2])] + } + + def clip_point_cloud(self, polygon_coords: List[List[float]], crs: str = "EPSG:4326") -> pv.PolyData: + """ + 裁剪点云数据 + + Args: + polygon_coords: 多边形坐标列表 + crs: 坐标系 + + Returns: + pyvista.PolyData: 裁剪后的点云数据 + """ + polygon_str = self.create_polygon_string(polygon_coords) + + # PDAL管道配置 + pipeline_config = { + "pipeline": [ + { + "type": "readers.las", + "filename": self.point_cloud_path, + "spatialreference": crs + }, + { + "type": "filters.crop", + "polygon": polygon_str + } + ] + } + + # 执行PDAL管道 + pipeline = pdal.Pipeline(json.dumps(pipeline_config)) + + try: + pipeline.execute() + + if len(pipeline.arrays) == 0: + raise ValueError("多边形区域内没有找到点云数据") + + # 获取裁剪后的点云数据 + points = pipeline.arrays[0] + x = points["X"] + y = points["Y"] + z = points["Z"] + + return pv.PolyData(np.column_stack((x, y, z))) + + except RuntimeError as e: + print(f"PDAL执行失败: {str(e)}") + # 如果没有PDAL数据，生成模拟数据用于测试 + return self.generate_mock_point_cloud(polygon_coords) + + def generate_mock_point_cloud(self, polygon_coords: List[List[float]]) -> pv.PolyData: + """ + 生成模拟点云数据（仅用于测试） + + Args: + polygon_coords: 多边形坐标列表 + + Returns: + pyvista.PolyData: 模拟点云数据 + """ + print("使用模拟数据进行测试...") + + n_points = 1000 + + # 获取多边形边界 + x_coords = [c[0] for c in polygon_coords] + y_coords = [c[1] for c in polygon_coords] + + x_min, x_max = min(x_coords), max(x_coords) + y_min, y_max = min(y_coords), max(y_coords) + + # 生成随机点 + x = np.random.uniform(x_min, x_max, n_points) + y = np.random.uniform(y_min, y_max, n_points) + z = np.random.uniform(100, 120, n_points) # 模拟高程在100-120米之间 + + return pv.PolyData(np.column_stack((x, y, z))) + + def create_tin_mesh(self, point_cloud: pv.PolyData) -> pv.PolyData: + """ + 创建三角网（TIN算法使用） + + Args: + point_cloud: 点云数据 + + Returns: + pyvista.PolyData: 三角网格 + """ + if len(point_cloud.points) < 3: + raise ValueError("点云数据不足，无法构网") + + try: + return point_cloud.delaunay_2d() + except Exception as e: + raise ValueError(f"三角网构网失败: {str(e)}") + + def calculate_volumes_by_tin(self, mesh: pv.PolyData, design_elev: float) -> Dict[str, Any]: + """ + TIN算法计算土方量 + + Args: + mesh: 三角网格 + design_elev: 设计高程 + + Returns: + Dict: 包含体积计算结果和额外信息 + """ + points = mesh.points + elev_diff = points[:, 2] - design_elev + + cut_volume = 0.0 + fill_volume = 0.0 + triangle_count = 0 + + # 遍历所有三角形面片计算体积 + cells = mesh.cells.reshape(-1, 4) + + if len(cells) == 0: + raise ValueError("无法生成有效的三角网") + + for cell in cells: + if cell[0] == 3: # 三角形（VTK格式：3个顶点） + triangle_count += 1 + vertex_indices = cell[1:] + pts = points[vertex_indices] + + # 计算三角形面积 + v1 = pts[1] - pts[0] + v2 = pts[2] - pts[0] + area = 0.5 * np.linalg.norm(np.cross(v1, v2)) + + # 计算平均高程差 + avg_diff = np.mean(elev_diff[vertex_indices]) + vol = area * avg_diff + + if vol > 0: + cut_volume += vol + else: + fill_volume += abs(vol) + + return { + "cut_volume": cut_volume, + "fill_volume": fill_volume, + "net_volume": cut_volume - fill_volume, + "triangle_count": triangle_count + } + + def calculate_volumes_by_grid(self, point_cloud: pv.PolyData, design_elev: float, + grid_size: float = 1.0) -> Dict[str, Any]: + """ + GRID算法计算土方量 + + Args: + point_cloud: 点云数据 + design_elev: 设计高程 + grid_size: 网格尺寸 + + Returns: + Dict: 包含体积计算结果和额外信息 + """ + points = point_cloud.points + + if len(points) == 0: + raise ValueError("点云数据为空") + + # 计算边界 + x_min, y_min = np.min(points[:, :2], axis=0) + x_max, y_max = np.max(points[:, :2], axis=0) + + # 创建网格 + x_edges = np.arange(x_min, x_max + grid_size, grid_size) + y_edges = np.arange(y_min, y_max + grid_size, grid_size) + + grid_count = (len(x_edges) - 1) * (len(y_edges) - 1) + + cut_volume = 0.0 + fill_volume = 0.0 + + # 对每个网格计算土方量 + for i in range(len(x_edges) - 1): + for j in range(len(y_edges) - 1): + # 获取当前网格内的点 + mask = (points[:, 0] >= x_edges[i]) & (points[:, 0] < x_edges[i+1]) & \ + (points[:, 1] >= y_edges[j]) & (points[:, 1] < y_edges[j+1]) + + grid_points = points[mask] + + if len(grid_points) > 0: + # 计算网格内点的平均高程 + avg_elevation = np.mean(grid_points[:, 2]) + + # 计算高程差 + elev_diff = avg_elevation - design_elev + + # 计算体积 + cell_area = grid_size * grid_size + vol = cell_area * elev_diff + + if vol > 0: + cut_volume += vol + else: + fill_volume += abs(vol) + + return { + "cut_volume": cut_volume, + "fill_volume": fill_volume, + "net_volume": cut_volume - fill_volume, + "grid_count": grid_count + } + + def calculate_volumes_by_prism(self, point_cloud: pv.PolyData, design_elev: float, + influence_radius: float = 0.5) -> Dict[str, Any]: + """ + PRISM算法计算土方量 + + Args: + point_cloud: 点云数据 + design_elev: 设计高程 + influence_radius: 影响半径 + + Returns: + Dict: 包含体积计算结果和额外信息 + """ + points = point_cloud.points + + if len(points) == 0: + raise ValueError("点云数据为空") + + cut_volume = 0.0 + fill_volume = 0.0 + + # 每个点的影响面积 + influence_area = np.pi * influence_radius ** 2 + prism_count = len(points) + + for point in points: + # 计算高程差 + elev_diff = point[2] - design_elev + + # 计算体积 + vol = influence_area * elev_diff + + if vol > 0: + cut_volume += vol + else: + fill_volume += abs(vol) + + return { + "cut_volume": cut_volume, + "fill_volume": fill_volume, + "net_volume": cut_volume - fill_volume, + "prism_count": prism_count + } + + def calculate_statistics(self, point_cloud: pv.PolyData, mesh: pv.PolyData = None) -> Dict[str, float]: + """ + 计算统计数据 + + Args: + point_cloud: 点云数据 + mesh: 三角网格（仅TIN算法需要） + + Returns: + Dict: 统计结果 + """ + elevations = point_cloud.points[:, 2] + + stats = { + "area": 0.0, + "max_elevation": np.max(elevations) if len(elevations) > 0 else 0.0, + "min_elevation": np.min(elevations) if len(elevations) > 0 else 0.0, + "avg_elevation": np.mean(elevations) if len(elevations) > 0 else 0.0, + "points_count": len(point_cloud.points) + } + + # 计算面积 + if mesh is not None: + stats["area"] = mesh.area + else: + # 对于非TIN算法，使用多边形面积近似 + if len(point_cloud.points) > 0: + # 使用点云的凸包面积 + try: + hull = point_cloud.delaunay_2d() + stats["area"] = hull.area + except: + # 如果无法计算凸包，使用边界框面积 + x_min, x_max = np.min(point_cloud.points[:, 0]), np.max(point_cloud.points[:, 0]) + y_min, y_max = np.min(point_cloud.points[:, 1]), np.max(point_cloud.points[:, 1]) + stats["area"] = (x_max - x_min) * (y_max - y_min) + + return stats + + def calculate_earthwork(self, + polygon_coords: List[List[float]], + design_elev: float, + algorithm: str = EarthworkAlgorithm.TIN.value, + algorithm_params: Optional[Dict[str, Any]] = None, + crs: str = "EPSG:4326", + volume_accuracy: Optional[float] = None, + resolution: Optional[float] = None) -> EarthworkResultPointCloud: + """ + 主计算方法：执行完整的土方量计算流程 + + Args: + polygon_coords: 多边形坐标列表 + design_elev: 设计高程 + algorithm: 计算算法，可选值：'grid', 'tin', 'prism' + algorithm_params: 算法特定参数 + crs: 坐标系 + volume_accuracy: 计算精度（0-1之间） + resolution: 计算分辨率 + + Returns: + EarthworkResultPointCloud: 计算结果 + """ + try: + # 1. 验证输入 + is_valid, message = self.validate_inputs(polygon_coords, design_elev) + if not is_valid: + raise ValueError(message) + + design_elev = float(design_elev) + + # 2. 验证算法参数 + if algorithm not in [a.value for a in EarthworkAlgorithm]: + raise ValueError(f"不支持的算法: {algorithm}。支持的算法: {[a.value for a in EarthworkAlgorithm]}") + + # 3. 设置默认参数 + if algorithm_params is None: + algorithm_params = {} + + # 4. 裁剪点云 + point_cloud = self.clip_point_cloud(polygon_coords, crs) + + # 5. 根据算法选择计算方法 + algorithm_name = "" + mesh = None + + if algorithm == EarthworkAlgorithm.TIN.value: + algorithm_name = "TIN三角网法" + mesh = self.create_tin_mesh(point_cloud) + volumes = self.calculate_volumes_by_tin(mesh, design_elev) + algorithm_params = { + "grid_size": algorithm_params.get("grid_size", 1.0) + } + + elif algorithm == EarthworkAlgorithm.GRID.value: + algorithm_name = "GRID格网法" + grid_size = algorithm_params.get("grid_size", 1.0) + algorithm_name = f"GRID格网法(网格尺寸={grid_size}m)" + volumes = self.calculate_volumes_by_grid(point_cloud, design_elev, grid_size) + algorithm_params = { + "grid_size": grid_size + } + + elif algorithm == EarthworkAlgorithm.PRISM.value: + algorithm_name = "PRISM棱柱体法" + influence_radius = algorithm_params.get("influence_radius", 0.5) + algorithm_name = f"PRISM棱柱体法(影响半径={influence_radius}m)" + volumes = self.calculate_volumes_by_prism(point_cloud, design_elev, influence_radius) + algorithm_params = { + "influence_radius": influence_radius + } + + # 6. 计算统计数据 + stats = self.calculate_statistics(point_cloud, mesh) + + # 7. 计算边界框 + bounding_box = self.calculate_bounding_box(point_cloud.points) + + # 8. 计算或使用传入的精度和分辨率 + if volume_accuracy is None: + # 根据算法和点云密度自动估算精度 + volume_accuracy = self.estimate_accuracy(algorithm, point_cloud) + + if resolution is None: + # 根据点云密度自动估算分辨率 + resolution = self.estimate_resolution(point_cloud) + + # 9. 创建EarthworkResultPointCloud对象 + result = EarthworkResultPointCloud( + cut_volume=volumes["cut_volume"], + fill_volume=volumes["fill_volume"], + net_volume=volumes["net_volume"], + area=stats["area"], + avg_elevation=stats["avg_elevation"], + min_elevation=stats["min_elevation"], + max_elevation=stats["max_elevation"], + points_count=stats["points_count"], + triangle_count=volumes.get("triangle_count", 0), + grid_count=volumes.get("grid_count", 0), + prism_count=volumes.get("prism_count", 0), + bounding_box=bounding_box, + volume_accuracy=volume_accuracy, + algorithm=algorithm_name, + resolution=resolution, + algorithm_params=algorithm_params + ) + + return result + + except Exception as e: + print(f"计算错误: {str(e)}") + print(traceback.format_exc()) + raise + + def estimate_accuracy(self, algorithm: str, point_cloud: pv.PolyData) -> float: + """根据算法和点云密度估计计算精度""" + point_density = len(point_cloud.points) / max(point_cloud.area, 0.1) + + # 基础精度 + base_accuracy = { + EarthworkAlgorithm.TIN.value: 0.95, + EarthworkAlgorithm.GRID.value: 0.90, + EarthworkAlgorithm.PRISM.value: 0.85 + }.get(algorithm, 0.90) + + # 根据点云密度调整精度 + if point_density > 10: # 高密度点云 + accuracy_boost = min(0.05, point_density * 0.002) + elif point_density < 1: # 低密度点云 + accuracy_boost = -0.05 + else: + accuracy_boost = 0 + + estimated_accuracy = base_accuracy + accuracy_boost + + # 确保精度在合理范围内 + return max(0.7, min(0.99, estimated_accuracy)) + + def estimate_resolution(self, point_cloud: pv.PolyData) -> float: + """根据点云密度估计分辨率""" + if len(point_cloud.points) < 2: + return 1.0 + + area = point_cloud.area + if area > 0: + point_count = len(point_cloud.points) + avg_spacing = np.sqrt(area / point_count) + return float(round(avg_spacing, 2)) + + return 1.0 + + def get_algorithm_info(self) -> Dict[str, Any]: + """ + 获取支持的算法信息 + + Returns: + Dict: 算法信息 + """ + return { + "supported_algorithms": [ + { + "id": EarthworkAlgorithm.TIN.value, + "name": "TIN三角网法", + "description": "通过构建不规则三角网计算土方量，精度高，适合复杂地形", + "parameters": [ + { + "name": "grid_size", + "type": "float", + "default": 1.0, + "description": "网格尺寸(m)，用于点云预处理", + "required": False + } + ] + }, + { + "id": EarthworkAlgorithm.GRID.value, + "name": "GRID格网法", + "description": "将区域划分为规则网格计算土方量，计算速度快，适合大规模区域", + "parameters": [ + { + "name": "grid_size", + "type": "float", + "default": 1.0, + "description": "网格尺寸(m)", + "required": True, + "min": 0.1, + "max": 10.0 + } + ] + }, + { + "id": EarthworkAlgorithm.PRISM.value, + "name": "PRISM棱柱体法", + "description": "将每个点视为一个棱柱体计算土方量，计算简单快速", + "parameters": [ + { + "name": "influence_radius", + "type": "float", + "default": 0.5, + "description": "点的影响半径(m)", + "required": True, + "min": 0.1, + "max": 5.0 + } + ] + } + ], + "default_algorithm": EarthworkAlgorithm.TIN.value + } + + +# 使用示例 +if __name__ == "__main__": + # 创建计算器实例 + calculator = EarthworkCalculatorPointCloud("./data/sample_point_cloud.laz") + + # 获取支持的算法信息 + algorithm_info = calculator.get_algorithm_info() + print("支持的算法:") + for algo in algorithm_info["supported_algorithms"]: + print(f" {algo['id']}: {algo['name']} - {algo['description']}") + + # 定义多边形区域 + polygon = [ + [116.3974, 39.9093], + [116.4084, 39.9093], + [116.4084, 39.9193], + [116.3974, 39.9193] + ] + + # 设计高程 + design_elevation = 100.0 + + # 测试不同算法 + algorithms = [ + (EarthworkAlgorithm.TIN.value, {}, "TIN算法"), + (EarthworkAlgorithm.GRID.value, {"grid_size": 2.0}, "GRID算法(2米网格)"), + (EarthworkAlgorithm.PRISM.value, {"influence_radius": 1.0}, "PRISM算法(1米影响半径)") + ] + + for algo_id, params, description in algorithms: + print(f"\n使用{description}计算:") + try: + result = calculator.calculate_earthwork( + polygon_coords=polygon, + design_elev=design_elevation, + algorithm=algo_id, + algorithm_params=params + ) + + print(f" 挖方量: {result.cut_volume:.3f} m³") + print(f" 填方量: {result.fill_volume:.3f} m³") + print(f" 净方量: {result.net_volume:.3f} m³") + print(f" 计算面积: {result.area:.3f} m²") + print(f" 计算精度: {result.volume_accuracy:.3%}") + print(f" 分辨率: {result.resolution:.2f} m") + + except Exception as e: + print(f" 计算失败: {str(e)}") \ No newline at end of file diff --git a/b3dm/glb_with_draco.py b/b3dm/glb_with_draco.py new file mode 100644 index 0000000..460528b --- /dev/null +++ b/b3dm/glb_with_draco.py @@ -0,0 +1,469 @@ +import json +import struct +import numpy as np +import DracoPy + +class DracoGLBParser: + """使用 DracoPy 解析包含 Draco 压缩的 GLB 文件""" + + def __init__(self, glb_file_path): + self.glb_file_path = glb_file_path + self.json_data = None + self.binary_data = None + self.decoded_meshes = [] # 缓存解码后的网格数据 + + def parse_glb_structure(self): + """解析 GLB 文件结构""" + with open(self.glb_file_path, 'rb') as f: + # 读取 GLB 头部 + magic = f.read(4) + version = struct.unpack(' 0: + if isinstance(faces_data[0], list) or isinstance(faces_data[0], tuple): + # 如果是列表的列表 + result['faces'] = np.array(faces_data, dtype=np.uint32) + else: + # 如果是扁平化的数组 + result['faces'] = np.array(faces_data, dtype=np.uint32).reshape(-1, 3) + + print(f" 面数量: {len(result['faces']) if result['faces'] is not None else 0}") + + # 获取属性数据 + if hasattr(draco_decoder, 'attributes'): + attrs = draco_decoder.attributes + + # 根据属性映射查找数据 + for gltf_attr_name, draco_attr_id in attributes.items(): + if draco_attr_id in attrs: + attr_data = attrs[draco_attr_id] + + if gltf_attr_name == 'POSITION': + result['vertices'] = np.array(attr_data, dtype=np.float32) + elif gltf_attr_name == 'TEXCOORD_0': + result['texcoords'] = np.array(attr_data, dtype=np.float32) + elif gltf_attr_name == '_BATCHID': + result['batch_ids'] = np.array(attr_data, dtype=np.uint32) + elif gltf_attr_name == 'NORMAL': + result['normals'] = np.array(attr_data, dtype=np.float32) + elif gltf_attr_name == 'COLOR_0': + result['colors'] = np.array(attr_data, dtype=np.float32) + + print(f" 已提取属性: {gltf_attr_name} (ID: {draco_attr_id})") + + # 如果没有通过attributes获取到顶点，尝试其他方式 + if result['vertices'] is None and hasattr(draco_decoder, 'get_points'): + try: + result['vertices'] = np.array(draco_decoder.get_points(), dtype=np.float32) + except: + pass + + return result + + def save_decoded_meshes(self, meshes, output_format='obj'): + """保存解码后的网格""" + import os + + base_name = os.path.splitext(os.path.basename(self.glb_file_path))[0] + output_dir = f"{base_name}_decoded" + + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + for mesh in meshes: + filename = f"{output_dir}/{mesh['name']}_p{mesh['primitive_idx']}.{output_format}" + + if output_format == 'obj': + self._save_as_obj(mesh, filename) + elif output_format == 'ply': + self._save_as_ply(mesh, filename) + elif output_format == 'npz': + self._save_as_npz(mesh, filename) + else: + print(f"不支持的格式: {output_format}") + continue + + print(f"已保存: {filename}") + + def _save_as_obj(self, mesh, filename): + """保存为 OBJ 格式""" + with open(filename, 'w') as f: + # 写入顶点 + if mesh['vertices'] is not None: + for v in mesh['vertices']: + f.write(f"v {v[0]} {v[1]} {v[2]}\n") + + # 写入纹理坐标 + if mesh['texcoords'] is not None: + for uv in mesh['texcoords']: + f.write(f"vt {uv[0]} {uv[1]}\n") + + # 写入法线 + if mesh['normals'] is not None: + for n in mesh['normals']: + f.write(f"vn {n[0]} {n[1]} {n[2]}\n") + + # 写入面 + if mesh['faces'] is not None: + for face in mesh['faces']: + # OBJ 索引从1开始 + face_indices = [str(idx + 1) for idx in face] + f.write(f"f {' '.join(face_indices)}\n") + + def _save_as_ply(self, mesh, filename): + """保存为 PLY 格式""" + from plyfile import PlyData, PlyElement + import numpy as np + + vertices = mesh['vertices'] + faces = mesh['faces'] + + if vertices is None: + return + + # 创建顶点数据 + vertex_data = np.zeros(len(vertices), dtype=[ + ('x', 'f4'), ('y', 'f4'), ('z', 'f4') + ]) + + vertex_data['x'] = vertices[:, 0] + vertex_data['y'] = vertices[:, 1] + vertex_data['z'] = vertices[:, 2] + + # 创建面数据 + if faces is not None: + face_data = np.zeros(len(faces), dtype=[('vertex_indices', 'i4', (3,))]) + face_data['vertex_indices'] = faces + + # 写入文件 + vertex_element = PlyElement.describe(vertex_data, 'vertex') + + if faces is not None: + face_element = PlyElement.describe(face_data, 'face') + PlyData([vertex_element, face_element], text=False).write(filename) + else: + PlyData([vertex_element], text=False).write(filename) + + def _save_as_npz(self, mesh, filename): + """保存为 NPZ 格式""" + np.savez( + filename, + vertices=mesh['vertices'], + faces=mesh['faces'], + texcoords=mesh['texcoords'], + batch_ids=mesh['batch_ids'], + normals=mesh['normals'], + colors=mesh['colors'] + ) + +# 使用示例 +def main(): + # 初始化解析器 + parser = DracoGLBParser(r"D:\devForBdzlWork\ai_project_v1\b3dm\test\temp_glb\temp_6e895637.glb") + + # 解析 GLB 结构 + parser.parse_glb_structure() + + # 分析结构 + parser.analyze_structure() + + # 解码 Draco 网格 + meshes = parser.decode_draco_meshes() + + # 使用新增的顶点获取方法 + print("\n" + "=" * 60) + print("顶点获取方法演示:") + print("=" * 60) + + # 1. 获取第一个网格的第一个图元的顶点 + vertices = parser.get_vertices(mesh_idx=0, primitive_idx=0) + if vertices is not None: + print(f"1. 获取第一个网格顶点:") + print(f" 形状: {vertices.shape}") + print(f" 数据类型: {vertices.dtype}") + print(f" 前5个顶点: \n{vertices[:5]}") + + # 2. 获取所有顶点（合并） + all_vertices = parser.get_all_vertices() + if all_vertices is not None: + print(f"\n2. 获取所有网格顶点（合并）:") + print(f" 总顶点数: {len(all_vertices)}") + print(f" 形状: {all_vertices.shape}") + + # 3. 获取总顶点数 + total_vertices = parser.get_vertex_count() + print(f"\n3. 总顶点数: {total_vertices}") + + # 4. 根据网格名称获取顶点 + if meshes: + mesh_name = meshes[0]['name'] + vertices_list = parser.get_vertices_by_mesh_name(mesh_name) + print(f"\n4. 根据名称 '{mesh_name}' 获取的顶点:") + for i, verts in enumerate(vertices_list): + print(f" 图元 {i}: {verts.shape if verts is not None else 'None'}") + + # 保存解码后的网格 + parser.save_decoded_meshes(meshes, output_format='obj') + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/b3dm/terrain3d_analyzer_color.py b/b3dm/terrain3d_analyzer_color.py new file mode 100644 index 0000000..6dfa19c --- /dev/null +++ b/b3dm/terrain3d_analyzer_color.py @@ -0,0 +1,380 @@ +# terrain3d_analyzer_color13 +import numpy as np +import matplotlib.pyplot as plt +from scipy.ndimage import sobel +import matplotlib as mpl +import rasterio +from osgeo import gdal +from matplotlib.patches import FancyArrowPatch +from mpl_toolkits.mplot3d import proj3d +import os + +# 自定义3D箭头类 +class Arrow3D(FancyArrowPatch): + def __init__(self, xs, ys, zs, *args, **kwargs): + super().__init__((0,0), (0,0), *args, **kwargs) + self._verts3d = xs, ys, zs + + def do_3d_projection(self, renderer=None): + xs3d, ys3d, zs3d = self._verts3d + xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, self.axes.M) + self.set_positions((xs[0], ys[0]), (xs[1], ys[1])) + return min(zs) + +def read_dem_rasterio(filepath): + """使用rasterio读取DEM数据""" + try: + with rasterio.open(filepath) as src: + dem_data = src.read(1) + transform = src.transform + bounds = src.bounds + crs = src.crs + + print(f"DEM信息:") + print(f" 尺寸: {dem_data.shape}") + print(f" 范围: {bounds}") + print(f" 坐标系: {crs}") + print(f" 高程范围: {np.nanmin(dem_data):.2f} - {np.nanmax(dem_data):.2f}") + + if np.isnan(dem_data).any(): + print(" 检测到NaN值") + dem_data = np.nan_to_num(dem_data, nan=np.nanmean(dem_data)) + + rows, cols = dem_data.shape + x = np.linspace(bounds.left, bounds.right, cols) + y = np.linspace(bounds.bottom, bounds.top, rows) + X, Y = np.meshgrid(x, y) + + return X, Y, dem_data + + except Exception as e: + print(f"使用rasterio读取失败: {e}") + return None + +def read_slope_aspect_by_dem(dem_path, overall_3d_output_path=None) : + # ---------------------- 核心：解决中文乱码配置 ---------------------- + setup_chinese_font() + + # 尝试读取DEM数据 + print("正在读取DEM文件...") + X, Y, Z = read_dem_rasterio(dem_path) + + if X is None: + raise FileNotFoundError(f"无法读取DEM文件: {dem_path}") + + # 检查数据有效性 + if Z.size == 0 or np.all(Z == 0): + raise ValueError("DEM数据无效或全部为0") + + # 重采样 + if Z.shape[0] > 200 or Z.shape[1] > 200: + print(f"原始DEM尺寸较大 ({Z.shape})，进行重采样...") + from scipy.ndimage import zoom + scale_factor = min(200/Z.shape[0], 200/Z.shape[1]) + Z = zoom(Z, scale_factor, order=1) + X = zoom(X, scale_factor, order=1) + Y = zoom(Y, scale_factor, order=1) + print(f"重采样后尺寸: {Z.shape}") + + # ---------------------- 2. 计算坡度和坡向 ---------------------- + print("计算坡度和坡向...") + + dx_pixel = np.abs(X[0,1] - X[0,0]) * 111000 + dy_pixel = np.abs(Y[1,0] - Y[0,0]) * 111000 + + dx = sobel(Z, axis=1) / (2 * dx_pixel) + dy = sobel(Z, axis=0) / (2 * dy_pixel) + + slope_rad = np.arctan(np.sqrt(dx**2 + dy**2)) + slope_deg = slope_rad * (180 / np.pi) + + aspect_rad = np.arctan2(-dx, dy) + aspect_deg = aspect_rad * (180 / np.pi) + aspect_deg[aspect_deg < 0] += 360 + + print(f"坡度范围: {np.min(slope_deg):.2f}° - {np.max(slope_deg):.2f}°") + print(f"坡向范围: {np.min(aspect_deg):.2f}° - {np.max(aspect_deg):.2f}°") + + # ---------------------- 3. 计算整体坡向 ---------------------- + print("\n计算整体坡向...") + + def calculate_overall_aspect(aspect_deg, slope_deg, method='weighted_mean'): + """计算整体坡向""" + aspect_rad = np.deg2rad(aspect_deg) + + if method == 'weighted_mean': + # 坡度加权平均法 + u = np.sin(aspect_rad) + v = np.cos(aspect_rad) + + weights = slope_deg.flatten() + weighted_u = np.nansum(u.flatten() * weights) / np.nansum(weights) + weighted_v = np.nansum(v.flatten() * weights) / np.nansum(weights) + + weighted_aspect_rad = np.arctan2(weighted_u, weighted_v) + weighted_aspect_deg = np.rad2deg(weighted_aspect_rad) + + if weighted_aspect_deg < 0: + weighted_aspect_deg += 360 + + weighted_strength = np.sqrt(weighted_u**2 + weighted_v**2) + return weighted_aspect_deg, weighted_strength, "坡度加权平均法" + + elif method == 'vector_mean': + # 向量平均法 + u = np.sin(aspect_rad) + v = np.cos(aspect_rad) + + mean_u = np.nanmean(u) + mean_v = np.nanmean(v) + + mean_aspect_rad = np.arctan2(mean_u, mean_v) + mean_aspect_deg = np.rad2deg(mean_aspect_rad) + + if mean_aspect_deg < 0: + mean_aspect_deg += 360 + + vector_strength = np.sqrt(mean_u**2 + mean_v**2) + return mean_aspect_deg, vector_strength, "向量平均法" + + # 使用坡度加权平均法计算整体坡向 + overall_aspect, overall_strength, method_name = calculate_overall_aspect(aspect_deg, slope_deg, 'weighted_mean') + + # 将整体坡向转换为方向描述 + if overall_aspect < 22.5 or overall_aspect >= 337.5: + overall_direction = "北" + elif 22.5 <= overall_aspect < 67.5: + overall_direction = "东北" + elif 67.5 <= overall_aspect < 112.5: + overall_direction = "东" + elif 112.5 <= overall_aspect < 157.5: + overall_direction = "东南" + elif 157.5 <= overall_aspect < 202.5: + overall_direction = "南" + elif 202.5 <= overall_aspect < 247.5: + overall_direction = "西南" + elif 247.5 <= overall_aspect < 292.5: + overall_direction = "西" + else: + overall_direction = "西北" + + print(f"整体坡向 ({method_name}):") + print(f" 角度: {overall_aspect:.1f}°") + print(f" 方向: {overall_direction}") + # print(f" 一致性: {overall_strength:.3f}") + + # ---------------------- 5. 3D可视化（俯视图，包含整体坡向和关键点坡向） ---------------------- + print("\n生成3D俯视图可视化...") + fig3d, ax3d = plt.subplots(figsize=(16, 12), subplot_kw={"projection": "3d"}) + + # 绘制地形曲面 - 俯视图需要更清晰的地形表现 + norm = mpl.colors.Normalize(vmin=np.percentile(slope_deg, 5), + vmax=np.percentile(slope_deg, 95)) + + plot_skip = max(2, Z.shape[0] // 60) # 增加采样密度，使俯视图更清晰 + X_plot = X[::plot_skip, ::plot_skip] + Y_plot = Y[::plot_skip, ::plot_skip] + Z_plot = Z[::plot_skip, ::plot_skip] + slope_plot = slope_deg[::plot_skip, ::plot_skip] + + surf = ax3d.plot_surface( + X_plot, Y_plot, Z_plot, + cmap="viridis_r", + alpha=0.85, # 增加透明度，使箭头更明显 + linewidth=0.1, # 很细的网格线 + facecolors=plt.cm.viridis_r(norm(slope_plot)), + zorder=1 + ) + + # ---------------------- 绘制整体坡向箭头（中心位置，红色粗箭头） ---------------------- + center_x = (X.min() + X.max()) / 2 + center_y = (Y.min() + Y.max()) / 2 + + # 整体坡向箭头长度 + arrow_length_overall = 0.15 * min(X.max()-X.min(), Y.max()-Y.min()) + + # 计算整体坡向箭头方向 + scale_factor = 1.5 + overall_aspect_rad = np.deg2rad(overall_aspect) + dx_overall = np.sin(overall_aspect_rad) * arrow_length_overall * scale_factor + dy_overall = np.cos(overall_aspect_rad) * arrow_length_overall * scale_factor + + # 方案2：如果希望箭头在地形上方一定高度 + terrain_max_z = Z.max() # 地形最高点 + float_height = 0.2 * terrain_max_z # 在地形最高点上方20%的高度 + + # 绘制整体坡向箭头（红色，粗） + arrow_overall = Arrow3D( + [center_x, center_x + dx_overall], + [center_y, center_y + dy_overall], + [terrain_max_z + float_height, terrain_max_z + float_height], # 在地形上方 + mutation_scale=25, # 稍大一些的箭头 + lw=5, # 更粗的线 + arrowstyle='-|>', + color='red', + alpha=0.98, # 更高的透明度 + zorder=20 # 更高的绘制顺序 + ) + ax3d.add_artist(arrow_overall) + + # 在整体坡向箭头起点添加标记 + ax3d.scatter(center_x, center_y, terrain_max_z + float_height, + color='red', s=120, edgecolor='white', linewidth=2, zorder=21) + + # ---------------------- 整体坡向面板放置在左上角 ---------------------- + # 计算左上角位置 + panel_x = X.min() + 0.02 * (X.max() - X.min()) # 左边留2%的边距 + panel_y = Y.max() - 0.02 * (Y.max() - Y.min()) # 上边留2%的边距 + panel_z = Z.max() + 0.5 * (Z.max() - Z.min()) # 进一步提高Z坐标，确保在视野内 + + # 整体坡向面板内容 + overall_info = ( + f"整体坡向分析\n" + f"角度: {overall_aspect:.1f}°\n" + f"方向: {overall_direction}" + # f"一致性: {overall_strength:.3f}" + ) + + # 绘制整体坡向面板（左上角） + ax3d.text(panel_x, panel_y, panel_z, + overall_info, + fontsize=12, fontweight='bold', + bbox=dict(facecolor='white', alpha=0.5, boxstyle="round,pad=0.5", + edgecolor='red', linewidth=2.5), + ha='left', va='top', zorder=30) + + # ---------------------- 设置3D图参数 ---------------------- + ax3d.set_xlabel("经度 (X)", fontsize=12) + ax3d.set_ylabel("纬度 (Y)", fontsize=12) + ax3d.set_zlabel("高程 (m)", fontsize=12) + + # 获取文件名用于标题 + filename = os.path.basename(dem_path) + ax3d.set_title(f"DEM三维俯视图 - 坡向分析\n文件: {filename}", + fontsize=14, fontweight='bold', pad=20) + + # 添加颜色条 + cbar = plt.colorbar( + mpl.cm.ScalarMappable(norm=norm, cmap='viridis_r'), + ax=ax3d, shrink=0.6, aspect=25, pad=0.15 + ) + cbar.set_label("坡度 (°)", fontsize=12) + + # ---------------------- 设置为俯视图（高仰角） ---------------------- + view_elev = 85 # 接近90度的仰角，俯视效果 + view_azim = overall_aspect + 180 # 从坡向的相反方向观看，可以看到坡面 + + # 确保方位角在0-360度范围内 + view_azim = view_azim % 360 + + print(f"\n设置3D俯视图视角:") + print(f" 整体坡向: {overall_aspect:.1f}° ({overall_direction})") + print(f" 视角方位角: {view_azim:.1f}°") + print(f" 视角仰角: {view_elev:.1f}°") + + # 应用俯视图视角设置 + ax3d.view_init(elev=view_elev, azim=view_azim) + + # 调整相机距离，使视角更广 + ax3d.dist = 9.0 # 增加相机距离 + + # 设置坐标轴范围，确保所有元素都显示 + z_min, z_max = Z.min(), Z.max() + z_padding = 0.6 * (z_max - z_min) # 适当增加Z轴范围 + ax3d.set_zlim(z_min - 0.1*z_padding, z_max + z_padding) + + # 设置XY轴范围 + x_margin = 0.1 * (X.max() - X.min()) + y_margin = 0.1 * (Y.max() - Y.min()) + ax3d.set_xlim(X.min() - x_margin, X.max() + x_margin) + ax3d.set_ylim(Y.min() - y_margin, Y.max() + y_margin) + + plt.tight_layout() + + # 保存图片 + if not overall_3d_output_path : + overall_3d_output_path = dem_path.replace('.tif', '_slopeAspect_3D_overlook.png') + plt.savefig(overall_3d_output_path, dpi=250, bbox_inches='tight', facecolor='white') + # plt.show() + + os.remove(dem_path) + print(f"\n3D俯视图已保存: {overall_3d_output_path}") + print("分析完成！") + print("="*60) + + return overall_3d_output_path + +# ---------------------- 跨平台中文字体配置 ---------------------- +def setup_chinese_font(): + """设置中文字体支持，兼容Windows、Linux、macOS""" + import platform + + plt.rcParams['axes.unicode_minus'] = False # 正确显示负号 + + system = platform.system() + + # 尝试添加系统中文字体路径 + font_paths = [] + + if system == 'Windows': + # Windows 字体路径 + font_paths = [ + 'C:/Windows/Fonts/simhei.ttf', # 黑体 + 'C:/Windows/Fonts/simkai.ttf', # 楷体 + 'C:/Windows/Fonts/simsun.ttc', # 宋体 + 'C:/Windows/Fonts/microsoftyahei.ttf', # 微软雅黑 + ] + elif system == 'Darwin': # macOS + # macOS 字体路径 + font_paths = [ + '/System/Library/Fonts/PingFang.ttc', # 苹方 + '/System/Library/Fonts/STHeiti Light.ttc', # 华文黑体 + '/System/Library/Fonts/STHeiti Medium.ttc', + '/Library/Fonts/Arial Unicode.ttf', # Arial Unicode + ] + elif system == 'Linux': + # Linux 字体路径 + font_paths = [ + '/usr/share/fonts/truetype/wqy/wqy-microhei.ttc', # 文泉驿微米黑 + '/usr/share/fonts/truetype/wqy/wqy-zenhei.ttc', # 文泉驿正黑 + '/usr/share/fonts/truetype/arphic/uming.ttc', # AR PL UMing + '/usr/share/fonts/opentype/noto/NotoSansCJK-Regular.ttc', # Noto + ] + + # 尝试找到并添加第一个可用的中文字体 + font_added = False + for font_path in font_paths: + try: + if os.path.exists(font_path): + font_prop = mpl.font_manager.FontProperties(fname=font_path) + font_name = font_prop.get_name() + mpl.font_manager.fontManager.addfont(font_path) + plt.rcParams['font.sans-serif'] = [font_name] + plt.rcParams['font.sans-serif'] + font_added = True + print(f"已添加中文字体: {font_name} ({font_path})") + break + except Exception as e: + print(f"添加字体 {font_path} 失败: {e}") + continue + + # 如果以上方法都失败，使用通用备选方案 + if not font_added: + if system == 'Windows': + fallback_fonts = ['SimHei', 'Microsoft YaHei', 'KaiTi', 'FangSong'] + elif system == 'Darwin': + fallback_fonts = ['PingFang SC', 'Hiragino Sans GB', 'Apple LiGothic Medium'] + else: # Linux and others + fallback_fonts = ['DejaVu Sans', 'WenQuanYi Micro Hei', + 'Noto Sans CJK SC', 'Heiti TC', 'AR PL UMing CN'] + + current_fonts = plt.rcParams.get('font.sans-serif', []) + plt.rcParams['font.sans-serif'] = fallback_fonts + current_fonts + print(f"使用备选字体方案: {fallback_fonts[:2]}...") + + # 设置字体族 + plt.rcParams['font.family'] = 'sans-serif' + + +if __name__ == "__main__": + dem_path = r'D:\devForBdzlWork\ai_project_v1\b3dm\test\o_dem_df67f1cc.tif' + read_slope_aspect_by_dem(dem_path) \ No newline at end of file diff --git a/b3dm/terrain_api.py b/b3dm/terrain_api.py new file mode 100644 index 0000000..b22e906 --- /dev/null +++ b/b3dm/terrain_api.py @@ -0,0 +1,342 @@ +from sanic import Blueprint, Request, json +from pydantic import BaseModel, Field, field_validator +from typing import List, Optional, Dict, Any +import numpy as np +import logging +import time +import uuid +import urllib.parse +import threading +from terrain_calculator import TerrainCalculator + +terrain_bp = Blueprint("terrain", url_prefix="") +MINIO_SUB_PATH = "slopeAspectPng" + +# 配置日志 +logging.basicConfig(level=logging.INFO) +logger = logging.getLogger(__name__) + +# 请求模型 +class NormalVector(BaseModel): + """法向量模型""" + nx: float = Field(..., description="法向量X分量") + ny: float = Field(..., description="法向量Y分量") + nz: float = Field(..., description="法向量Z分量") + + @field_validator('nx', 'ny', 'nz') + def check_finite(cls, v): + if not np.isfinite(v): + raise ValueError(f"值必须是有限数字，得到: {v}") + return v + + def to_list(self): + return [self.nx, self.ny, self.nz] + +class BatchRequest(BaseModel): + """批量请求模型""" + vectors: List[List[float]] = Field(..., description="法向量列表") + + @field_validator('vectors') + def validate_vectors(cls, v): + if len(v) > 1000: + raise ValueError("批量处理最多支持1000个向量") + for vec in v: + if len(vec) != 3: + raise ValueError("每个向量必须是长度为3的列表") + if not all(isinstance(x, (int, float)) for x in vec): + raise ValueError("向量元素必须是数字") + return v + +class PointItem(BaseModel): + """单个点模型""" + x: float = Field(..., description="x坐标") + y: float = Field(..., description="y坐标") + z: float = Field(..., description="z坐标") + +class PointRequest(BaseModel): + """使用Pydantic HttpUrl类型的严格点批量请求模型""" + points: List[PointItem] = Field(..., description="点列表") + url: str = Field(..., description="URL地址") + + @field_validator('points') + def validate_points_count(cls, v): + if len(v) > 10: + raise ValueError("批量处理最多支持10个点") + return v + +class AnalysisConfig(BaseModel): + """分析配置""" + classify: bool = Field(default=True, description="是否进行分类") + include_percent: bool = Field(default=True, description="是否包含坡度百分比") + include_direction: bool = Field(default=True, description="是否包含方向描述") + +# 中间件：请求计时 +@terrain_bp.middleware("request") +async def add_start_time(request: Request): + request.ctx.start_time = time.time() + +@terrain_bp.middleware("response") +async def add_response_time(request: Request, response): + if hasattr(request.ctx, "start_time"): + process_time = (time.time() - request.ctx.start_time) * 1000 + response.headers["X-Process-Time"] = f"{process_time:.2f}ms" + +@terrain_bp.post("/api/v1/calculate/slope") +async def calculate_slope(request: Request): + """计算坡度""" + try: + data = request.json + if not data: + return json({"error": "请求体不能为空"}, status=400) + + # 验证输入 + try: + vector = NormalVector(**data) + except Exception as e: + return json({"error": f"输入验证失败: {str(e)}"}, status=400) + + # 计算坡度 + result = TerrainCalculator.calculate_slope(vector.to_list()) + + # 检查是否有错误 + if "error" in result and result["error"]: + return json(result, status=400) + + return json({ + "success": True, + "data": result, + "request": { + "input_vector": vector.to_list(), + "timestamp": time.time() + } + }) + + except Exception as e: + logger.error(f"坡度计算API错误: {e}") + return json({ + "error": f"服务器内部错误: {str(e)}", + "success": False + }, status=500) + +@terrain_bp.post("/api/v1/calculate/aspect") +async def calculate_aspect1(request: Request): + """计算坡向""" + try: + data = request.json + if not data: + return json({"error": "请求体不能为空"}, status=400) + + # 验证输入 + try: + vector = NormalVector(**data) + except Exception as e: + return json({"error": f"输入验证失败: {str(e)}"}, status=400) + + # 计算坡向 + result = TerrainCalculator.calculate_aspect(vector.to_list()) + + # 检查是否有错误 + if "error" in result and result["error"]: + return json(result, status=400) + + return json({ + "success": True, + "data": result, + "request": { + "input_vector": vector.to_list(), + "timestamp": time.time() + } + }) + + except Exception as e: + logger.error(f"坡向计算API错误: {e}") + return json({ + "error": f"服务器内部错误: {str(e)}", + "success": False + }, status=500) + +@terrain_bp.post("/api/v1/calculate/slopeAspect") +async def calculate_slopeAspect(request: Request): + """生成坡向坡度俯视图""" + try: + data = request.json + if not data: + return json({"error": "请求体不能为空"}, status=400) + + # 验证输入 + try: + vector = PointRequest(**data) + except Exception as e: + return json({"error": f"输入验证失败: {str(e)}"}, status=400) + + # 生成坡向坡度俯视图 + region_coords = [(point.x, point.y, point.z) for point in vector.points] + overall_3d_png_name = f"o_dem_{uuid.uuid4().hex[:8]}_slopeAspect.png" + # 创建并启动线程 + thread1 = threading.Thread(target=TerrainCalculator.generate_slopeAspect_3d_overlook, args=(region_coords, vector.url, overall_3d_png_name, MINIO_SUB_PATH)) + # 启动线程 + thread1.start() + url_prefix = extract_and_rebuild_url(vector.url) + return json({ + "success": True, + "data": f"{url_prefix}/{MINIO_SUB_PATH}/{overall_3d_png_name}", + "request": { + "input_vector": vector.model_dump(), + "timestamp": time.time() + } + }) + + except Exception as e: + logger.error(f"坡向计算API错误: {e}") + return json({ + "error": f"服务器内部错误: {str(e)}", + "success": False + }, status=500) + +@terrain_bp.post("/api/v1/calculate/both") +async def calculate_both(request: Request): + """同时计算坡度和坡向""" + try: + data = request.json + if not data: + return json({"error": "请求体不能为空"}, status=400) + + # 验证输入 + try: + vector = NormalVector(**data) + except Exception as e: + return json({"error": f"输入验证失败: {str(e)}"}, status=400) + + # 计算坡度和坡向 + result = TerrainCalculator.calculate_slope_aspect(vector.to_list()) + + # 检查是否有错误 + if "error" in result and result["error"]: + return json(result, status=400) + + return json({ + "success": True, + "data": result, + "request": { + "input_vector": vector.to_list(), + "timestamp": time.time() + } + }) + + except Exception as e: + logger.error(f"综合计算API错误: {e}") + return json({ + "error": f"服务器内部错误: {str(e)}", + "success": False + }, status=500) + +@terrain_bp.post("/api/v1/calculate/batch") +async def batch_calculate(request: Request): + """批量计算""" + try: + data = request.json + if not data: + return json({"error": "请求体不能为空"}, status=400) + + # 验证输入 + try: + batch_request = BatchRequest(**data) + except Exception as e: + return json({"error": f"输入验证失败: {str(e)}"}, status=400) + + # 批量计算 + start_time = time.time() + result = TerrainCalculator.batch_calculate(batch_request.vectors) + process_time = (time.time() - start_time) * 1000 + + if "error" in result and result["error"]: + return json(result, status=400) + + return json({ + "success": True, + "data": result, + "performance": { + "process_time_ms": process_time, + "vectors_per_second": len(batch_request.vectors) / (process_time / 1000) if process_time > 0 else 0 + }, + "request": { + "vector_count": len(batch_request.vectors), + "timestamp": time.time() + } + }) + + except Exception as e: + logger.error(f"批量计算API错误: {e}") + return json({ + "error": f"服务器内部错误: {str(e)}", + "success": False + }, status=500) + +@terrain_bp.get("/api/v1/example") +async def get_examples(request: Request): + """获取示例数据""" + examples = { + "flat": { + "nx": 0.0, + "ny": 0.0, + "nz": 1.0, + "expected_slope": 0.0, + "description": "完全水平面" + }, + "north_slope_30": { + "nx": 0.0, + "ny": -0.5, + "nz": 0.8660254, + "expected_slope": 30.0, + "expected_aspect": 0.0, + "description": "朝北30度斜坡" + }, + "east_slope_45": { + "nx": 0.7071068, + "ny": 0.0, + "nz": 0.7071068, + "expected_slope": 45.0, + "expected_aspect": 90.0, + "description": "朝东45度斜坡" + }, + "vertical": { + "nx": 1.0, + "ny": 0.0, + "nz": 0.0, + "expected_slope": 90.0, + "description": "垂直面" + } + } + + return json({ + "examples": examples, + "count": len(examples) + }) + +def extract_and_rebuild_url(url): + """提取URL的三部分并重建""" + # 解析URL + parsed = urllib.parse.urlparse(url) + + # 1. 提取协议部分 (http/https) + scheme = parsed.scheme or "http" # 如果没有协议，默认用http + + # 2. 提取IP端口/主机部分 + netloc = parsed.netloc + + # 3. 提取第一个路径分段 + path = parsed.path.strip("/") # 去掉首尾的斜杠 + path_parts = path.split("/") + + if path_parts and path_parts[0]: + first_segment = path_parts[0] + else: + first_segment = "" + + # 重建URL + if first_segment: + rebuilt_url = f"{scheme}://{netloc}/{first_segment}" + else: + rebuilt_url = f"{scheme}://{netloc}" + + return rebuilt_url \ No newline at end of file diff --git a/b3dm/terrain_calculator.py b/b3dm/terrain_calculator.py new file mode 100644 index 0000000..5f45082 --- /dev/null +++ b/b3dm/terrain_calculator.py @@ -0,0 +1,361 @@ +import numpy as np +from typing import List, Tuple, Dict, Any +import logging +import os +import uuid +from data_3dtiles_manager import MinIO3DTilesManager +import data_3dtiles_to_dem +import terrain3d_analyzer_color + +logger = logging.getLogger(__name__) + +class TerrainCalculator: + """地形坡度和坡向计算器""" + + def generate_slopeAspect_3d_overlook(region_coords, url, overall_3d_png_name, minio_sub_path) : + # 下载3dtiles地图数据 + ENDPOINT_URL = "222.212.85.86:9000" + ACCESS_KEY = "WuRenJi" + SECRET_KEY = "WRJ@2024" + manager = MinIO3DTilesManager( + endpoint_url=ENDPOINT_URL, + access_key=ACCESS_KEY, + secret_key=SECRET_KEY, + secure=False + ) + script_dir = os.path.dirname(os.path.abspath(__file__)) + success, entry_local_path = manager.download_full_tileset( + tileset_url=url, + save_dir=f"{script_dir}/data_3dtiles", + region_filter=None + ) + if not success : + logger.info(f"下载地图数据失败: {url},{region_coords}") + return "下载地图数据失败", None + + tileset_path = entry_local_path + dem_path = os.path.join(script_dir, f"o_dem_{uuid.uuid4().hex[:8]}.tif") + data_3dtiles_to_dem.generate_dem(region_coords, tileset_path, dem_path) + + if not os.path.exists(dem_path) : + logger.info(f"生成坡度坡向俯视图失败: {url},{region_coords}") + return "生成坡度坡向俯视图失败", None + + overall_3d_png_path = os.path.join(script_dir, overall_3d_png_name) + terrain3d_analyzer_color.read_slope_aspect_by_dem(dem_path, overall_3d_png_path) + logger.info(f"生成成功: {url},{region_coords},{overall_3d_png_path}") + + entry_bucket, _ = manager.parse_minio_url(url); + success, minio_path = manager.upload_file(entry_bucket, f"{minio_sub_path}/{overall_3d_png_name}", overall_3d_png_path) + if success : + return "生成成功", minio_path + else : + return "生成失败", None + + @staticmethod + def validate_vector(vector: List[float]) -> bool: + """验证输入向量是否有效""" + if len(vector) != 3: + return False + if not all(isinstance(v, (int, float)) for v in vector): + return False + norm = np.linalg.norm(vector) + return norm > 1e-10 # 避免零向量 + + @staticmethod + def normalize_vector(vector: List[float]) -> np.ndarray: + """向量归一化""" + arr = np.array(vector, dtype=np.float64) + norm = np.linalg.norm(arr) + return arr / norm if norm > 0 else arr + + @staticmethod + def calculate_slope(normal_vector: List[float]) -> Dict[str, Any]: + """ + 计算坡度 + + Args: + normal_vector: 法向量 [nx, ny, nz] + + Returns: + dict: 包含坡度(度)和相关信息 + """ + try: + # 验证输入 + if not TerrainCalculator.validate_vector(normal_vector): + return { + "error": "无效的法向量,必须是长度为3的数值列表且不能为零向量", + "slope_deg": None + } + + # 归一化 + n = TerrainCalculator.normalize_vector(normal_vector) + + # 计算坡度（使用arccos法） + nz_abs = abs(n[2]) + + # 处理数值误差 + if nz_abs > 1.0: + nz_abs = 1.0 + elif nz_abs < 0.0: + nz_abs = 0.0 + + # 计算坡度（弧度） + if abs(nz_abs - 1.0) < 1e-10: # 完全水平 + slope_rad = 0.0 + elif abs(nz_abs) < 1e-10: # 完全垂直 + slope_rad = np.pi / 2 + else: + slope_rad = np.arccos(nz_abs) + + # 转换为度 + slope_deg = np.degrees(slope_rad) + + # 计算坡度百分比 + slope_percent = np.tan(slope_rad) * 100 if slope_rad < np.pi/2 else float('inf') + + return { + "slope_deg": float(slope_deg), + "slope_rad": float(slope_rad), + "slope_percent": float(slope_percent), + "normalized_vector": n.tolist(), + "classification": TerrainCalculator.classify_slope(slope_deg) + } + + except Exception as e: + logger.error(f"坡度计算错误: {e}") + return { + "error": f"计算失败: {str(e)}", + "slope_deg": None + } + + @staticmethod + def calculate_aspect(normal_vector: List[float]) -> Dict[str, Any]: + """ + 计算坡向 + + Args: + normal_vector: 法向量 [nx, ny, nz] + + Returns: + dict: 包含坡向(度)和相关信息 + """ + try: + # 验证输入 + if not TerrainCalculator.validate_vector(normal_vector): + return { + "error": "无效的法向量,必须是长度为3的数值列表且不能为零向量", + "aspect_deg": None + } + + # 归一化 + n = TerrainCalculator.normalize_vector(normal_vector) + + # 检查是否为水平面 + nx, ny, nz = n + horizontal_magnitude = np.sqrt(nx*nx + ny*ny) + + if horizontal_magnitude < 1e-10: # 水平面,坡向无定义 + return { + "aspect_deg": None, + "aspect_rad": None, + "is_flat": True, + "message": "水平面,坡向无定义", + "normalized_vector": n.tolist() + } + + # 计算原始坡向（四象限反正切） + # 注意：arctan2(nx, ny) 不是 arctan2(ny, nx) + raw_angle_rad = np.arctan2(nx, ny) + + # 转换为坡向（下坡方向 = 法向量方向 + 180°） + aspect_rad = raw_angle_rad + np.pi + + # 转换为度 + aspect_deg = np.degrees(aspect_rad) + + # 归一化到 [0, 360) 范围 + aspect_deg = aspect_deg % 360.0 + + # 转换为八方向 + direction = TerrainCalculator.aspect_to_direction(aspect_deg) + + return { + "aspect_deg": float(aspect_deg), + "aspect_rad": float(aspect_rad % (2*np.pi)), + "direction": direction, + "is_flat": False, + "normalized_vector": n.tolist(), + "raw_angle_deg": float(np.degrees(raw_angle_rad)) + } + + except Exception as e: + logger.error(f"坡向计算错误: {e}") + return { + "error": f"计算失败: {str(e)}", + "aspect_deg": None + } + + @staticmethod + def calculate_slope_aspect(normal_vector: List[float]) -> Dict[str, Any]: + """ + 同时计算坡度和坡向 + + Args: + normal_vector: 法向量 [nx, ny, nz] + + Returns: + dict: 包含坡度和坡向的综合结果 + """ + try: + slope_result = TerrainCalculator.calculate_slope(normal_vector) + aspect_result = TerrainCalculator.calculate_aspect(normal_vector) + + result = { + "slope": slope_result, + "aspect": aspect_result, + "input_vector": normal_vector, + "calculation_time": None # 可在调用处添加时间戳 + } + + # 如果有错误,合并错误信息 + errors = [] + if "error" in slope_result and slope_result["error"]: + errors.append(f"坡度: {slope_result['error']}") + if "error" in aspect_result and aspect_result["error"]: + errors.append(f"坡向: {aspect_result['error']}") + + if errors: + result["errors"] = errors + + return result + + except Exception as e: + logger.error(f"综合计算错误: {e}") + return { + "error": f"综合计算失败: {str(e)}" + } + + @staticmethod + def classify_slope(slope_deg: float) -> Dict[str, Any]: + """坡度分类""" + if slope_deg < 2: + return {"category": "平坦", "level": 0, "description": "基本平坦"} + elif slope_deg < 5: + return {"category": "缓坡", "level": 1, "description": "适合农业"} + elif slope_deg < 15: + return {"category": "斜坡", "level": 2, "description": "适合建设"} + elif slope_deg < 30: + return {"category": "陡坡", "level": 3, "description": "需要工程措施"} + elif slope_deg < 45: + return {"category": "急陡坡", "level": 4, "description": "高风险区域"} + else: + return {"category": "峭壁", "level": 5, "description": "危险区域"} + + @staticmethod + def aspect_to_direction(aspect_deg: float) -> Dict[str, Any]: + """将坡向转换为八方向""" + directions = ["北", "东北", "东", "东南", "南", "西南", "西", "西北"] + + # 计算方向索引 (45°一个区间) + index = int((aspect_deg + 22.5) % 360 / 45) + + return { + "chinese": directions[index], + "english": ["N", "NE", "E", "SE", "S", "SW", "W", "NW"][index], + "degree_range": { + "min": (index * 45 - 22.5) % 360, + "max": (index * 45 + 22.5) % 360 + } + } + + @staticmethod + def batch_calculate(vectors: List[List[float]]) -> Dict[str, Any]: + """批量计算多个法向量""" + try: + results = [] + errors = [] + + for i, vec in enumerate(vectors): + try: + result = TerrainCalculator.calculate_slope_aspect(vec) + result["index"] = i + results.append(result) + except Exception as e: + errors.append({ + "index": i, + "vector": vec, + "error": str(e) + }) + + return { + "total": len(vectors), + "successful": len(results), + "failed": len(errors), + "results": results, + "errors": errors, + "statistics": TerrainCalculator.calculate_statistics(results) + } + + except Exception as e: + logger.error(f"批量计算错误: {e}") + return {"error": f"批量计算失败: {str(e)}"} + + @staticmethod + def calculate_statistics(results: List[Dict]) -> Dict[str, Any]: + """计算统计信息""" + if not results: + return {} + + slope_values = [] + aspect_values = [] + + for r in results: + if "slope" in r and "slope_deg" in r["slope"] and r["slope"]["slope_deg"] is not None: + slope_values.append(r["slope"]["slope_deg"]) + if "aspect" in r and "aspect_deg" in r["aspect"] and r["aspect"]["aspect_deg"] is not None: + aspect_values.append(r["aspect"]["aspect_deg"]) + + if slope_values: + slope_arr = np.array(slope_values) + aspect_arr = np.array(aspect_values) if aspect_values else np.array([]) + + stats = { + "slope": { + "count": len(slope_values), + "mean": float(np.mean(slope_arr)), + "std": float(np.std(slope_arr)), + "min": float(np.min(slope_arr)), + "max": float(np.max(slope_arr)), + "median": float(np.median(slope_arr)) + } + } + + if aspect_values: + # 坡向统计需要循环统计 + stats["aspect"] = { + "count": len(aspect_values), + "mean_vector": TerrainCalculator.circular_mean(aspect_arr), + "concentration": TerrainCalculator.circular_concentration(aspect_arr) + } + + return stats + + return {} + + @staticmethod + def circular_mean(angles_deg: np.ndarray) -> float: + """计算循环数据的平均值（角度）""" + angles_rad = np.radians(angles_deg) + x = np.mean(np.cos(angles_rad)) + y = np.mean(np.sin(angles_rad)) + mean_rad = np.arctan2(y, x) + return np.degrees(mean_rad) % 360 + + @staticmethod + def circular_concentration(angles_deg: np.ndarray) -> float: + """计算角度数据的集中度 (0-1)""" + angles_rad = np.radians(angles_deg) + x = np.mean(np.cos(angles_rad)) + y = np.mean(np.sin(angles_rad)) + return np.sqrt(x*x + y*y) \ No newline at end of file diff --git a/b3dm/tileset_data_source.py b/b3dm/tileset_data_source.py new file mode 100644 index 0000000..bcfcb19 --- /dev/null +++ b/b3dm/tileset_data_source.py @@ -0,0 +1,167 @@ +# tileset_data_source_py +import numpy as np +from typing import List, Dict, Optional, Tuple +import asyncio +from concurrent.futures import ThreadPoolExecutor +import logging +import os + +logger = logging.getLogger(__name__) + +class TilesetDataSource: + """使用py3dtiles库的数据源""" + + def __init__(self, tileset_path: str, cache_size: int = 1000): + self.tileset_path = os.path.abspath(tileset_path) + self.tileset_dir = os.path.dirname(self.tileset_path) + self.cache_size = cache_size + self._tileset = None + self._point_cache = {} + self._executor = ThreadPoolExecutor(max_workers=4) + self._crs = "EPSG:4979" + + async def initialize(self): + """初始化""" + try: + # 尝试导入py3dtiles + try: + import py3dtiles + from py3dtiles.tileset import TileSet + except ImportError: + logger.warning("py3dtiles未安装，将使用简化数据源") + raise ImportError("请安装py3dtiles: pip install py3dtiles") + + loop = asyncio.get_event_loop() + self._tileset = await loop.run_in_executor( + self._executor, + TileSet.from_file, + self.tileset_path + ) + + logger.info(f"py3dtiles数据源初始化完成: {self.tileset_path}") + + except Exception as e: + logger.error(f"py3dtiles初始化失败: {str(e)}") + # 回退到简化数据源 + self._tileset = None + + async def get_points_in_polygon(self, + polygon_coords: List[List[float]], + z_range: Optional[Tuple[float, float]] = None) -> np.ndarray: + """获取点数据""" + if self._tileset is None: + # 使用简化数据源 + return await self._get_simulated_points(polygon_coords, z_range) + + try: + # 使用py3dtiles API获取数据 + points = [] + polygon_np = np.array(polygon_coords) + + # 遍历tileset中的所有tile + for tile in self._tileset.root_tile.traverse(): + tile_points = self._extract_tile_points(tile) + if tile_points.size > 0: + # 筛选多边形内的点 + points_in_polygon = self._filter_points_by_polygon(tile_points, polygon_np) + if points_in_polygon.size > 0: + points.append(points_in_polygon) + + if points: + all_points = np.vstack(points) + if z_range: + mask = (all_points[:, 2] >= z_range[0]) & (all_points[:, 2] <= z_range[1]) + all_points = all_points[mask] + return all_points + + return np.array([]) + + except Exception as e: + logger.error(f"py3dtiles获取数据失败: {str(e)}") + return await self._get_simulated_points(polygon_coords, z_range) + + def _extract_tile_points(self, tile) -> np.ndarray: + """从tile提取点数据""" + try: + if hasattr(tile, 'content') and tile.content: + # 尝试获取点数据 + if hasattr(tile.content, 'points'): + return tile.content.points.positions + elif hasattr(tile.content, 'body'): + # 处理其他格式 + return np.array([]) + return np.array([]) + except: + return np.array([]) + + def _filter_points_by_polygon(self, points: np.ndarray, polygon: np.ndarray) -> np.ndarray: + """筛选多边形内的点""" + from matplotlib.path import Path + if points.size == 0: + return points + + path = Path(polygon[:, :2]) + mask = path.contains_points(points[:, :2]) + return points[mask] + + async def _get_simulated_points(self, + polygon_coords: List[List[float]], + z_range: Optional[Tuple[float, float]] = None) -> np.ndarray: + """获取模拟点数据""" + # 同上一个类的生成模拟数据方法 + polygon_np = np.array(polygon_coords) + + if polygon_np.shape[0] < 3: + return np.array([]) + + x_min, y_min = polygon_np.min(axis=0) + x_max, y_max = polygon_np.max(axis=0) + + grid_size = 100 + x = np.linspace(x_min, x_max, grid_size) + y = np.linspace(y_min, y_max, grid_size) + xx, yy = np.meshgrid(x, y) + + points = np.column_stack([xx.ravel(), yy.ravel()]) + + np.random.seed(42) + base_elevation = 50.0 + terrain_variation = 10.0 + z = base_elevation + terrain_variation * np.sin(0.1 * xx).ravel() * np.cos(0.1 * yy).ravel() + + points = np.column_stack([points[:, 0], points[:, 1], z]) + + from matplotlib.path import Path + path = Path(polygon_np[:, :2]) + mask = path.contains_points(points[:, :2]) + points = points[mask] + + if z_range: + mask = (points[:, 2] >= z_range[0]) & (points[:, 2] <= z_range[1]) + points = points[mask] + + logger.info(f"生成 {points.shape[0]} 个模拟点") + return points + + async def get_data_bounds(self) -> Dict[str, List[float]]: + """获取数据边界""" + if self._tileset is None: + await self.initialize() + + bounds = { + "min": [float('inf'), float('inf'), float('inf')], + "max": [-float('inf'), -float('inf'), -float('inf')] + } + + if self._tileset and hasattr(self._tileset.root_tile, 'bounding_volume'): + bv = self._tileset.root_tile.bounding_volume + if hasattr(bv, 'get_corners'): + corners = bv.get_corners() + if corners is not None: + bounds["min"] = corners.min(axis=0).tolist() + bounds["max"] = corners.max(axis=0).tolist() + + return bounds + + def get_crs(self) -> str: + return self._crs \ No newline at end of file diff --git a/b3dm/tileset_to_ply.py b/b3dm/tileset_to_ply.py new file mode 100644 index 0000000..aae94ce --- /dev/null +++ b/b3dm/tileset_to_ply.py @@ -0,0 +1,395 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +3D Tiles Tileset to PLY Converter +将整个3D Tiles tileset转换为单个PLY文件 +""" + +import json +import struct +import os +import sys +from pathlib import Path +import numpy as np + +try: + import DracoPy +except ImportError: + print("警告: DracoPy库未安装，无法处理Draco压缩的数据") + print("请运行: pip install DracoPy") + DracoPy = None + +class TilesetToPLYConverter: + def __init__(self): + self.all_vertices = [] + self.vertex_count = 0 + + def multiply_matrix_vector(self, matrix, vector): + """4x4矩阵与4D向量相乘""" + # matrix是16个元素的列表，按列主序排列 + # 转换为4x4矩阵（行主序） + m = [ + [matrix[0], matrix[4], matrix[8], matrix[12]], + [matrix[1], matrix[5], matrix[9], matrix[13]], + [matrix[2], matrix[6], matrix[10], matrix[14]], + [matrix[3], matrix[7], matrix[11], matrix[15]] + ] + + # 向量扩展为齐次坐标 [x, y, z, 1] + v = [vector[0], vector[1], vector[2], 1.0] + + # 矩阵乘法 + result = [ + m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2] + m[0][3]*v[3], + m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2] + m[1][3]*v[3], + m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2] + m[2][3]*v[3] + ] + + return result + + def multiply_matrices(self, m1, m2): + """两个4x4矩阵相乘""" + # 将16元素列表转换为4x4矩阵 + def list_to_matrix(lst): + return [ + [lst[0], lst[4], lst[8], lst[12]], + [lst[1], lst[5], lst[9], lst[13]], + [lst[2], lst[6], lst[10], lst[14]], + [lst[3], lst[7], lst[11], lst[15]] + ] + + def matrix_to_list(mat): + return [ + mat[0][0], mat[1][0], mat[2][0], mat[3][0], + mat[0][1], mat[1][1], mat[2][1], mat[3][1], + mat[0][2], mat[1][2], mat[2][2], mat[3][2], + mat[0][3], mat[1][3], mat[2][3], mat[3][3] + ] + + a = list_to_matrix(m1) + b = list_to_matrix(m2) + + result = [[0 for _ in range(4)] for _ in range(4)] + for i in range(4): + for j in range(4): + for k in range(4): + result[i][j] += a[i][k] * b[k][j] + + return matrix_to_list(result) + + def apply_transform_to_vertices(self, vertices, transform_matrix): + """对顶点应用变换矩阵""" + if not transform_matrix: + return vertices + + transformed_vertices = [] + for vertex in vertices: + transformed = self.multiply_matrix_vector(transform_matrix, vertex) + transformed_vertices.append(transformed) + + return transformed_vertices + + def parse_tileset_json(self, tileset_path, parent_transform=None): + """解析tileset.json文件，收集B3DM文件和变换矩阵""" + try: + with open(tileset_path, 'r', encoding='utf-8') as f: + tileset_data = json.load(f) + + b3dm_files = [] + + def process_node(node, base_path, accumulated_transform): + # 获取当前节点的变换矩阵 + current_transform = node.get('transform') + + # 计算累积变换矩阵 + if current_transform and accumulated_transform: + # 矩阵相乘：accumulated_transform * current_transform + final_transform = self.multiply_matrices(accumulated_transform, current_transform) + elif current_transform: + final_transform = current_transform + else: + final_transform = accumulated_transform + + if 'content' in node and 'uri' in node['content']: + uri = node['content']['uri'] + if uri.endswith('.b3dm'): + full_path = os.path.join(base_path, uri) + if os.path.exists(full_path): + b3dm_files.append((full_path, final_transform)) + elif uri.endswith('.json'): + # 递归处理子tileset + sub_tileset_path = os.path.join(base_path, uri) + if os.path.exists(sub_tileset_path): + sub_files = self.parse_tileset_json(sub_tileset_path, final_transform) + b3dm_files.extend(sub_files) + + if 'children' in node: + for child in node['children']: + process_node(child, base_path, final_transform) + + base_path = os.path.dirname(tileset_path) + if 'root' in tileset_data: + process_node(tileset_data['root'], base_path, parent_transform) + + return b3dm_files + + except Exception as e: + print(f"解析tileset.json时出错: {e}") + return [] + + def parse_b3dm_file(self, file_path): + """解析B3DM文件""" + try: + with open(file_path, 'rb') as f: + # 读取B3DM头部 + magic = f.read(4) + if magic != b'b3dm': + print(f"警告: {file_path} 不是有效的B3DM文件") + return None + + version = struct.unpack(' len(glb_data): + break + + chunk_length = struct.unpack(' [输出PLY文件路径]") + print("示例: python tileset_to_ply.py tileset.json output.ply") + return + + tileset_path = sys.argv[1] + output_path = sys.argv[2] if len(sys.argv) > 2 else "merged_tileset.ply" + + if not os.path.exists(tileset_path): + print(f"错误: 文件不存在 {tileset_path}") + return + + converter = TilesetToPLYConverter() + success = converter.convert_tileset_to_ply(tileset_path, output_path) + + if success: + print("\n转换完成!") + else: + print("\n转换失败!") + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/b3dm/volume_calculator.py b/b3dm/volume_calculator.py new file mode 100644 index 0000000..3719bfb --- /dev/null +++ b/b3dm/volume_calculator.py @@ -0,0 +1,684 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +三维模型体积计算器 +基于指定地理范围内的三维模型数据，使用三角构网方法计算体积 +""" + +import json +import struct +import os +import sys +from pathlib import Path +import numpy as np +from scipy.spatial import Delaunay +import math +from tqdm import tqdm + +try: + import DracoPy +except ImportError: + print("警告: DracoPy库未安装，无法处理Draco压缩的数据") + print("请运行: pip install DracoPy") + DracoPy = None + +class VolumeCalculator: + def __init__(self, location_file): + self.location_bounds = self.load_location_bounds(location_file) + self.all_vertices = [] + self.filtered_vertices = [] + + def load_location_bounds(self, location_file): + """加载地理范围边界""" + try: + with open(location_file, 'r', encoding='utf-8') as f: + coords = json.load(f) + + # 提取经纬度范围 + lons = [coord[0] for coord in coords] + lats = [coord[1] for coord in coords] + elevs = [coord[2] for coord in coords] + + bounds = { + 'min_lon': min(lons), + 'max_lon': max(lons), + 'min_lat': min(lats), + 'max_lat': max(lats), + 'min_elev': min(elevs), + 'max_elev': max(elevs), + 'coords': coords + } + + print(f"地理范围边界:") + print(f" 经度: {bounds['min_lon']:.8f} ~ {bounds['max_lon']:.8f}") + print(f" 纬度: {bounds['min_lat']:.8f} ~ {bounds['max_lat']:.8f}") + print(f" 高程: {bounds['min_elev']:.2f} ~ {bounds['max_elev']:.2f}") + + return bounds + + except Exception as e: + print(f"加载地理范围文件失败: {e}") + return None + + def wgs84_to_cartesian(self, lon, lat, elev): + """WGS84坐标转换为笛卡尔坐标（高精度算法）""" + # WGS84椭球参数（EPSG:4326） + a = 6378137.0 # 长半轴 (米) + f = 1/298.257223563 # 扁率 + e2 = 2*f - f*f # 第一偏心率平方 + + # 角度转弧度 + lon_rad = math.radians(lon) + lat_rad = math.radians(lat) + + # 计算卯酉圈曲率半径 + sin_lat = math.sin(lat_rad) + cos_lat = math.cos(lat_rad) + sin_lon = math.sin(lon_rad) + cos_lon = math.cos(lon_rad) + + N = a / math.sqrt(1 - e2 * sin_lat * sin_lat) + + # 高精度笛卡尔坐标计算 + x = (N + elev) * cos_lat * cos_lon + y = (N + elev) * cos_lat * sin_lon + z = (N * (1 - e2) + elev) * sin_lat + + return [x, y, z] + + def cartesian_to_wgs84(self, x, y, z): + """笛卡尔坐标转换为WGS84坐标（高精度迭代法）""" + # WGS84椭球参数 + a = 6378137.0 # 长半轴 + f = 1/298.257223563 # 扁率 + e2 = 2*f - f*f # 第一偏心率平方 + ep2 = e2 / (1 - e2) # 第二偏心率平方 + + # 计算经度（精确值） + lon = math.atan2(y, x) + + # 计算纬度和高程（使用Bowring迭代法） + p = math.sqrt(x*x + y*y) + + if p == 0: + # 极点情况 + lat = math.pi/2 if z > 0 else -math.pi/2 + elev = abs(z) - a * math.sqrt(1 - e2) + return [math.degrees(lon), math.degrees(lat), elev] + + # 初始估计 + theta = math.atan2(z, p * (1 - f)) + lat_prev = math.atan2(z + ep2 * a * (1 - f) * math.sin(theta)**3, + p - e2 * a * math.cos(theta)**3) + + # 迭代求解纬度 + max_iterations = 10 + tolerance = 1e-12 + + for i in range(max_iterations): + N = a / math.sqrt(1 - e2 * math.sin(lat_prev)**2) + elev = p / math.cos(lat_prev) - N + + # 更新纬度估计 + lat_new = math.atan2(z + e2 * N * math.sin(lat_prev), p) + + # 检查收敛性 + if abs(lat_new - lat_prev) < tolerance: + break + + lat_prev = lat_new + + # 最终计算高程 + N = a / math.sqrt(1 - e2 * math.sin(lat_prev)**2) + elev = p / math.cos(lat_prev) - N + + return [math.degrees(lon), math.degrees(lat_prev), elev] + + def is_point_in_bounds(self, vertex): + """检查点是否在指定的地理范围内""" + if not self.location_bounds: + return True + + # 将笛卡尔坐标转换为WGS84 + try: + lon, lat, elev = self.cartesian_to_wgs84(vertex[0], vertex[1], vertex[2]) + + # 检查是否在边界范围内 + return (self.location_bounds['min_lon'] <= lon <= self.location_bounds['max_lon'] and + self.location_bounds['min_lat'] <= lat <= self.location_bounds['max_lat']) + except: + return False + + def multiply_matrix_vector(self, matrix, vector): + """4x4矩阵与4D向量相乘""" + m = [ + [matrix[0], matrix[4], matrix[8], matrix[12]], + [matrix[1], matrix[5], matrix[9], matrix[13]], + [matrix[2], matrix[6], matrix[10], matrix[14]], + [matrix[3], matrix[7], matrix[11], matrix[15]] + ] + + v = [vector[0], vector[1], vector[2], 1.0] + + result = [ + m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2] + m[0][3]*v[3], + m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2] + m[1][3]*v[3], + m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2] + m[2][3]*v[3] + ] + + return result + + def multiply_matrices(self, m1, m2): + """两个4x4矩阵相乘""" + def list_to_matrix(lst): + return [ + [lst[0], lst[4], lst[8], lst[12]], + [lst[1], lst[5], lst[9], lst[13]], + [lst[2], lst[6], lst[10], lst[14]], + [lst[3], lst[7], lst[11], lst[15]] + ] + + def matrix_to_list(mat): + return [ + mat[0][0], mat[1][0], mat[2][0], mat[3][0], + mat[0][1], mat[1][1], mat[2][1], mat[3][1], + mat[0][2], mat[1][2], mat[2][2], mat[3][2], + mat[0][3], mat[1][3], mat[2][3], mat[3][3] + ] + + a = list_to_matrix(m1) + b = list_to_matrix(m2) + + result = [[0 for _ in range(4)] for _ in range(4)] + for i in range(4): + for j in range(4): + for k in range(4): + result[i][j] += a[i][k] * b[k][j] + + return matrix_to_list(result) + + def apply_transform_to_vertices(self, vertices, transform_matrix): + """对顶点应用变换矩阵""" + if not transform_matrix: + return vertices + + transformed_vertices = [] + for vertex in vertices: + transformed = self.multiply_matrix_vector(transform_matrix, vertex) + transformed_vertices.append(transformed) + + return transformed_vertices + + def parse_tileset_json(self, tileset_path, parent_transform=None): + """解析tileset.json文件，收集B3DM文件和变换矩阵""" + try: + with open(tileset_path, 'r', encoding='utf-8') as f: + tileset_data = json.load(f) + + b3dm_files = [] + + def process_node(node, base_path, accumulated_transform): + current_transform = node.get('transform') + + if current_transform and accumulated_transform: + final_transform = self.multiply_matrices(accumulated_transform, current_transform) + elif current_transform: + final_transform = current_transform + else: + final_transform = accumulated_transform + + if 'content' in node and 'uri' in node['content']: + uri = node['content']['uri'] + if uri.endswith('.b3dm'): + full_path = os.path.join(base_path, uri) + if os.path.exists(full_path): + b3dm_files.append((full_path, final_transform)) + elif uri.endswith('.json'): + sub_tileset_path = os.path.join(base_path, uri) + if os.path.exists(sub_tileset_path): + sub_files = self.parse_tileset_json(sub_tileset_path, final_transform) + b3dm_files.extend(sub_files) + + if 'children' in node: + for child in node['children']: + process_node(child, base_path, final_transform) + + base_path = os.path.dirname(tileset_path) + if 'root' in tileset_data: + process_node(tileset_data['root'], base_path, parent_transform) + + return b3dm_files + + except Exception as e: + print(f"解析tileset.json时出错: {e}") + return [] + + def parse_b3dm_file(self, file_path): + """解析B3DM文件""" + try: + with open(file_path, 'rb') as f: + magic = f.read(4) + if magic != b'b3dm': + return None + + version = struct.unpack(' len(glb_data): + break + + chunk_length = struct.unpack(' 0 and max_side / min_side > max_aspect_ratio: + return False + + return True + + def calculate_circumcenter_and_radius(self, p1, p2, p3): + """计算三角形外接圆圆心和半径（高精度算法）""" + try: + x1, y1 = p1[0], p1[1] + x2, y2 = p2[0], p2[1] + x3, y3 = p3[0], p3[1] + + # 使用C#代码中的高精度外接圆计算公式 + d = 2 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) + + if abs(d) < 1e-10: # 三点共线 + return None, float('inf') + + ux = ((x1*x1 + y1*y1) * (y2 - y3) + (x2*x2 + y2*y2) * (y3 - y1) + (x3*x3 + y3*y3) * (y1 - y2)) / d + uy = ((x1*x1 + y1*y1) * (x3 - x2) + (x2*x2 + y2*y2) * (x1 - y3) + (x3*x3 + y3*y3) * (x2 - x1)) / d + + # 计算半径 + radius = np.sqrt((ux - x1)**2 + (uy - y1)**2) + + return np.array([ux, uy]), radius + except: + return None, float('inf') + + def calculate_volume_delaunay(self, vertices, base_elevation=None, min_angle=10.0, use_quality_filter=True): + """使用Delaunay三角剖分计算体积""" + if len(vertices) < 4: + print("顶点数量不足，无法进行三角剖分") + return 0.0 + + try: + points = np.array(vertices) + + if base_elevation is None: + base_elevation = np.min(points[:, 2]) + + print(f"Delaunay方法使用基准面高程: {base_elevation:.2f} 米") + print(f"质量控制参数: 最小角度={min_angle}°, 质量过滤={'开启' if use_quality_filter else '关闭'}") + + adjusted_points = points.copy() + adjusted_points[:, 2] = np.maximum(points[:, 2] - base_elevation, 0) + + print("正在进行Delaunay三角剖分...") + tri = Delaunay(adjusted_points) + + valid_simplices = [] + total_simplices = len(tri.simplices) + + if use_quality_filter: + print("正在进行三角形质量检查...") + for simplex in tqdm(tri.simplices, desc="质量检查", unit="个"): + p0 = adjusted_points[simplex[0]] + p1 = adjusted_points[simplex[1]] + p2 = adjusted_points[simplex[2]] + p3 = adjusted_points[simplex[3]] + + faces = [(p0, p1, p2), (p0, p1, p3), (p0, p2, p3), (p1, p2, p3)] + valid_faces = 0 + + for face in faces: + if self.is_valid_triangle(face[0], face[1], face[2], min_angle): + valid_faces += 1 + + if valid_faces >= 3: + valid_simplices.append(simplex) + + print(f"质量过滤: {len(valid_simplices)}/{total_simplices} 个四面体通过质量检查") + else: + valid_simplices = tri.simplices + + total_volume = 0.0 + valid_volume_count = 0 + + print(f"计算 {len(valid_simplices)} 个四面体的体积...") + + for simplex in tqdm(valid_simplices, desc="计算四面体体积", unit="个"): + p0 = adjusted_points[simplex[0]] + p1 = adjusted_points[simplex[1]] + p2 = adjusted_points[simplex[2]] + p3 = adjusted_points[simplex[3]] + + v1 = p1 - p0 + v2 = p2 - p0 + v3 = p3 - p0 + + det = np.linalg.det(np.array([v1, v2, v3])) + volume = abs(det) / 6.0 + + if volume > 1e-12: + total_volume += volume + valid_volume_count += 1 + + print(f"有效体积计算: {valid_volume_count}/{len(valid_simplices)} 个四面体") + + return total_volume + + except Exception as e: + print(f"Delaunay三角剖分计算体积失败: {e}") + return 0.0 + + + def load_and_filter_vertices(self, tileset_path): + """加载并过滤指定范围内的顶点""" + print(f"开始处理tileset: {tileset_path}") + + # 解析tileset获取所有B3DM文件 + b3dm_files = self.parse_tileset_json(tileset_path) + print(f"找到 {len(b3dm_files)} 个B3DM文件") + + if not b3dm_files: + print("未找到任何B3DM文件") + return False + + # 处理每个B3DM文件 + processed_count = 0 + total_vertices = 0 + filtered_count = 0 + + for i, (b3dm_file, transform_matrix) in enumerate(tqdm(b3dm_files, desc="处理B3DM文件", unit="文件")): + # print(f"处理文件 {i+1}/{len(b3dm_files)}: {os.path.basename(b3dm_file)}") + + vertices = self.parse_b3dm_file(b3dm_file) + if vertices: + # 应用变换矩阵 + if transform_matrix: + vertices = self.apply_transform_to_vertices(vertices, transform_matrix) + + # 过滤范围内的顶点 + for vertex in vertices: + total_vertices += 1 + if self.is_point_in_bounds(vertex): + self.filtered_vertices.append(vertex) + filtered_count += 1 + + self.all_vertices.extend(vertices) + processed_count += 1 + # tqdm.write(f" {os.path.basename(b3dm_file)}: 提取到 {len(vertices)} 个顶点") + + print(f"\n成功处理 {processed_count}/{len(b3dm_files)} 个文件") + print(f"总计提取 {total_vertices} 个顶点") + print(f"范围内顶点 {filtered_count} 个") + + return len(self.filtered_vertices) > 0 + + + + def calculate_volume(self, tileset_path, base_elevation=None, min_angle=10.0, use_quality_filter=True): + """计算指定范围内三维模型的体积 + + Args: + tileset_path: tileset.json文件路径 + base_elevation: 基准面高程 + min_angle: 最小角度约束（度） + use_quality_filter: 是否启用质量过滤 + """ + if not self.load_and_filter_vertices(tileset_path): + print("未找到范围内的顶点数据") + return 0.0 + + print(f"\n开始计算体积，使用Delaunay三角剖分方法") + print(f"参与计算的顶点数: {len(self.filtered_vertices)}") + + points = np.array(self.filtered_vertices) + if base_elevation is None: + base_elevation = np.min(points[:, 2]) + + print(f"统一基准面高程: {base_elevation:.2f} 米") + + volume = self.calculate_volume_delaunay(self.filtered_vertices, base_elevation, min_angle, use_quality_filter) + + print(f"\n计算结果:") + print(f"体积: {volume:.6f} 立方米") + print(f"体积: {volume/1000000:.6f} 立方千米") + + return volume + +def main(): + if len(sys.argv) < 3: + print("用法: python volume_calculator.py [基准面高程] [最小角度] [质量过滤]") + print("基准面高程: 基准面高程(米)，默认使用最低点") + print("最小角度: 三角形最小角度约束(度)，默认10.0") + print("质量过滤: 是否启用质量过滤(true/false)，默认true") + print("示例: python volume_calculator.py boundary.json tileset.json 100.0 15.0 true") + return + + location_file = sys.argv[1] + tileset_path = sys.argv[2] + + # 解析可选参数 + base_elevation = None + min_angle = 10.0 + use_quality_filter = True + + if len(sys.argv) > 3: + try: + base_elevation = float(sys.argv[3]) + except ValueError: + print("错误：基准面高程必须是数字") + return + + if len(sys.argv) > 4: + try: + min_angle = float(sys.argv[4]) + except ValueError: + print("错误：最小角度必须是数字") + return + + if len(sys.argv) > 5: + use_quality_filter = sys.argv[5].lower() in ['true', '1', 'yes', 'on'] + + if not os.path.exists(location_file): + print(f"错误: 位置文件不存在 {location_file}") + return + + if not os.path.exists(tileset_path): + print(f"错误: tileset文件不存在 {tileset_path}") + return + + calculator = VolumeCalculator(location_file) + volume = calculator.calculate_volume(tileset_path, base_elevation, min_angle, use_quality_filter) + + if volume > 0: + print("\n体积计算完成!") + else: + print("\n体积计算失败!") + +if __name__ == "__main__": + main() \ No newline at end of file