# Sequentially distilling several raster data sets into final map based on first time pixel’s value was greater than zero?

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
rasters.append(rast)

# 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)
else:
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)
``````