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])
mlateb
  • 139
  • 1
  • 6

1 Answers1

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])
user2856
  • 65,736
  • 6
  • 115
  • 196