# 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 (简单高效)                                   │
│     ┌───┬───┬───┬───┐                                      │
│     │   │ ● │   │   │  规则网格分桶                        │
│     ├───┼───┼───┼───┤  适合：均匀分布数据                  │
│     │   │   │ ● │   │                                      │
│     └───┴───┴───┴───┘                                      │
│                                                             │
└─────────────────────────────────────────────────────────────┘
```

---

## 设计原理

### 表征选择决策树

```
                        开始
                         │
                         ▼
                  数据是什么类型？
                  ┌──────┴──────┐
                  │             │
             连续场/表面      离散要素
             (高程、温度)    (边界、道路)
                  │             │
                  ▼             ▼
              栅格优先      需要精确边界？
                  │          ┌──┴──┐
                  │          │     │
                  │         是    否
                  │          │     │
                  │          ▼     ▼
                  │       矢量   栅格也可
                  │          │     │
                  │          ▼     │
                  │    需要拓扑分析？
                  │          │     │
                  │      ┌───┴───┐ │
                  │      │       │ │
                  │     是      否 │
                  │      │       │ │
                  ▼      ▼       ▼ ▼
              栅格   拓扑矢量  简单矢量
```

### 混合表征策略

实践中，最佳方案往往是混合使用：

```python
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：栅格数据处理

```python
"""
栅格数据处理示例
"""
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：矢量数据处理

```python
"""
矢量数据处理示例
"""
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：图构建与分析

```python
"""
空间图表示示例
"""
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（生态网络分析智能体）在处理生态网络时采用了混合表征策略：

```python
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)
```

**关键设计决策**：

1. **源地用矢量**：需要精确边界和面积计算
2. **阻力面用栅格**：需要邻域分析和成本距离计算
3. **网络用图**：需要连通性分析和路径优化

这种混合策略充分发挥了各种表征的优势。

---

## 反思与延伸

### 思考问题

1. **尺度效应**：在不同分析尺度下，同一空间现象的表征会有什么变化？

2. **不确定性传播**：从一种表征转换到另一种（如矢量转栅格）时，不确定性如何传播？

3. **动态数据**：时变的空间数据应该如何表征？

4. **三维扩展**：如何将这些二维表征扩展到三维空间？

5. **存储与效率**：当数据量达到TB级别时，表征策略需要做什么调整？

### 延伸阅读

- **"Fundamentals of Geographic Information Systems"** (Demers) - 空间数据模型基础
- **"Geographic Information Systems and Science"** (Longley) - 第3-4章
- **"Spatial Databases"** (Rigaux) - 空间索引原理
- Shapely Documentation - Python几何操作
- Rasterio Documentation - Python栅格处理

---

## 关键要点

1. **空间表征是空间智能的基础**：选择合适的表征方式直接影响后续分析的效率和准确性

2. **栅格和矢量各有优势**：栅格适合连续场和邻域分析，矢量适合离散要素和精确边界

3. **图表示连接空间与算法**：将空间问题转化为图问题，可以应用丰富的图算法

4. **多尺度是现实需求**：空间智能系统必须能处理不同尺度的数据和分析

5. **混合策略往往是最佳选择**：ENAgent的实践表明，根据数据类型和操作需求选择表征方式是最有效的