topic Re: Trend Analysis Through Time Series of Raster Data in Python Questions

Trend Analysis Through Time Series of Raster Data

ShouvikJha — Thu, 27 Jul 2017 18:14:17 GMT

Hi All,

I have total 13 raster files for different years (2001 - 2013). All data are in tiff format and same spatial extent (row and column). I am trying to perform regression line slope/trend analysis between each grid points for 13 raster data sets. I have tried in ArcGIS raster calculator but I won't run the task due to complexity.

The outcome would measure the net change between pixels through my time series data.

So my question is, can it be performed using raster calculator or ArcPy Script? For such type of statistical analysis using ArcGIS, really I am facing one of the big limitation with ArcGIS.

Using excel, it can be calculated using SLOPE command where X (Time) and Y (value) has been used. But for creating a Spatial map of Slope/trend analysis, we need to perform using raster data.

The download link is given below for Sample data sets (Due to size limit, I could not attach the data with here )

Data Link: Dropbox - Water_Use_Efficiency.zip

Please share your experience.

Thanks in advance.

Re: Trend Analysis Through Time Series of Raster Data

ShouvikJha — Fri, 28 Jul 2017 05:48:41 GMT

Xander Bakker Patterson‌ , A request to you peoples, We need your kind cooperation to solve the issue. Since long i am trying but not solved yet. Thanks

Re: Trend Analysis Through Time Series of Raster Data

XanderBakker — Fri, 10 Dec 2021 22:33:01 GMT

Hishouvik jha ,

Recently I needed to detect strange behavior in monthly consumption by water consumers. I created a data structure that had the 12 months of data for each client as attributes of a point file. The code below is a snippet that takes a list of data (this could be the 13 values for each pixel in your case) and calculates a regression line based on the data.

The slope will tell you if there is an increase or decrease in the range of data. In my case to detect this abnormal behavior I determined the average value, and the mean distance for each value to the regression line divided by mean value. The higher the value the less constant the consumption.

I think you could use this code (or parts) to create the numpy arrays for each raster and loop through the individual pixels (create a list of values for the 13 years you have for each pixel, do the calculation and fill an output array to write the result of the analysis to an output raster. I'm pretty sure there must be an easier way to do this, or even a standard way to calculate this on a time series of data. I will investigate a little more and let you know.

from numpy import arange,array,ones,linalg

def main():

    # some test input data (lists of 12 elements with possible NoData)
    lst_data = [[2.0, None, None, None, None, 1.0, 4.0, 6.0, 5.0, 7.0, 9.0, 9.0],
           [None, None, None, None, None, None, None, None, None, None, 15.0, 16.0],
           [14.0, 11.0, 25.0, 4.0, 4.0, 4.0, 3.0, 4.0, 2.0, 3.0, 4.0, 4.0],
           [7.0, 9.0, 7.0, 3.0, 2.0, 2.0, None, 13.0, 5.0, 60.0, 25.0, 25.0],
           [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0],
           [1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 1.0, 2.0],
           [11.0, 12.0, 11.0, 12.0, 11.0, 12.0, 11.0, 12.0, 11.0, 12.0, 11.0, 12.0]]

    # loop through input data
    for lst in lst_data:
        print "\ndata:", lst
        # calculate regression line: y = ax + b
        a, b, sum_abs_dif, avg_abs_dif = RegressionLine(lst)
        if not a is None:
            cons_prom = getAvg(lst)
            print " - mean value:", cons_prom
            print " - slope     (a):", a
            print " - intercept (b):", b
            print " - total distance from regression line:", sum_abs_dif
            print " - mean distance from regression line :", avg_abs_dif
            # print a, b, sum_abs_dif, avg_abs_dif
            print " - score:", avg_abs_dif / cons_prom
        else:
            print " - Insufficient data"


def RegressionLine(lst):
    from numpy import arange, array, ones, linalg
    months = range(1, 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]

        # 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
    else:
        return None, None, None, None

def getAvg(lst):
    corr_lst = [a for a in lst if not a is None]
    cnt = len(corr_lst)
    if cnt == 0:
        return None
    else:
        return sum(corr_lst) / float(cnt)

if __name__ == '__main__':
    main()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Results:

data: [2.0, None, None, None, None, 1.0, 4.0, 6.0, 5.0, 7.0, 9.0, 9.0]
 - mean value: 5.375
 - slope     (a): 0.738095238095
 - intercept (b): -0.529761904762
 - total distance from regression line: 9.29761904762
 - mean distance from regression line : 1.16220238095
 - score: 0.216223698782

data: [None, None, None, None, None, None, None, None, None, None, 15.0, 16.0]
 - mean value: 15.5
 - slope     (a): 1.0
 - intercept (b): 4.0
 - total distance from regression line: 3.5527136788e-15
 - mean distance from regression line : 1.7763568394e-15
 - score: 1.14603667058e-16

data: [14.0, 11.0, 25.0, 4.0, 4.0, 4.0, 3.0, 4.0, 2.0, 3.0, 4.0, 4.0]
 - mean value: 6.83333333333
 - slope     (a): -1.18181818182
 - intercept (b): 14.5151515152
 - total distance from regression line: 42.303030303
 - mean distance from regression line : 3.52525252525
 - score: 0.515890613452

data: [7.0, 9.0, 7.0, 3.0, 2.0, 2.0, None, 13.0, 5.0, 60.0, 25.0, 25.0]
 - mean value: 14.3636363636
 - slope     (a): 2.69171974522
 - intercept (b): -3.0101910828
 - total distance from regression line: 103.946496815
 - mean distance from regression line : 9.44968152866
 - score: 0.65788922035

data: [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0]
 - mean value: 6.5
 - slope     (a): 1.0
 - intercept (b): -1.46316339947e-15
 - total distance from regression line: 3.17523785043e-14
 - mean distance from regression line : 2.64603154202e-15
 - score: 4.07081775696e-16

data: [1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 1.0, 2.0]
 - mean value: 1.5
 - slope     (a): 0.020979020979
 - intercept (b): 1.36363636364
 - total distance from regression line: 5.87412587413
 - mean distance from regression line : 0.48951048951
 - score: 0.32634032634

data: [11.0, 12.0, 11.0, 12.0, 11.0, 12.0, 11.0, 12.0, 11.0, 12.0, 11.0, 12.0]
 - mean value: 11.5
 - slope     (a): 0.020979020979
 - intercept (b): 11.3636363636
 - total distance from regression line: 5.87412587413
 - mean distance from regression line : 0.48951048951
 - score: 0.042566129522‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Trend Analysis Through Time Series of Raster Data

XanderBakker — Fri, 28 Jul 2017 12:55:17 GMT

I also came across this tool: Curve Fit : A Pixel Level Raster Regression Tool - Desktop Decision Support Tools by the USGS

Re: Trend Analysis Through Time Series of Raster Data

ShouvikJha — Sat, 29 Jul 2017 07:28:35 GMT

Hi Xander Bakker‌, Thank you.

The script is quite interesting. in above script, there are seven client's data has been used for the analysis. Instead of using a numeric value which you assigned as an input data, it can be possible to assign raster file or ASC file in there? If is it possible to use raster or ASC file as an input data then output also can be written in raster format. (Plz give your opinion also)

Somewhere I have read that if we stack all the year's data and make it a single raster. Might be it will help us to perform a loop on all pixels (number of the pixels at same grid point for all years, so in my case at the same grid point, the number of pixels is 13) at same grid point for the raster based statistical analysis such as trend, correlation etc.

basically, Since long the geospatial research communities are also looking for the same analysis on Python or ArcGIS raster calculator platform, as following

1) Raster based Trend Analysis through time series data.

2) Raster based correlation analysis between two variables. Just assume the correlation analysis between Temperature and precipitation raster data and output also be in raster format.

If the above issues can be solved on Python platform so it would be a very good news for all researchers who are looking for the solution since long.

The tool which you suggested, It is very helpful and the number of functions is inbuilt with the tool but the only issue is, it is only functioning on ArcGIS 10.1 version. It is not compatible with the later version of ArcGIS 10.1 and taking too long time for the analysis. However, I am going through with that tool. In addition, I am trying to further explore the script which you provided.

It is our kind request behalf of the geospatial research community, India, please show a path to overcome the issues. If you need the sample data, I can provide you some datasets for further investigation on the matter.

Once again thank you very much for your cooperation.

With warm regards

Shouvik Jha

Re: Trend Analysis Through Time Series of Raster Data

XanderBakker — Sat, 29 Jul 2017 15:19:02 GMT

Hi shouvik jha ,

For your first question, is it possible to connect to a datasource, the answer is yes. In my case the data was available as point data and the 12 values as 12 fields for each feature. I read them into a list of values and process the results and write the output statistic values to additional field in the point feature class.

In your case you have a list of 13 raster where each raster aligns with the other raster and the number of rows and columns match. You will have to loop through each pixel location in a numpy array that was created from converting the each raster to a numpy, and for each pixel location you will read the 13 values out (from the 13 numpy arrays), construct the list, perform the analysis and populate a new numpy array as the output result.

It would be interesting to have access to the 13 rasters to see what can be done. So if you can share them, please attach them at the thread.

You could also explore the How Band Collection Statistics works—Help | ArcGIS Desktop tool, although this will not create a raster as a result, but a report using the pixels in a diagonal of the bands (you will need to check the "compute covariance and correlation matrices" option).

Re: Trend Analysis Through Time Series of Raster Data

ShouvikJha — Sat, 29 Jul 2017 18:26:18 GMT

Hi Xander Bakker‌, Thank you.

I have attached the Dropbox link at the thread for accessing the data. Due to the size limit, I could not attach the data with here, however, you can download the data by that link and please have a look.

Re: Trend Analysis Through Time Series of Raster Data

DanPatterson_Retired — Fri, 10 Dec 2021 22:33:04 GMT

You shouldn't have to cycle through the cells of a stacked raster .... consider

# ---- create a 'raster' consisting of 4 rows and 5 columns with 12 
#      months of data... ie. 12 months x 4 rows x 5 columns

a = np.arange(12*4*5).reshape(12, 4, 5)

a[0]  # ---- the first month
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19]])

a[-1] # ---- the last month ..... you get the drift I hope
array([[220, 221, 222, 223, 224],
       [225, 226, 227, 228, 229],
       [230, 231, 232, 233, 234],
       [235, 236, 237, 238, 239]])

top_left_12months = a[:, 0, 0]  # --- let's get the top-left cell's data
array([  0,  20,  40, ..., 180, 200, 220])  # ---- makes sense doesn't it?

np.mean(top_left_12months)  # ---- the 12 month mean for the top left cell
110.0

np.mean(a, axis=0)  # ---- Since I am lazy, I will get every cell's average

array([[ 110.00,  111.00,  112.00,  113.00,  114.00],
       [ 115.00,  116.00,  117.00,  118.00,  119.00],
       [ 120.00,  121.00,  122.00,  123.00,  124.00],
       [ 125.00,  126.00,  127.00,  128.00,  129.00]])

# ---- Notice the top left cell value is correct
#‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

The point being,... once you snag the data you can do anything you want to it. In the above example, I simply used a builtin function to calculate the mean... that doesn't mean that I am restricted to using builtin functions. The possibilities are endless since calculations can be 'vectorized' in any dimension... in this case, by row, by column or by depth/band/month/whatever is the 3rd dimension.

PS handling nodata values is trivial when working with arrays... in fact, many are builtin... instead of np.mean... there is amazingly ... np.nanmean()

Re: Trend Analysis Through Time Series of Raster Data

ShouvikJha — Mon, 31 Jul 2017 05:06:40 GMT

Thanks, Dan_Patterson‌, I am trying to further explore the script.

Re: Trend Analysis Through Time Series of Raster Data

XanderBakker — Fri, 10 Dec 2021 22:33:07 GMT

See below an example which I created from the slope of the regression line based on the 13 values. This means that if the values increment in time they will have a positive slope (red) and in case of a decreasing regression line a negative slope (green):

The code used to create this raster is:

#-------------------------------------------------------------------------------
# Name:        stats_numpy_series.py
# Purpose:
#
# Author:      Xander
#
# Created:     30-07-2017
#-------------------------------------------------------------------------------

def main():
    print "load modules..."
    import arcpy
    import numpy as np
    # template = r'C:\GeoNet\LoopMultiFolderCalcMeanSameFolder\Test_Data\2004A\{0}162004.tif'
    template = r'C:\GeoNet\WUE\r{0}_WUE.TIF'
    nodata = -3.4028235e+38
    out_ras = r'C:\GeoNet\WUE\WUE_slope.tif'

    print "create numpy list..."
    lst_np_ras = []
    for i in range(1, 14):
        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 "create output numpy array..."
    ras_path = template.format("%03d" % (1,))
    raster = arcpy.Raster(ras_path)
    ras_np_res = np.ndarray((rows, cols)) # np.ndarray((raster.height, raster.width))

    print "loop through pixels..."
    pix_cnt = 0
    for row in range(rows):
        for col in range(cols):
            pix_cnt += 1
            if pix_cnt % 10000 == 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)
            if cnt_elem > 1:
                val = a
##                print " - ", lst_vals
##                print " - slope     (a):", a
##                print " - intercept (b):", b
##                print " - total distance from regression line:", sum_abs_dif
##                print " - mean distance from regression line :", avg_abs_dif
            else:
                val = nodata
            try:
                ras_np_res[row, col] = val
            except:
                print "err"

    pnt = arcpy.Point(raster.extent.XMin, raster.extent.YMin)
    xcellsize = raster.meanCellWidth
    ycellsize = raster.meanCellHeight
    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]]")


def ReplaceNoData(lst, nodata):
    res = []
    for a in lst:
        if a < nodata / 2.0:
            res.append(None)
        else:
            res.append(a)
    return res


def RegressionLine(lst):
    from numpy import arange, array, ones, linalg
    months = range(1, 14)
    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]

        # 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()
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

I have attached the raster for you to examine.

Kind regards, Xander

Re: Trend Analysis Through Time Series of Raster Data

XanderBakker — Mon, 31 Jul 2017 12:02:18 GMT

Sorry Dan Patterson for not implementing a three dimensional ndarray as you suggested. I was not sure if this or my method could run into memory problems, but if the rasters are too big for this, there is a nice example of using blocks in the help page: RasterToNumPyArray—Help | ArcGIS Desktop (scroll down to example 2).

Re: Trend Analysis Through Time Series of Raster Data

DanPatterson_Retired — Mon, 31 Jul 2017 17:21:23 GMT

xander_bakker‌ blocks are good... what I would call strides with or without overlap... they can also reduce processing times if one is memory starved as well.

Re: Trend Analysis Through Time Series of Raster Data

DanPatterson_Retired — Fri, 10 Dec 2021 22:33:10 GMT

Well I wouldn't even worry about memory blocks ... 300, 000, 000 values

# ---- 5000 x 5000 cell raster ----

a = np.random.randint(0, 10, size=(12, 5000, 5000))

np.mean(a, axis=0)  # ---- let's try 5k x 5k raster
Out[5]: 
array([[ 3.17,  5.50,  4.75, ...,  4.00,  4.33,  5.25],
       [ 5.25,  6.50,  6.00, ...,  4.42,  5.00,  3.42],
       [ 3.25,  5.33,  4.58, ...,  3.67,  4.83,  4.58],
       ..., 
       [ 4.17,  2.83,  3.67, ...,  2.58,  4.08,  5.42],
       [ 4.83,  4.17,  4.83, ...,  6.58,  3.50,  3.67],
       [ 5.42,  4.33,  5.75, ...,  4.67,  3.00,  2.92]])
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

on the 8 Gb Surface Book...

even double (float64) 10,000 x 10,000 x 12 months caused a minor murmur

a = np.random.randint(0, 10, size=(12, 10000, 10000)) * 1.0

np.mean(a, axis=0)  # ---- let's try 10k x 10k raster
Out[12]: 
array([[ 5.00,  4.25,  5.17, ...,  3.67,  2.08,  5.08],
       [ 4.67,  5.50,  3.08, ...,  4.17,  3.50,  4.75],
       [ 3.83,  4.17,  4.00, ...,  3.50,  4.58,  4.42],
       ..., 
       [ 5.75,  5.00,  4.42, ...,  3.92,  5.25,  4.50],
       [ 4.17,  3.83,  4.50, ...,  5.17,  4.67,  4.50],
       [ 4.75,  4.00,  3.08, ...,  4.17,  5.33,  5.83]])

a.size
Out[13]: 1200000000  # 1,200,000,000 1.2 Billion values‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Trend Analysis Through Time Series of Raster Data

ShouvikJha — Mon, 31 Jul 2017 18:57:05 GMT

Hi Xander Bakker Really a speechless gratitude for such wonderful work. And it is appreciated by our team.

I ran the script, after that, i cross check the output pixel values with manual calculation in excel corresponding to the same pixel. The output values are not matching with manual pixel wise calculation. It is slightly different.

Our kind request please have a look on the matter. I have attached my Excel calculation sheet, where both values are assigned, from manual and script output value.

Once again thanks for your such cooperation.

Thanks & regards

Shouvik Jha

Re: Trend Analysis Through Time Series of Raster Data

XanderBakker — Mon, 31 Jul 2017 19:20:19 GMT

You are right. Good catch! It seems there is something going wrong in the calculation. I just tried a manual list from pixel 1 into the function and a slope of -1.8379670329670326 (same as you Excel result) was returned. I will need to investigate what is going wrong here.

Re: Trend Analysis Through Time Series of Raster Data

XanderBakker — Tue, 01 Aug 2017 12:36:22 GMT

Hi shouvik jha ,

Apparently the number of lines of the output raster was 1 more than the input rasters and that caused a shift of 1 pixel. If have attached a corrected version of the output raster to the thread. I will look at the code to see what is wrong and attach the corrected code at a later moment. Could you revise the new raster?

Kind regards, Xander

Re: Trend Analysis Through Time Series of Raster Data

ShouvikJha — Tue, 01 Aug 2017 13:46:49 GMT

Hi Xander Bakker‌, Thank you very much. really a wonderful work has been done by you.

I cross checked the output and it matched exactly with pixel wise manual calculation.

Re: Trend Analysis Through Time Series of Raster Data

XanderBakker — Tue, 01 Aug 2017 14:18:57 GMT

That is great news. I will get back with the corrected code as soon as I detect the error.

Re: Trend Analysis Through Time Series of Raster Data

XanderBakker — Fri, 10 Dec 2021 22:33:13 GMT

Hi shouvik jha ,

To be able to debug a little faster I extracted a small part of the larger rasters and ran the script again. This time using 10.5 and not 10.3 which I have at home and there was no pixel shift. Or something changed between versions or something very weird is happening.

See the code below which in this case includes additional lines to generate 3 other rasters too (intercept, total distance to regression line and mean distance to regression line)

#-------------------------------------------------------------------------------
# Name:        stats_numpy_series.py
# Purpose:
#
# Author:      Xander
#
# Created:     30-07-2017
#-------------------------------------------------------------------------------

def main():
    print "load modules..."
    import arcpy
    arcpy.env.overwriteOutput = True

    import numpy as np
    template = r'C:\GeoNet\WUE2\r{0}_WUE.TIF'
    nodata = -3.4028235e+38
    out_template = r'C:\GeoNet\WUE2\WUE2_{0}.TIF'
    out_names = ['slope', 'intersept', 'tot_dist', 'mean_dist']

    print "create numpy list..."
    lst_np_ras = []
    for i in range(1, 14):
        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" % (1,))
    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)
            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 RegressionLine(lst):
    from numpy import arange, array, ones, linalg
    months = range(1, 14)
    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]

        # 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()
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

See attached ZIP with the data used, script and resulting rasters.

Re: Trend Analysis Through Time Series of Raster Data

ShouvikJha — Wed, 02 Aug 2017 12:59:35 GMT

Hi Xander Bakker‌ Thank you very much. Really this code helps us in many ways to do the statistical calculation.

When I use the whole area that time also produced the same error (Increasing the row number) as earlier error, then i check it using another raster and I got the exact right result. Then I realized the Problem was not with the code, the problem was with data. After that, i just use the same data (WUE) but this time i used some tricky, First, i fill no data value as a Zero then extract the area of interest using my shape file, then I use the raster as an input of trend line analysis and got the exact result.

What I realized, If the raster data contain any no data value that time the problem may occur.

And you and Dan_Patterson‌ grateful thanks to you people behalf of our team for your kind cooperation.