topic Re: Mann Kendall Test in Python Questions

Mann Kendall Test

AhsanAbbas — Sun, 12 Dec 2021 02:39:37 GMT

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

Re: Mann Kendall Test

DanPatterson_Retired — Sat, 27 Dec 2014 09:39:58 GMT

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

Re: Mann Kendall Test

DanPatterson_Retired — Sun, 12 Dec 2021 02:39:40 GMT

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)

Re: Mann Kendall Test

AhsanAbbas — Sat, 27 Dec 2014 11:43:45 GMT

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

Re: Mann Kendall Test

DanPatterson_Retired — Sun, 12 Dec 2021 02:39:43 GMT

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

Re: Mann Kendall Test

AhsanAbbas — Sun, 28 Dec 2014 07:05:53 GMT

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

Re: Mann Kendall Test

DanPatterson_Retired — Sun, 28 Dec 2014 10:21:49 GMT

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

Re: Mann Kendall Test

AhsanAbbas — Sun, 28 Dec 2014 10:27:54 GMT

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

Re: Mann Kendall Test

XanderBakker — Sun, 12 Dec 2021 02:39:46 GMT

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.

Re: Mann Kendall Test

AhsanAbbas — Tue, 30 Dec 2014 06:50:28 GMT

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 ?

Re: Mann Kendall Test

XanderBakker — Tue, 30 Dec 2014 23:14:33 GMT

I made some minor changes to the code:

I put the program flow in the def main
changed the code to create the list of pairs
used str.format to create the output raster name (please note that if the inputs are .img, you can better strip the extension to get a valid output name)
the numpy.sign is executed on the entire array, not for each cell individually
I used the os.path.join to create the output path, name for the raster

To assign an output projection, you will probably have to set the output coordinate system in the environments. The numpy to raster conversion seems to honor this setting (arcpy.env.outputCoordinateSystem):

ArcGIS Help (10.2, 10.2.1, and 10.2.2)

Re: Mann Kendall Test

AhsanAbbas — Wed, 31 Dec 2014 07:11:24 GMT

I tried this method but i get error, one problem more, the spatial extent of output raster is different, mean height is different but mean cell width is same, i don't understand ?

Re: Mann Kendall Test

XanderBakker — Thu, 01 Jan 2015 06:57:29 GMT

Do all the input rasters have the same dimension (width, height, resolution and location)? If not, you will have to account for these differences.