Bfast interface file

.py

import numpy
from PIL import Image
import os, fnmatch, time, copy
from datetime import datetime
from bfast import BFASTMonitor
from bfast.utils import crop_data_dates
import matplotlib.pyplot as plt
import matplotlib

# path = '../Images/Tile_1'
path = '../Images/subset_1'

files = []
# data = numpy.empty((124,1584,1354),"int16") # for Tile_1
data = numpy.empty((124,17,17),"int16") # for subset_1
i = 0 # iterator

#% Open Tif files data
files = fnmatch.filter(os.listdir(path), '*.tif')
for file in files:
    if '.tif' in file:
        print(file)
        im = Image.open(path + "//" + file)
        marray = numpy.array(im,"int16")
        data[i,:,:] = marray
        i +=1

#% Get dates
dates = []
for file in files:
    dates.append( datetime.strptime(file, 'X%Y%m%d.tif'))

#% parameters
k = 3
freq = 365
trend = False
hfrac = 0.25
level = 0.05
start_hist = dates[0]
start_monitor = datetime(2018, 6, 18)
end_monitor = dates[-1]

# # fit BFASTMontiro model
model = BFASTMonitor(
            start_monitor,
            freq=freq,
            k=k,
            hfrac=hfrac,
            trend=trend,
            level=level,
            backend='python',
            verbose=1,
            detailed_results=False,
            # device_id=0,
            )

# only apply on a small subset
# data = data[:,120:422,:] # for test Tile_1 error

# data, dates = crop_data_dates(data, dates, start_hist, end_monitor)
print("First date: {}".format(dates[0]))
print("Last date: {}".format(dates[-1]))
print("Start monitor date: {}".format(start_monitor))
print("Shape of data array: {}".format(data.shape))

start_time = time.time()
model.fit(data, dates, n_chunks=None, nan_value=0)
end_time = time.time()
print("All computations have taken {} seconds.".format(end_time - start_time))

# visualize results
breaks = model.breaks
means = model.means

no_breaks_indices = (breaks == -1)
means[no_breaks_indices] = 0
means[means > 0] = 0
breaks_plot = breaks.astype(numpy.float)
breaks_plot[breaks == -2] = numpy.nan
breaks_plot[breaks == -1] = numpy.nan
breaks_plot[means >= 0] = numpy.nan

dates_monitor = []

# collect dates for monitor period
for i in range(len(dates)):
    if start_monitor <= dates[i]:
        dates_monitor.append(dates[i])
 
dates_array = numpy.array(dates_monitor)
idxs = []
idxs.append(numpy.argmax((dates_array >= datetime(2013, 1, 1)) > False))
idxs.append(numpy.argmax((dates_array >= datetime(2014, 1, 1)) > False))
idxs.append(numpy.argmax((dates_array >= datetime(2015, 1, 1)) > False))
idxs.append(numpy.argmax((dates_array >= datetime(2016, 1, 1)) > False))
idxs.append(numpy.argmax((dates_array >= datetime(2017, 1, 1)) > False))
idxs.append(numpy.argmax((dates_array >= datetime(2018, 1, 1)) > False))
idxs.append(numpy.argmax((dates_array >= datetime(2019, 1, 1)) > False))
idxs.append(numpy.argmax((dates_array >= datetime(2020, 1, 1)) > False))
idxs.append(numpy.argmax((dates_array >= datetime(2021, 1, 1)) > False))

i = 1
ticks = []
ticklabels = []

breaks_plot_years = copy.deepcopy(breaks_plot)
breaks_plot_years[numpy.where(breaks_plot <= idxs[start_monitor.year-2013])] = i
ticks.append(i)
ticklabels.append(str(start_monitor.year-1))

for x in range(start_monitor.year, end_monitor.year):
    i += 1
    breaks_plot_years[numpy.where(numpy.logical_and(idxs[x-2013] < breaks_plot, breaks_plot <= idxs[x-2013+1]))] = i
    ticks.append(i)
    ticklabels.append(str(x))
    
i += 1
breaks_plot_years[numpy.where(idxs[end_monitor.year-2013] < breaks_plot)] = i
ticks.append(i)
ticklabels.append(str(end_monitor.year))
ticks.append(i+1)
ticklabels.append(str(end_monitor.year+1))

print(ticklabels)

bounds = numpy.linspace(1, i+1, i+1)
cmap = plt.get_cmap("gist_rainbow")
cmaplist = [cmap(i) for i in range(cmap.N)]
cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)
norm = matplotlib.colors.BoundaryNorm(bounds, cmap.N)

fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(20, 20))
im = axes.imshow(breaks_plot_years, cmap=cmap, vmin=0, vmax=7, norm=norm)
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
fig.colorbar(im, cax=cbar_ax, ticks=ticks)
labels = cbar_ax.set_yticklabels(ticklabels)

plt.show()