Mann Kendall Test

[AhsanAbbas]

Hello Guys, I need to perform Mann Kendall test on 5 series raster data yearly basis, from 2007 to 2011.

if need to iterate pixel by pixel on raster cell values for calculating the sign. How can i do it... Below is code but i have problem with numpy...

import arcpy
import numpy
from arcpy import env


env.workspace = r"D:\PythonProject\Data"
rasterlist = arcpy.ListRasters()


for raster in rasterlist:
    print raster


n = len(rasterlist)
pairs = (n * ( n - 1)) / 2
print "Pairs to be formed: ", pairs


def pairsRasters(raster1, raster2):
    print raster1, " - ", raster2
    
    inRas1 = arcpy.Raster(raster1)
    lowerLeft1 = arcpy.Point(inRas1.extent.XMin,inRas1.extent.YMin)
    #cellSize = inRas.meanCellWidth
    my_array1 = arcpy.RasterToNumPyArray(inRas1,nodata_to_value=0)


    inRas2 = arcpy.Raster(raster2)
    lowerLeft2 = arcpy.Point(inRas2.extent.XMin,inRas2.extent.YMin)
    #cellSize = inRas.meanCellWidth
    my_array2 = arcpy.RasterToNumPyArray(inRas2,nodata_to_value=0)
    


for current in reversed(rasterlist):
    for previous in rasterlist:
        if current == previous:
            pass
        elif current > previous:
            pairsRasters(current, previous)
        else:
            pass

Suppose the two raster of similar extent, i need to calculate difference, if it's greater than 1 the value store in output raster would be 1, if less than 1 the value would be -1 and else zero. How can i iterate on pixel by pixel of two rasters. so i can save it into new raster. That's it..

Column 1	Column 2	Column 3	Column 4	Column 5
50	25	98	102	98
122	98	10	210	100
57	66	92	21	10
98	45	255	9	59
58	98	17	31	54

Message was edited by: Ahsan Abbas

[DanPatterson_Retired]

can you show a sample of your inputs and your output...It seems you are comparing the whole rasters once without any just checking to see whether they are the same.

Also see the Mann-Kendall implementation here

[DanPatterson_Retired]

Not familiar with the test that much, but if you are only interested in the signs then you can look at this example

>>> a
array([[ 1.,  2.],
       [ 3.,  4.],
       [ 5.,  6.]])
>>> b
array([[ 1., -2.],
       [ 3.,  4.],
       [ 5., -6.]])
>>> sign_ab = np.sign(a-b)
>>> sign_ab
array([[ 0.,  1.],
       [ 0.,  0.],
       [ 0.,  1.]])
>>>

You can iterate from there is you like to see if the sign changes through a sequence of rasters (oh yes, 0 means the sign didn't change 1, it did)

[AhsanAbbas]

what would be value if it's (5 - 9), is it -1 or -4 ?

[DanPatterson_Retired · Accepted solution]

try the following to see if this is what you want

3 arrays with differing values and sign

>>> a
array([ 6.,  5.,  5.,  3.,  1.,  1.])
>>> b
array([ 1., -2.,  3.,  4.,  5., -6.])
>>> c
array([-3., -2., -1.,  4., -5., -6.])
>>> np.sign(b-a)
array([-1., -1., -1.,  1.,  1., -1.])
>>> np.sign(c-b)
array([-1.,  0., -1.,  0., -1.,  0.])
>>> np.sign(a-c)
array([ 1.,  1.,  1., -1.,  1.,  1.])
>>>

[AhsanAbbas]

Thanks Dan Patterson‌, for help, i have bug in my program, but i find it successfully...

[DanPatterson_Retired]

Ahsan... Good! ...you should post the result the result so others may find it

[AhsanAbbas]

I will post it when it will be complete, half done yet...

The half code....Dan Patterson‌, please check it...

import arcpy

import numpy

from arcpy import env

import numpy as np

arcpy.env.overwriteOutput = True

env.workspace = r"D:\PythonProject\Data"

rasterlist = arcpy.ListRasters()

for raster in rasterlist:

    print raster

n = len(rasterlist)

pairs = (n * ( n - 1)) / 2

print "Pairs to be formed: ", pairs

def pairsRasters(raster1, raster2):

    print raster1, " - ", raster2

    desc1 = arcpy.Describe(raster1)

    desc2 = arcpy.Describe(raster2)

    nameRaster = desc1.basename + " - " + desc2.basename

    #print nameRaster   

    array1 = arcpy.RasterToNumPyArray(raster1,nodata_to_value=0)

    array2 = arcpy.RasterToNumPyArray(raster2,nodata_to_value=0)

    #print "\n ", array1.shape

    # to see how many dimension the array has

    (height, width) = array1.shape

    newRaster = numpy.zeros([height, width])

    for row in range(height):

        for col in range(width):

            #print "Row, Column: ", row, col

            #print float(array2[row][col]), float(array1[row][col])

            #print (float(array2[row][col]) - float(array1[row][col]))

            newRaster[row][col] = np.sign(float(array2[row][col]) - float(array1[row][col]))

    finalRaster = arcpy.NumPyArrayToRaster(newRaster)

    outname = "D:/PythonProject/Outputs/" + (nameRaster + ".img") 

    #arcpy.CopyRaster_management(name, outname ,"DEFAULTS","0","9","","","8_BIT_UNSIGNED")

    finalRaster.save(outname)

for current in reversed(rasterlist):

    for previous in rasterlist:

        if current == previous:

            pass

        elif current > previous:

            pairsRasters(current, previous)

        else:

            pass

[XanderBakker]

Although I have no idea on what the Mann Kendall Test will get you, it seems that the code below works.

import arcpy
import numpy
import os

def main():
    ws_in = r"D:\PythonProject\Data"
    ws_out = r"D:\PythonProject\Outputs"

    arcpy.env.overwriteOutput = True
    arcpy.env.workspace = ws_in
    rasterlist = arcpy.ListRasters()

    n = len(rasterlist)
    pairs = (n * ( n - 1)) / 2
    print "Pairs to be formed: {0}".format(pairs)

    for i in rasterlist:
        for j in rasterlist[rasterlist.index(i)+1:]:
            print i, j
            pairsRasters(i, j, ws_out)


def pairsRasters(raster1, raster2, ws):
    desc1 = arcpy.Describe(raster1)
    desc2 = arcpy.Describe(raster2)
    nameRaster = "{0}_{1}.img".format(desc1.basename, desc2.basename)
    print nameRaster

    array1 = arcpy.RasterToNumPyArray(raster1,nodata_to_value=0)
    array2 = arcpy.RasterToNumPyArray(raster2,nodata_to_value=0)

    newRaster = numpy.sign(array2 - array1)
    finalRaster = arcpy.NumPyArrayToRaster(newRaster)

    outname = os.path.join(ws, nameRaster)
    finalRaster.save(outname)


if __name__ == '__main__':
    main()

You may want to check the names of the input rasters. If they have extensions, you may want to strip them to create a valid output name.

[AhsanAbbas]

Xander Bakker ,The Menn Kendall test is non-parametric test, use to find the trend of data either increasing, decreasing or remains same. I just want to find the trend, And what is wrong with my code is it not understand able ?

I also wants to add these save rasters to one raster, and assign the projection of original raster, Where should i change the code ?