I am trying to use the Theil-Sen analysis on the Landsat images to compute vegetation indices, help me figure out how to come up with the median from the Theil-Sen equation in ArcGIS, ERDAS IMAGINE or Excel. The equation is as follows: where:
This notation is terrible in communicating application to a time-series. Keep in mind that the subscripts indicate periods of time so, it is the median of all pairwise comparisons. Because the median slope is derived from all pairwise comparisons in the time-series
n(n-1)/2, without scripting in Python, this will be very difficult to implement in ArcGIS. With just 8 time-periods you would have to take the median across 28 pairwise rasters.
It is also important to note that the Theil-Sen median slope can absorb ~30% noise without undue influence on the slope estimate. Over long time-series data, it is important to smooth the data to reduce outlier effects on noise. Additionally, without a test, such as the Mann-Kendall, you have no way to evaluate the significance of any monotonic trend through the time-series. It is uncommon that the Theil-Sen is applied with no additional test statistics.
Honestly, you would be better off using Kendall's Tau and its derived Theil-Sen slope. Application of this statistic also returns the p-values and +/- confidence intervals. I just added a function
raster.kendall in the R package spatialEco, that implements this statistic on raster time series data.
You compute the formula and as a result a new raster is generated. In ArcGIS you just go to Raster statistics. In ERDAS it's pretty much the same (Right click on the image and see the statistics of the raster). In Excel you go to an empty cell and use the MEDIAN function having all of the cells selected.