使用 Rasterio 和 Fiona 在 Python 中使用多个 shapefile 子集/裁剪多个栅格

问题描述

我正在尝试采用现有代码根据此链接在 Python 3.6 中使用多个多边形裁剪多个栅格：

https://gis.stackexchange.com/questions/382783/subsetting-multiple-rasters-with-multiple-shapefiles-using-python

我使用以下代码创建了 shp_list 和 raster_list。如何定义栅格、shapefile 名称和外路径？ 另外，如何修改代码将每个shapefile的剪辑保存到单独的文件夹中？

import os
import numpy as np
from PIL import Image
import shapefile
from Rasterio.mask import mask

path= 'C:/dataset/test'
os.chdir(path)

imageFiles = list(range(1,6))           # images = [1,2,3,4,5]
shpFiles= list(range(1,3))           # shapefiles = [1,2]

# Make empty list of images
raster_list = []

# Load all images,primary and ancillary
for i in [*imageFiles]:
    raster_names = f'ras{i}.tif'
    print(raster_names)
    im = np.array(Image.open(raster_names))
    raster_list.insert(i,im)

result:
ras1.tif
ras2.tif
ras3.tif
ras4.tif
ras5.tif

# Make empty list of shapefiles
shp_list = []

# Load all shapefiles
for i in [*shpFiles]:
    shp_names = f'no{i}.shp'
    print(shp_names)
    shp = shapefile.Reader(shp_names)
    shp_list.insert(i,shp)

result:
no1.shp
no2.shp

def subs(shp_list,shp_names,raster_list,raster_names,outpath):
    for i,shp in enumerate(shp_list):
        with fiona.open(shp,"r") as shapefile:
            shapes = [feature["geometry"] for feature in shapefile]

            for j,ras in enumerate(raster_list):
                for scene in raster_list:
                    with Rasterio.open(ras,"r") as src:
                        out_image,out_transform = mask(src,shapes,crop=True)
                        out_Meta = src.Meta

                        out_Meta.update({"driver": "GTiff","height": out_image.shape[1],"width": out_image.shape[2],"transform": out_transform})
                                         
                        ras_path = f"{outpath}{raster_names[j]}{shp_names[i]}{ras_extension[1:]}"

                        with Rasterio.open(ras_path,"w",**out_Meta) as dest:
                            dest.write(out_image)

    # number of created subsets
    subset_count = (i+1) * (j+1)
    if len(shp_list) * len(raster_list) == subset_count:
        print(f"Done. \n{subset_count} subsets created")

# calling the function
subs(shp_list,outpath)

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）