4

I have a stack of 10 Landsat satellite images over the same area that I have calculated NDVI for. I want to calculate NDVI trend over time on a per-pixel basis. This means that for each pixel in the images I want a value for the slope of the regression NDVI vs. Time.

From a previous question asked here I came across the super cool Curve Fit ArcGIS add on. I input my NDVI rasters into the tool, and asked for an output being linear regression slope and p-value. I got those results, but only for a small subset of my study area. I found that any pixel set as NoData (aka areas I had masked due to cloud cover) led to that entire pixel stack being out put as NoData. So even if 1 out of 10 values in a pixel stack was NoData the linear regression would fail and the output would be NoData.

My goal is to find a way to run linear regressions on all the pixels in my study area, even if one or more of the values in a pixel stack are masked due to cloud cover (meaning if 1 value is set to NoData then the linear regression will simply be out of 9 values out of 10 in a stack, rather than just being output as NoData).

  • 1
    The best would be to script it, rather easy by using Raster To Numpy Array function. If this is not an option perhaps try replacing NO DATA pixels by cell statistics MEAN – FelixIP Dec 6 '16 at 0:57
  • I don't want to overstate the obvious but are you sure there is no 'ignore value' or similar option that would allow you to skip 0 values? – Nathan Thomas Dec 6 '16 at 5:03
2

Using technique described here I populated the table of integer grid by values from 3 rasters of interest.

I used IS NULL query to populate records with no match (NO DATA) by -999 during data transfer from zonal statistics table. Finally I added field “SLOPE” and computed it using Python field calculator expression: import numpy

def getSlope(yList):
    x,y=[],[]
    years=[1984, 1990, 1991, 1992, 1996, 2002, 2006, 2011, 2013,2016]
    for i,v in enumerate(yList):
        if v==-999:continue
        x.append(years[i]);y.append(v)
    ab=numpy.polyfit(x, y, 1)
    return ab[0]
#------------------------------
getSlope([ !GRD_01!, !GRD_2!, !GRD_03!])

Assuming time interval between rasters is the same.

RESULT:

enter image description here

NOTE:

  • Values transfer and slope calculations both based on field calculator, which makes it a bit slow
  • Use LOOKUP in spatial analyst Reclass to convert SLOPE field to raster.
  • Thanks for the response! A colleague in my department is proficient with R and is giving me a hand with my code but I am sure that the correct work flow is somewhat similar to your script so it should be helpful. My rasters are not from subsequent years (the 10 years are 1984, 1990, 1991, 1992, 1996, 2002, 2006, 2011, 2013 and 2016) so that is a wrinkle I will have to figure out as well. – Mitchell Dec 7 '16 at 0:42
  • Easy to fix, see update – FelixIP Dec 7 '16 at 1:04

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.