python矢量裁剪栅格(按掩模提取)

前言

矢量裁剪栅格,也就是将栅格中对应矢量区域的像素裁剪出来,形成一个新栅格的过程,对应Arcgis中空间分析工具里的按掩模提取功能。

过程

gdal方法实现

其中一种方法是利用gdal.Warp工具进行裁剪。

gdal.Warp各个参数

一般情况下,只需要用到几个参数,具体使用方式如下:

try:
    import gdal
except:
    from osgeo import gdal
def extract_by_shp_gdal(in_shp_path, in_raster_path, out_raster_path):
    input_raster = gdal.Open(in_raster_path)
    # 利用gdal.Warp进行裁剪
    # https://gdal.org/api/python/osgeo.gdal.html#osgeo.gdal.Warp
    result = gdal.Warp(
        out_raster_path,
        input_raster,
        format = 'GTiff',
        cutlineDSName = in_shp_path, # 用于裁剪的矢量
        cropToCutline = True, # 是否使用cutlineDSName的extent作为输出的界线
        dstNodata = 0 # 输出数据的nodata值
    result.FlushCache()
    del result

rasterio方法实现

gdal方法好处就是安装了gdal就行,但是偶尔可能会出现小bug。rasterio则比较稳健,非常好用。

import fiona
import rasterio
from rasterio.mask import mask
    import gdal
except:
    from osgeo import gdal
def extract_by_shp_rasterio(in_shp_path, in_raster_path, out_raster_path):
    with fiona.open(in_shp_path, "r", encoding='utf-8') as shapefile:
        # 获取所有要素feature的形状geometry
        geoms = [feature["geometry"] for feature in shapefile]
    with rasterio.open(in_raster_path) as src:
        out_image, out_transform = mask(src, geoms, crop=True)
        out_meta = src.meta.copy()
    # 更新输出文件的元数据
    out_meta.update({"driver": "GTiff",
                     "height": out_image.shape[1],