I have a 9 year weekly time series (~500 raster grids of equal cell size and extent). I'm interested in obtaining the regression line slope between pixels (Imagine stacking all 500 grids on top of one another and running a linear regression between each individual pixel). Essentially, this would measure the net change between pixels through my time series.
My main question is about the coefficients defined in the first answer to Making linear regression across multiple raster layers using ArcGIS Desktop?
Is the coefficient defined as 12 due to the temporal scale of the original question being at the monthly level and, if so, should my code's coefficient be set to 7 due to the weekly temporal scale between grids?
I have been referencing these Q&As that are similar:
I wrote this python code after reading the links referenced above:
from __future__ import division
import arcpy
import os
from arcpy import env
from arcpy.sa import*
# define workspace
arcpy.env.workspace=r"C:\Users\mbs7038\Documents\Remove_Gauge_Outliers\test\forSlope.gdb"
# enable overwriting
arcpy.env.overwriteOutput=True
# check spatial analyst extension
arcpy.CheckOutExtension('Spatial')
# define output paths
slopePath=r"C:\Users\mbs7038\Documents\Remove_Gauge_Outliers\SpatialTrends\SlopeTEST.gdb"
# list all rasters in the workspace
rasters=arcpy.ListRasters('*')
# sort rasters numerically
rasters.sort()
# get the number of rasters
n=len(rasters)
print n
# setup index
i=1
# define division
seed=(n*n*n)-n
print 'the global seed is {0}'.format(seed)
for raster in rasters:
print i
coef=(12*(i)-((6*n)+(6*1)))/seed
print 'Raster {0} is {1}:'.format(i,raster)
print 'the coef for raster {0} is {1}'.format(i,coef)
# Multiple raster by coefficient
if i==1:
outSlope=(Raster(raster)*coef)
i+=1 # same as saying i=i+1
else:
print 'adding {0} to outSlope'.format(raster)
outSlope=outSlope+(Raster(raster)*coef)
i+=1
if i==6:
break
# Save final slope grid
print 'saving final slope grid'
outSlope.save(slopePath + "/" + "floodChange_py")
print 'script is complete'
I created this code as a test which appeared to work on a 5 week subset of my data. If what I have referenced above is correct then this code will work on a time series (of equal cell size and extent) of any length. Although there were no bugs in the code, I have a feeling I may have made a mistake someplace.
On a side note I'm using PyScripter, hence the "from future import division" in the first line. This disables the automatic floor function which rounds floating point numbers to the nearest integer away from zero.
No comments:
Post a Comment