[GIS] Performing per-pixel linear regression over multiple rasters with NoData

landsatndviregressiontime series

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 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).

Best Answer

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.
Related Question