【Python&RS】基于Rasterio库裁剪栅格数据&压缩栅格影像

        之前分享过【Python&RS】Rasterio库安装+基础函数使用教程，大家有兴趣的可以去看看。由于最近有涉及到栅格裁剪和压缩的问题，所以研究了一下今天和大家分享分享。

原创作者：RS迷途小书童

博客地址：https://blog.csdn.net/m0_56729804?type=blog

1 需要的库

import os
import rasterio
import geopandas as gpd
from rasterio.mask import mask
from rasterio.plot import show
from rasterio.windows import Window, from_bounds
from shapely.geometry import Point, Polygon, box, mapping, LineString

2 自定义矩形框裁剪

        这里的矩形可以自己选择大小，当然也可以通过读取矢量的空间范围去裁剪。同时代码中还加入了压缩参数。

def compress_tif(tif_path, output_file):
    with rasterio.open(tif_path) as src:
        transform = src.transform
        # minx, miny, maxx, maxy = geometry1.bounds(要素的地理字段)
        expanded_bbox = box(minx=0, miny=0, maxx=100, maxy=100)  # 创建一个新的扩展后的边界框（Shapely box）
        window = from_bounds(expanded_bbox.bounds[0], expanded_bbox.bounds[1],
                             expanded_bbox.bounds[2], expanded_bbox.bounds[3],
                             transform=transform)  # 将Shapely几何对象转换为rasterio可以理解的边界
        clipped = src.read(window=window)  # 读取并裁剪TIFF数据

        out_meta = src.meta.copy()
        out_meta.update({"driver": "GTiff", "compress": 'lzw'})  # 更新元数据,rle，lzw等
        del src
    with rasterio.open(output_file, 'w', **out_meta) as dest:  # 写入裁剪后的TIFF文件
        dest.write(clipped)

3 使用矢量要素裁剪

        这里使用到了geopandas库用来读取每个要素的空间范围。

def clip_raster_from_features(vector_file=r"彭俊喜/1.shp", raster_file=r"彭俊喜/1.tif"):
    """
    :param vector_file: 输入需要裁剪的面矢量（多面）
    :param raster_file: 输入需要裁剪的影像
    :return: None
    """
    # 读取面矢量数据
    gdf = gpd.read_file(vector_file)
    # 循环遍历面矢量中的每个面要素
    for index, row in gdf.iterrows():
        # 获取当前面要素的几何形状
        geometry1 = row.geometry
        # 确保面要素不是空的
        if geometry1.is_empty:
            print(f"Skipping empty geometry for feature {index}")
            continue
            # 打开影像文件
        with rasterio.open(raster_file) as src:
            # 将面要素的边界转换为shapely的box对象
            geometry_bounds = box(*geometry1.bounds)
            # 将rasterio的bounds转换为shapely的box对象
            src_bounds = box(*src.bounds)
            # 检查面要素的边界是否与影像的边界相交
            if not geometry_bounds.intersects(src_bounds):
                print(f"Skipping feature {index} as it does not intersect with the raster.")
                continue
                # 转换几何形状为Rasterio可以理解的格式
            geom_for_rasterio = mapping(geometry1)
            # 使用面要素裁剪影像
            try:
                out_image, out_transform = mask(src, [geom_for_rasterio], crop=True)
            except ValueError as e:
                print(f"Error clipping feature {index}: {e}")
                continue
                # 检查是否成功裁剪出影像
            if out_image is None:
                print(f"No data was clipped for feature {index}. Skipping.")
                continue
                # 为裁剪后的影像设置输出路径和文件名
            output_file = f'clipped_image_{index}.tif'
            output_path = os.path.join(r'Z:\Shanghai Metro Automatic Identification System\Data Source\weigui\1',
                                       output_file)
            # 创建输出文件，并写入裁剪后的影像数据
            dtypes = [src.dtypes[i] for i in range(src.count)]  # 确保数据类型列表与波段数量匹配
            # 假设所有波段的数据类型都是相同的，并且你想要保持与输入影像相同的数据类型
            dtype = src.dtypes[0]  # 获取第一个波段的数据类型
            # 使用这个数据类型打开输出文件
            with rasterio.open(output_path, 'w', driver='GTiff', height=out_image.shape[1],
                               width=out_image.shape[2], count=src.count, dtype=dtype,
                               crs=src.crs, transform=out_transform) as dest:
                dest.write(out_image)
            print(f'Clipped image {output_file} has been saved.')
    print('All features have been processed.')