6

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])

1 Answer 1

1

I would try something like (untested) which skips the first mask as you've already done that and uses numpy.squeeze which returns a view instead of slicing so no extra memory is used.

import numpy as np

# etc...

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

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

# Update meta due to stacking
# etc...

#Read each layer and write it to stack
with rasterio.open(output_stack, 'w', **meta) as dst:
    
    # Write the 1st image you already masked
    dst.write(np.squeeze(out_image).astype('int16'), 1)  # squeeze returns a view so no extra memory used
    dst.set_band_description(1, desc)

    for id, layer in enumerate(list_to_stack[1:], start=2):  # Start from the 2nd element (list_to_stack[1])
        with rasterio.open(layer) as src1:
            out_image, _ = rasterio.mask.mask(src1, features, crop=True)

            dst.write(np.squeeze(out_image).astype('int16'), id)  # squeeze returns a view so no extra memory used
            # add source dataset description
            dst.set_band_description(id, src1.descriptions[0])

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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