pengxiao
219232de74
以讲义内容为骨架迁移到标准目录格式： - 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>
983 lines
30 KiB
Python
983 lines
30 KiB
Python
"""
|
|
空间表征示例 (Spatial Representation Example)
|
|
============================================
|
|
|
|
本示例展示空间智能系统中的空间表征方法。
|
|
空间表征是对地理空间现象的抽象和建模，是空间推理的基础。
|
|
|
|
核心概念:
|
|
1. 空间对象模型 - 点、线、面等几何对象
|
|
2. 空间关系模型 - 拓扑、距离、方向关系
|
|
3. 空间场模型 - 连续表面的表示
|
|
4. 空间索引结构 - 加速空间查询
|
|
5. 多尺度表征 - 不同详细程度的表示
|
|
|
|
应用场景:
|
|
- 地理信息系统 (GIS)
|
|
- 空间数据库
|
|
- 空间推理引擎
|
|
- 地理可视化
|
|
|
|
作者: CC4SI 项目组
|
|
"""
|
|
|
|
import math
|
|
import json
|
|
from typing import List, Dict, Tuple, Optional, Any, Set
|
|
from dataclasses import dataclass, field
|
|
from enum import Enum
|
|
from abc import ABC, abstractmethod
|
|
import random
|
|
|
|
|
|
# ============================================================================
|
|
# 几何类型与空间对象
|
|
# ============================================================================
|
|
|
|
class GeometryType(Enum):
|
|
"""几何类型枚举"""
|
|
POINT = "Point"
|
|
LINESTRING = "LineString"
|
|
POLYGON = "Polygon"
|
|
MULTIPOINT = "MultiPoint"
|
|
MULTILINESTRING = "MultiLineString"
|
|
MULTIPOLYGON = "MultiPolygon"
|
|
GEOMETRYCOLLECTION = "GeometryCollection"
|
|
|
|
|
|
class SpatialReference(Enum):
|
|
"""空间参考系统"""
|
|
WGS84 = "EPSG:4326" # 经纬度
|
|
WEB_MERCATOR = "EPSG:3857" # Web墨卡托
|
|
CGCS2000 = "EPSG:4490" # 中国大地坐标系统
|
|
|
|
|
|
@dataclass
|
|
class Point:
|
|
"""点几何对象"""
|
|
x: float
|
|
y: float
|
|
srid: int = 4326 # 空间参考ID
|
|
|
|
def __iter__(self):
|
|
"""支持解包"""
|
|
return iter((self.x, self.y))
|
|
|
|
def __iter__(self):
|
|
return iter((self.x, self.y))
|
|
|
|
def to_tuple(self) -> Tuple[float, float]:
|
|
"""转换为元组"""
|
|
return (self.x, self.y)
|
|
|
|
def to_geojson(self) -> Dict:
|
|
"""转换为GeoJSON"""
|
|
return {
|
|
"type": "Point",
|
|
"coordinates": [self.x, self.y]
|
|
}
|
|
|
|
def distance_to(self, other: 'Point') -> float:
|
|
"""计算到另一个点的欧氏距离"""
|
|
return math.sqrt((self.x - other.x)**2 + (self.y - other.y)**2)
|
|
|
|
def __repr__(self) -> str:
|
|
return f"Point({self.x:.4f}, {self.y:.4f})"
|
|
|
|
|
|
@dataclass
|
|
class LineString:
|
|
"""线几何对象"""
|
|
coordinates: List[Point]
|
|
srid: int = 4326
|
|
|
|
@property
|
|
def length(self) -> float:
|
|
"""计算线的长度"""
|
|
if len(self.coordinates) < 2:
|
|
return 0.0
|
|
total = 0.0
|
|
for i in range(len(self.coordinates) - 1):
|
|
total += self.coordinates[i].distance_to(self.coordinates[i + 1])
|
|
return total
|
|
|
|
def to_geojson(self) -> Dict:
|
|
"""转换为GeoJSON"""
|
|
return {
|
|
"type": "LineString",
|
|
"coordinates": [[p.x, p.y] for p in self.coordinates]
|
|
}
|
|
|
|
def get_point_at(self, ratio: float) -> Point:
|
|
"""
|
|
获取线上指定比例位置的点
|
|
|
|
Args:
|
|
ratio: 0到1之间的比例值
|
|
|
|
Returns:
|
|
该位置的点
|
|
"""
|
|
if ratio <= 0:
|
|
return self.coordinates[0]
|
|
if ratio >= 1:
|
|
return self.coordinates[-1]
|
|
|
|
target_length = self.length * ratio
|
|
accumulated = 0.0
|
|
|
|
for i in range(len(self.coordinates) - 1):
|
|
segment_length = self.coordinates[i].distance_to(self.coordinates[i + 1])
|
|
if accumulated + segment_length >= target_length:
|
|
# 在此段上
|
|
segment_ratio = (target_length - accumulated) / segment_length
|
|
p1 = self.coordinates[i]
|
|
p2 = self.coordinates[i + 1]
|
|
return Point(
|
|
p1.x + (p2.x - p1.x) * segment_ratio,
|
|
p1.y + (p2.y - p1.y) * segment_ratio
|
|
)
|
|
accumulated += segment_length
|
|
|
|
return self.coordinates[-1]
|
|
|
|
def __repr__(self) -> str:
|
|
return f"LineString({len(self.coordinates)} points)"
|
|
|
|
|
|
@dataclass
|
|
class Polygon:
|
|
"""面几何对象"""
|
|
exterior: List[Point] # 外环
|
|
interiors: List[List[Point]] = field(default_factory=list) # 内环（空洞）
|
|
srid: int = 4326
|
|
|
|
@property
|
|
def area(self) -> float:
|
|
"""使用鞋带公式计算多边形面积"""
|
|
return self._ring_area(self.exterior) - sum(
|
|
self._ring_area(interior) for interior in self.interiors
|
|
)
|
|
|
|
def _ring_area(self, ring: List[Point]) -> float:
|
|
"""计算环的面积（绝对值）"""
|
|
if len(ring) < 3:
|
|
return 0.0
|
|
area = 0.0
|
|
for i in range(len(ring)):
|
|
j = (i + 1) % len(ring)
|
|
area += ring[i].x * ring[j].y
|
|
area -= ring[j].x * ring[i].y
|
|
return abs(area) / 2
|
|
|
|
@property
|
|
def centroid(self) -> Point:
|
|
"""计算多边形质心"""
|
|
if len(self.exterior) < 3:
|
|
return self.exterior[0] if self.exterior else Point(0, 0)
|
|
|
|
# 简化：使用外环的顶点平均值
|
|
avg_x = sum(p.x for p in self.exterior) / len(self.exterior)
|
|
avg_y = sum(p.y for p in self.exterior) / len(self.exterior)
|
|
return Point(avg_x, avg_y)
|
|
|
|
def contains_point(self, point: Point) -> bool:
|
|
"""
|
|
判断点是否在多边形内（射线法）
|
|
|
|
Args:
|
|
point: 待判断的点
|
|
|
|
Returns:
|
|
是否在多边形内
|
|
"""
|
|
return self._point_in_ring(point, self.exterior) and \
|
|
all(not self._point_in_ring(point, interior)
|
|
for interior in self.interiors)
|
|
|
|
def _point_in_ring(self, point: Point, ring: List[Point]) -> bool:
|
|
"""射线法判断点是否在环内"""
|
|
if len(ring) < 3:
|
|
return False
|
|
|
|
x, y = point.x, point.y
|
|
inside = False
|
|
|
|
for i in range(len(ring)):
|
|
j = (i + 1) % len(ring)
|
|
xi, yi = ring[i].x, ring[i].y
|
|
xj, yj = ring[j].x, ring[j].y
|
|
|
|
# 检查射线与边的交点
|
|
if ((yi > y) != (yj > y)) and \
|
|
(x < (xj - xi) * (y - yi) / (yj - yi + 1e-10) + xi):
|
|
inside = not inside
|
|
|
|
return inside
|
|
|
|
def to_geojson(self) -> Dict:
|
|
"""转换为GeoJSON"""
|
|
coords = [[[p.x, p.y] for p in self.exterior]]
|
|
coords.extend([[[p.x, p.y] for p in interior] for interior in self.interiors])
|
|
|
|
return {
|
|
"type": "Polygon",
|
|
"coordinates": coords
|
|
}
|
|
|
|
def __repr__(self) -> str:
|
|
holes = len(self.interiors)
|
|
return f"Polygon({len(self.exterior)} vertices{f', {holes} holes' if holes > 0 else ''})"
|
|
|
|
|
|
@dataclass
|
|
class Envelope:
|
|
"""包围盒"""
|
|
min_x: float
|
|
max_x: float
|
|
min_y: 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 area(self) -> float:
|
|
return self.width * self.height
|
|
|
|
@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 intersects(self, other: 'Envelope') -> bool:
|
|
"""判断是否与另一个包围盒相交"""
|
|
return not (self.max_x < other.min_x or self.min_x > other.max_x or
|
|
self.max_y < other.min_y or self.min_y > other.max_y)
|
|
|
|
def union(self, other: 'Envelope') -> 'Envelope':
|
|
"""计算与另一个包围盒的并集"""
|
|
return Envelope(
|
|
min_x=min(self.min_x, other.min_x),
|
|
max_x=max(self.max_x, other.max_x),
|
|
min_y=min(self.min_y, other.min_y),
|
|
max_y=max(self.max_y, other.max_y)
|
|
)
|
|
|
|
def __repr__(self) -> str:
|
|
return f"Envelope([{self.min_x:.2f}, {self.min_y:.2f}] -> [{self.max_x:.2f}, {self.max_y:.2f}])"
|
|
|
|
|
|
# ============================================================================
|
|
# 空间特征对象
|
|
# ============================================================================
|
|
|
|
@dataclass
|
|
class SpatialFeature:
|
|
"""
|
|
空间特征
|
|
|
|
包含几何和属性信息的完整空间对象。
|
|
"""
|
|
id: str
|
|
geometry: Any # Point, LineString, Polygon等
|
|
properties: Dict[str, Any] = field(default_factory=dict)
|
|
|
|
def to_geojson(self) -> Dict:
|
|
"""转换为GeoJSON Feature"""
|
|
geom_data = None
|
|
if isinstance(self.geometry, Point):
|
|
geom_data = self.geometry.to_geojson()
|
|
elif isinstance(self.geometry, LineString):
|
|
geom_data = self.geometry.to_geojson()
|
|
elif isinstance(self.geometry, Polygon):
|
|
geom_data = self.geometry.to_geojson()
|
|
|
|
return {
|
|
"type": "Feature",
|
|
"id": self.id,
|
|
"geometry": geom_data,
|
|
"properties": self.properties
|
|
}
|
|
|
|
def to_geojson_collection(self) -> Dict:
|
|
"""转换为GeoJSON FeatureCollection"""
|
|
return {
|
|
"type": "FeatureCollection",
|
|
"features": [self.to_geojson()]
|
|
}
|
|
|
|
def envelope(self) -> Envelope:
|
|
"""计算包围盒"""
|
|
if isinstance(self.geometry, Point):
|
|
return Envelope(
|
|
self.geometry.x, self.geometry.x,
|
|
self.geometry.y, self.geometry.y
|
|
)
|
|
elif isinstance(self.geometry, (LineString, list)):
|
|
coords = self.geometry.coordinates if isinstance(self.geometry, LineString) else self.geometry
|
|
xs = [p.x for p in coords]
|
|
ys = [p.y for p in coords]
|
|
return Envelope(min(xs), max(xs), min(ys), max(ys))
|
|
elif isinstance(self.geometry, Polygon):
|
|
xs = [p.x for p in self.geometry.exterior]
|
|
ys = [p.y for p in self.geometry.exterior]
|
|
return Envelope(min(xs), max(xs), min(ys), max(ys))
|
|
else:
|
|
raise ValueError(f"Unsupported geometry type: {type(self.geometry)}")
|
|
|
|
def __repr__(self) -> str:
|
|
return f"SpatialFeature(id={self.id}, geom={type(self.geometry).__name__})"
|
|
|
|
|
|
# ============================================================================
|
|
# 空间关系
|
|
# ============================================================================
|
|
|
|
class SpatialRelation(Enum):
|
|
"""空间关系类型"""
|
|
EQUALS = "equals"
|
|
DISJOINT = "disjoint"
|
|
INTERSECTS = "intersects"
|
|
TOUCHES = "touches"
|
|
CROSSES = "crosses"
|
|
WITHIN = "within"
|
|
CONTAINS = "contains"
|
|
OVERLAPS = "overlaps"
|
|
|
|
|
|
@dataclass
|
|
class TopologyRelation:
|
|
"""
|
|
拓扑关系描述
|
|
|
|
基于DE-9IM (Dimensionally Extended nine-Intersection Model)模型。
|
|
"""
|
|
relation: SpatialRelation
|
|
confidence: float = 1.0
|
|
|
|
def __repr__(self) -> str:
|
|
return f"TopologyRelation({self.relation.value}, conf={self.confidence:.2f})"
|
|
|
|
|
|
def calculate_relation(geom1: Any, geom2: Any) -> TopologyRelation:
|
|
"""
|
|
计算两个几何对象的空间关系
|
|
|
|
Args:
|
|
geom1: 第一个几何对象
|
|
geom2: 第二个几何对象
|
|
|
|
Returns:
|
|
拓扑关系
|
|
"""
|
|
# 点-点关系
|
|
if isinstance(geom1, Point) and isinstance(geom2, Point):
|
|
if geom1.x == geom2.x and geom1.y == geom2.y:
|
|
return TopologyRelation(SpatialRelation.EQUALS)
|
|
return TopologyRelation(SpatialRelation.DISJOINT)
|
|
|
|
# 点-多边形关系
|
|
if isinstance(geom1, Point) and isinstance(geom2, Polygon):
|
|
if geom2.contains_point(geom1):
|
|
return TopologyRelation(SpatialRelation.WITHIN)
|
|
# 检查是否在边界上
|
|
for i in range(len(geom2.exterior)):
|
|
p1, p2 = geom2.exterior[i], geom2.exterior[(i + 1) % len(geom2.exterior)]
|
|
if point_on_segment(geom1, p1, p2):
|
|
return TopologyRelation(SpatialRelation.TOUCHES)
|
|
return TopologyRelation(SpatialRelation.DISJOINT)
|
|
|
|
# 多边形-点关系
|
|
if isinstance(geom1, Polygon) and isinstance(geom2, Point):
|
|
relation = calculate_relation(geom2, geom1)
|
|
if relation.relation == SpatialRelation.WITHIN:
|
|
return TopologyRelation(SpatialRelation.CONTAINS)
|
|
return relation
|
|
|
|
# 多边形-多边形关系（简化版：只检查相交）
|
|
if isinstance(geom1, Polygon) and isinstance(geom2, Polygon):
|
|
env1 = envelope_of_polygon(geom1)
|
|
env2 = envelope_of_polygon(geom2)
|
|
|
|
if not env1.intersects(env2):
|
|
return TopologyRelation(SpatialRelation.DISJOINT)
|
|
|
|
# 简化：检查是否有交点
|
|
if geom1.contains_point(geom2.exterior[0]):
|
|
return TopologyRelation(SpatialRelation.CONTAINS)
|
|
if geom2.contains_point(geom1.exterior[0]):
|
|
return TopologyRelation(SpatialRelation.WITHIN)
|
|
|
|
return TopologyRelation(SpatialRelation.INTERSECTS)
|
|
|
|
return TopologyRelation(SpatialRelation.DISJOINT)
|
|
|
|
|
|
def point_on_segment(point: Point, seg_start: Point, seg_end: Point,
|
|
tolerance: float = 1e-6) -> bool:
|
|
"""判断点是否在线段上"""
|
|
# 检查点是否在线段的包围盒内
|
|
if not (min(seg_start.x, seg_end.x) - tolerance <= point.x <=
|
|
max(seg_start.x, seg_end.x) + tolerance and
|
|
min(seg_start.y, seg_end.y) - tolerance <= point.y <=
|
|
max(seg_start.y, seg_end.y) + tolerance):
|
|
return False
|
|
|
|
# 检查三点共线
|
|
cross = (seg_end.x - seg_start.x) * (point.y - seg_start.y) - \
|
|
(seg_end.y - seg_start.y) * (point.x - seg_start.x)
|
|
return abs(cross) < tolerance
|
|
|
|
|
|
def envelope_of_polygon(polygon: Polygon) -> Envelope:
|
|
"""计算多边形的包围盒"""
|
|
xs = [p.x for p in polygon.exterior]
|
|
ys = [p.y for p in polygon.exterior]
|
|
return Envelope(min(xs), max(xs), min(ys), max(ys))
|
|
|
|
|
|
# ============================================================================
|
|
# 空间索引
|
|
# ============================================================================
|
|
|
|
class RTreeNode:
|
|
"""R树节点"""
|
|
def __init__(self, is_leaf: bool = False):
|
|
self.is_leaf = is_leaf
|
|
self.envelope: Optional[Envelope] = None
|
|
self.children: List['RTreeNode'] = []
|
|
self.features: List[SpatialFeature] = []
|
|
|
|
def update_envelope(self):
|
|
"""更新节点的包围盒"""
|
|
if self.is_leaf and self.features:
|
|
envelopes = [f.envelope() for f in self.features]
|
|
self.envelope = envelopes[0]
|
|
for env in envelopes[1:]:
|
|
self.envelope = self.envelope.union(env)
|
|
elif not self.is_leaf and self.children:
|
|
self.envelope = self.children[0].envelope
|
|
for child in self.children[1:]:
|
|
if child.envelope:
|
|
self.envelope = self.envelope.union(child.envelope)
|
|
|
|
|
|
class RTree:
|
|
"""
|
|
R树空间索引
|
|
|
|
用于加速空间查询的树状索引结构。
|
|
"""
|
|
|
|
def __init__(self, max_children: int = 4):
|
|
"""
|
|
初始化R树
|
|
|
|
Args:
|
|
max_children: 每个节点的最大子节点数
|
|
"""
|
|
self.max_children = max_children
|
|
self.root = RTreeNode(is_leaf=True)
|
|
self.size = 0
|
|
|
|
def insert(self, feature: SpatialFeature) -> None:
|
|
"""插入空间特征"""
|
|
self._insert(self.root, feature)
|
|
self.size += 1
|
|
|
|
def _insert(self, node: RTreeNode, feature: SpatialFeature) -> None:
|
|
"""递归插入"""
|
|
feature_env = feature.envelope()
|
|
|
|
if node.is_leaf:
|
|
node.features.append(feature)
|
|
node.update_envelope()
|
|
|
|
# 如果超过容量，分裂节点
|
|
if len(node.features) > self.max_children:
|
|
self._split(node)
|
|
else:
|
|
# 选择最佳子节点
|
|
best_child = self._choose_best_child(node, feature_env)
|
|
self._insert(best_child, feature)
|
|
node.update_envelope()
|
|
|
|
def _choose_best_child(self, node: RTreeNode, env: Envelope) -> RTreeNode:
|
|
"""选择插入代价最小的子节点"""
|
|
best = None
|
|
best_increase = float('inf')
|
|
|
|
for child in node.children:
|
|
if child.envelope is None:
|
|
continue
|
|
union_env = child.envelope.union(env)
|
|
increase = union_env.area - child.envelope.area
|
|
|
|
if increase < best_increase:
|
|
best_increase = increase
|
|
best = child
|
|
|
|
return best or node.children[0]
|
|
|
|
def _split(self, node: RTreeNode) -> None:
|
|
"""分裂节点 (简化版)"""
|
|
if node.is_leaf:
|
|
# 简单分裂：将特征分成两组
|
|
mid = len(node.features) // 2
|
|
group1 = node.features[:mid]
|
|
group2 = node.features[mid:]
|
|
|
|
node.features = group1
|
|
|
|
new_leaf = RTreeNode(is_leaf=True)
|
|
new_leaf.features = group2
|
|
new_leaf.update_envelope()
|
|
|
|
# 更新父节点
|
|
if node == self.root and not node.children:
|
|
# 根节点分裂
|
|
new_root = RTreeNode(is_leaf=False)
|
|
new_root.children = [node, new_leaf]
|
|
new_root.update_envelope()
|
|
self.root = new_root
|
|
else:
|
|
# 简化：不处理非根节点的分裂
|
|
pass
|
|
|
|
def query(self, envelope: Envelope) -> List[SpatialFeature]:
|
|
"""查询与包围盒相交的所有特征"""
|
|
results = []
|
|
self._query(self.root, envelope, results)
|
|
return results
|
|
|
|
def _query(self, node: RTreeNode, envelope: Envelope,
|
|
results: List[SpatialFeature]) -> None:
|
|
"""递归查询"""
|
|
if node.envelope and not node.envelope.intersects(envelope):
|
|
return
|
|
|
|
if node.is_leaf:
|
|
for feature in node.features:
|
|
if feature.envelope().intersects(envelope):
|
|
results.append(feature)
|
|
else:
|
|
for child in node.children:
|
|
self._query(child, envelope, results)
|
|
|
|
def nearest_neighbor(self, point: Point, k: int = 1) -> List[Tuple[SpatialFeature, float]]:
|
|
"""
|
|
K近邻查询
|
|
|
|
Args:
|
|
point: 查询点
|
|
k: 返回的最近邻数量
|
|
|
|
Returns:
|
|
(特征, 距离) 列表
|
|
"""
|
|
candidates = []
|
|
self._collect_candidates(self.root, candidates)
|
|
|
|
# 计算距离并排序
|
|
distances = []
|
|
for feature in candidates:
|
|
if isinstance(feature.geometry, Point):
|
|
dist = feature.geometry.distance_to(point)
|
|
distances.append((feature, dist))
|
|
else:
|
|
# 非点几何：使用包围盒中心距离
|
|
env = feature.envelope()
|
|
center = env.center
|
|
dist = math.sqrt((center.x - point.x)**2 + (center.y - point.y)**2)
|
|
distances.append((feature, dist))
|
|
|
|
distances.sort(key=lambda x: x[1])
|
|
return distances[:k]
|
|
|
|
def _collect_candidates(self, node: RTreeNode, results: List[SpatialFeature]) -> None:
|
|
"""收集所有候选特征"""
|
|
if node.is_leaf:
|
|
results.extend(node.features)
|
|
else:
|
|
for child in node.children:
|
|
self._collect_candidates(child, results)
|
|
|
|
|
|
# ============================================================================
|
|
# 空间场模型
|
|
# ============================================================================
|
|
|
|
class GridField:
|
|
"""
|
|
栅格场模型
|
|
|
|
使用规则网格表示连续空间现象。
|
|
"""
|
|
|
|
def __init__(self, bounds: Envelope, rows: int, cols: int,
|
|
nodata: float = -9999):
|
|
"""
|
|
初始化栅格场
|
|
|
|
Args:
|
|
bounds: 空间范围
|
|
rows: 行数
|
|
cols: 列数
|
|
nodata: 无数据值
|
|
"""
|
|
self.bounds = bounds
|
|
self.rows = rows
|
|
self.cols = cols
|
|
self.nodata = nodata
|
|
self.data = [[nodata for _ in range(cols)] for _ in range(rows)]
|
|
|
|
@property
|
|
def cell_width(self) -> float:
|
|
"""获取单元格宽度"""
|
|
return self.bounds.width / self.cols
|
|
|
|
@property
|
|
def cell_height(self) -> float:
|
|
"""获取单元格高度"""
|
|
return self.bounds.height / self.rows
|
|
|
|
def get_cell_index(self, point: Point) -> Optional[Tuple[int, int]]:
|
|
"""
|
|
获取点对应的栅格索引
|
|
|
|
Args:
|
|
point: 空间点
|
|
|
|
Returns:
|
|
(行索引, 列索引) 或 None
|
|
"""
|
|
if not self.bounds.contains(point):
|
|
return None
|
|
|
|
col = int((point.x - self.bounds.min_x) / self.cell_width)
|
|
row = int((self.bounds.max_y - point.y) / self.cell_height)
|
|
|
|
col = max(0, min(col, self.cols - 1))
|
|
row = max(0, min(row, self.rows - 1))
|
|
|
|
return (row, col)
|
|
|
|
def set_value(self, row: int, col: int, value: float) -> None:
|
|
"""设置栅格值"""
|
|
if 0 <= row < self.rows and 0 <= col < self.cols:
|
|
self.data[row][col] = value
|
|
|
|
def get_value(self, row: int, col: int) -> float:
|
|
"""获取栅格值"""
|
|
if 0 <= row < self.rows and 0 <= col < self.cols:
|
|
return self.data[row][col]
|
|
return self.nodata
|
|
|
|
def get_value_at_point(self, point: Point) -> float:
|
|
"""获取点位置的值（最近邻）"""
|
|
idx = self.get_cell_index(point)
|
|
if idx:
|
|
return self.get_value(idx[0], idx[1])
|
|
return self.nodata
|
|
|
|
def interpolate_at_point(self, point: Point) -> float:
|
|
"""双线性插值获取点位置的值"""
|
|
idx = self.get_cell_index(point)
|
|
if not idx:
|
|
return self.nodata
|
|
|
|
row, col = idx
|
|
|
|
# 获取四个角点的值
|
|
values = []
|
|
for r in range(row, min(row + 2, self.rows)):
|
|
for c in range(col, min(col + 2, self.cols)):
|
|
values.append((r, c, self.data[r][c]))
|
|
|
|
if len(values) < 4 or any(v[2] == self.nodata for v in values):
|
|
return self.get_value(row, col) # 回退到最近邻
|
|
|
|
# 双线性插值
|
|
# 简化实现
|
|
return self.get_value(row, col)
|
|
|
|
def get_statistics(self) -> Dict[str, float]:
|
|
"""获取统计信息"""
|
|
values = [v for row in self.data for v in row if v != self.nodata]
|
|
|
|
if not values:
|
|
return {"count": 0, "min": self.nodata, "max": self.nodata,
|
|
"mean": self.nodata, "std": 0}
|
|
|
|
import statistics
|
|
return {
|
|
"count": len(values),
|
|
"min": min(values),
|
|
"max": max(values),
|
|
"mean": statistics.mean(values),
|
|
"std": statistics.stdev(values) if len(values) > 1 else 0
|
|
}
|
|
|
|
def to_geojson(self) -> Dict:
|
|
"""转换为GeoJSON（简化：输出为点集）"""
|
|
features = []
|
|
for r in range(self.rows):
|
|
for c in range(self.cols):
|
|
val = self.data[r][c]
|
|
if val != self.nodata:
|
|
# 计算点坐标
|
|
x = self.bounds.min_x + (c + 0.5) * self.cell_width
|
|
y = self.bounds.max_y - (r + 0.5) * self.cell_height
|
|
|
|
features.append({
|
|
"type": "Feature",
|
|
"geometry": {"type": "Point", "coordinates": [x, y]},
|
|
"properties": {"value": val}
|
|
})
|
|
|
|
return {
|
|
"type": "FeatureCollection",
|
|
"features": features
|
|
}
|
|
|
|
|
|
# ============================================================================
|
|
# 主程序
|
|
# ============================================================================
|
|
|
|
def main():
|
|
"""主程序 - 演示空间表征的使用"""
|
|
|
|
print("="*70)
|
|
print("空间表征示例演示")
|
|
print("="*70)
|
|
|
|
# ========================================================================
|
|
# 1. 基础几何对象
|
|
# ========================================================================
|
|
print("\n[部分 1] 基础几何对象")
|
|
print("-" * 50)
|
|
|
|
# 创建点
|
|
point1 = Point(120.5, 30.2)
|
|
point2 = Point(121.0, 30.5)
|
|
print(f"\n点对象:")
|
|
print(f" point1 = {point1}")
|
|
print(f" point2 = {point2}")
|
|
print(f" 距离 = {point1.distance_to(point2):.4f}")
|
|
|
|
# 创建线
|
|
linestring = LineString([
|
|
Point(120.0, 30.0),
|
|
Point(120.5, 30.2),
|
|
Point(121.0, 30.5),
|
|
Point(121.5, 30.3)
|
|
])
|
|
print(f"\n线对象:")
|
|
print(f" {linestring}")
|
|
print(f" 长度 = {linestring.length:.4f}")
|
|
print(f" 中点 = {linestring.get_point_at(0.5)}")
|
|
|
|
# 创建多边形
|
|
polygon = Polygon(exterior=[
|
|
Point(120.0, 30.0),
|
|
Point(121.0, 30.0),
|
|
Point(121.0, 31.0),
|
|
Point(120.0, 31.0),
|
|
Point(120.0, 30.0)
|
|
])
|
|
print(f"\n多边形对象:")
|
|
print(f" {polygon}")
|
|
print(f" 面积 = {polygon.area:.4f}")
|
|
print(f" 质心 = {polygon.centroid}")
|
|
|
|
# 点包含测试
|
|
test_inside = Point(120.5, 30.5)
|
|
test_outside = Point(121.5, 30.5)
|
|
print(f" {test_inside} 在多边形内: {polygon.contains_point(test_inside)}")
|
|
print(f" {test_outside} 在多边形内: {polygon.contains_point(test_outside)}")
|
|
|
|
# ========================================================================
|
|
# 2. 空间特征
|
|
# ========================================================================
|
|
print("\n\n[部分 2] 空间特征")
|
|
print("-" * 50)
|
|
|
|
feature1 = SpatialFeature(
|
|
id="poi_001",
|
|
geometry=Point(120.5, 30.2),
|
|
properties={
|
|
"name": "杭州西湖",
|
|
"type": "scenic_spot",
|
|
"rating": 4.8
|
|
}
|
|
)
|
|
|
|
feature2 = SpatialFeature(
|
|
id="zone_001",
|
|
geometry=polygon,
|
|
properties={
|
|
"name": "开发区A",
|
|
"type": "industrial_zone",
|
|
"area_ha": 10000
|
|
}
|
|
)
|
|
|
|
print(f"\n特征1: {feature1}")
|
|
print(f" 属性: {feature1.properties}")
|
|
print(f" 包围盒: {feature1.envelope()}")
|
|
|
|
print(f"\n特征2: {feature2}")
|
|
print(f" 属性: {feature2.properties}")
|
|
print(f" 包围盒: {feature2.envelope()}")
|
|
|
|
# GeoJSON输出
|
|
print(f"\nGeoJSON输出:")
|
|
print(json.dumps(feature1.to_geojson(), ensure_ascii=False, indent=2))
|
|
|
|
# ========================================================================
|
|
# 3. 空间关系
|
|
# ========================================================================
|
|
print("\n\n[部分 3] 空间关系")
|
|
print("-" * 50)
|
|
|
|
# 点与多边形的关系
|
|
point_a = Point(120.5, 30.5) # 在多边形内
|
|
point_b = Point(121.5, 30.5) # 在多边形外
|
|
point_c = Point(120.0, 30.5) # 在边界上
|
|
|
|
print(f"\n点-多边形关系测试:")
|
|
print(f" {point_a} 与多边形: {calculate_relation(point_a, polygon)}")
|
|
print(f" {point_b} 与多边形: {calculate_relation(point_b, polygon)}")
|
|
print(f" {point_c} 与多边形: {calculate_relation(point_c, polygon)}")
|
|
|
|
# 多边形-多边形关系
|
|
poly2 = Polygon(exterior=[
|
|
Point(120.5, 29.5),
|
|
Point(121.5, 29.5),
|
|
Point(121.5, 30.5),
|
|
Point(120.5, 30.5),
|
|
Point(120.5, 29.5)
|
|
])
|
|
print(f"\n多边形-多边形关系:")
|
|
print(f" poly1 与 poly2: {calculate_relation(polygon, poly2)}")
|
|
|
|
# ========================================================================
|
|
# 4. 空间索引
|
|
# ========================================================================
|
|
print("\n\n[部分 4] R树空间索引")
|
|
print("-" * 50)
|
|
|
|
# 创建R树
|
|
rtree = RTree(max_children=4)
|
|
|
|
# 插入一些特征
|
|
features = []
|
|
for i in range(20):
|
|
x = random.uniform(119, 122)
|
|
y = random.uniform(29, 32)
|
|
feature = SpatialFeature(
|
|
id=f"feature_{i:03d}",
|
|
geometry=Point(x, y),
|
|
properties={"value": random.uniform(0, 100)}
|
|
)
|
|
features.append(feature)
|
|
rtree.insert(feature)
|
|
|
|
print(f"\n已插入 {rtree.size} 个特征到R树")
|
|
|
|
# 范围查询
|
|
query_env = Envelope(120, 121, 30, 31)
|
|
results = rtree.query(query_env)
|
|
print(f"\n范围查询 {query_env}:")
|
|
print(f" 找到 {len(results)} 个特征")
|
|
for f in results[:5]:
|
|
print(f" - {f.id}: {f.geometry}")
|
|
|
|
# K近邻查询
|
|
query_point = Point(120.5, 30.5)
|
|
neighbors = rtree.nearest_neighbor(query_point, k=5)
|
|
print(f"\nK近邻查询 (中心点: {query_point}):")
|
|
for i, (f, dist) in enumerate(neighbors, 1):
|
|
print(f" {i}. {f.id}: 距离 = {dist:.4f}")
|
|
|
|
# ========================================================================
|
|
# 5. 栅格场模型
|
|
# ========================================================================
|
|
print("\n\n[部分 5] 栅格场模型")
|
|
print("-" * 50)
|
|
|
|
# 创建栅格场
|
|
grid = GridField(
|
|
bounds=Envelope(119, 122, 29, 32),
|
|
rows=30,
|
|
cols=30,
|
|
nodata=-9999
|
|
)
|
|
|
|
# 填充一些模拟数据 (温度场)
|
|
for r in range(grid.rows):
|
|
for c in range(grid.cols):
|
|
# 计算坐标
|
|
x = grid.bounds.min_x + (c + 0.5) * grid.cell_width
|
|
y = grid.bounds.max_y - (r + 0.5) * grid.cell_height
|
|
|
|
# 模拟温度场 (简单的径向基函数)
|
|
center_x, center_y = 120.5, 30.5
|
|
dist = math.sqrt((x - center_x)**2 + (y - center_y)**2)
|
|
temp = 25 - dist * 2 # 中心25度，向外递减
|
|
|
|
grid.set_value(r, c, round(temp, 2))
|
|
|
|
print(f"\n栅格场信息:")
|
|
print(f" 范围: {grid.bounds}")
|
|
print(f" 尺寸: {grid.rows} x {grid.cols}")
|
|
print(f" 单元格大小: {grid.cell_width:.4f} x {grid.cell_height:.4f}")
|
|
|
|
stats = grid.get_statistics()
|
|
print(f" 统计: {stats}")
|
|
|
|
# 点查询
|
|
sample_point = Point(120.5, 30.5)
|
|
value = grid.get_value_at_point(sample_point)
|
|
print(f"\n点位置 {sample_point} 的值: {value}")
|
|
|
|
# ========================================================================
|
|
# 6. GeoJSON导出
|
|
# ========================================================================
|
|
print("\n\n[部分 6] GeoJSON导出")
|
|
print("-" * 50)
|
|
|
|
# 创建特征集合
|
|
feature_collection = {
|
|
"type": "FeatureCollection",
|
|
"features": [
|
|
feature1.to_geojson(),
|
|
feature2.to_geojson()
|
|
]
|
|
}
|
|
|
|
print("\n特征集合:")
|
|
print(json.dumps(feature_collection, ensure_ascii=False, indent=2))
|
|
|
|
print("\n" + "="*70)
|
|
print("演示完成!")
|
|
print("="*70)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|