# -*- coding: utf-8 -*-

#==========================================================#
# Draw reference points on satellite image
#==========================================================#

# Preamble
import os
import ee
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import pandas as pd
from datetime import datetime
import pickle
import pdb
import pytz
from pylab import ginput

# image processing modules
import skimage.filters as filters 
import skimage.exposure as exposure
import skimage.transform as transform
import sklearn.decomposition as decomposition
import skimage.morphology as morphology
import skimage.measure as measure

# my functions
import functions.utils as utils
import functions.sds as sds

np.seterr(all='ignore') # raise/ignore divisions by 0 and nans
ee.Initialize()

# collection
input_col = ee.ImageCollection('LANDSAT/LC08/C01/T1_RT_TOA')

# location (Narrabeen-Collaroy beach)
#polygon = [[[151.301454, -33.700754],
#               [151.311453, -33.702075], 
#               [151.307237, -33.739761],
#               [151.294220, -33.736329],
#               [151.301454, -33.700754]]];
# location (Oldbar shoreline)
#polygon = [[[152.664508, -31.896163],
#               [152.665827, -31.897112], 
#               [152.631516, -31.924846],
#               [152.629285, -31.923362],
#               [152.664508, -31.896163]]];
# location (Oldbar inlet)
polygon = [[[152.676283, -31.866784],
   [152.709174, -31.869993], 
   [152.678229, -31.892082],
   [152.670366, -31.886360],
   [152.676283, -31.866784]]];

# dates
start_date = '2017-01-30'
end_date = '2017-02-02'
#start_date = '2017-01-30'
#end_date = '2018-02-02'
# filter by location
flt_col = input_col.filterBounds(ee.Geometry.Polygon(polygon)).filterDate(start_date, end_date)

n_img = flt_col.size().getInfo()
print('Number of images covering the area:', n_img)
im_all = flt_col.getInfo().get('features')

satname = 'L8'
#sitename = 'NARRA'
sitename = 'OLDBAR_inlet'

# parameters
plot_bool = False # if you want the plots
prob_high = 99.9 # upper probability to clip and rescale pixel intensity
min_contour_points = 100 # minimum number of points contained in each water line
output_epsg = 28356 # GDA94 / MGA Zone 56
cloud_threshold = 0.8

# find image in ee database
im = ee.Image(im_all[0].get('id')) 
# load image as np.array
im_pan, im_ms, im_cloud, crs, meta = sds.read_eeimage(im, polygon, satname, plot_bool)

# rescale intensities
im_ms = sds.rescale_image_intensity(im_ms, im_cloud, prob_high, plot_bool)
im_pan = sds.rescale_image_intensity(im_pan, im_cloud, prob_high, plot_bool)

# pansharpen rgb image
im_ms_ps = sds.pansharpen(im_ms[:,:,[0,1,2]], im_pan, im_cloud, plot_bool)

plt.figure()
plt.imshow(im_ms_ps[:,:,[2,1,0]])
plt.show()

pts = ginput(n=50, timeout=1000, show_clicks=True)
points = np.array(pts)
plt.plot(points[:,0], points[:,1], 'ko')
plt.show()

pts_coords = sds.convert_pix2world(points[:,[1,0]], crs['crs_15m'])
pts = sds.convert_epsg(pts_coords, crs['epsg_code'], output_epsg)

filepath = os.path.join(os.getcwd(), 'data', satname, sitename)
with open(os.path.join(filepath, sitename + '_refpoints2.pkl'), 'wb') as f:
    pickle.dump(pts, f)