4

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"))
5

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)
1
  • 1
    Thank you very much. It works perfect! – FelixB. Dec 23 '20 at 18:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.