import skimage.measure as skmeasure
import skimage.io as skio
import scipy.stats as stats
import scipy.ndimage as ndimage
import matplotlib.pyplot as plt
import os


def threshold(im, threshmin=1, threshmax=20, display=False):

    image_border = (im == 0) / 2
    water = stats.threshold(im, threshmin, threshmax)

    #select only the largest water area, removes all noise and inland lakes and rivers
    labeled_areas, _ = ndimage.label(water)
    mainregion = max(skmeasure.regionprops(labeled_areas), key=lambda r: r.filled_area)

    water = labeled_areas == mainregion.label
    water = water + image_border

    if display:
        fig, (ax1, ax2) = plt.subplots(1, 2)

        ax1.imshow(im, cmap=plt.cm.gray)
        ax1.axis('off')
        ax1.set_title('raw')

        ax2.imshow(water, cmap=plt.cm.gray)
        ax2.axis('off')
        ax2.set_title('cleaned')

        plt.show()

    return water



if __name__ == "__main__":
    #Takes a landsat band 5 image
    filename = os.path.join('testdata', 'LT50400372011313PAC01_B5.TIF')
    im = skio.imread(filename)
    data = threshold(im, display=True)