by
1 Solution
Accepted Solutions
by
by
by
by
by
XanderBakker
See below an example which I created from the slope of the regression line based on the 13 values. This means that if the values increment in time they will have a positive slope (red) and in case of a decreasing regression line a negative slope (green):
The code used to create this raster is:
#-------------------------------------------------------------------------------
# Name: stats_numpy_series.py
# Purpose:
#
# Author: Xander
#
# Created: 30-07-2017
#-------------------------------------------------------------------------------
def main():
print "load modules..."
import arcpy
import numpy as np
# template = r'C:\GeoNet\LoopMultiFolderCalcMeanSameFolder\Test_Data\2004A\{0}162004.tif'
template = r'C:\GeoNet\WUE\r{0}_WUE.TIF'
nodata = -3.4028235e+38
out_ras = r'C:\GeoNet\WUE\WUE_slope.tif'
print "create numpy list..."
lst_np_ras = []
for i in range(1, 14):
ras_path = template.format("%03d" % (i,))
print " - ", ras_path
ras_np = arcpy.RasterToNumPyArray(ras_path)
lst_np_ras.append(ras_np)
print "read props numpy raster..."
ras_np = lst_np_ras[0]
rows = ras_np.shape[0]
cols = ras_np.shape[1]
print "create output numpy array..."
ras_path = template.format("%03d" % (1,))
raster = arcpy.Raster(ras_path)
ras_np_res = np.ndarray((rows, cols)) # np.ndarray((raster.height, raster.width))
print "loop through pixels..."
pix_cnt = 0
for row in range(rows):
for col in range(cols):
pix_cnt += 1
if pix_cnt % 10000 == 0:
print "row:", row, " col:", col, " pixel:", pix_cnt
lst_vals = []
for ras_np in lst_np_ras:
lst_vals.append(ras_np[row, col])
lst_vals = ReplaceNoData(lst_vals, nodata)
# perform calculation on list
a, b, sum_abs_dif, avg_abs_dif, cnt_elem = RegressionLine(lst_vals)
if cnt_elem > 1:
val = a
## print " - ", lst_vals
## print " - slope (a):", a
## print " - intercept (b):", b
## print " - total distance from regression line:", sum_abs_dif
## print " - mean distance from regression line :", avg_abs_dif
else:
val = nodata
try:
ras_np_res[row, col] = val
except:
print "err"
pnt = arcpy.Point(raster.extent.XMin, raster.extent.YMin)
xcellsize = raster.meanCellWidth
ycellsize = raster.meanCellHeight
ras_res = arcpy.NumPyArrayToRaster(ras_np_res, lower_left_corner=pnt, x_cell_size=xcellsize,
y_cell_size=ycellsize, value_to_nodata=nodata)
ras_res.save(out_ras)
arcpy.DefineProjection_management(in_dataset=out_ras, coor_system="GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
def ReplaceNoData(lst, nodata):
res = []
for a in lst:
if a < nodata / 2.0:
res.append(None)
else:
res.append(a)
return res
def RegressionLine(lst):
from numpy import arange, array, ones, linalg
months = range(1, 14)
lst_x = []
lst_y = []
cnt = 0
for a in lst:
if not a is None:
lst_x.append(months[cnt])
lst_y.append(a)
cnt += 1
cnt_elem = len(lst_x)
if cnt_elem >= 2:
A = array([ lst_x, ones(len(lst_x))])
# linearly generated sequence
w = linalg.lstsq(A.T, lst_y)[0]
# obtaining the parameters of the regression line
a = w[0] # slope
b = w[1] # offset
# calculate acumulated abs difference from trend line
zip_lst = zip(lst_x, lst_y)
sum_abs_dif = 0
for xy in zip_lst:
x = xy[0]
y = xy[1]
calc = x * a + b
abs_dif = abs(y - calc)
sum_abs_dif += abs_dif
avg_abs_dif = sum_abs_dif / cnt_elem
return a, b, sum_abs_dif, avg_abs_dif, cnt_elem
else:
return None, None, None, None, None
if __name__ == '__main__':
main()
I have attached the raster for you to examine.
Kind regards, Xander
