""" 空间分析助手类 (Spatial Analysis Assistant) ======================================== 这是一个完整的空间分析助手实现,展示了如何构建一个基础的 空间智能系统。该系统可以帮助用户进行空间数据处理、分析和可视化。 主要功能: 1. 空间数据加载与管理 2. 空间关系计算 (距离、方位、包含关系等) 3. 空间统计分析 4. 空间插值与预测 5. 多准则决策分析 作者: CC4SI 项目组 日期: 2025-01 """ import json import math from typing import List, Dict, Tuple, Optional, Any, Union from dataclasses import dataclass, field from enum import Enum import random # ============================================================================ # 数据结构定义 # ============================================================================ class GeometryType(Enum): """几何类型枚举""" POINT = "Point" LINESTRING = "LineString" POLYGON = "Polygon" MULTIPOINT = "MultiPoint" @dataclass class Point: """点几何类""" x: float y: float z: Optional[float] = None properties: Dict[str, Any] = field(default_factory=dict) def __repr__(self) -> str: if self.z is not None: return f"Point({self.x:.2f}, {self.y:.2f}, {self.z:.2f})" return f"Point({self.x:.2f}, {self.y:.2f})" @dataclass class BoundingBox: """边界框类""" min_x: float min_y: float max_x: float max_y: float @property def width(self) -> float: """获取边界框宽度""" return self.max_x - self.min_x @property def height(self) -> float: """获取边界框高度""" return self.max_y - self.min_y @property def center(self) -> Point: """获取边界框中心点""" return Point( (self.min_x + self.max_x) / 2, (self.min_y + self.max_y) / 2 ) def contains(self, point: Point) -> bool: """检查点是否在边界框内""" return (self.min_x <= point.x <= self.max_x and self.min_y <= point.y <= self.max_y) def __repr__(self) -> str: return f"BoundingBox[{self.min_x:.2f},{self.min_y:.2f} -> {self.max_x:.2f},{self.max_y:.2f}]" @dataclass class SpatialFeature: """空间要素类""" id: str geometry: Union[Point, List[Point]] # 简化: 点或点列表 properties: Dict[str, Any] = field(default_factory=dict) def __repr__(self) -> str: return f"SpatialFeature(id={self.id}, geometry={type(self.geometry).__name__})" # ============================================================================ # 空间分析助手主类 # ============================================================================ class SpatialHelper: """ 空间分析助手类 这是系统的核心类,提供空间数据处理和分析的主要功能。 """ def __init__(self, name: str = "空间分析助手", version: str = "1.0.0"): """ 初始化空间分析助手 Args: name: 助手名称 version: 版本号 """ self.name = name self.version = version self.features: List[SpatialFeature] = [] self.coordinate_system: str = "EPSG:4326" # 默认 WGS84 self.metadata: Dict[str, Any] = {} print(f"[{self.name}] v{self.version} 初始化完成") # ------------------------------------------------------------------------ # 数据加载与管理 # ------------------------------------------------------------------------ def add_feature(self, feature: SpatialFeature) -> bool: """ 添加空间要素 Args: feature: 要添加的空间要素 Returns: 是否添加成功 """ try: self.features.append(feature) return True except Exception as e: print(f"添加要素失败: {e}") return False def add_point(self, x: float, y: float, feature_id: str = None, properties: Dict[str, Any] = None) -> SpatialFeature: """ 添加点要素 Args: x: X坐标 y: Y坐标 feature_id: 要素ID properties: 属性字典 Returns: 创建的空间要素 """ if feature_id is None: feature_id = f"point_{len(self.features)}" point = Point(x, y) feature = SpatialFeature( id=feature_id, geometry=point, properties=properties or {} ) self.add_feature(feature) return feature def load_from_geojson(self, geojson_str: str) -> int: """ 从 GeoJSON 字符串加载数据 Args: geojson_str: GeoJSON 格式字符串 Returns: 加载的要素数量 """ try: data = json.loads(geojson_str) count = 0 if data.get("type") == "FeatureCollection": for feature_data in data.get("features", []): feature = self._parse_geojson_feature(feature_data) if feature: self.add_feature(feature) count += 1 print(f"从 GeoJSON 加载了 {count} 个要素") return count except json.JSONDecodeError as e: print(f"GeoJSON 解析失败: {e}") return 0 def _parse_geojson_feature(self, feature_data: Dict) -> Optional[SpatialFeature]: """解析 GeoJSON 要素""" try: feature_id = feature_data.get("id", f"feature_{len(self.features)}") properties = feature_data.get("properties", {}) geometry = feature_data.get("geometry", {}) geom_type = geometry.get("type") if geom_type == "Point": coordinates = geometry.get("coordinates", []) point = Point(coordinates[0], coordinates[1]) return SpatialFeature(id=feature_id, geometry=point, properties=properties) return None except Exception: return None def get_feature_by_id(self, feature_id: str) -> Optional[SpatialFeature]: """根据ID获取要素""" for feature in self.features: if feature.id == feature_id: return feature return None def get_feature_count(self) -> int: """获取要素数量""" return len(self.features) def clear_features(self) -> None: """清空所有要素""" self.features.clear() print("已清空所有要素") # ------------------------------------------------------------------------ # 空间关系计算 # ------------------------------------------------------------------------ @staticmethod def calculate_distance(point1: Point, point2: Point) -> float: """ 计算两点间的欧氏距离 Args: point1: 第一个点 point2: 第二个点 Returns: 距离值 """ dx = point2.x - point1.x dy = point2.y - point1.y return math.sqrt(dx * dx + dy * dy) @staticmethod def calculate_bearing(point1: Point, point2: Point) -> float: """ 计算从 point1 到 point2 的方位角 (度) Args: point1: 起始点 point2: 目标点 Returns: 方位角 (0-360度) """ dx = point2.x - point1.x dy = point2.y - point1.y radians = math.atan2(dy, dx) degrees = math.degrees(radians) return (degrees + 360) % 360 @staticmethod def calculate_midpoint(point1: Point, point2: Point) -> Point: """ 计算两点间的中点 Args: point1: 第一个点 point2: 第二个点 Returns: 中点 """ return Point( (point1.x + point2.x) / 2, (point1.y + point2.y) / 2 ) def find_nearest_neighbor(self, target_point: Point, max_distance: float = float('inf')) -> Optional[Tuple[SpatialFeature, float]]: """ 查找最近邻要素 Args: target_point: 目标点 max_distance: 最大搜索距离 Returns: (最近的要素, 距离) 或 None """ nearest_feature = None min_distance = float('inf') for feature in self.features: if isinstance(feature.geometry, Point): dist = self.calculate_distance(target_point, feature.geometry) if dist < min_distance and dist <= max_distance: min_distance = dist nearest_feature = feature if nearest_feature: return nearest_feature, min_distance return None def find_neighbors_within_distance(self, target_point: Point, distance: float) -> List[Tuple[SpatialFeature, float]]: """ 查找指定距离内的所有要素 Args: target_point: 目标点 distance: 搜索半径 Returns: (要素, 距离) 列表,按距离排序 """ neighbors = [] for feature in self.features: if isinstance(feature.geometry, Point): dist = self.calculate_distance(target_point, feature.geometry) if dist <= distance: neighbors.append((feature, dist)) neighbors.sort(key=lambda x: x[1]) return neighbors def calculate_bounding_box(self) -> Optional[BoundingBox]: """ 计算所有要素的边界框 Returns: 边界框对象,如果没有要素则返回 None """ if not self.features: return None points = [f.geometry for f in self.features if isinstance(f.geometry, Point)] if not points: return None min_x = min(p.x for p in points) max_x = max(p.x for p in points) min_y = min(p.y for p in points) max_y = max(p.y for p in points) return BoundingBox(min_x, min_y, max_x, max_y) # ------------------------------------------------------------------------ # 空间统计分析 # ------------------------------------------------------------------------ def calculate_centroid(self) -> Optional[Point]: """ 计算所有点要素的质心 Returns: 质心点 """ points = [f.geometry for f in self.features if isinstance(f.geometry, Point)] if not points: return None avg_x = sum(p.x for p in points) / len(points) avg_y = sum(p.y for p in points) / len(points) return Point(avg_x, avg_y) def calculate_mean_center(self, weight_field: str = None) -> Optional[Point]: """ 计算加权或未加权的平均中心 Args: weight_field: 权重字段名 Returns: 平均中心点 """ points = [] weights = [] for feature in self.features: if isinstance(feature.geometry, Point): points.append(feature.geometry) if weight_field: weights.append(feature.properties.get(weight_field, 1.0)) else: weights.append(1.0) if not points: return None total_weight = sum(weights) avg_x = sum(p.x * w for p, w in zip(points, weights)) / total_weight avg_y = sum(p.y * w for p, w in zip(points, weights)) / total_weight return Point(avg_x, avg_y) def calculate_standard_distance(self) -> Optional[float]: """ 计算标准距离 (标准差圆) Returns: 标准距离值 """ centroid = self.calculate_centroid() if not centroid: return None points = [f.geometry for f in self.features if isinstance(f.geometry, Point)] if not points: return None n = len(points) squared_distances = [(p.x - centroid.x)**2 + (p.y - centroid.y)**2 for p in points] return math.sqrt(sum(squared_distances) / n) def calculate_spatial_autocorrelation(self, field: str) -> Optional[float]: """ 计算 Moran's I 空间自相关指数 Args: field: 要分析的属性字段 Returns: Moran's I 值 """ # 简化实现: 使用距离权重 points_data = [] for feature in self.features: if isinstance(feature.geometry, Point) and field in feature.properties: points_data.append((feature.geometry, feature.properties[field])) n = len(points_data) if n < 2: return None mean_value = sum(v for _, v in points_data) / n # 计算权重矩阵 (距离倒数) weights = {} total_weight = 0 for i, (p1, v1) in enumerate(points_data): for j, (p2, v2) in enumerate(points_data): if i != j: dist = self.calculate_distance(p1, p2) w = 1 / (dist + 0.001) # 避免除零 weights[(i, j)] = w total_weight += w # 计算 Moran's I numerator = 0 denominator = 0 for i, (p1, v1) in enumerate(points_data): for j, (p2, v2) in enumerate(points_data): if i != j: w = weights.get((i, j), 0) numerator += w * (v1 - mean_value) * (v2 - mean_value) denominator += (v1 - mean_value) ** 2 if denominator == 0: return None morans_i = (n / total_weight) * (numerator / denominator) return morans_i # ------------------------------------------------------------------------ # 空间插值与预测 # ------------------------------------------------------------------------ def inverse_distance_weighting(self, target_point: Point, power: float = 2.0, field: str = "value", max_distance: float = None) -> Optional[float]: """ 反距离加权插值 (IDW) Args: target_point: 目标点 power: 距离幂次 field: 插值字段 max_distance: 最大搜索距离 Returns: 插值结果 """ points_data = [] for feature in self.features: if isinstance(feature.geometry, Point) and field in feature.properties: points_data.append((feature.geometry, feature.properties[field])) if not points_data: return None numerator = 0.0 denominator = 0.0 for point, value in points_data: dist = self.calculate_distance(target_point, point) if max_distance and dist > max_distance: continue if dist < 1e-10: # 几乎重合 return value weight = 1.0 / (dist ** power) numerator += weight * value denominator += weight if denominator == 0: return None return numerator / denominator def simple_trend_prediction(self, field: str, target_x: float, target_y: float) -> Optional[float]: """ 基于简单趋势的预测 (线性回归) Args: field: 预测字段 target_x: 目标X坐标 target_y: 目标Y坐标 Returns: 预测值 """ points_data = [] for feature in self.features: if isinstance(feature.geometry, Point) and field in feature.properties: points_data.append({ 'x': feature.geometry.x, 'y': feature.geometry.y, 'z': feature.properties[field] }) if len(points_data) < 3: return None # 简单的线性趋势: z = a + b*x + c*y n = len(points_data) sum_x = sum(p['x'] for p in points_data) sum_y = sum(p['y'] for p in points_data) sum_z = sum(p['z'] for p in points_data) sum_xx = sum(p['x']**2 for p in points_data) sum_yy = sum(p['y']**2 for p in points_data) sum_xy = sum(p['x'] * p['y'] for p in points_data) sum_xz = sum(p['x'] * p['z'] for p in points_data) sum_yz = sum(p['y'] * p['z'] for p in points_data) # 简化: 只使用x方向趋势 try: # z = a + b*x b = (n * sum_xz - sum_x * sum_z) / (n * sum_xx - sum_x**2) a = (sum_z - b * sum_x) / n return a + b * target_x except ZeroDivisionError: return None # ------------------------------------------------------------------------ # 多准则决策分析 (MCDA) # ------------------------------------------------------------------------ def weighted_sum_model(self, criteria: List[str], weights: List[float], feature_ids: List[str] = None) -> List[Tuple[SpatialFeature, float]]: """ 加权求和模型 (WSM) Args: criteria: 评价准则列表 weights: 各准则权重 feature_ids: 参与评价的要素ID列表 Returns: (要素, 得分) 列表,按得分降序排列 """ if len(criteria) != len(weights): print("错误: 准则数量与权重数量不匹配") return [] if abs(sum(weights) - 1.0) > 0.001: print(f"警告: 权重总和为 {sum(weights)}, 建议归一化为 1.0") # 确定评价范围 features_to_eval = self.features if feature_ids: features_to_eval = [f for f in self.features if f.id in feature_ids] results = [] # 归一化参数 min_max = {} for criterion in criteria: values = [] for f in features_to_eval: if criterion in f.properties: values.append(f.properties[criterion]) if values: min_max[criterion] = (min(values), max(values)) for feature in features_to_eval: score = 0.0 valid = True for criterion, weight in zip(criteria, weights): if criterion not in feature.properties: valid = False break value = feature.properties[criterion] cmin, cmax = min_max.get(criterion, (0, 1)) # 归一化 (假设值越大越好) if cmax - cmin > 0: normalized = (value - cmin) / (cmax - cmin) else: normalized = 0.5 score += weight * normalized if valid: results.append((feature, score)) results.sort(key=lambda x: x[1], reverse=True) return results def weighted_product_model(self, criteria: List[str], weights: List[float], feature_ids: List[str] = None) -> List[Tuple[SpatialFeature, float]]: """ 加权乘积模型 (WPM) Args: criteria: 评价准则列表 weights: 各准则权重 feature_ids: 参与评价的要素ID列表 Returns: (要素, 得分) 列表,按得分降序排列 """ if len(criteria) != len(weights): print("错误: 准则数量与权重数量不匹配") return [] features_to_eval = self.features if feature_ids: features_to_eval = [f for f in self.features if f.id in feature_ids] results = [] # 归一化参数 min_max = {} for criterion in criteria: values = [] for f in features_to_eval: if criterion in f.properties: values.append(f.properties[criterion]) if values: min_max[criterion] = (min(values), max(values)) for feature in features_to_eval: product = 1.0 valid = True for criterion, weight in zip(criteria, weights): if criterion not in feature.properties: valid = False break value = feature.properties[criterion] cmin, cmax = min_max.get(criterion, (0, 1)) if cmax - cmin > 0: normalized = (value - cmin) / (cmax - cmin) else: normalized = 1.0 product *= normalized ** weight if valid: results.append((feature, product)) results.sort(key=lambda x: x[1], reverse=True) return results def topsis(self, criteria: List[str], weights: List[float], benefit_criteria: List[bool] = None) -> List[Tuple[SpatialFeature, float]]: """ TOPSIS (逼近理想解排序法) Args: criteria: 评价准则列表 weights: 各准则权重 benefit_criteria: 是否为效益型准则 (True=越大越好, False=越小越好) Returns: (要素, 相对贴近度) 列表,按贴近度降序排列 """ if benefit_criteria is None: benefit_criteria = [True] * len(criteria) features_to_eval = [f for f in self.features if all(c in f.properties for c in criteria)] if len(features_to_eval) == 0: return [] n = len(features_to_eval) m = len(criteria) # 构建决策矩阵 decision_matrix = [] for feature in features_to_eval: row = [feature.properties[c] for c in criteria] decision_matrix.append(row) # 归一化决策矩阵 normalized_matrix = [] for j in range(m): column = [decision_matrix[i][j] for i in range(n)] norm = math.sqrt(sum(x**2 for x in column)) for i in range(n): if j == 0: normalized_matrix.append([]) normalized_matrix[i].append(decision_matrix[i][j] / norm if norm > 0 else 0) # 加权归一化矩阵 weighted_matrix = [] for i in range(n): weighted_matrix.append([normalized_matrix[i][j] * weights[j] for j in range(m)]) # 确定理想解和负理想解 ideal_positive = [] ideal_negative = [] for j in range(m): column = [weighted_matrix[i][j] for i in range(n)] if benefit_criteria[j]: ideal_positive.append(max(column)) ideal_negative.append(min(column)) else: ideal_positive.append(min(column)) ideal_negative.append(max(column)) # 计算距离和相对贴近度 results = [] for i, feature in enumerate(features_to_eval): dist_positive = math.sqrt( sum((weighted_matrix[i][j] - ideal_positive[j])**2 for j in range(m)) ) dist_negative = math.sqrt( sum((weighted_matrix[i][j] - ideal_negative[j])**2 for j in range(m)) ) closeness = dist_negative / (dist_positive + dist_negative) if (dist_positive + dist_negative) > 0 else 0 results.append((feature, closeness)) results.sort(key=lambda x: x[1], reverse=True) return results # ------------------------------------------------------------------------ # 缓冲区分析 # ------------------------------------------------------------------------ def create_buffer_analysis(self, feature_id: str, buffer_distance: float) -> Dict[str, Any]: """ 缓冲区分析 Args: feature_id: 中心要素ID buffer_distance: 缓冲距离 Returns: 分析结果字典 """ target_feature = self.get_feature_by_id(feature_id) if not target_feature or not isinstance(target_feature.geometry, Point): return {"error": "找不到指定的点要素"} center = target_feature.geometry # 查找缓冲区内的要素 features_in_buffer = self.find_neighbors_within_distance(center, buffer_distance) # 计算统计信息 values_in_buffer = [f.properties for f, _ in features_in_buffer] return { "center_feature": feature_id, "buffer_distance": buffer_distance, "count": len(features_in_buffer), "features": [(f.id, dist) for f, dist in features_in_buffer], "statistics": { "avg_distance": sum(dist for _, dist in features_in_buffer) / len(features_in_buffer) if features_in_buffer else 0 } } # ------------------------------------------------------------------------ # 可视化辅助 (文本形式) # ------------------------------------------------------------------------ def print_summary(self) -> None: """打印数据摘要""" print(f"\n{'='*50}") print(f"空间分析助手摘要: {self.name} v{self.version}") print(f"{'='*50}") print(f"要素数量: {len(self.features)}") print(f"坐标系: {self.coordinate_system}") bbox = self.calculate_bounding_box() if bbox: print(f"边界范围: {bbox}") centroid = self.calculate_centroid() if centroid: print(f"质心位置: {centroid}") print(f"{'='*50}\n") def print_features(self, limit: int = 10) -> None: """打印要素列表""" print(f"\n要素列表 (显示前 {min(limit, len(self.features))} 个):") print("-" * 60) for i, feature in enumerate(self.features[:limit]): if isinstance(feature.geometry, Point): print(f"{i+1}. ID: {feature.id:15s} 位置: {feature.geometry} 属性: {feature.properties}") if len(self.features) > limit: print(f"... 还有 {len(self.features) - limit} 个要素") print("-" * 60 + "\n") # ============================================================================ # 辅助函数 # ============================================================================ def create_sample_data(helper: SpatialHelper, n_points: int = 20) -> None: """创建示例数据""" print(f"生成 {n_points} 个随机样本点...") random.seed(42) # 可重现的随机数 # 生成随机点 for i in range(n_points): x = random.uniform(0, 100) y = random.uniform(0, 100) value = random.uniform(0, 100) population = random.randint(100, 10000) accessibility = random.uniform(0.3, 0.95) helper.add_point( x=x, y=y, feature_id=f"point_{i:03d}", properties={ "value": value, "population": population, "accessibility": accessibility, "name": f"位置_{i+1}" } ) print(f"已生成 {helper.get_feature_count()} 个样本点") # ============================================================================ # 主程序 # ============================================================================ def main(): """主程序 - 演示空间分析助手的使用""" print("="*70) print("空间分析助手 - 完整示例演示") print("="*70) # 1. 创建助手实例 print("\n[步骤 1] 创建空间分析助手") helper = SpatialHelper(name="城市空间分析助手", version="1.0.0") # 2. 添加示例数据 print("\n[步骤 2] 添加示例数据") create_sample_data(helper, n_points=15) # 3. 打印数据摘要 print("\n[步骤 3] 数据摘要") helper.print_summary() helper.print_features() # 4. 空间关系计算 print("\n[步骤 4] 空间关系计算") print("-" * 50) test_point = Point(50, 50) print(f"测试点: {test_point}") nearest = helper.find_nearest_neighbor(test_point) if nearest: feature, dist = nearest print(f"最近邻: {feature.id}, 距离: {dist:.2f}") neighbors = helper.find_neighbors_within_distance(test_point, 25) print(f"半径25内的邻居数量: {len(neighbors)}") # 5. 空间统计分析 print("\n[步骤 5] 空间统计分析") print("-" * 50) centroid = helper.calculate_centroid() print(f"质心: {centroid}") std_dist = helper.calculate_standard_distance() print(f"标准距离: {std_dist:.2f}") # 6. 多准则决策分析 print("\n[步骤 6] 多准则决策分析 (TOPSIS)") print("-" * 50) criteria = ["accessibility", "population"] weights = [0.6, 0.4] # 可达性权重更高 results = helper.topsis(criteria, weights) print("选址优先级排序 (基于可达性和人口):") for i, (feature, score) in enumerate(results[:5]): print(f" {i+1}. {feature.properties.get('name', feature.id)}: 得分={score:.4f}") # 7. 空间插值 print("\n[步骤 7] 空间插值预测") print("-" * 50) predict_point = Point(45, 55) predicted = helper.inverse_distance_weighting(predict_point, field="value") print(f"在 {predict_point} 处的插值预测: {predicted:.2f}") # 8. 缓冲区分析 print("\n[步骤 8] 缓冲区分析") print("-" * 50) buffer_result = helper.create_buffer_analysis("point_000", 30) print(f"以 point_000 为中心,半径30的缓冲区:") print(f" 包含要素数: {buffer_result['count']}") print(f" 平均距离: {buffer_result['statistics']['avg_distance']:.2f}") print("\n" + "="*70) print("演示完成!") print("="*70) if __name__ == "__main__": main()