F-test for trend analysis

Posted on Edited on In
Symbols count in article: 0 Reading time ≈ 1 mins. 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
# -*- coding: utf-8 -*-
import arcpy
import math
from arcpy.sa import *
arcpy.CheckOutExtension("spatial")

def F_TestForTrendAnalysis(dataFolder,resultSaveFolder):
    workSpace1 = dataFolder
    arcpy.env.workspace = unicode(workSpace1, "utf8")
    ys = arcpy.ListRasters("1*", "tif")
    # Generate time series layers
    i = 0
    for y in ys:
        i = i + 1
        out = dataFolder+"\\"
        if i < 10:
            out1 = unicode(out, "utf8") + "20" + str(i) + ".tif"
        else:
            out1 = unicode(out, "utf8") + "2" + str(i) + ".tif"
        outx = Con(1, i, y)
        outx.save(out1)
    # Read the time series layer
    arcpy.env.workspace = unicode(workSpace1, "utf8")
    xs = arcpy.ListRasters("2*", "tif")
    # Calculate
    x_ = CellStatistics(xs, "MEAN", "DATA")
    y_ = CellStatistics(ys, "MEAN", "DATA")
    n = len(xs)
    xy = 0
    x2 = 0
    St = 0
    for i in range(0, len(xs)):
        St = (Raster(ys[i]) - y_) ** 2 + St
        xy = Raster(xs[i]) * Raster(ys[i]) + xy
        x2 = Raster(xs[i]) ** 2 + x2
    k = (xy - n * x_ * y_) / (x2 - n * (x_ ** 2))
    outpath1 = resultSaveFolder+"\\"
    sty1 = "Trend" + ".tif"
    out1 = outpath1 + sty1
    k.save(out1)
    # F-test
    Sr = (x2 - n * (x_ ** 2)) * (k ** 2)
    Se = St - Sr
    F = (n - 2) * Sr / Se
    outpath2 = resultSaveFolder+"\\"
    sty2 = "F" + ".tif"
    out2 = outpath2 + sty2
    F.save(out2)