Based on answers from previous questions 1,2 and the Masking a raster using a shapefile guide from rasterio cookbook, I wrote the code below.

The code runs correctly. I would like to improve it.

1- Notice that I'm using rasterio.mask.mask(src0, features,crop=True) twice. is there a way to keep only one?

2- Before the last line of code I'm slicing out_image[0,:,:] array in order to get a 2D array. Can avoid this operation?

import rasterio
import rasterio.mask
import fiona

# rasterio.__version__ =='1.0b1'

#open a shapefile, read geometries
with fiona.open(shapefile, "r") as shapefile:
    features = [feature["geometry"] for feature in shapefile]

# Read metadata of first file
with rasterio.open(list_to_stack[0]) as src0:
    meta = src0.meta

    # extract new mask raster info
    out_image, out_transform = rasterio.mask.mask(src0, features,crop=True)

# Update meta due to stacking
meta.update({"count": len(list_to_stack), #set the number of stack layers
             "driver": "GTiff"}) #set the input is ENVI format

# Update meta with masking info
meta.update({"height": out_image.shape[1], 
             "width": out_image.shape[2], 
             "transform": out_transform})

#Read each layer and write it to stack 
with rasterio.open(output_stack, 'w', **meta) as dst:
    for id, layer in enumerate(list_to_stack):
        with rasterio.open(layer) as src1:
            out_image, _ = rasterio.mask.mask(src1, features,crop=True)

            dst.write_band(id + 1,out_image[0,:,:].astype(rasterio.int16))
            # add source dataset description
            dst.set_band_description(id+1 , src1.descriptions[0])

Your Answer

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

Browse other questions tagged or ask your own question.