第10章：几何属性与度量

几何属性与度量是空间数据分析的基础。GeoPandas 通过 Shapely 库提供了丰富的几何属性查询和度量计算功能，可以方便地获取几何对象的类型、面积、长度、边界框、质心等信息。本章将全面介绍 GeoPandas 中的几何属性与度量方法。

10.1 几何类型属性

10.1.1 geom_type 属性

geom_type 属性返回每个几何对象的类型名称（字符串）：

import geopandas as gpd
from shapely.geometry import Point, LineString, Polygon, MultiPoint

# 创建包含不同几何类型的 GeoDataFrame
gdf = gpd.GeoDataFrame(
    {"name": ["点", "线", "面", "多点"]},
    geometry=[
        Point(0, 0),
        LineString([(0, 0), (1, 1), (2, 0)]),
        Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
        MultiPoint([(0, 0), (1, 1)])
    ]
)

# 查看每个要素的几何类型
print(gdf.geom_type)

输出：

0         Point
1    LineString
2       Polygon
3    MultiPoint
dtype: object

10.1.2 所有几何类型

GeoPandas/Shapely 支持的几何类型：

类型	说明	geom_type 值
点	单个坐标点	`Point`
线	有序点序列	`LineString`
线环	首尾相连的线	`LinearRing`
多边形	由外环和可选内环组成	`Polygon`
多点	点的集合	`MultiPoint`
多线	线的集合	`MultiLineString`
多面	多边形的集合	`MultiPolygon`
几何集合	任意几何的集合	`GeometryCollection`

10.1.3 几何类型判断

# 判断是否为特定类型
is_point = gdf.geom_type == "Point"
is_polygon = gdf.geom_type.isin(["Polygon", "MultiPolygon"])
is_line = gdf.geom_type.isin(["LineString", "MultiLineString"])

print(f"点要素: {is_point.sum()} 个")
print(f"面要素: {is_polygon.sum()} 个")
print(f"线要素: {is_line.sum()} 个")

# 按类型筛选
gdf_points = gdf[gdf.geom_type == "Point"]
gdf_polygons = gdf[gdf.geom_type.isin(["Polygon", "MultiPolygon"])]

10.1.4 统计几何类型分布

# 统计各类型数量
type_counts = gdf.geom_type.value_counts()
print("几何类型分布:")
print(type_counts)

# 检查是否为单一类型
if len(gdf.geom_type.unique()) == 1:
    print(f"数据集为单一类型: {gdf.geom_type.unique()[0]}")
else:
    print("数据集包含混合类型")

10.2 面积计算

10.2.1 area 属性

area 属性返回每个几何对象的面积：

import geopandas as gpd
from shapely.geometry import Polygon, Point, LineString

gdf = gpd.GeoDataFrame(
    {"name": ["正方形", "三角形", "点", "线"]},
    geometry=[
        Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
        Polygon([(0, 0), (2, 0), (1, 1)]),
        Point(0, 0),
        LineString([(0, 0), (1, 1)])
    ]
)

print(gdf.area)

输出：

0    1.0
1    1.0
2    0.0
3    0.0
dtype: float64

注意： 点和线的面积始终为 0。

10.2.2 投影坐标系的重要性

这是最容易出错的地方之一！ 在地理坐标系（如 EPSG:4326）下，.area 返回的值的单位是平方度，这个数值没有实际物理意义。

import geopandas as gpd
from shapely.geometry import box

# 创建一个近似北京市区范围的矩形
beijing_approx = gpd.GeoDataFrame(
    {"name": ["北京近似"]},
    geometry=[box(116.0, 39.6, 116.8, 40.2)],
    crs="EPSG:4326"
)

# ❌ 错误：在地理坐标系下计算面积
area_degrees = beijing_approx.area.iloc[0]
print(f"地理坐标系面积: {area_degrees:.6f} 平方度（无意义）")

# ✅ 正确：转换到投影坐标系后计算
utm_crs = beijing_approx.estimate_utm_crs()
beijing_utm = beijing_approx.to_crs(utm_crs)
area_sqm = beijing_utm.area.iloc[0]
print(f"投影坐标系面积: {area_sqm:,.0f} 平方米")
print(f"投影坐标系面积: {area_sqm / 1e6:,.2f} 平方公里")

10.2.3 面积计算的最佳实践

def calculate_area(gdf, unit="sqm"):
    """
    计算面积，自动处理坐标系转换

    参数:
        gdf: GeoDataFrame
        unit: 面积单位 - "sqm"(平方米), "sqkm"(平方公里),
              "ha"(公顷), "mu"(亩)
    """
    # 如果是地理坐标系，先转换到投影坐标系
    if gdf.crs and gdf.crs.is_geographic:
        utm_crs = gdf.estimate_utm_crs()
        gdf_proj = gdf.to_crs(utm_crs)
    else:
        gdf_proj = gdf

    area_sqm = gdf_proj.area

    # 单位转换
    conversions = {
        "sqm": 1,
        "sqkm": 1e-6,
        "ha": 1e-4,
        "mu": 1 / 666.667,  # 1亩 ≈ 666.667 平方米
        "acre": 1 / 4046.86
    }

    factor = conversions.get(unit, 1)
    return area_sqm * factor

# 使用示例
gdf["area_sqkm"] = calculate_area(gdf, unit="sqkm")
gdf["area_mu"] = calculate_area(gdf, unit="mu")
print(gdf[["name", "area_sqkm", "area_mu"]])

10.3 长度计算

10.3.1 length 属性

length 属性返回几何对象的长度或周长：

import geopandas as gpd
from shapely.geometry import LineString, Polygon, Point

gdf = gpd.GeoDataFrame(
    {"name": ["线段", "折线", "多边形", "点"]},
    geometry=[
        LineString([(0, 0), (3, 4)]),              # 直线段
        LineString([(0, 0), (1, 0), (1, 1)]),      # 折线
        Polygon([(0, 0), (4, 0), (4, 3), (0, 3)]), # 矩形
        Point(0, 0)                                  # 点
    ]
)

print(gdf.length)

输出：

0     5.0      # sqrt(3² + 4²) = 5
1     2.0      # 1 + 1 = 2
2    14.0      # 矩形周长 = 2*(4+3) = 14
3     0.0      # 点的长度为 0
dtype: float64

10.3.2 线要素与多边形周长

几何类型	length 含义
Point	0
LineString	线的总长度
LinearRing	环的周长
Polygon	外环周长（不含内环）
MultiLineString	所有线段长度之和
MultiPolygon	所有多边形外环周长之和

from shapely.geometry import Polygon

# 带洞的多边形
polygon_with_hole = Polygon(
    [(0, 0), (10, 0), (10, 10), (0, 10)],  # 外环
    [[(2, 2), (4, 2), (4, 4), (2, 4)]]     # 内环（洞）
)

# 注意：Polygon 的 length 返回的是外环的周长
print(f"多边形周长（外环）: {polygon_with_hole.length}")
# 40.0

# 如果需要包括内环的总长度
exterior_length = polygon_with_hole.exterior.length
interior_length = sum(ring.length for ring in polygon_with_hole.interiors)
total_length = exterior_length + interior_length
print(f"外环周长: {exterior_length}")
print(f"内环周长: {interior_length}")
print(f"总周长: {total_length}")

10.3.3 长度计算同样需要投影坐标系

# ❌ 地理坐标系下的长度单位是"度"，无意义
gdf_geo = gpd.GeoDataFrame(
    geometry=[LineString([(116.0, 39.0), (117.0, 39.0)])],
    crs="EPSG:4326"
)
print(f"地理坐标系长度: {gdf_geo.length.iloc[0]:.4f} 度")

# ✅ 转换到投影坐标系后计算
gdf_utm = gdf_geo.to_crs(gdf_geo.estimate_utm_crs())
print(f"投影坐标系长度: {gdf_utm.length.iloc[0]:,.0f} 米")
# 约 85,394 米（在纬度 39° 处，1° 经度约 85.4 公里）

10.4 边界框

10.4.1 bounds — 每个要素的边界框

bounds 属性返回每个几何对象的最小边界矩形（MBR）：

import geopandas as gpd
from shapely.geometry import Polygon, LineString, Point

gdf = gpd.GeoDataFrame(
    {"name": ["三角形", "折线", "点"]},
    geometry=[
        Polygon([(0, 0), (4, 0), (2, 3)]),
        LineString([(1, 1), (5, 3), (3, 5)]),
        Point(2, 2)
    ]
)

print(gdf.bounds)

输出：

   minx  miny  maxx  maxy
0   0.0   0.0   4.0   3.0
1   1.0   1.0   5.0   5.0
2   2.0   2.0   2.0   2.0

10.4.2 bounds 的返回格式

bounds 返回一个 DataFrame，包含四列：

列名	说明
`minx`	最小 X 坐标（最西端经度）
`miny`	最小 Y 坐标（最南端纬度）
`maxx`	最大 X 坐标（最东端经度）
`maxy`	最大 Y 坐标（最北端纬度）

# 将边界框信息添加到 GeoDataFrame
bounds = gdf.bounds
gdf["bbox_width"] = bounds["maxx"] - bounds["minx"]
gdf["bbox_height"] = bounds["maxy"] - bounds["miny"]
gdf["bbox_area"] = gdf["bbox_width"] * gdf["bbox_height"]

print(gdf[["name", "bbox_width", "bbox_height", "bbox_area"]])

10.4.3 total_bounds — 整体范围

total_bounds 属性返回整个 GeoDataFrame 的总边界框（numpy 数组）：

# 所有要素的总范围
total = gdf.total_bounds
print(f"总范围: {total}")
# [minx, miny, maxx, maxy]

print(f"X 范围: {total[0]:.4f} - {total[2]:.4f}")
print(f"Y 范围: {total[1]:.4f} - {total[3]:.4f}")

# 用途示例：创建覆盖所有数据的边界矩形
from shapely.geometry import box
bbox_polygon = box(*total)
print(f"边界矩形: {bbox_polygon}")

10.4.4 实际应用示例

import geopandas as gpd

# 读取中国省级数据
# gdf = gpd.read_file("china_provinces.shp")

# 获取每个省份的边界范围
# province_bounds = gdf.bounds
# province_bounds["province"] = gdf["name"]
# print(province_bounds.head())

# 获取中国的整体范围
# china_extent = gdf.total_bounds
# print(f"中国范围:")
# print(f"  经度: {china_extent[0]:.2f}° - {china_extent[2]:.2f}°")
# print(f"  纬度: {china_extent[1]:.2f}° - {china_extent[3]:.2f}°")
# 约为: 经度 73.5° - 135.1°, 纬度 18.2° - 53.6°

# 使用总范围进行空间查询
# other_gdf = gpd.read_file("data.shp", bbox=tuple(gdf.total_bounds))

10.5 质心

10.5.1 centroid 属性

centroid 属性返回每个几何对象的几何中心点：

import geopandas as gpd
from shapely.geometry import Polygon, LineString

gdf = gpd.GeoDataFrame(
    {"name": ["正方形", "三角形", "线段"]},
    geometry=[
        Polygon([(0, 0), (4, 0), (4, 4), (0, 4)]),
        Polygon([(0, 0), (6, 0), (3, 4)]),
        LineString([(0, 0), (10, 0)])
    ]
)

centroids = gdf.centroid
print("质心坐标:")
for name, centroid in zip(gdf["name"], centroids):
    print(f"  {name}: ({centroid.x:.2f}, {centroid.y:.2f})")

输出：

质心坐标:
  正方形: (2.00, 2.00)
  三角形: (3.00, 1.33)
  线段: (5.00, 0.00)

10.5.2 质心的特点与注意事项

from shapely.geometry import Polygon

# 注意：质心可能不在几何内部！
# 例如凹多边形（L 形）
l_shape = Polygon([
    (0, 0), (2, 0), (2, 1), (1, 1),
    (1, 3), (0, 3)
])

centroid = l_shape.centroid
print(f"L 形质心: ({centroid.x:.2f}, {centroid.y:.2f})")
print(f"质心在几何内部: {l_shape.contains(centroid)}")
# 质心可能落在 L 形的凹陷处，不在几何内部

10.5.3 为面要素添加质心标注

# 常见用途：在地图上为面要素添加标注点
gdf["centroid_x"] = gdf.centroid.x
gdf["centroid_y"] = gdf.centroid.y

# 创建质心点图层
gdf_centroids = gdf.copy()
gdf_centroids = gdf_centroids.set_geometry(gdf.centroid)
print(gdf_centroids.geom_type.unique())
# ['Point']

10.6 代表点

10.6.1 representative_point() 方法

representative_point() 返回一个保证在几何对象内部的点。与质心不同，代表点一定在几何内部。

import geopandas as gpd
from shapely.geometry import Polygon

# 创建 L 形（凹多边形）
l_shape = Polygon([
    (0, 0), (4, 0), (4, 2), (2, 2),
    (2, 4), (0, 4)
])

gdf = gpd.GeoDataFrame(
    {"name": ["L形"]},
    geometry=[l_shape]
)

# 对比质心和代表点
centroid = gdf.centroid.iloc[0]
rep_point = gdf.representative_point().iloc[0]

print(f"质心: ({centroid.x:.2f}, {centroid.y:.2f})")
print(f"代表点: ({rep_point.x:.2f}, {rep_point.y:.2f})")
print(f"质心在内部: {l_shape.contains(centroid)}")
print(f"代表点在内部: {l_shape.contains(rep_point)}")

10.6.2 质心 vs 代表点

特性	centroid	representative_point()
计算方式	几何重心	算法保证在内部
是否在内部	不一定	保证在内部
凸几何	一定在内部	在内部
凹几何	可能在外部	在内部
多边形有洞	可能在洞内	在实体部分
适用场景	统计分析	标注、定位

# 实际建议：
# - 统计分析（如加权平均中心）→ 使用 centroid
# - 地图标注 → 使用 representative_point()
# - 空间索引点 → 使用 representative_point()

# 为每个面要素生成标注点
gdf["label_point"] = gdf.representative_point()

10.7 空几何判断

10.7.1 is_empty 属性

is_empty 属性检查几何对象是否为空（没有坐标点的几何）：

import geopandas as gpd
from shapely.geometry import Point, LineString
from shapely import wkt

# 创建包含空几何的数据
gdf = gpd.GeoDataFrame(
    {"name": ["正常点", "空几何", "正常线"]},
    geometry=[
        Point(1, 2),
        wkt.loads("GEOMETRYCOLLECTION EMPTY"),
        LineString([(0, 0), (1, 1)])
    ]
)

print(gdf.is_empty)

输出：

0    False
1     True
2    False
dtype: bool

10.7.2 处理空几何

# 过滤掉空几何
gdf_valid = gdf[~gdf.is_empty]
print(f"移除空几何后: {len(gdf_valid)} 条记录")

# 统计空几何数量
empty_count = gdf.is_empty.sum()
print(f"空几何数量: {empty_count}")

# 替换空几何
# gdf.loc[gdf.is_empty, "geometry"] = Point(0, 0)  # 用默认点替代

10.8 有效性判断

10.8.1 is_valid 属性

is_valid 属性检查几何是否满足 OGC Simple Features 规范的有效性规则：

import geopandas as gpd
from shapely.geometry import Polygon

# 创建有效和无效的几何
valid_polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])

# "蝴蝶结"（自相交多边形）— 无效
invalid_polygon = Polygon([(0, 0), (1, 1), (1, 0), (0, 1)])

gdf = gpd.GeoDataFrame(
    {"name": ["有效", "无效"]},
    geometry=[valid_polygon, invalid_polygon]
)

print(gdf.is_valid)

输出：

0     True
1    False
dtype: bool

10.8.2 is_valid_reason — 获取无效原因

from shapely.validation import explain_validity

# 检查无效原因
for idx, row in gdf.iterrows():
    if not row.geometry.is_valid:
        reason = explain_validity(row.geometry)
        print(f"要素 {row['name']}: {reason}")

# 在 GeoSeries 级别
# GeoPandas 提供了 is_valid_reason 属性（取决于版本）
# print(gdf.geometry.is_valid_reason)

10.8.3 常见的无效几何类型

无效类型	说明	示例
自相交	多边形边界与自身相交	蝴蝶结形状
重复点	连续的相同坐标点	退化的线段
环方向	外环和内环方向不正确	内环应为顺时针
非闭合环	多边形的环首尾不相连	—
嵌套洞	内环在外环之外	—
零面积	退化为点或线的多边形	—

10.8.4 修复无效几何

from shapely.validation import make_valid

# 方法1：使用 make_valid()
gdf["geometry"] = gdf.geometry.make_valid()

# 方法2：使用 buffer(0) 修复（经典方法）
gdf["geometry"] = gdf.geometry.buffer(0)

# 方法3：逐个检查和修复
def fix_geometry(geom):
    """尝试修复无效几何"""
    if geom is None or geom.is_empty:
        return geom
    if not geom.is_valid:
        # 先尝试 make_valid
        fixed = make_valid(geom)
        if fixed.is_valid:
            return fixed
        # 再尝试 buffer(0)
        fixed = geom.buffer(0)
        return fixed
    return geom

gdf["geometry"] = gdf.geometry.apply(fix_geometry)

# 验证修复结果
print(f"所有几何有效: {gdf.is_valid.all()}")

10.9 简单性判断

10.9.1 is_simple 属性

is_simple 检查几何是否是"简单的"（不自相交）：

import geopandas as gpd
from shapely.geometry import LineString

# 简单线（不自相交）
simple_line = LineString([(0, 0), (1, 1), (2, 0)])

# 非简单线（自相交，形成8字形）
complex_line = LineString([(0, 0), (1, 1), (1, 0), (0, 1)])

gdf = gpd.GeoDataFrame(
    {"name": ["简单线", "自相交线"]},
    geometry=[simple_line, complex_line]
)

print(gdf.is_simple)

输出：

0     True
1    False
dtype: bool

10.9.2 简单性 vs 有效性

概念	适用类型	说明
简单性（is_simple）	主要针对线	线不自相交
有效性（is_valid）	主要针对面	符合 OGC 规范

# 对于多边形：有效则一定简单
# 对于线：简单不一定有效（线没有有效性概念）
# 点始终是简单的

10.10 环判断

10.10.1 is_ring 属性

is_ring 检查线要素是否构成闭合环（首尾相连且简单）：

import geopandas as gpd
from shapely.geometry import LineString

# 闭合环（首尾相连）
ring = LineString([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])

# 非闭合线
open_line = LineString([(0, 0), (1, 0), (1, 1)])

# 自相交的闭合线（不是环，因为不简单）
self_intersecting = LineString([(0, 0), (1, 1), (1, 0), (0, 1), (0, 0)])

gdf = gpd.GeoDataFrame(
    {"name": ["闭合环", "非闭合", "自相交闭合"]},
    geometry=[ring, open_line, self_intersecting]
)

print(gdf.is_ring)

输出：

0     True
1    False
2    False
dtype: bool

is_ring 要求同时满足：闭合（is_closed）且简单（is_simple）

10.10.2 is_closed 属性

is_closed 仅检查首尾是否相连，不要求简单性：

# is_closed 只检查首尾是否相连
print("is_closed:")
for name, geom in zip(gdf["name"], gdf.geometry):
    if hasattr(geom, 'is_closed'):
        print(f"  {name}: {geom.is_closed}")

# is_ring = is_closed AND is_simple

10.11 维度信息

10.11.1 has_z 属性

has_z 检查几何是否包含 Z（高程）坐标：

import geopandas as gpd
from shapely.geometry import Point

# 二维点
point_2d = Point(116.4, 39.9)

# 三维点（含高程）
point_3d = Point(116.4, 39.9, 50.0)

gdf = gpd.GeoDataFrame(
    {"name": ["2D点", "3D点"]},
    geometry=[point_2d, point_3d]
)

print(gdf.has_z)

输出：

0    False
1     True
dtype: bool

10.11.2 Z 值的使用场景

# Z 值常用于：
# - 地形高程数据
# - 建筑物高度
# - 地下管线深度
# - 三维城市模型

# 获取 Z 值
point_3d = Point(116.4, 39.9, 50.0)
print(f"X: {point_3d.x}")
print(f"Y: {point_3d.y}")
print(f"Z: {point_3d.z}")

# 注意：大多数二维空间分析会忽略 Z 值
# 面积、距离等计算通常只使用 X、Y 坐标

10.11.3 has_m 属性

M 值（Measure）是除 X、Y、Z 之外的第四个坐标维度，常用于线性参考：

# 注意：Shapely 对 M 值的支持有限
# has_m 属性在较新版本的 Shapely 中可用

# M 值常用于：
# - 道路里程桩
# - 管线测量
# - 时间序列轨迹

# 检查是否支持
# print(gdf.has_m)  # 取决于 Shapely 版本

10.12 坐标计数

10.12.1 count_coordinates()

count_coordinates() 返回每个几何的坐标点总数：

import geopandas as gpd
from shapely.geometry import Point, LineString, Polygon
import shapely

# 创建不同复杂度的几何
gdf = gpd.GeoDataFrame(
    {"name": ["点", "线(3点)", "三角形", "矩形"]},
    geometry=[
        Point(0, 0),                                    # 1 个坐标
        LineString([(0, 0), (1, 1), (2, 0)]),           # 3 个坐标
        Polygon([(0, 0), (1, 0), (0, 1)]),              # 4 个坐标（闭合）
        Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])       # 5 个坐标（闭合）
    ]
)

# 使用 shapely 的 count_coordinates 函数
coord_counts = shapely.count_coordinates(gdf.geometry)
print("坐标数量:")
for name, count in zip(gdf["name"], coord_counts):
    print(f"  {name}: {count}")

10.12.2 count_geometries()

count_geometries() 返回多几何中子几何的数量：

from shapely.geometry import MultiPoint, MultiPolygon, Point, Polygon
import shapely

gdf = gpd.GeoDataFrame(
    {"name": ["单点", "3个多点", "2个多面"]},
    geometry=[
        Point(0, 0),
        MultiPoint([(0, 0), (1, 1), (2, 2)]),
        MultiPolygon([
            Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
            Polygon([(2, 2), (3, 2), (3, 3), (2, 3)])
        ])
    ]
)

geom_counts = shapely.count_geometries(gdf.geometry)
print("子几何数量:")
for name, count in zip(gdf["name"], geom_counts):
    print(f"  {name}: {count}")

输出：

子几何数量:
  单点: 1
  3个多点: 3
  2个多面: 2

10.12.3 count_interior_rings()

count_interior_rings() 返回多边形中内环（洞）的数量：

from shapely.geometry import Polygon
import shapely

# 无洞多边形
simple_poly = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])

# 有 1 个洞的多边形
one_hole = Polygon(
    [(0, 0), (10, 0), (10, 10), (0, 10)],
    [[(2, 2), (4, 2), (4, 4), (2, 4)]]
)

# 有 2 个洞的多边形
two_holes = Polygon(
    [(0, 0), (10, 0), (10, 10), (0, 10)],
    [
        [(1, 1), (3, 1), (3, 3), (1, 3)],
        [(5, 5), (7, 5), (7, 7), (5, 7)]
    ]
)

gdf = gpd.GeoDataFrame(
    {"name": ["无洞", "1个洞", "2个洞"]},
    geometry=[simple_poly, one_hole, two_holes]
)

ring_counts = shapely.count_interior_rings(gdf.geometry)
print("内环数量:")
for name, count in zip(gdf["name"], ring_counts):
    print(f"  {name}: {count}")

输出：

内环数量:
  无洞: 0
  1个洞: 1
  2个洞: 2

10.13 精度控制

10.13.1 get_precision() — 获取精度

import shapely
from shapely.geometry import Point

# 默认精度为 0（无限精度/浮点精度）
point = Point(116.40723456789, 39.90456789123)
precision = shapely.get_precision(point)
print(f"默认精度: {precision}")
# 0.0

10.13.2 set_precision() — 设置精度

import shapely
from shapely.geometry import Point, Polygon
import geopandas as gpd

# 降低精度（四舍五入坐标）
point = Point(116.40723456789, 39.90456789123)

# 精度为 0.001（约 100 米级别）
point_low = shapely.set_precision(point, grid_size=0.001)
print(f"原始: ({point.x}, {point.y})")
print(f"降精: ({point_low.x}, {point_low.y})")

输出：

原始: (116.40723456789, 39.90456789123)
降精: (116.407, 39.905)

10.13.3 对 GeoSeries 设置精度

# 对整个 GeoSeries 设置精度
gdf = gpd.GeoDataFrame(
    {"name": ["A", "B"]},
    geometry=[
        Point(116.40723456789, 39.90456789123),
        Point(121.47456789012, 31.23056789012)
    ],
    crs="EPSG:4326"
)

# 将坐标精度降低到小数点后 4 位
gdf["geometry"] = shapely.set_precision(gdf.geometry.values, grid_size=0.0001)
print(gdf.geometry)

10.13.4 精度设置的应用场景

# 1. 减小文件大小（较低精度 → 更少的存储空间）
# 地理坐标系下，不同精度对应的近似空间精度：
# 0.1      → 约 11 km
# 0.01     → 约 1.1 km
# 0.001    → 约 110 m
# 0.0001   → 约 11 m
# 0.00001  → 约 1.1 m
# 0.000001 → 约 0.11 m

# 2. 避免浮点数精度问题导致的拓扑错误
# 3. 匹配不同数据源的精度级别

10.14 面积与长度的单位问题

10.14.1 单位取决于坐标系

这是空间分析中最重要的概念之一：面积和长度的单位由坐标参考系统决定。

CRS 类型	area 单位	length 单位	是否有实际意义
地理坐标系（度）	平方度	度	❌ 无实际意义
投影坐标系（米）	平方米	米	✅ 有实际意义
投影坐标系（英尺）	平方英尺	英尺	✅ 有实际意义

10.14.2 完整的单位处理流程

import geopandas as gpd
from shapely.geometry import Polygon

# 创建数据（WGS84 地理坐标系）
gdf = gpd.GeoDataFrame(
    {"name": ["区域A", "区域B"]},
    geometry=[
        Polygon([(116.0, 39.5), (116.5, 39.5), (116.5, 40.0), (116.0, 40.0)]),
        Polygon([(121.0, 31.0), (121.5, 31.0), (121.5, 31.5), (121.0, 31.5)])
    ],
    crs="EPSG:4326"
)

# 步骤1：检查当前 CRS
print(f"CRS: {gdf.crs}")
print(f"是地理坐标系: {gdf.crs.is_geographic}")

# 步骤2：转换到投影坐标系
# 方法A：使用 estimate_utm_crs()
utm_crs = gdf.estimate_utm_crs()
gdf_utm = gdf.to_crs(utm_crs)

# 方法B：使用等面积投影（更适合面积计算）
# gdf_aea = gdf.to_crs("+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +datum=WGS84 +units=m")

# 步骤3：计算面积和周长
gdf["area_sqm"] = gdf_utm.area
gdf["area_sqkm"] = gdf_utm.area / 1e6
gdf["perimeter_m"] = gdf_utm.length
gdf["perimeter_km"] = gdf_utm.length / 1000

print(gdf[["name", "area_sqkm", "perimeter_km"]])

10.14.3 不同投影下面积的差异

import geopandas as gpd
from shapely.geometry import Polygon

# 创建一个近似北京市区的矩形
beijing = gpd.GeoDataFrame(
    {"name": ["北京"]},
    geometry=[Polygon([
        (116.0, 39.6), (116.8, 39.6),
        (116.8, 40.2), (116.0, 40.2)
    ])],
    crs="EPSG:4326"
)

# 使用不同投影计算面积
projections = {
    "UTM 50N": "EPSG:32650",
    "Web Mercator": "EPSG:3857",
    "Albers 等面积": "+proj=aea +lat_1=25 +lat_2=47 +lon_0=105 +datum=WGS84 +units=m"
}

print("不同投影下的面积对比:")
for name, crs in projections.items():
    area = beijing.to_crs(crs).area.iloc[0] / 1e6
    print(f"  {name}: {area:,.2f} 平方公里")

# 建议：
# - 面积计算优先使用等面积投影
# - UTM 在小范围内也足够准确
# - Web Mercator 不适合面积计算（面积变形严重）

10.14.4 单位转换参考表

转换	系数
平方米 → 平方公里	× 0.000001 (÷ 1,000,000)
平方米 → 公顷	× 0.0001 (÷ 10,000)
平方米 → 亩	÷ 666.667
平方米 → 英亩	÷ 4,046.86
平方米 → 平方英里	÷ 2,589,988
米 → 公里	÷ 1,000
米 → 英里	÷ 1,609.34
米 → 海里	÷ 1,852

# 单位转换工具函数
class AreaConverter:
    """面积单位转换器"""
    FACTORS = {
        "sqm": 1,
        "sqkm": 1e-6,
        "ha": 1e-4,
        "mu": 1 / 666.667,
        "acre": 1 / 4046.86,
        "sqmi": 1 / 2589988
    }

    @staticmethod
    def convert(value_sqm, to_unit):
        """将平方米转换为目标单位"""
        factor = AreaConverter.FACTORS.get(to_unit)
        if factor is None:
            raise ValueError(f"不支持的单位: {to_unit}")
        return value_sqm * factor

class LengthConverter:
    """长度单位转换器"""
    FACTORS = {
        "m": 1,
        "km": 0.001,
        "mi": 1 / 1609.34,
        "nmi": 1 / 1852,
        "ft": 3.28084
    }

    @staticmethod
    def convert(value_m, to_unit):
        factor = LengthConverter.FACTORS.get(to_unit)
        if factor is None:
            raise ValueError(f"不支持的单位: {to_unit}")
        return value_m * factor

# 使用示例
area_sqm = 16400000  # 16.4 平方公里
print(f"面积: {AreaConverter.convert(area_sqm, 'sqkm'):.2f} 平方公里")
print(f"面积: {AreaConverter.convert(area_sqm, 'ha'):.0f} 公顷")
print(f"面积: {AreaConverter.convert(area_sqm, 'mu'):.0f} 亩")

10.15 本章小结

本章全面介绍了 GeoPandas 中的几何属性与度量功能，涵盖了以下核心内容：

主题	属性/方法	关键要点
几何类型	`geom_type`	返回字符串类型名，支持类型过滤
面积	`area`	单位取决于 CRS，必须使用投影坐标系
长度	`length`	线的长度或面的周长，同样需要投影坐标系
边界框	`bounds` / `total_bounds`	单个要素或整体的最小外接矩形
质心	`centroid`	几何中心，可能不在几何内部
代表点	`representative_point()`	保证在几何内部，适合标注
空几何	`is_empty`	检查几何是否为空
有效性	`is_valid`	OGC 规范有效性，可用 make_valid() 修复
简单性	`is_simple`	检查线是否自相交
环	`is_ring`	闭合且简单的线
Z/M 维度	`has_z` / `has_m`	是否包含高程或度量维度
坐标计数	`count_coordinates()`	统计坐标点数量
精度控制	`set_precision()`	控制坐标精度，减小存储

核心要诀：

面积和长度计算必须在投影坐标系下进行，地理坐标系下的值无意义
面积计算推荐使用等面积投影或 UTM 投影
标注定位使用 representative_point() 而非 centroid
处理数据前检查有效性（is_valid），必要时修复（make_valid()）
不同投影下的面积/长度结果会有差异，选择合适的投影很重要

下一章我们将学习 GeoPandas 中的几何操作方法，包括缓冲区分析、叠加分析等空间分析功能。

此内容由惯性聚合(RSS阅读器)自动聚合整理，仅供阅读参考。原文来自 — 版权归原作者所有。

推荐订阅源

博客园 - 我才是银古