pengxiao/2026_DesignAI
1
0
Fork 0
Code Issues 1 Pull Requests Actions Packages Projects 1 Releases Wiki Activity
Files
main
2026_DesignAI/dofile/examples/00-introduction/setup-first-assistant/spatial_helper.py
T
pengxiao 219232de74 refactor: 重组项目目录结构 
以讲义内容为骨架迁移到标准目录格式：
- officefile/ 主内容（12章 + 附录 + CC4SI补充）
- dofile/ 代码示例（11个Python脚本）
- data/ 图片资源
- output/ 生成输出（忽略）
- Archive/ 归档旧目录（忽略）
- .claude/skills/ 保留markdown-to-docx工具链
- .pandoc/ 保留CSL和本地化配置

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
2026-05-25 14:00:56 +08:00

963 lines
30 KiB
Python
Raw Permalink Blame History

"""
空间分析助手类 (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()
Reference in New Issue View Git Blame Copy Permalink