import glob import datetime import os import subprocess import warnings 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): lst = sorted(glob.glob(fld + '/' + template + '*.' + ext)) cmd = ['otbcli_TemporalSmoothing','-il'] with open(fld + '/dates.txt','wb') as df: for f in lst: dt = os.path.basename(f)[date_pos:date_pos+date_len] if date_jul: y,jd = int(dt[0:4]),int(dt[4:7]) dt = (datetime.datetime(y,1,1) + datetime.timedelta(jd)).strftime('%Y%m%d') df.write(dt+'\n') cmd += [f] cmd += ['-dates',fld+'/dates.txt','-out',fld+'/SmoothedSeries.tif','uint16','-interp','sg','-interp.sg.deg',str(sgdeg),'-interp.sg.rad',str(sgrad)] subprocess.call(cmd) if clip_valid > 0: cmd = ['otbcli_BandMathX', '-il', fld+'/SmoothedSeries.tif', '-out', fld+'ClippedSmoothedSeries.tif','uint16','-exp'] 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 = [] for n in range(clip_valid,clip_end): expr.append('im1b'+str(n+1)) cdf.write(dts[n]) cmd += ['{' + ';'.join(expr) + '}'] 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,N = 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)