Latest Contributions by ShouvikJha

Hi Xander Bakker Thank you very much for such wonderful script. However I am running the script on my system and will inform you shortly.   ... View more

Hi Xander Bakker‌, I am facing same problem here also.  we eagerly await your response Thanks & regards  Shouvik Jha  ... View more

Hi xander_bakker‌, I am still stuck into this matter to modify the code based on SciPy.  My kind request please cooperate to solve the issue.   ... View more

Dan_Patterson‌ Thanks.  Yes, Scipy is more efficient to produce the statistical result and P value tells the how significant the trend. The significance of trend is more important while you interpret any output of trend or correlation. Everybody wants to know the How significance relation.   ... View more

Hi Xander Bakker‌, Thank you for the suggestion.  I trying to modify the code as you suggested but i don't understand, where I have to call the p Value raster.  From our earlier script, We got the output only Correlation raster, but this time I want Correlation raster as well as P value raster.  Snap from our earlier script on Input and output:   import numpy as np     template1 = r'D:\Correlation\Parameter1\r{0}_NPP.TIF'     template2 = r'D:\Correlation\Parameter2\r{0}_WUE.TIF'     nodata = -3.4028235e+38     out_ras_Corr = r'D:\Correlation\correlation.TIF' # Already exist in our script     out_ras_p    = r'D:\Correlation\Significance.TIF' # this one i want to add     ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ ... View more

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:  #------------------------------------------------------------------------------- # Name:        Regression Analysis.py  # Purpose:     Regression Analysis through time series data # Author:      Xander Bakker # Source:      Geonet # Created:     30-07-2017 #------------------------------------------------------------------------------- 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', 'tot_dist', 'mean_dist']     out_names = ['slope', 'intersept', 'r_value', 'p_value', 'std_err']     print "create numpy list..."     lst_np_ras = []     for i in range(2, 13): # 2 number will start of raster name e.g. 2 = r002_WUE, 3 = r003_WUE         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,)) # 2 number will take the number of raster name e.g. 2 = r002_WUE, 3 = r003_WUE     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)             # perform calculation on list             # a, b, sum_abs_dif, avg_abs_dif, cnt_elem = RegressionLine(lst_vals)             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) #  - raster.meanCellHeight     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 # def ReplaceNoData(lst, nodata):     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))])         # linearly generated sequence         # w = linalg.lstsq(A.T, lst_y)[0]         w = stats.linregress(A.T, lst_y)[0]         # obtaining the parameters of the regression line         a = w[0]  # slope         b = w[1]  # offset         # calculate acumulated abs difference from trend line         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‍‍‍‍‍‍ ... View more

Hi Xander Bakker‌, Recently we got a very useful script from you regarding regression analysis.  Here, I would like to add one more parameter in this script which is P Value (calculate the significance). In the above script, we have been used the Numpy library numpy.linalg.lstsq which gives the output as Slope, intercept etc. However, P value is not calculated. P value also needed because it tells the how significant increasing and decreasing the trend while we interpret the output statistically in front of the scientific community. So while reading the python library, I found the Scipy library scipy.stats.linregress and i just ran one  script on my system and i got the slope, P value. Script is attached below  from numpy import arange,array,ones#,random,linalg from scipy import stats xi = arange(0,11) A = array([ xi, ones(11)]) # linearly generated sequence y = [154.96, 151.39, 143.87, 152.24, 142.97, 154.23, 141.15, 135.37, 151.50, 143.82, 143.77] slope, intercept, r_value, p_value, std_err = stats.linregress(xi,y) print 'slope', slope print  'p_value', p_value ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ which gives the output as  Returns: slope -0.899090909091 slope of the regression line p-value : 0.142283969084 two-sided p-value for a hypothesis test whose null hypothesis is that the slope is zero.  So instead of Numpy, can we use here Numpy and Scipy both library to get the addition p value or in the existing script can we add this function to obtain the p Value as a raster.  Please guide us to solve this.  Thanks & regards  Shouvik Jha  ... View more

Nicely done. Great Work  ... View more

Hi xander_bakker‌, thank you very much. The script is very useful to generate statistics from the raster.  ... View more

Hi Xander Bakker   Yes, now it's appearing perfectly and got exact correct pixel calculation. thanks for such wonderful script and it will save the time and also reduce the chances of human introduce error.  And the last column (acz_raster) did you type it manually into excel, based on the zone sequence in the zone raster or it's automatically come through script? Because for me it's not appearing. it looks like (below attached) after paste into excel, after "cnt_nodata", there is no other column. The table which you attached there is complete information. but when i am running not getting the same, is there anything need to change?  Thank you   regards   Shouvik Jha   ... View more

Hi Xander Bakker‌, Thanks.  I ran the script using same data sets which you used, after running the script I am getting such output which is not matching with your output file which you attached and the Zone column 1 is not appearing. the last column acz_raster also not appearing in my output. I copied the output from PyScripter after running and paste it into Excel.  Below I have attached the output which I received after running the script.  I am running the script on PyScripter and ArcGIS 10.3.  Zone     sum     cnt     mean     min     max     cnt_neg     cnt_pos     cnt_nodata 2      568906,785936     800059     0,711081040194     -5,48433862166     9,41357835423     91714     626752     774 3      304342,163442     389143     0,782083099122     -6,14446017525     4,65285074061     32275     327596     407 4      75973,486289     254708     0,298276796524     -9,89411033284     16,0140284267     64581     146669     1436 5      138842,463549     182611     0,760318182087     -6,15788269043     11,1494489421     26080     147168     157 6      505731,281974     691091     0,731786815303     -8,45667485757     9,71765823364     3174     273474     16118 7      60633,5277998     113542     0,534018493595     -16,3086192868     20,4710965796     18460     85396     495 255      10,6494542035     62     0,171765390379     0,0     2,51514136574     0     8     2371928‍‍‍‍‍‍‍‍ ... View more

Hi Xander Bakker, thanks.   I am mentioning below the input path which I am using and error message also. I tried both with gdb and without gdb in the input path but the error message is same.   However, I have attached my input file both raster and Zone vector file with the previous thread. location Just after the script which you posted.       Input path name:  import numpy as np     zones_ras = r'D:\Raster_Statistics\Count_Pixel\Zones'     values_ras = r'D:\Raster_Statistics\Count_Pixel\NPP' ‍‍‍‍ Error massage:  load modules... create numpy rasters... Traceback (most recent call last):   File "<module1>", line 70, in <module>   File "<module1>", line 12, in main   File "C:\Program Files\ArcGIS\Desktop10.3\ArcPy\arcpy\__init__.py", line 2244, in RasterToNumPyArray     return _RasterToNumPyArray(*args, **kwargs) RuntimeError: ERROR 000732: Input Raster: Dataset D:\Raster_Statistics\Count_Pi ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ ... View more

Hi jayanta.poddar‌, Thanks. I tried the same, but still same error. i think problem with Name of zone field in feature class.  ... View more

Hi Xander Bakker‌ Thank you.  I ran the code but the same error I am getting which I mentioned above. Please, could you check, where is the problem.  I have made the gdb file using below code, import os import sys from arcpy import env import arcpy env.workspace = "D:\Raster_Statistics\Count_Pixel" # Set local variables out_folder_path = "D:\Raster_Statistics\Count_Pixel" out_name = "Count.gdb" # Execute CreateFileGDB arcpy.CreateFileGDB_management(out_folder_path, out_name)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ While I using the Count_Pixel.gdb as an input path, that time it's producing the error. I think somewhere name of Zone is not matching.   if you can attach the data (Zone_ras and value_ras) which you used for generating the result, so I can compare and check with my data. Thank you   ... View more