Cross-year date interval selection and trend computation (pre).

NDVITrend.py

@@ -2,6 +2,7 @@ import glob
 import datetime
 import os
 import subprocess
+import warnings
 
 def temporalSmoothing(fld,
                       template='MOD13Q1.006__250m_16_days_NDVI_',
@@ -43,4 +44,70 @@ def temporalSmoothing(fld,
                 expr.append('im1b'+str(n+1))
                 cdf.write(dts[n])
             cmd += ['{' + ';'.join(expr) + '}']
-        subprocess.call(cmd)
\ No newline at end of file
+        subprocess.call(cmd)
+
+def getCoherentReferencePeriod(dates_fn,start_date=(1,1),end_date=(12,31),dates_per_year=23):
+    with open(dates_fn,'rb') as df:
+        lst = df.readlines()
+    date_list = [datetime.datetime.strptime(x.strip(),'%Y%m%d') for x in lst]
+    # Get closest date to start_date
+    gaps = []
+    for dt in date_list:
+        d = 366
+        for y in range(dt.year-1,dt.year+2):
+            tmp = (dt-datetime.datetime(y,start_date[0],start_date[1])).days
+            # Allow only existing starting dates after start_date
+            if tmp > 0 and tmp < d:
+                d = tmp
+            """
+            # Allow search for closest date before start_date
+            if abs(tmp) < d:
+                d = abs(tmp)
+            """
+        gaps.append(d)
+
+    start_index = gaps.index(min(gaps))
+
+    # Find the next date closest to (and lower or equal to) end_date
+    i = start_index + 1
+    ey = date_list[i].year
+    if (date_list[i]-datetime.datetime(ey,end_date[0],end_date[1])).days > 0:
+        ey += 1
+    ed = datetime.datetime(ey,end_date[0],end_date[1])
+    while date_list[i] <= ed:
+        i += 1
+    # Always get (month,day) below the limit to avoid overlaps on yearly series
+    end_index = i-1
+    """
+    # Allow closest date search beyond the limit (dangerous on yearly series)
+    end_index = i
+    if abs((date_list[end_index - 1] - ed).days) < abs((date_list[end_index] - ed).days):
+        end_index -= 1
+    """
+
+    # Get relative indices
+    idx_delay = end_index-start_index
+    start_index %= dates_per_year
+    end_index = start_index + idx_delay
+
+    # Overlap windows and get largest time interval
+    min_start = (date_list[start_index].month,date_list[start_index].day)
+    max_end = (date_list[end_index].month, date_list[end_index].day)
+    for i in range(start_index+dates_per_year,len(date_list),dates_per_year):
+        if (date_list[i].month < min_start[0]) or (date_list[i].month == min_start[0] and date_list[i].day < min_start[1]):
+            min_start = (date_list[i].month, date_list[i].day)
+        if (i+idx_delay) < len(lst):
+            if (date_list[i+idx_delay].month > max_end[0]) or (date_list[i+idx_delay].month == max_end[0] and date_list[i+idx_delay].day > max_end[1]):
+                max_end = (date_list[i+idx_delay].month, date_list[i+idx_delay].day)
+        else:
+            warnings.warn('Time series does not cover last period entirely!')
+            break
+
+    return min_start,max_end,idx_delay+1
+
+def computeNDVITrends(series_fn,dates_fn,start_date=(1,1),end_date=(12,31),dates_per_year=23):
+    cmd = ['otbcli_TimeSeriesIndexTrend', '-ndvits', series_fn, '-ndvidates', dates_fn, '-ndvi.reduce', 'cumul']
+    sd,ed = getCoherentReferencePeriod(dates_fn,start_date,end_date,dates_per_year)
+    cmd += ['-ndvi.reduce.cumul.month1', str(sd[0]),'-ndvi.reduce.cumul.day1', str(sd[1])]
+    cmd += ['-ndvi.reduce.cumul.month2', str(ed[0]), '-ndvi.reduce.cumul.day2', str(ed[1])]
+    subprocess.call(cmd)