I'm working with 52 rasters of a specific region over a period of one year in python. I also have a mask raster of the same area that is a binary color (250 black pixels for disparate areas of interest, nothing for everything else).

I'm trying to use a black pixel from the mask raster as bounds to grab a single pixel over time (52 rasters/weeks) as a numpy array, and I'd like to do this for all 250 black pixels, grabbing their corresponding pixels over time into separate arrays.

Unfortunately I cannot seem to find a method to pull this off. The closest thing I've found is using a shapefile to clip a raster. Any suggestions on how to achieve this in python would be golden.

  • In particular GIS?
    – FelixIP
    Aug 22 '18 at 18:59
  • Yeah trying to do temporal analysis on the pixels. Aug 22 '18 at 19:05
  • Do the 52 rasters and the mask have the same size/boundaries? Sep 29 '18 at 10:31
  • They have the exact same dimensions. Sep 29 '18 at 10:47

Hope this helps!

import rasterio
import numpy as np
import os

mask_path = r"mask.tif"
tiff_paths = [r"tiff1.tif", r"tiff2.tif", ...]
# tiff_paths = os.listdir(folderpath_with_tiffs)

with rasterio.open(mask_path) as src:
    mask = src.read(1)

stack = None
for tiff_path in tiff_paths:
    with rasterio.open(tiff_path) as src:
        array = src.read(1)
        extract = array[mask==0]  # [3 4]
        if stack is None:
            stack = extract.copy()
            stack = np.vstack((stack, extract))


If this is our binary mask

[[0 1 1]
 [0 1 1]
 [1 1 1]] 

and this is the array of the first tiff raster

[[3 5 2]
 [4 7 6]
 [8 8 1]]

and we want to extract the pixel values at the positions where the mask is 0

[3 4]

we would get a stack like this after looping through all tiffs.

[[3 4]    # extract of first tiff
 [1 6]    # second
 [2 1]]   # third and so on

We can get the time series of the first mask pixel by

# [3 1 2]
  • Question wouldn't np.dstack be a better choice than np.vstack? Oct 1 '18 at 10:29
  • Really depends on how you want to arrange and later access your results array. np.dstack would give you a third array dimension to deal with though. Give np.column_stack() a try, this would give you one pixel time series per row. Oct 1 '18 at 12:18

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.