Hi Xander Bakker, Thank you for the suggestions.
I tried the same, I just changed the function in line no 53 and 102 as using Numpy and Scipy both. but I getting the error.
Whatever line I have updated based on the source line which is marked by #
Error regarding 'tuple index out of range'. I have attached below updated code and Error massage. KIndly have a look
I just want to modify our earlier script to add one more parameter as p value raster. Earlier script we have calculated the parameters such ['slope', 'intersept', 'tot_dist', 'mean_dist']. I just want to remove the 'tot_dist', 'mean_dist' and it will replace by ['r_value', 'p_value', 'std_err'] remain will be the same. So the final output raster would be ['slope', 'intersept', 'r_value', 'p_value', 'std_err']
Updated script:
def main():
print "load modules..."
import arcpy
from scipy import stats
from numpy import arange,array,ones
arcpy.env.overwriteOutput = True
import numpy as np
template = r'D:\NPP_GUJARAT_MONSOON\MODIS_ET\WUE4\r{0}_WUE.TIF'
nodata = -3.4028235e+38
out_template = r'D:\NPP_GUJARAT_MONSOON\MODIS_ET\WUE4\WUE_{0}.TIF'
out_names = ['slope', 'intersept', 'r_value', 'p_value', 'std_err']
print "create numpy list..."
lst_np_ras = []
for i in range(2, 13):
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 " - rows:", rows
print " - cols:", cols
print "create output numpy array..."
ras_path = template.format("%03d" % (2,))
raster = arcpy.Raster(ras_path)
lst_ras_np_res = []
for out_name in out_names:
print "out_name", out_name
ras_np_res = np.ndarray((rows, cols))
lst_ras_np_res.append(ras_np_res)
print " - out rows:", ras_np_res.shape[0]
print " - out cols:", ras_np_res.shape[1]
print "loop through pixels..."
pix_cnt = 0
for row in range(rows):
for col in range(cols):
pix_cnt += 1
if pix_cnt % 5000 == 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)
a, b, sum_abs_dif, avg_abs_dif, cnt_elem = stats.linregress(lst_vals)
if cnt_elem > 1:
slope = a
intersept = b
lst_ras_np_res[0][row, col] = slope
lst_ras_np_res[1][row, col] = intersept
lst_ras_np_res[2][row, col] = sum_abs_dif
lst_ras_np_res[3][row, col] = avg_abs_dif
else:
lst_ras_np_res[0][row, col] = nodata
lst_ras_np_res[1][row, col] = nodata
lst_ras_np_res[2][row, col] = nodata
lst_ras_np_res[3][row, col] = nodata
pnt = arcpy.Point(raster.extent.XMin, raster.extent.YMin)
xcellsize = raster.meanCellWidth
ycellsize = raster.meanCellHeight
print "Write output rasters..."
i = 0
for out_name in out_names:
out_ras = out_template.format(out_name)
print " - ", out_ras
ras_np_res = lst_ras_np_res[i]
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]]")
i += 1
def ReplaceNoData(lst, nodata):
res = []
for a in lst:
if a < nodata / 2.0:
res.append(None)
else:
res.append(a)
return res
stats.linregress(lst_x, lst_y)
from numpy import arange, array, ones, linalg
months = range(2, 13)
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))])
w = stats.linregress(A.T, lst_y)[0]
a = w[0]
b = w[1]
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()
Error message:
Traceback (most recent call last):
File "<module1>", line 142, in <module>
File "<module1>", line 61, in main
File "C:\Python27\ArcGIS10.3\lib\site-packages\scipy\stats\stats.py", line 2997, in linregress
elif x.shape[1] == 2:
IndexError: tuple index out of range
