Rasterio掩膜提取?Mask函数怎么用?

wowwwai GIS研习社 · 工具流程与项目排障

引言：很多同学在做 Python GIS 栅格处理时，会遇到“Rasterio掩膜提取?Mask函数怎么用?”这个问题：手里有一幅 GeoTIFF 影像和一个矢量边界，想把边界内的栅格裁剪出来，但不知道 rasterio.mask.mask 的参数该怎么写，输出结果为什么会黑边、空值或坐标不对。本文用一个可复现的流程，讲清楚 Rasterio 掩膜提取的核心用法、常见错误和检查方法。

背景：Rasterio掩膜提取通常解决什么问题

Rasterio 是 Python 中常用的栅格数据读写库，适合处理 GeoTIFF、影像裁剪、重采样、波段读取、元数据更新等任务。所谓 Rasterio掩膜提取，通常指使用一个矢量范围，例如行政区边界、研究区面、多边形 GeoJSON，把栅格中范围内的像元保留下来，范围外的像元设为 NoData，或者直接裁剪到边界外包矩形。

典型场景包括：

• 用县界裁剪 DEM，只保留县域范围内的高程数据。
• 用研究区边界裁剪 Landsat、Sentinel 或无人机影像。
• 批量按多个行政区输出独立的栅格文件。
• 在建模前先把大范围栅格裁剪成项目需要的小范围数据。

在 Rasterio 中，最常用的函数是 rasterio.mask.mask。很多问题并不出在函数本身，而是出在坐标系不一致、矢量几何格式不对、NoData 没设置、输出元数据没有更新等细节上。

原理：mask函数到底做了什么

rasterio.mask.mask 的核心作用是：根据输入的几何范围，对栅格像元进行空间筛选，并返回两个结果：

• out_image：裁剪或掩膜后的栅格数组。
• out_transform：输出栅格新的仿射变换参数，用来描述像元和地理坐标之间的关系。

一个常见调用方式如下：

out_image, out_transform = mask(
    dataset,
    shapes,
    crop=True,
    nodata=dataset.nodata
)

这里的 dataset 是 Rasterio 打开的栅格对象，shapes 是 GeoJSON-like 格式的几何对象列表。所谓 GeoJSON-like，就是类似下面这样的 Python 字典结构：

{
    "type": "Polygon",
    "coordinates": [...]
}

crop=True 表示输出结果会被裁剪到矢量几何的外包矩形范围。如果设置为 False，输出栅格尺寸会保持原图大小，只是范围外像元被掩膜。

理解这一点很重要：Rasterio mask函数不是简单地按影像行列号截取，而是根据栅格坐标参考系中的几何范围进行空间计算。所以，栅格和矢量必须在同一个 CRS，即坐标参考系统中。

步骤：Rasterio mask函数怎么用

1. 安装需要的 Python 库

建议在独立环境中安装 Rasterio、GeoPandas 和 Fiona。GeoPandas 用来读取 Shapefile、GeoJSON 等矢量数据，Rasterio 用来读取和写出栅格。

pip install rasterio geopandas

如果你在 Windows 上安装 Rasterio 遇到 GDAL 相关错误，建议优先使用 Conda 环境：

conda create -n gis-raster python=3.11
conda activate gis-raster
conda install -c conda-forge rasterio geopandas

2. 准备输入数据

假设你有两个文件：

• input.tif：待裁剪的 GeoTIFF 栅格。
• boundary.shp：研究区边界面数据。

注意，掩膜边界最好是面要素。如果是线或点，通常不能直接作为栅格裁剪边界，除非先做缓冲区或转换为面。

3. 读取矢量边界并检查坐标系

下面代码先读取矢量边界，再打开栅格，并检查二者坐标系是否一致。如果不一致，就把矢量重投影到栅格坐标系。

import geopandas as gpd
import rasterio
from rasterio.mask import mask

raster_path = "input.tif"
vector_path = "boundary.shp"
output_path = "output_mask.tif"

gdf = gpd.read_file(vector_path)

with rasterio.open(raster_path) as src:
    print("Raster CRS:", src.crs)
    print("Vector CRS:", gdf.crs)

    if gdf.crs != src.crs:
        gdf = gdf.to_crs(src.crs)

这是 Rasterio 掩膜提取中最关键的一步。很多“输出为空”“裁剪位置偏移”“结果全是 NoData”的问题，都来自矢量和栅格坐标系不一致。

4. 转换为 mask函数需要的几何格式

rasterio.mask.mask 不能直接接收 GeoDataFrame，它需要的是几何对象列表。可以这样转换：

geometries = [geom for geom in gdf.geometry if geom is not None and not geom.is_empty]

如果你的边界文件包含多个面，以上写法会把所有面一起作为掩膜范围。如果你只想裁剪某一个行政区，可以先筛选属性：

gdf_one = gdf[gdf["name"] == "某某区"]
geometries = [geom for geom in gdf_one.geometry if geom is not None and not geom.is_empty]

5. 执行 Rasterio 掩膜提取

完整裁剪代码如下：

import geopandas as gpd
import rasterio
from rasterio.mask import mask

raster_path = "input.tif"
vector_path = "boundary.shp"
output_path = "output_mask.tif"

gdf = gpd.read_file(vector_path)

with rasterio.open(raster_path) as src:
    if gdf.crs != src.crs:
        gdf = gdf.to_crs(src.crs)

    geometries = [geom for geom in gdf.geometry if geom is not None and not geom.is_empty]

    out_image, out_transform = mask(
        src,
        geometries,
        crop=True,
        nodata=src.nodata
    )

    out_meta = src.meta.copy()
    out_meta.update({
        "driver": "GTiff",
        "height": out_image.shape[1],
        "width": out_image.shape[2],
        "transform": out_transform,
        "nodata": src.nodata
    })

    with rasterio.open(output_path, "w", **out_meta) as dst:
        dst.write(out_image)

这段代码就是最基础、最常用的 Rasterio mask函数用法。它会根据矢量边界裁剪栅格，并输出一个新的 GeoTIFF 文件。

6. 如果原始栅格没有 NoData 怎么办

有些 GeoTIFF 没有设置 src.nodata，这时 src.nodata 可能是 None。建议根据数据类型手动指定一个合理的 NoData 值。例如，DEM 可以使用 -9999，整数分类栅格要避免和真实类别值冲突。

nodata_value = src.nodata
if nodata_value is None:
    nodata_value = -9999

out_image, out_transform = mask(
    src,
    geometries,
    crop=True,
    nodata=nodata_value
)

out_meta.update({
    "nodata": nodata_value
})

如果是无符号整型数据，例如 uint8，不要随便使用 -9999，因为数据类型无法表示负数。可以考虑使用 0 或 255，但前提是这个值不代表有效分类。

7. 验证输出结果是否正确

写出结果后，建议至少检查四项：

• 输出文件能否在 QGIS 或 ArcGIS Pro 中正常打开。
• 输出栅格的 CRS 是否和原始栅格一致。
• 输出栅格范围是否覆盖研究区。
• NoData 是否正确显示为透明或无值区域。

也可以用 Rasterio 快速查看输出元数据：

with rasterio.open(output_path) as dst:
    print(dst.crs)
    print(dst.bounds)
    print(dst.width, dst.height)
    print(dst.nodata)
    print(dst.transform)

常见坑：Rasterio掩膜提取为什么会失败

1. 坐标系不一致导致输出为空

这是最常见的问题。比如栅格是 UTM 投影坐标，矢量是 WGS84 经纬度坐标，如果不先 to_crs(src.crs)，mask函数会认为二者不在同一个空间位置，结果可能为空或报错。

排查方法：

• 打印 src.crs 和 gdf.crs。
• 在 QGIS 中叠加查看是否真正重合。
• 不要只看图层“看起来在同一个地方”，要确认图层 CRS 和数据本身坐标正确。

2. 矢量边界没有和栅格相交

如果矢量范围和栅格范围完全不相交，mask 可能会抛出类似 Input shapes do not overlap raster 的错误。

可以先打印二者范围：

print("Raster bounds:", src.bounds)
print("Vector bounds:", gdf.total_bounds)

如果范围差异很大，优先检查坐标系、单位和数据来源。

3. 忘记更新 height、width 和 transform

当 crop=True 时，输出数组尺寸已经变化。如果写文件时仍然使用原始栅格的 height、width 和 transform，就会出现空间位置错误或写入失败。

正确做法是更新：

out_meta.update({
    "height": out_image.shape[1],
    "width": out_image.shape[2],
    "transform": out_transform
})

4. NoData 值和真实像元值冲突

如果你把 0 设置为 NoData，但原始数据中 0 本来就是有效值，那么输出结果会误把真实数据当成空值。分类栅格、土地利用数据尤其要注意这一点。

建议先查看数据统计或类别表，再选择 NoData 值。

5. 多波段影像输出维度理解错误

Rasterio 读取的数组通常是 (bands, rows, cols)，不是 (rows, cols, bands)。所以更新高度和宽度时要使用：

• out_image.shape[1] 作为高度。
• out_image.shape[2] 作为宽度。

方法比较：mask、窗口裁剪和GDAL裁剪怎么选

方法	适用场景	优点	注意事项
rasterio.mask.mask	按矢量边界裁剪栅格	代码简洁，适合 Python 自动化流程	需要保证栅格和矢量 CRS 一致
Rasterio window	按矩形范围或行列号裁剪	速度快，适合大栅格局部读取	不能直接按不规则多边形边界裁剪
GDAL Warp	命令行批处理、重投影加裁剪	功能强，适合生产环境脚本	参数较多，新手排错成本略高
QGIS 按掩膜图层裁剪栅格	可视化操作、少量数据处理	不用写代码，方便检查效果	不适合大量文件自动化处理

如果你的目标是写 Python 脚本自动裁剪多个栅格，Rasterio 掩膜提取是很合适的选择。如果只是偶尔裁剪一幅影像，QGIS 的“按掩膜图层裁剪栅格”工具更直观。如果还要同时做重投影、分辨率调整和格式转换，GDAL 的 gdalwarp 也值得考虑。

检查清单：写 Rasterio mask 脚本前先确认这些项

• 输入栅格是否是带地理参考的 GeoTIFF。
• 矢量边界是否为面要素，而不是点或线。
• 矢量 CRS 是否存在，不是 None。
• 矢量 CRS 是否已经转换到栅格 CRS。
• 矢量范围是否和栅格范围相交。
• 是否过滤了空几何和无效几何。
• crop=True 后是否更新了 height、width、transform。
• NoData 值是否适合当前栅格数据类型。
• 输出文件是否能在 QGIS 或 ArcGIS Pro 中正确叠加。

FAQ：Rasterio掩膜提取常见问题

Rasterio mask函数可以直接读取 Shapefile 吗？

rasterio.mask.mask 本身不直接读取 Shapefile。通常做法是用 GeoPandas 或 Fiona 读取 Shapefile，再把几何传给 mask 函数。也就是说，读取矢量和裁剪栅格是两个步骤。

为什么 Rasterio 掩膜提取结果全是黑色？

常见原因有三个：第一，NoData 没有正确设置；第二，输出数据范围很窄，但显示软件没有正确拉伸；第三，矢量和栅格坐标系不一致导致有效像元很少。建议先在 QGIS 中查看像元值，再检查 dst.nodata 和图层渲染方式。

mask函数的 crop=True 和 crop=False 有什么区别？

crop=True 会把输出栅格裁剪到几何范围的外包矩形，文件尺寸通常会变小。crop=False 会保留原始栅格尺寸，只把几何范围外的像元设为 NoData。大多数研究区裁剪任务会使用 crop=True。

Rasterio mask函数能处理多个多边形吗？

可以。只要把多个几何对象放进列表传给 mask 函数即可。多个多边形会共同作为保留范围。如果你想分别输出多个区域，需要循环每一个要素，并为每个要素单独执行 Rasterio 掩膜提取。

矢量边界有自相交或无效几何怎么办？

无效几何可能导致裁剪失败或输出异常。可以在 GeoPandas 中先检查：

print(gdf.is_valid.value_counts())

如果存在无效几何，可以尝试使用：

gdf["geometry"] = gdf.geometry.buffer(0)

但 buffer(0) 不是万能修复方法。重要边界数据建议在 QGIS 中用“修复几何”工具检查后再用于批处理。

输出 GeoTIFF 为什么叠加位置不对？

通常是输出元数据没有更新，尤其是 transform 仍然使用原始栅格的值。使用 mask 返回的 out_transform 更新输出文件，是保证空间位置正确的关键。

结论：掌握 Rasterio mask 的关键不是一行函数，而是完整流程

Rasterio掩膜提取的基本代码并不复杂，核心就是用 rasterio.mask.mask 读取矢量几何并裁剪栅格。但在真实项目中，结果是否正确取决于完整流程：坐标系是否一致、几何是否有效、NoData 是否合理、输出元数据是否更新。

如果你只记住一个原则，就是：先让矢量边界进入栅格的坐标系，再执行 mask，并用返回的 out_transform 更新输出 GeoTIFF。这样写出来的 Rasterio mask函数脚本，才更适合用于 GIS 项目中的批量裁剪和自动化处理。