"""
空间表征示例 (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()