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
Solved! Go to 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.]) >>>
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
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)
what would be value if it's (5 - 9), is it -1 or -4 ?
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.]) >>>
Thanks Dan Patterson, for help, i have bug in my program, but i find it successfully...
Ahsan... Good! ...you should post the result the result so others may find it
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
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.
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 ?