Subsetting multiple rasters with multiple shapefiles using Python

Question

I want to subset multiple raster files with multiple shapefiles and the following code is working (it is not the whole code, only the part for the subsetting). I am using Python 3.6.8

i_shp = 0  # iterate over .shp-files
while i_shp < len(shp_list):
    with fiona.open(shp_list[i_shp], "r") as shapefile:
        shapes = [feature["geometry"] for feature in shapefile]

        i_ras = 0  # iterate over raster files
        while i_ras < len(raster_list):
            for scene in raster_list:
                with rasterio.open(raster_list[i_ras], "r") as src:
                    out_image, out_transform = rasterio.mask.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})

                    with rasterio.open(
                            str(outpath) + raster_names[i_ras] + str(shp_names[i_shp]) + str(ras_extension[1:])
                            , "w", **out_meta) as dest:
                        dest.write(out_image)
                i_ras = i_ras + 1
    i_shp = i_shp + 1

I`ve got 7 rasters and 2 shapefiles and after I run the code I have 14 subsets.

But I want to place this code inside a function but then it's not working correctly. I just get 7 subsets (only from the first shapefile in my shp_list). I think that the while-loop isn't working correctly, but I cannot figure out why.

def subs(shp_list, shp_names, raster_list, raster_names, outpath):
    # subsetting all raster files with all .shp-files
    i_shp = 0  # iterate over .shp-files
    while i_shp < len(shp_list):
        with fiona.open(shp_list[i_shp], "r") as shapefile:
            shapes = [feature["geometry"] for feature in shapefile]

            i_ras = 0  # iterate over raster files
            while i_ras < len(raster_list):
                for scene in raster_list:
                    with rasterio.open(raster_list[i_ras], "r") as src:
                        out_image, out_transform = rasterio.mask.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})

                        with rasterio.open(
                                str(outpath) + raster_names[i_ras] + str(shp_names[i_shp]) + str(ras_extension[1:])
                                , "w", **out_meta) as dest:
                            dest.write(out_image)
                    i_ras = i_ras + 1
        i_shp = i_shp + 1

        # number of created subsets
        subset_count = i_shp * i_ras
        if len(shp_list) * len(raster_list) == subset_count:
            (str("Done. \n ") + str(subset_count) + str(" subsets created"))

        return print(str("Done. \n ") + str(subset_count) + str(" subsets created"))

Answer 1

You can use the script below. You don't need while here. Instead, use for .... in enumerate(...): for a cleaner script. Also I've changed string definitions to f-strings which work in Python 3.6 and provide a clenaer string definition.

# other imports

from rasterio.mask import mask

# other stuffs -> shp_list, shp_names etc.

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, shp_names, raster_list, raster_names, outpath)