I have 120 daily data sets that contain a value only when a threshold temperature is exceeded. If the threshold has been reached that pixel was assigned a value indicating the day it occurred (the Julian day), all the rest of the cells for that day were assigned a value of zero. I’d like to combine these raster data sets to create a final map indicating the Julian day when the threshold was first reached or exceeded.

I could use the raster calculator (in QGIS it would look something like this):

(“julian1” = 1) * “julian1” + (“julian2” = 2) * “julian2” + (“julian3” = 3)* “julian3” 

up to 120 but I have 30 years of data parsed into 120 day blocks and it would take a really long time. I’ll bet there is a way better solution in either QGIS or ArcMap or a python script.

  • ArcGIS cell statistics - maximum will do if I understand data structure correctly – FelixIP Jan 19 '16 at 21:48

You could build a raster stack from the 120 layers with numpy and then apply a function to each cell to determine the first non-zero position. Don't know about the speed for large rasters though.

import os
import numpy as np

# create some random rasters
ncol = 20
nrow = 10
for x in range(1,6):
    rast = np.random.randint(0,2,(nrow,ncol)) * x
    rast[rast<x] = 0

# build raster stack
stackRast = np.dstack(rasters)

# define function to find first non-zero
def my_func(a):
    if max(a)>0:
        b = (a>0).tolist().index(1)
        b = -1
    print a, b
    return b

# apply that function to the temporal axis of the raster stack
firstNonZero = np.apply_along_axis(my_func, 2, stackRast)

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.