znlgis 博客

GIS开发与技术分享

第18章：裁剪与掩膜

裁剪（Clip）与掩膜（Mask）是地理空间数据处理中常用的操作，用于将数据限制在特定的空间范围内。裁剪保留落在指定区域内的要素，掩膜则保留指定区域外的要素。本章将详细介绍 GeoPandas 中的 clip() 函数及掩膜操作的实现方法。

18.1 裁剪与掩膜概述

18.1.1 基本概念

裁剪（Clip）与掩膜（Mask）是一对互补操作：

  • 裁剪：保留落在裁剪区域内部的要素或要素部分
  • 掩膜：保留落在掩膜区域外部的要素或要素部分
原始数据          裁剪区域          裁剪结果          掩膜结果
+--------+       +----+          +----+           +--+  +--+
|* *  * *|       |    |          |* * |           |* |  |* |
| *  *   |       |    |          | *  |           |  |  |  |
|  * * * |       +----+          +----+           +--+  +--+
|*   *  *|                                        |*   *  *|
+--------+                                        +--------+

18.1.2 裁剪与叠加分析的关系

裁剪实际上是叠加分析的一种特殊情况：

操作 等价的叠加分析 属性处理
裁剪（clip） overlay(how='intersection') 只保留被裁剪层的属性
掩膜（mask） overlay(how='difference') 只保留被掩膜层的属性

主要区别在于：clip() 更简洁，不合并裁剪区域的属性到结果中。

18.1.3 适用场景

场景 推荐操作
提取某行政区内的所有数据 裁剪
排除水体区域的陆地数据 掩膜
限定研究范围 裁剪
去除云层覆盖区域 掩膜
将全国数据限定到省级范围 裁剪

18.2 clip() 函数基础用法

18.2.1 基本语法

geopandas.clip(gdf, mask, keep_geom_type=False)

或使用 GeoDataFrame 的方法形式：

gdf.clip(mask, keep_geom_type=False)

18.2.2 参数说明

参数 类型 说明
gdf GeoDataFrame / GeoSeries 要裁剪的数据
mask GeoDataFrame / GeoSeries / Polygon 裁剪区域
keep_geom_type bool 是否只保留与输入相同的几何类型（默认 False）

18.2.3 第一个裁剪操作

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

np.random.seed(42)

# 创建点数据
n = 100
points = gpd.GeoDataFrame({
    'id': range(n),
    'value': np.random.randint(1, 100, n),
    'geometry': [Point(x, y) for x, y in zip(
        np.random.uniform(0, 10, n),
        np.random.uniform(0, 10, n)
    )]
})

# 定义裁剪区域
clip_area = gpd.GeoDataFrame({
    'geometry': [Polygon([(2, 2), (7, 2), (7, 8), (2, 8)])]
})

# 执行裁剪
clipped = gpd.clip(points, clip_area)
print(f"裁剪前点数: {len(points)}")
print(f"裁剪后点数: {len(clipped)}")
print(f"\n裁剪结果前5行:")
print(clipped.head())

18.2.4 使用 Shapely 几何对象作为 mask

mask 参数不仅可以接受 GeoDataFrame，也可以直接使用 Shapely 几何对象：

from shapely.geometry import Polygon, box

# 使用 Polygon 对象
polygon_mask = Polygon([(2, 2), (7, 2), (7, 8), (2, 8)])
clipped_1 = gpd.clip(points, polygon_mask)

# 使用 box 创建矩形
box_mask = box(3, 3, 6, 6)
clipped_2 = gpd.clip(points, box_mask)

print(f"多边形裁剪: {len(clipped_1)} 个点")
print(f"矩形裁剪: {len(clipped_2)} 个点")

18.2.5 使用 GeoSeries 作为 mask

# GeoSeries 作为 mask（多个裁剪区域的并集）
mask_series = gpd.GeoSeries([
    Polygon([(1, 1), (3, 1), (3, 3), (1, 3)]),
    Polygon([(6, 6), (9, 6), (9, 9), (6, 9)])
])

clipped = gpd.clip(points, mask_series)
print(f"多区域裁剪结果: {len(clipped)} 个点")

18.3 按多边形裁剪点数据

18.3.1 基本裁剪

点数据的裁剪最为直观——落在裁剪区域内的点被保留，区域外的点被移除：

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

np.random.seed(42)

# 气象站点数据
stations = gpd.GeoDataFrame({
    'station': [f'站点{i}' for i in range(20)],
    'temperature': np.random.uniform(15, 35, 20),
    'rainfall': np.random.uniform(0, 100, 20),
    'geometry': [Point(x, y) for x, y in zip(
        np.random.uniform(110, 120, 20),
        np.random.uniform(30, 40, 20)
    )]
}, crs="EPSG:4326")

# 省级边界
province = gpd.GeoDataFrame({
    'name': ['目标省份'],
    'geometry': [Polygon([
        (113, 32), (117, 32), (118, 35),
        (116, 37), (113, 36)
    ])]
}, crs="EPSG:4326")

# 裁剪 - 获取省内气象站
provincial_stations = gpd.clip(stations, province)

print(f"全部气象站: {len(stations)}")
print(f"省内气象站: {len(provincial_stations)}")
print(f"\n省内站点详情:")
print(provincial_stations[['station', 'temperature', 'rainfall']])

18.3.2 裁剪后的统计分析

# 省内气温统计
print(f"省内平均气温: {provincial_stations['temperature'].mean():.1f}°C")
print(f"省内最高气温: {provincial_stations['temperature'].max():.1f}°C")
print(f"省内总降雨量: {provincial_stations['rainfall'].sum():.1f}mm")

18.3.3 边界上的点

默认情况下，恰好位于裁剪区域边界上的点也会被保留：

# 创建恰好在边界上的点
boundary_points = gpd.GeoDataFrame({
    'name': ['内部', '边界上', '外部'],
    'geometry': [
        Point(1, 1),    # 在内部
        Point(2, 2),    # 在边界上
        Point(5, 5)     # 在外部
    ]
})

mask = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
result = gpd.clip(boundary_points, mask)
print("裁剪结果:")
print(result['name'].tolist())

18.4 按多边形裁剪线数据

18.4.1 线数据裁剪的特点

裁剪线数据时，穿越裁剪区域边界的线段会被切割，只保留位于裁剪区域内的部分：

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

# 道路数据
roads = gpd.GeoDataFrame({
    'road_name': ['主干道', '环路', '跨区公路', '区内道路'],
    'road_class': ['一级', '二级', '一级', '三级'],
    'geometry': [
        LineString([(0, 5), (10, 5)]),       # 横穿全域
        LineString([(3, 2), (7, 2), (7, 8), (3, 8), (3, 2)]),  # 部分在内
        LineString([(1, 1), (9, 9)]),         # 斜向穿越
        LineString([(4, 4), (6, 6)])          # 完全在内
    ]
})

# 裁剪区域
clip_zone = gpd.GeoDataFrame({
    'geometry': [Polygon([(3, 3), (8, 3), (8, 7), (3, 7)])]
})

# 裁剪
clipped_roads = gpd.clip(roads, clip_zone)

print("裁剪前后对比:")
for idx, row in clipped_roads.iterrows():
    original_length = roads.loc[idx, 'geometry'].length
    clipped_length = row.geometry.length
    print(f"  {row['road_name']}: 原长 {original_length:.2f} -> 裁剪后 {clipped_length:.2f}")

18.4.2 保持几何类型

裁剪可能将线段切割成多段，产生 MultiLineString。使用 keep_geom_type=True 可以确保结果保持线类型：

# 不限制几何类型
result_all = gpd.clip(roads, clip_zone, keep_geom_type=False)
print("不限制类型:", [g.geom_type for g in result_all.geometry])

# 仅保留线类型
result_lines = gpd.clip(roads, clip_zone, keep_geom_type=True)
print("仅保留线:", [g.geom_type for g in result_lines.geometry])

18.4.3 计算裁剪后的长度

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

# 假设已投影到平面坐标系
roads_proj = gpd.GeoDataFrame({
    'name': ['路段A', '路段B', '路段C'],
    'geometry': [
        LineString([(0, 500), (1000, 500)]),
        LineString([(200, 0), (200, 1000)]),
        LineString([(800, 0), (800, 1000)])
    ]
}, crs="EPSG:32650")

zone = gpd.GeoDataFrame({
    'geometry': [Polygon([(100, 100), (600, 100), (600, 800), (100, 800)])]
}, crs="EPSG:32650")

clipped = gpd.clip(roads_proj, zone)
clipped['length_m'] = clipped.geometry.length

print("裁剪后道路长度（米）:")
print(clipped[['name', 'length_m']])
print(f"\n区域内道路总长度: {clipped['length_m'].sum():.1f} 米")

18.5 按多边形裁剪面数据

18.5.1 面数据裁剪的特点

裁剪面数据时，跨越裁剪边界的多边形会被切割，只保留位于裁剪区域内的部分：

import geopandas as gpd
from shapely.geometry import Polygon

# 土地利用数据
landuse = gpd.GeoDataFrame({
    'type': ['耕地', '林地', '建设用地', '水域'],
    'area_ha': [500, 300, 200, 100],
    'geometry': [
        Polygon([(0, 0), (5, 0), (5, 5), (0, 5)]),
        Polygon([(5, 0), (10, 0), (10, 5), (5, 5)]),
        Polygon([(0, 5), (5, 5), (5, 10), (0, 10)]),
        Polygon([(5, 5), (10, 5), (10, 10), (5, 10)])
    ]
})

# 研究区边界
study_area = gpd.GeoDataFrame({
    'geometry': [Polygon([(2, 2), (8, 2), (8, 8), (2, 8)])]
})

# 裁剪
clipped_landuse = gpd.clip(landuse, study_area)

# 对比裁剪前后面积
print("裁剪前后面积对比:")
for idx, row in clipped_landuse.iterrows():
    original_area = landuse.loc[idx, 'geometry'].area
    clipped_area = row.geometry.area
    ratio = clipped_area / original_area * 100
    print(f"  {row['type']}: {original_area:.1f} -> {clipped_area:.1f} ({ratio:.1f}%)")

18.5.2 面积重新计算

裁剪后，原始面积属性不再准确，需要重新计算：

# 重新计算面积
clipped_landuse['new_area'] = clipped_landuse.geometry.area

# 按面积比例重新分配属性
for idx, row in clipped_landuse.iterrows():
    original = landuse.loc[idx]
    ratio = row['new_area'] / original.geometry.area
    clipped_landuse.loc[idx, 'adjusted_ha'] = original['area_ha'] * ratio

print("调整后的面积（公顷）:")
print(clipped_landuse[['type', 'area_ha', 'adjusted_ha']])

18.5.3 处理复杂边界

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

# 不规则裁剪区域（如自然边界）
irregular_boundary = gpd.GeoDataFrame({
    'geometry': [Polygon([
        (1, 1), (4, 0.5), (7, 2),
        (9, 5), (7, 8), (4, 9),
        (1, 7), (0, 4)
    ])]
})

# 裁剪
result = gpd.clip(landuse, irregular_boundary)
print(f"不规则裁剪结果: {len(result)} 个要素")

# 检查是否产生了 MultiPolygon
for idx, row in result.iterrows():
    print(f"  {row['type']}: {row.geometry.geom_type}")

18.5.4 使用 keep_geom_type 参数

# 裁剪可能产生 GeometryCollection（包含点或线碎片）
result_all = gpd.clip(landuse, study_area, keep_geom_type=False)
result_poly = gpd.clip(landuse, study_area, keep_geom_type=True)

print(f"不限类型: {len(result_all)} 个要素")
print(f"仅多边形: {len(result_poly)} 个要素")

18.6 使用 mask 进行掩膜操作

18.6.1 掩膜的实现方式

GeoPandas 没有直接的 mask() 函数，但可以通过叠加分析的差集操作实现：

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

np.random.seed(42)

# 原始数据
data = gpd.GeoDataFrame({
    'id': range(50),
    'value': np.random.randint(1, 100, 50),
    'geometry': [Point(x, y) for x, y in zip(
        np.random.uniform(0, 10, 50),
        np.random.uniform(0, 10, 50)
    )]
})

# 掩膜区域（需要排除的区域）
mask_area = gpd.GeoDataFrame({
    'geometry': [Polygon([(3, 3), (7, 3), (7, 7), (3, 7)])]
})

# 方法一：使用 overlay 差集
masked = gpd.overlay(data, mask_area, how='difference')
print(f"原始数据: {len(data)} 个点")
print(f"掩膜后: {len(masked)} 个点")

18.6.2 基于空间关系的掩膜

对于点数据，可以使用空间关系判断实现更灵活的掩膜：

# 方法二：使用空间谓词
from shapely.geometry import Polygon

mask_polygon = Polygon([(3, 3), (7, 3), (7, 7), (3, 7)])

# 找出不在掩膜区域内的点
outside_mask = data[~data.geometry.within(mask_polygon)]
print(f"掩膜外的点: {len(outside_mask)}")

# 方法三：使用 sjoin 反选
joined = gpd.sjoin(data, mask_area, predicate='within')
mask_indices = joined.index
outside = data[~data.index.isin(mask_indices)]
print(f"sjoin 反选结果: {len(outside)}")

18.6.3 多区域掩膜

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

np.random.seed(42)

# 数据点
n = 200
all_points = gpd.GeoDataFrame({
    'id': range(n),
    'geometry': [Point(x, y) for x, y in zip(
        np.random.uniform(0, 10, n),
        np.random.uniform(0, 10, n)
    )]
})

# 多个需要排除的区域
exclusion_zones = gpd.GeoDataFrame({
    'zone_name': ['水域', '保护区', '军事区'],
    'geometry': [
        Polygon([(1, 1), (3, 1), (3, 3), (1, 3)]),
        Polygon([(5, 5), (8, 5), (8, 8), (5, 8)]),
        Polygon([(7, 0), (9, 0), (9, 2), (7, 2)])
    ]
})

# 合并所有排除区域后掩膜
merged_exclusion = exclusion_zones.union_all()
valid_points = all_points[~all_points.geometry.within(merged_exclusion)]

print(f"总数据点: {len(all_points)}")
print(f"排除后: {len(valid_points)}")
print(f"排除了: {len(all_points) - len(valid_points)} 个点")

18.6.4 面数据的掩膜

import geopandas as gpd
from shapely.geometry import Polygon

# 研究区域
study = gpd.GeoDataFrame({
    'name': ['研究区'],
    'geometry': [Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])]
})

# 需要排除的区域
water_bodies = gpd.GeoDataFrame({
    'type': ['湖泊', '河流缓冲区'],
    'geometry': [
        Polygon([(2, 3), (4, 3), (4, 5), (2, 5)]),
        Polygon([(6, 0), (7, 0), (7, 10), (6, 10)])
    ]
})

# 差集操作实现掩膜
land_only = gpd.overlay(study, water_bodies, how='difference')
print(f"原始面积: {study.geometry.area.sum():.1f}")
print(f"排除水体后面积: {land_only.geometry.area.sum():.1f}")
print(f"水体面积: {study.geometry.area.sum() - land_only.geometry.area.sum():.1f}")

18.7 裁剪与叠加分析的区别

18.7.1 功能对比

特征 clip() overlay(how=’intersection’)
结果几何 裁剪后的几何 交集几何（相同）
结果属性 仅被裁剪层的属性 两层的属性合并
裁剪区域属性 不包含在结果中 包含在结果中
性能 通常更快 需要更多计算
适用几何 所有类型 主要用于面数据

18.7.2 代码对比

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

np.random.seed(42)

points = gpd.GeoDataFrame({
    'poi': ['A', 'B', 'C'],
    'geometry': [Point(1, 1), Point(3, 3), Point(5, 5)]
})

zone = gpd.GeoDataFrame({
    'zone_name': ['核心区'],
    'zone_level': [1],
    'geometry': [Polygon([(0, 0), (4, 0), (4, 4), (0, 4)])]
})

# 方法一：clip
clip_result = gpd.clip(points, zone)
print("clip 结果列:", clip_result.columns.tolist())
# 输出: ['poi', 'geometry'] — 不包含 zone 属性

# 方法二：overlay
overlay_result = gpd.overlay(points, zone, how='intersection')
print("overlay 结果列:", overlay_result.columns.tolist())
# 输出: ['poi', 'zone_name', 'zone_level', 'geometry'] — 包含 zone 属性

# 方法三：sjoin 后筛选
sjoin_result = gpd.sjoin(points, zone, predicate='within')
print("sjoin 结果列:", sjoin_result.columns.tolist())
# 输出: ['poi', 'geometry', 'index_right', 'zone_name', 'zone_level']

18.7.3 选择建议

需求 推荐方法
仅提取区域内的数据 clip()
需要裁剪区域的属性 overlay() 或 sjoin()
点数据在区域内外的判断 sjoin()
线/面数据的精确裁剪 clip()
面数据的交叉分析 overlay()

18.8 实际应用案例

18.8.1 案例一：提取城市范围内的道路网络

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

# 模拟道路网络
roads = gpd.GeoDataFrame({
    'road_id': [f'R{i}' for i in range(8)],
    'road_type': ['高速', '主干道', '主干道', '次干道',
                  '次干道', '支路', '支路', '高速'],
    'geometry': [
        LineString([(0, 5), (10, 5)]),
        LineString([(2, 0), (2, 10)]),
        LineString([(8, 0), (8, 10)]),
        LineString([(0, 3), (10, 3)]),
        LineString([(0, 7), (10, 7)]),
        LineString([(4, 0), (4, 10)]),
        LineString([(6, 0), (6, 10)]),
        LineString([(0, 0), (10, 10)])
    ]
})

# 城市边界
city_boundary = gpd.GeoDataFrame({
    'geometry': [Polygon([(1, 1), (9, 1), (9, 9), (1, 9)])]
})

# 裁剪
city_roads = gpd.clip(roads, city_boundary)
city_roads['length'] = city_roads.geometry.length

# 按道路类型统计
road_stats = city_roads.groupby('road_type').agg(
    条数=('road_id', 'count'),
    总长度=('length', 'sum'),
    平均长度=('length', 'mean')
).round(2)

print("城市道路统计:")
print(road_stats)
print(f"\n城市内道路总长度: {city_roads['length'].sum():.2f}")

18.8.2 案例二：提取流域内的土地利用

import geopandas as gpd
from shapely.geometry import Polygon

# 土地利用数据
landuse_data = gpd.GeoDataFrame({
    'lu_type': ['耕地', '林地', '草地', '建设用地', '水域', '未利用地'],
    'area_km2': [120, 80, 50, 30, 15, 25],
    'geometry': [
        Polygon([(0, 0), (4, 0), (4, 4), (0, 4)]),
        Polygon([(4, 0), (8, 0), (8, 4), (4, 4)]),
        Polygon([(8, 0), (12, 0), (12, 4), (8, 4)]),
        Polygon([(0, 4), (4, 4), (4, 8), (0, 8)]),
        Polygon([(4, 4), (8, 4), (8, 8), (4, 8)]),
        Polygon([(8, 4), (12, 4), (12, 8), (8, 8)])
    ]
})

# 流域边界
watershed = gpd.GeoDataFrame({
    'name': ['小流域'],
    'geometry': [Polygon([
        (2, 1), (7, 0.5), (10, 3),
        (9, 6), (5, 7), (1, 5)
    ])]
})

# 裁剪
watershed_lu = gpd.clip(landuse_data, watershed)

# 重新计算面积和比例
watershed_lu['clipped_area'] = watershed_lu.geometry.area
total_area = watershed_lu['clipped_area'].sum()
watershed_lu['proportion'] = watershed_lu['clipped_area'] / total_area * 100

print("流域内土地利用结构:")
print(watershed_lu[['lu_type', 'clipped_area', 'proportion']].to_string(
    formatters={
        'clipped_area': '{:.2f}'.format,
        'proportion': '{:.1f}%'.format
    }
))
print(f"\n流域总面积: {total_area:.2f}")

18.8.3 案例三：分区裁剪与批量处理

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

np.random.seed(42)

# 全局数据
n = 500
global_data = gpd.GeoDataFrame({
    'id': range(n),
    'category': np.random.choice(['A', 'B', 'C'], n),
    'value': np.random.uniform(0, 100, n),
    'geometry': [Point(x, y) for x, y in zip(
        np.random.uniform(0, 10, n),
        np.random.uniform(0, 10, n)
    )]
})

# 多个分区
zones = gpd.GeoDataFrame({
    'zone': ['东区', '西区', '南区', '北区'],
    'geometry': [
        Polygon([(5, 5), (10, 5), (10, 10), (5, 10)]),
        Polygon([(0, 5), (5, 5), (5, 10), (0, 10)]),
        Polygon([(0, 0), (5, 0), (5, 5), (0, 5)]),
        Polygon([(5, 0), (10, 0), (10, 5), (5, 5)])
    ]
})

# 批量裁剪并统计
results = {}
for _, zone_row in zones.iterrows():
    zone_gdf = gpd.GeoDataFrame([zone_row], columns=zones.columns)
    clipped = gpd.clip(global_data, zone_gdf)
    results[zone_row['zone']] = {
        '数据量': len(clipped),
        '平均值': clipped['value'].mean(),
        'A类占比': (clipped['category'] == 'A').mean() * 100
    }

import pandas as pd
summary = pd.DataFrame(results).T
print("分区统计:")
print(summary.round(1))

18.8.4 案例四：排除特定区域的空间分析

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

np.random.seed(42)

# 采样点
samples = gpd.GeoDataFrame({
    'sample_id': range(30),
    'concentration': np.random.uniform(0, 50, 30),
    'geometry': [Point(x, y) for x, y in zip(
        np.random.uniform(0, 10, 30),
        np.random.uniform(0, 10, 30)
    )]
})

# 需要排除的异常区域
anomaly_zones = gpd.GeoDataFrame({
    'reason': ['工厂影响', '历史污染'],
    'geometry': [
        Polygon([(1, 1), (3, 1), (3, 3), (1, 3)]),
        Polygon([(7, 7), (9, 7), (9, 9), (7, 9)])
    ]
})

# 裁剪获取异常区域内的点
anomaly_points = gpd.sjoin(samples, anomaly_zones, predicate='within')
clean_samples = samples[~samples.index.isin(anomaly_points.index)]

print(f"总采样点: {len(samples)}")
print(f"异常区域内: {len(anomaly_points)}")
print(f"有效采样点: {len(clean_samples)}")
print(f"\n有效采样点平均浓度: {clean_samples['concentration'].mean():.2f}")
print(f"异常区域平均浓度: {anomaly_points['concentration'].mean():.2f}")

18.9 本章小结

本章详细介绍了 GeoPandas 中的裁剪与掩膜操作。主要内容回顾：

核心函数

操作 方法 说明
裁剪 gpd.clip(gdf, mask) 保留 mask 区域内的数据
掩膜 gpd.overlay(gdf, mask, how='difference') 排除 mask 区域的数据
点掩膜 gdf[~gdf.within(mask_geom)] 基于空间谓词的快速掩膜

不同几何类型的裁剪行为

几何类型 裁剪行为
点 保留在裁剪区域内的点
线 切割并保留区域内的线段
面 切割并保留区域内的面片

关键参数

参数 说明 默认值
mask 裁剪区域 必需
keep_geom_type 保留相同几何类型 False

使用建议

  1. 选择合适的方法：仅需限定范围时用 clip()，需要属性关联时用 sjoin() 或 overlay()
  2. 注意 CRS 一致性：裁剪前确保数据与裁剪区域的坐标参考系统一致
  3. 重新计算属性：裁剪后面积、长度等属性需要重新计算
  4. 处理碎片几何：使用 keep_geom_type=True 避免产生不期望的几何类型
  5. 大数据优化：大数据集建议先用边界框预筛选，再执行精确裁剪

在下一章中，我们将学习聚合与融合（Dissolve）操作，它将多个要素按属性合并为更大的区域。

↑