代码收藏家技术教程 2024-05-28

Python批量裁剪（掩膜）tif文件使用shp面数据

文章目录

前言

一、基于哪些库？

二、使用步骤

1.引入库

2.定义裁剪函数

3.引用函数裁剪

总结

前言

不知道各位在日常处理遥感影像时，是否会遇到这种情况：就是假设有一个shp文件，里面有很多面，需要基于这些面对tif分别进行裁剪。这种裁剪其实arcgis也可以，但是就比较麻烦，于是乎我想到了利用python代码直接批量裁剪，请看下文。

一、基于那些库？

主要利用的库有geopandas、rasterio、os 。

1.geopandas主要用于处理地理空间数据。提供了一个GeoDataFrame对象，它是pandas的DataFrame对象的扩展，能够存储和操作地理空间数据。

2.rasterio用于读取、写入和处理栅格数据集。它建立在GDAL（Geospatial Data Abstraction Library）之上，提供了一种简单而高效的方式来处理遥感和地理信息系统数据。

3.os不做介绍，懂得都懂。

二、使用步骤

1.引入库

代码如下：

import os
import geopandas as gpd
import rasterio
from rasterio.mask import mask
from rasterio.crs import CRS

2.定义裁剪函数

代码如下：

def clip_raster_by_shp(raster_path, shp_path, output_folder):
    # 读取shp文件
    gdf = gpd.read_file(shp_path)

    # 读取栅格数据
    with rasterio.open(raster_path) as src:
        # 获取栅格数据的元数据
        src_meta = src.meta

        if src.crs:
            gdf_crs = CRS.from_user_input(gdf.crs) if gdf.crs else None
            src_crs = CRS.from_user_input(src.crs) if src.crs else None

            if gdf_crs != src_crs:
                gdf_projected = gdf.to_crs(src.crs)
            else:
                gdf_projected = gdf
        else:
            gdf_projected = gdf

        # 遍历shp文件中的每个区域
        for index, row in gdf_projected.iterrows():
            # 获取当前区域的几何形状
            geometry_projected = row['geometry']

            # 使用区域的几何形状裁剪栅格数据
            out_image, out_transform = mask(src, [geometry_projected], crop=True)

            # 更新元数据
            out_meta = src_meta.copy()
            out_meta.update({
                "driver": "GTiff",
                "height": out_image.shape[1],
                "width": out_image.shape[2],
                "transform": out_transform
            })

            # 保存裁剪后的栅格数据
            out_filename = os.path.join(output_folder, f"clip_{index}.tif")
            with rasterio.open(out_filename, "w", **out_meta) as dest:
                dest.write(out_image)

3.引用函数裁剪

代码如下：

raster_path =r'demo.tif'      # 设置需要裁剪的tif路径

shp_path =   r'fanwei.shp'    # 设置shp路径

output_folder = r'cut_result' # 设置结果保存路径

# 调用函数实现裁剪
clip_raster_by_shp(raster_path, shp_path, output_folder)