NDVITrend.py 4.56 KB
Newer Older
1
2
import glob
import datetime
3
import os
4
import subprocess
5
import warnings
6

7
8
9
10
11
12
13
14
15
16
def temporalSmoothing(fld,
                      template='MOD13Q1.006__250m_16_days_NDVI_',
                      date_pos=34,
                      date_len=7,
                      date_jul=True,
                      ext='tif',
                      sgdeg=4,
                      sgrad=6,
                      clip_valid=23,
                      clip_both_sides=True):
17
18
19
20
21
22
23

    lst = sorted(glob.glob(fld + '/' + template + '*.' + ext))

    cmd = ['otbcli_TemporalSmoothing','-il']

    with open(fld + '/dates.txt','wb') as df:
        for f in lst:
24
            dt = os.path.basename(f)[date_pos:date_pos+date_len]
25
            if date_jul:
26
                y,jd = int(dt[0:4]),int(dt[4:7])
27
28
                dt = (datetime.datetime(y,1,1) + datetime.timedelta(jd)).strftime('%Y%m%d')
            df.write(dt+'\n')
29
            cmd += [f]
30

31
    cmd += ['-dates',fld+'/dates.txt','-out',fld+'/SmoothedSeries.tif','uint16','-interp','sg','-interp.sg.deg',str(sgdeg),'-interp.sg.rad',str(sgrad)]
32
    subprocess.call(cmd)
33
34

    if clip_valid > 0:
35
        cmd = ['otbcli_BandMathX', '-il', fld+'/SmoothedSeries.tif', '-out', fld+'ClippedSmoothedSeries.tif','uint16','-exp']
36
37
38
39
40
41
42
        clip_end = len(lst)
        if clip_both_sides:
            clip_end -= clip_valid
        with open(fld + '/dates.clipped.txt','wb') as cdf:
            with open(fld + '/dates.txt','rb') as df:
                dts = df.readlines()
            expr = []
43
            for n in range(clip_valid,clip_end):
44
45
46
                expr.append('im1b'+str(n+1))
                cdf.write(dts[n])
            cmd += ['{' + ';'.join(expr) + '}']
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
        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']
110
    sd,ed,N = getCoherentReferencePeriod(dates_fn,start_date,end_date,dates_per_year)
111
112
113
    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)