pengxiao
a90f7adfa1
将 Markdown 源文件移入 md/，LaTeX 工作目录保留在 latex/， Word 导出移入 word/；删除临时脚本、调试截图和空 stub。 Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
35 KiB
35 KiB
02.1 空间表征
核心问题
机器如何"看懂"地理空间？ 栅格和矢量，谁更智能？还是"小孩子才做选择"？
概念讲解
空间表征的本质
空间表征是空间智能的基石——它解决的是**"如何用计算机能理解的方式描述空间"**这一根本问题。
现实世界 → 空间表征 → 计算操作
─────────────────────────────────────────────
真实景观 计算机表示 算法处理
🏔️ → 栅格/矢量 → 分析/推理
🌲 ↘
🏙️ 图表示 决策/优化
🛣️
表征类型对比
| 维度 | 栅格 (Raster) | 矢量 (Vector) |
|---|---|---|
| 基本单元 | 像元 (Pixel/Cell) | 点、线、面 (Point, Line, Polygon) |
| 数据结构 | 规则网格 | 坐标序列 |
| 适用场景 | 连续场、表面分析 | 离散要素、边界精确 |
| 典型操作 | 邻域分析、代数运算 | 拓扑分析、几何运算 |
| 存储效率 | 与分辨率强相关 | 与复杂度相关 |
| 代表数据 | DEM、遥感影像 | 行政边界、道路网 |
图表示 (Graph Representation)
许多空间问题可以抽象为图：
空间 → 图的转换
生态网络场景 图表示
───────────── ───────
源地A ──廊道──→ 源地B 节点A ──边(weight=5)──→ 节点B
│ │
└──廊道──→ 源地C └──边(weight=8)──→ 节点C
道路网络 图表示
交叉路口 节点
道路路段 加权边
图表示的优势：
- 将空间问题转化为成熟的图算法
- 天然支持连通性、路径分析
- 易于扩展（添加权重、方向）
多尺度表征
空间智能必须处理多尺度问题：
┌─────────────────────────────────────────────────────────────┐
│ 多尺度金字塔 │
├─────────────────────────────────────────────────────────────┤
│ │
│ Level 3 (1:1,000,000) ┌─────┐ │
│ 区域尺度 │ A │ │
│ └─────┘ │
│ 一个像元 = 1km² │
│ │
│ Level 2 (1:100,000) ┌─────┬─────┐ │
│ 景观尺度 │ A │ B │ │
│ └─────┴─────┘ │
│ 一个像元 = 100m × 100m │
│ │
│ Level 1 (1:10,000) ┌─┬─┬─┬─┬─┐ │
│ 局地尺度 │A│A│A│B│B│ │
│ └─┴─┴─┴─┴─┘ │
│ 一个像元 = 10m × 10m │
│ │
└─────────────────────────────────────────────────────────────┘
多尺度处理策略：
- 尺度金字塔：预先生成多分辨率版本
- 自适应网格：关键区域高分辨率，其他区域低分辨率
- 层次聚类：构建空间层次结构
空间索引 (Spatial Indexing)
当数据量大时，空间索引是效率的关键：
┌─────────────────────────────────────────────────────────────┐
│ 空间索引类型对比 │
├─────────────────────────────────────────────────────────────┤
│ │
│ 1. R-Tree (常用) │
│ ┌────────────────────┐ │
│ │ ┌───┐ ┌────┐ │ 层次包围盒 │
│ │ │ A│ │ B │ │ 适合：复杂几何体 │
│ │ └───┘ └────┘ │ │
│ └────────────────────┘ │
│ │
│ 2. Quadtree (栅格友好) │
│ ┌─────────┬─────────┐ │
│ │ ● │ │ 四叉递归分解 │
│ ├─────────┼─────────┤ 适合：点数据、栅格 │
│ │ ● │ ● │ │
│ └─────────┴─────────┘ │
│ │
│ 3. Grid Index (简单高效) │
│ ┌───┬───┬───┬───┐ │
│ │ │ ● │ │ │ 规则网格分桶 │
│ ├───┼───┼───┼───┤ 适合：均匀分布数据 │
│ │ │ │ ● │ │ │
│ └───┴───┴───┴───┘ │
│ │
└─────────────────────────────────────────────────────────────┘
设计原理
表征选择决策树
开始
│
▼
数据是什么类型？
┌──────┴──────┐
│ │
连续场/表面 离散要素
(高程、温度) (边界、道路)
│ │
▼ ▼
栅格优先 需要精确边界？
│ ┌──┴──┐
│ │ │
│ 是 否
│ │ │
│ ▼ ▼
│ 矢量 栅格也可
│ │ │
│ ▼ │
│ 需要拓扑分析？
│ │ │
│ ┌───┴───┐ │
│ │ │ │
│ 是 否 │
│ │ │ │
▼ ▼ ▼ ▼
栅格 拓扑矢量 简单矢量
混合表征策略
实践中，最佳方案往往是混合使用：
class HybridSpatialRepresentation:
"""
混合空间表征
设计理念：不同操作用最合适的表征
"""
def __init__(self, raster_resolution=30):
"""
Args:
raster_resolution: 栅格分辨率(米)
"""
self.raster_resolution = raster_resolution
self.vector_features = {} # 矢量要素
self.raster_surfaces = {} # 栅格表面
self.spatial_index = None # 空间索引
def add_vector_feature(self, feature_id, geometry, attributes):
"""添加矢量要素（适合精确边界）"""
import shapely.geometry as geom
self.vector_features[feature_id] = {
'geometry': geometry if isinstance(geometry, geom.base.BaseGeometry)
else geom.shape(geometry),
'attributes': attributes
}
# 更新空间索引
self._build_spatial_index()
def add_raster_surface(self, surface_id, array, transform=None):
"""添加栅格表面（适合连续场）"""
import numpy as np
self.raster_surfaces[surface_id] = {
'array': np.asarray(array),
'transform': transform,
'resolution': self.raster_resolution
}
def _build_spatial_index(self):
"""构建R-Tree空间索引"""
from rtree import index
idx = index.Index()
for i, (fid, feature) in enumerate(self.vector_features.items()):
idx.insert(i, feature['geometry'].bounds, fid)
self.spatial_index = idx
def query_by_location(self, point, buffer_distance=0):
"""
基于位置查询
使用矢量+索引：高效精确
"""
from shapely.geometry import Point
query_point = Point(point) if not isinstance(point, Point) else point
query_box = query_point.buffer(buffer_distance).bounds
# 使用空间索引快速筛选
candidates = []
for i in self.spatial_index.intersection(query_box):
fid = list(self.vector_features.keys())[i]
feature = self.vector_features[fid]
if feature['geometry'].intersects(query_point):
candidates.append({
'id': fid,
'geometry': feature['geometry'],
'attributes': feature['attributes']
})
return candidates
def extract_values_at_points(self, surface_id, points):
"""
在点上提取栅格值
使用栅格：快速采样
"""
import numpy as np
surface = self.raster_surfaces.get(surface_id)
if not surface:
raise ValueError(f"Surface {surface_id} not found")
array = surface['array']
values = []
for point in points:
# 坐标转换（如果有transform）
if surface['transform']:
# 这里简化处理，实际需要仿射变换
col = int(point[0] / self.raster_resolution)
row = int(point[1] / self.raster_resolution)
else:
col, row = int(point[0]), int(point[1])
# 边界检查
if 0 <= row < array.shape[0] and 0 <= col < array.shape[1]:
values.append(array[row, col])
else:
values.append(np.nan)
return values
def vector_to_raster(self, feature_id, value_field=None,
surface_id=None, default_value=1):
"""
矢量转栅格
适用场景：需要栅格分析（如邻域、成本距离）
"""
import numpy as np
from rasterio.features import rasterize
feature = self.vector_features.get(feature_id)
if not feature:
raise ValueError(f"Feature {feature_id} not found")
# 确定输出范围
geom = feature['geometry']
bounds = geom.bounds
width = int((bounds[2] - bounds[0]) / self.raster_resolution) + 1
height = int((bounds[3] - bounds[1]) / self.raster_resolution) + 1
# 确定栅格化值
if value_field and value_field in feature['attributes']:
value = feature['attributes'][value_field]
else:
value = default_value
# 栅格化
output_array = rasterize(
[(geom, value)],
out_shape=(height, width),
transform=None, # 简化处理
fill=0,
dtype=np.float32
)
# 保存栅格
new_surface_id = surface_id or f"{feature_id}_raster"
self.add_raster_surface(new_surface_id, output_array)
return new_surface_id
def to_graph(self, threshold_distance=None):
"""
转换为图表示
适用场景：连通性分析、路径优化
"""
import networkx as nx
from shapely.geometry import Point
G = nx.Graph()
# 添加节点（矢量要素）
for fid, feature in self.vector_features.items():
centroid = feature['geometry'].centroid
G.add_node(fid, pos=(centroid.x, centroid.y))
# 添加边（基于距离阈值）
if threshold_distance:
fids = list(self.vector_features.keys())
for i, fid1 in enumerate(fids):
for fid2 in fids[i+1:]:
geom1 = self.vector_features[fid1]['geometry']
geom2 = self.vector_features[fid2]['geometry']
dist = geom1.distance(geom2)
if dist <= threshold_distance:
G.add_edge(fid1, fid2, weight=dist)
return G
设计权衡
| 决策维度 | 选项A | 选项B | 权衡考量 |
|---|---|---|---|
| 数据结构 | 栅格 | 矢量 | 精度 vs 效率；连续 vs 离散 |
| 分辨率 | 高精度 | 低精度 | 存储成本 vs 信息保留 |
| 索引方式 | R-Tree | Quadtree | 数据分布、查询类型 |
| 单一 vs 混合 | 统一表征 | 按需选择 | 一致性 vs 灵活性 |
代码示例
示例1：栅格数据处理
"""
栅格数据处理示例
"""
import numpy as np
from typing import Tuple, List, Optional
from scipy.ndimage import convolve
class RasterProcessor:
"""栅格数据处理器"""
def __init__(self, array: np.ndarray, resolution: float = 1.0):
"""
Args:
array: 栅格数据数组
resolution: 分辨率(单位/像元)
"""
self.array = np.asarray(array)
self.resolution = resolution
self.nodata = -9999 # 无数据值
def slope(self) -> np.ndarray:
"""
计算坡度
Returns:
坡度数组(度)
"""
# Sobel算子
kernel_x = np.array([[-1, 0, 1],
[-2, 0, 2],
[-1, 0, 1]])
kernel_y = np.array([[-1, -2, -1],
[ 0, 0, 0],
[ 1, 2, 1]])
# 计算梯度
dz_dx = convolve(self.array, kernel_x) / (8 * self.resolution)
dz_dy = convolve(self.array, kernel_y) / (8 * self.resolution)
# 坡度(度)
slope = np.arctan(np.sqrt(dz_dx**2 + dz_dy**2)) * 180 / np.pi
return slope
def aspect(self) -> np.ndarray:
"""
计算坡向
Returns:
坡向数组(度, 0-360)
"""
kernel_x = np.array([[-1, 0, 1],
[-2, 0, 2],
[-1, 0, 1]])
kernel_y = np.array([[-1, -2, -1],
[ 0, 0, 0],
[ 1, 2, 1]])
dz_dx = convolve(self.array, kernel_x)
dz_dy = convolve(self.array, kernel_y)
aspect = np.arctan2(dz_dy, -dz_x) * 180 / np.pi
aspect = (90 - aspect) % 360
return aspect
def neighborhood_stats(self, radius: int = 1) -> dict:
"""
邻域统计
Args:
radius: 邻域半径(像元)
Returns:
统计指标字典
"""
from scipy.ndimage import uniform_filter
size = 2 * radius + 1
# 均值
mean = uniform_filter(self.array.astype(float), size=size, mode='reflect')
# 方差
squared_mean = uniform_filter((self.array ** 2).astype(float),
size=size, mode='reflect')
variance = squared_mean - mean ** 2
# 最值
from scipy.ndimage import maximum_filter, minimum_filter
maximum = maximum_filter(self.array, size=size, mode='reflect')
minimum = minimum_filter(self.array, size=size, mode='reflect')
return {
'mean': mean,
'std': np.sqrt(np.maximum(variance, 0)),
'max': maximum,
'min': minimum,
'range': maximum - minimum
}
def resample(self, target_resolution: float,
method: str = 'bilinear') -> 'RasterProcessor':
"""
重采样
Args:
target_resolution: 目标分辨率
method: 'nearest', 'bilinear', 'cubic'
Returns:
新的RasterProcessor
"""
from scipy.ndimage import zoom
scale_factor = self.resolution / target_resolution
if method == 'nearest':
order = 0
elif method == 'bilinear':
order = 1
elif method == 'cubic':
order = 3
else:
raise ValueError(f"Unknown method: {method}")
resampled = zoom(self.array, scale_factor, order=order)
return RasterProcessor(resampled, target_resolution)
# 使用示例
if __name__ == "__main__":
# 创建示例DEM
dem_data = np.array([
[100, 105, 110, 108, 102],
[102, 108, 115, 112, 105],
[105, 112, 120, 118, 110],
[108, 115, 122, 120, 112],
[106, 110, 115, 112, 108]
], dtype=float)
processor = RasterProcessor(dem_data, resolution=30)
# 计算坡度
slope = processor.slope()
print(f"Average slope: {np.mean(slope):.2f} degrees")
# 计算坡向
aspect = processor.aspect()
# 邻域统计
stats = processor.neighborhood_stats(radius=1)
print(f"Smoothed elevation mean: {np.mean(stats['mean']):.2f}")
示例2：矢量数据处理
"""
矢量数据处理示例
"""
import numpy as np
from typing import List, Dict, Any, Optional
from shapely.geometry import Point, LineString, Polygon, MultiPolygon
from shapely.ops import unary_union, voronoi_diagram
import geopandas as gpd
class VectorProcessor:
"""矢量数据处理器"""
def __init__(self, crs: str = "EPSG:4326"):
"""
Args:
crs: 坐标参考系统
"""
self.crs = crs
self.features = []
self.gdf = None
def add_feature(self, geometry: Any, attributes: Dict[str, Any]):
"""添加要素"""
# 确保是Shapely几何对象
if isinstance(geometry, dict):
from shapely.geometry import shape
geometry = shape(geometry)
self.features.append({
'geometry': geometry,
'attributes': attributes
})
def to_geodataframe(self) -> gpd.GeoDataFrame:
"""转换为GeoDataFrame"""
if not self.features:
return gpd.GeoDataFrame(geometry=[], crs=self.crs)
data = {
'geometry': [f['geometry'] for f in self.features],
**{k: [f['attributes'].get(k) for f in self.features]
for k in self.features[0]['attributes'].keys()}
}
self.gdf = gpd.GeoDataFrame(data, crs=self.crs)
return self.gdf
def buffer_all(self, distance: float,
resolution: int = 16) -> 'VectorProcessor':
"""
缓冲区分析
Args:
distance: 缓冲距离(与CRS单位一致)
resolution: 缓冲圆弧的分辨率
Returns:
新的VectorProcessor
"""
result = VectorProcessor(self.crs)
for feature in self.features:
buffered = feature['geometry'].buffer(
distance,
resolution=resolution
)
result.add_feature(buffered, feature['attributes'])
return result
def intersect_all(self, other: 'VectorProcessor') -> 'VectorProcessor':
"""
相交分析
Args:
other: 另一个VectorProcessor
Returns:
相交结果的新VectorProcessor
"""
result = VectorProcessor(self.crs)
for feat1 in self.features:
for feat2 in other.features:
intersection = feat1['geometry'].intersection(feat2['geometry'])
if not intersection.is_empty:
# 合并属性
merged_attrs = {
**{f"left_{k}": v for k, v in feat1['attributes'].items()},
**{f"right_{k}": v for k, v in feat2['attributes'].items()}
}
result.add_feature(intersection, merged_attrs)
return result
def centroid(self) -> List[Point]:
"""计算所有要素的质心"""
return [f['geometry'].centroid for f in self.features]
def area(self) -> List[float]:
"""计算所有面要素的面积"""
areas = []
for f in self.features:
geom = f['geometry']
if isinstance(geom, (Polygon, MultiPolygon)):
# 使用投影后的CRS计算更准确
areas.append(geom.area)
else:
areas.append(0.0)
return areas
def length(self) -> List[float]:
"""计算所有线要素的长度"""
lengths = []
for f in self.features:
geom = f['geometry']
if isinstance(geom, (LineString, Point)):
lengths.append(geom.length)
else:
lengths.append(0.0)
return lengths
def dissolve(self, by_attribute: Optional[str] = None) -> 'VectorProcessor':
"""
融合要素
Args:
by_attribute: 按此属性分组融合
Returns:
融合后的新VectorProcessor
"""
if not self.features:
return VectorProcessor(self.crs)
if by_attribute:
# 按属性分组
groups = {}
for f in self.features:
key = f['attributes'].get(by_attribute)
if key not in groups:
groups[key] = []
groups[key].append(f['geometry'])
result = VectorProcessor(self.crs)
for key, geometries in groups.items():
dissolved = unary_union(geometries)
result.add_feature(dissolved, {by_attribute: key})
return result
else:
# 全部融合
dissolved = unary_union([f['geometry'] for f in self.features])
result = VectorProcessor(self.crs)
result.add_feature(dissolved, {})
return result
# 使用示例
if __name__ == "__main__":
# 创建示例矢量数据
processor = VectorProcessor(crs="EPSG:3857") # 投影坐标系，单位米
# 添加一些面要素
processor.add_feature(
Polygon([(0, 0), (100, 0), (100, 100), (0, 100)]),
{'id': 1, 'type': 'forest'}
)
processor.add_feature(
Polygon([(120, 20), (200, 20), (200, 120), (120, 120)]),
{'id': 2, 'type': 'forest'}
)
processor.add_feature(
Polygon([(80, 80), (150, 80), (150, 150), (80, 150)]),
{'id': 3, 'type': 'wetland'}
)
# 转换为GeoDataFrame
gdf = processor.to_geodataframe()
print(f"Number of features: {len(gdf)}")
# 缓冲区分析
buffered = processor.buffer_all(distance=50)
print(f"Buffered features: {len(buffered.features)}")
# 融合分析
dissolved = processor.dissolve(by_attribute='type')
print(f"Dissolved groups: {len(dissolved.features)}")
示例3：图构建与分析
"""
空间图表示示例
"""
import numpy as np
import networkx as nx
from typing import List, Tuple, Dict, Optional
from shapely.geometry import Point, LineString
class SpatialGraphBuilder:
"""空间图构建器"""
@staticmethod
def from_points(points: List[Point],
distance_threshold: float,
distance_type: str = 'euclidean') -> nx.Graph:
"""
从点集构建图（基于距离阈值）
Args:
points: 点列表
distance_threshold: 连接阈值
distance_type: 'euclidean' 或 'manhattan'
Returns:
NetworkX图
"""
G = nx.Graph()
# 添加节点
for i, point in enumerate(points):
G.add_node(i, pos=(point.x, point.y))
# 添加边
for i in range(len(points)):
for j in range(i + 1, len(points)):
if distance_type == 'euclidean':
dist = points[i].distance(points[j])
else: # manhattan
dist = abs(points[i].x - points[j].x) + \
abs(points[i].y - points[j].y)
if dist <= distance_threshold:
G.add_edge(i, j, weight=dist)
return G
@staticmethod
def from_polygons(polygons: List[Polygon],
connectivity_type: str = 'shared_boundary') -> nx.Graph:
"""
从多边形构建图（基于拓扑关系）
Args:
polygons: 多边形列表
connectivity_type: 'shared_boundary' 或 'within_distance'
Returns:
NetworkX图
"""
G = nx.Graph()
# 添加节点
for i, poly in enumerate(polygons):
G.add_node(i, centroid=poly.centroid, area=poly.area)
# 添加边
for i in range(len(polygons)):
for j in range(i + 1, len(polygons)):
if connectivity_type == 'shared_boundary':
# 共享边界
if polygons[i].touches(polygons[j]):
# 计算共享边界长度
shared = polygons[i].intersection(polygons[j])
weight = shared.length if not shared.is_empty else 0
G.add_edge(i, j, weight=weight)
elif connectivity_type == 'within_distance':
# 距离阈值
dist = polygons[i].distance(polygons[j])
if dist >= 0 and dist <= 100: # 100米阈值
G.add_edge(i, j, weight=dist)
return G
@staticmethod
def from_network(lines: List[LineString]) -> nx.Graph:
"""
从线网络构建图（如道路网）
Args:
lines: 线列表
Returns:
NetworkX图
"""
G = nx.Graph()
# 收集所有端点
endpoints = []
for i, line in enumerate(lines):
coords = list(line.coords)
endpoints.append((i, Point(coords[0]), 'start'))
endpoints.append((i, Point(coords[-1]), 'end'))
# 构建节点（合并接近的端点）
tolerance = 1e-6
node_id = 0
point_to_node = {}
for line_idx, point, _ in endpoints:
# 检查是否已有接近的节点
matched = False
for existing_point, existing_node in point_to_node.items():
if existing_point.distance(point) < tolerance:
point_to_node[point] = existing_node
matched = True
break
if not matched:
point_to_node[point] = node_id
G.add_node(node_id, pos=(point.x, point.y))
node_id += 1
# 添加边
for i, line in enumerate(lines):
coords = list(line.coords)
start_point = Point(coords[0])
end_point = Point(coords[-1])
start_node = point_to_node[start_point]
end_node = point_to_node[end_point]
G.add_edge(start_node, end_node,
weight=line.length,
geometry=line,
edge_id=i)
return G
@staticmethod
def compute_connectivity_metrics(G: nx.Graph) -> Dict[str, float]:
"""
计算图的连通性指标
Args:
G: NetworkX图
Returns:
指标字典
"""
metrics = {}
# 基本指标
metrics['n_nodes'] = G.number_of_nodes()
metrics['n_edges'] = G.number_of_edges()
if G.number_of_nodes() == 0:
return metrics
# 连通分量
metrics['n_components'] = nx.number_connected_components(G)
# 最大连通分量
largest_cc = max(nx.connected_components(G), key=len) if G.nodes() else set()
metrics['largest_component_size'] = len(largest_cc)
metrics['largest_component_ratio'] = len(largest_cc) / G.number_of_nodes()
# 平均度
degrees = [d for n, d in G.degree()]
metrics['avg_degree'] = np.mean(degrees) if degrees else 0
# 网络密度
metrics['density'] = nx.density(G)
# 平均最短路径长度（仅当图连通时）
if nx.is_connected(G):
metrics['avg_path_length'] = nx.average_shortest_path_length(G)
metrics['diameter'] = nx.diameter(G)
else:
# 对最大连通分量计算
if largest_cc:
subgraph = G.subgraph(largest_cc)
metrics['avg_path_length_lcc'] = nx.average_shortest_path_length(subgraph)
metrics['diameter_lcc'] = nx.diameter(subgraph)
# 聚类系数
metrics['avg_clustering'] = nx.average_clustering(G)
return metrics
# 使用示例
if __name__ == "__main__":
# 示例1：从点构建图
points = [
Point(0, 0), Point(50, 0), Point(100, 0),
Point(0, 50), Point(50, 50), Point(100, 50),
Point(0, 100), Point(50, 100), Point(100, 100)
]
G_points = SpatialGraphBuilder.from_points(points, distance_threshold=80)
metrics = SpatialGraphBuilder.compute_connectivity_metrics(G_points)
print("Point Graph Metrics:", {k: round(v, 2) if isinstance(v, float) else v
for k, v in metrics.items()})
# 示例2：从多边形构建图
from shapely.geometry import box
polygons = [
box(0, 0, 50, 50),
box(50, 0, 100, 50),
box(0, 50, 50, 100),
box(50, 50, 100, 100)
]
G_poly = SpatialGraphBuilder.from_polygons(polygons)
metrics_poly = SpatialGraphBuilder.compute_connectivity_metrics(G_poly)
print("Polygon Graph Metrics:", {k: round(v, 2) if isinstance(v, float) else v
for k, v in metrics_poly.items()})
案例分析
ENAgent中的混合表征策略
ENAgent（生态网络分析智能体）在处理生态网络时采用了混合表征策略：
class ENAgentSpatialManager:
"""
ENAgent的空间管理模块
核心设计：不同数据类型用最合适的表征方式
"""
def __init__(self, raster_resolution=30):
# 矢量：生态源地（精确边界）
self.sources = VectorProcessor()
# 栅格：阻力面（连续场）
self.resistance_surface = None
# 图：生态网络（连通性分析）
self.network_graph = None
def load_sources_from_vector(self, vector_file):
"""从矢量文件加载源地"""
import geopandas as gpd
gdf = gpd.read_file(vector_file)
for _, row in gdf.iterrows():
self.sources.add_feature(
row.geometry,
{'id': row.get('id', len(self.sources.features)),
'name': row.get('name', ''),
'area': row.geometry.area}
)
def create_resistance_surface(self, land_use_raster, resistance_dict):
"""
创建阻力面（栅格）
Args:
land_use_raster: 土地利用栅格
resistance_dict: {土地类型: 阻力值}
"""
import numpy as np
# 栅格计算：矢量化重映射
land_use_array = self._read_raster(land_use_raster)
# 创建阻力面
resistance = np.zeros_like(land_use_array, dtype=float)
for land_type, resist_value in resistance_dict.items():
resistance[land_use_array == land_type] = resist_value
self.resistance_surface = resistance
return resistance
def build_network_graph(self, connectivity_threshold=5000):
"""
构建生态网络图
用于连通性分析和优化
"""
# 从源地质心构建点集
centroids = self.sources.centroid()
# 构建图
self.network_graph = SpatialGraphBuilder.from_points(
centroids,
distance_threshold=connectivity_threshold
)
return self.network_graph
def analyze_connectivity(self):
"""分析生态网络连通性"""
if self.network_graph is None:
self.build_network_graph()
return SpatialGraphBuilder.compute_connectivity_metrics(
self.network_graph
)
def _read_raster(self, raster_file):
"""读取栅格文件"""
import rasterio
with rasterio.open(raster_file) as src:
return src.read(1)
关键设计决策：
- 源地用矢量：需要精确边界和面积计算
- 阻力面用栅格：需要邻域分析和成本距离计算
- 网络用图：需要连通性分析和路径优化
这种混合策略充分发挥了各种表征的优势。
反思与延伸
思考问题
-
尺度效应：在不同分析尺度下，同一空间现象的表征会有什么变化？
-
不确定性传播：从一种表征转换到另一种（如矢量转栅格）时，不确定性如何传播？
-
动态数据：时变的空间数据应该如何表征？
-
三维扩展：如何将这些二维表征扩展到三维空间？
-
存储与效率：当数据量达到TB级别时，表征策略需要做什么调整？
延伸阅读
- "Fundamentals of Geographic Information Systems" (Demers) - 空间数据模型基础
- "Geographic Information Systems and Science" (Longley) - 第3-4章
- "Spatial Databases" (Rigaux) - 空间索引原理
- Shapely Documentation - Python几何操作
- Rasterio Documentation - Python栅格处理
关键要点
-
空间表征是空间智能的基础：选择合适的表征方式直接影响后续分析的效率和准确性
-
栅格和矢量各有优势：栅格适合连续场和邻域分析，矢量适合离散要素和精确边界
-
图表示连接空间与算法：将空间问题转化为图问题，可以应用丰富的图算法
-
多尺度是现实需求：空间智能系统必须能处理不同尺度的数据和分析
-
混合策略往往是最佳选择：ENAgent的实践表明，根据数据类型和操作需求选择表征方式是最有效的