Creating loop in raster algebra which compares each cell to values of other rasters in time series?

I have a list of rasters representing NDVI values (and plan to do similar thing with some other index values as well) and in that raster I also set all values representing clouds to NoData.

Now, for each raster I want to calculate difference in values between that raster and previous one BUT if the previous raster has NoData value (so it's a cloud) in a cell, I want to skip that raster and get values from next previous raster (and if that raster also has NoData then skip to the next and so on) and then calculate delta values.

So with ArcPy it would look something like this:

def DeltaNDVI(pre, post):
DNDVI = Raster(post) - Raster (pre)
return DNDVI
for i in range (1, len(NDVIList)):
preprepre = NDVIList[i-3]
prepre = NDVIList[i-2]
pre = NDVIList[i-1]
post = NDVIList[i]
DeltaRaster = Con(IsNull(pre), Con(IsNull(prepre), DeltaNDVI(preprepre, post),DeltaNDVI(prepre, post)), DeltaNDVI(pre, post))

This obviously isn't an optimal solution since the list I have contains a dozen rasters and there are multiple lists with different number of rasters in them. There should be some kind of a loop there but I don't think that I can use loop with ArcGIS map algebra.

The list is sorted based on a date and I have date of acquisition for each raster in date(day, month, year) format so calculating difference in days is not a problem(and I also know that the usual interval for Landsat 8 scenes is 16 days, but my list doesn't have all rasters from that year so the difference between 2 time adjacent rasters in a list can be for example, 32 instead of 16 days).

So with that in mind I'd like to get another band as an output where I would se the difference in days for each delta pixel. As you can see in the image, the target raster is "input4" (or in the code above "post"). If "input3"(or "pre") has actual value in its cell then output cell equals input4-input3, but if input3 has NoData (ND) in a cell then output cell equals input4-input2 (unless input2 also has NoData) and so on.

if the target raster in question was "input3" then it would only be compared to input2 and input1 and first raster in a list doesn't get compared to anything obviously. If there would be cloud on the same spot in every raster(quite unlikely) then final output would be NoData. And if the target raster itself has a cloud then output is also NoData.

Band 2 equals to difference in days for the calculated delta values (and in this case the time cycle between two Landsat 8 scenes is 16 days).

Can this be done with ArcPy script?

• If you wish to also investigate doing this using GDAL then please do that in a separate question that is focused on what you have tried and where you are stuck with GDAL. – PolyGeo Mar 16 at 21:12
• Looking at you image it seem to me you want to collapse rasters 1,2&3 into a single intermediate raster which is the value you want to take away from raster 4. Have a look such tools like Diff and Test in the Logical toolset in Spatial Analyst toolbox. I'm sure a combination of those will generate this intermediate raster – Hornbydd Mar 16 at 21:20