import datacube
import utils.data_cube_utilities.data_access_api as dc_api  

from datacube.utils.aws import configure_s3_access
configure_s3_access(requester_pays=True)

api = dc_api.DataAccessApi()
dc = datacube.Datacube(app = 'ardc_task_d')
api.dc = dc

# Enable plotting
%matplotlib inline

# Supress Warning 
import warnings
warnings.filterwarnings('ignore')

"""
Explanation: ARDC Training: Python Notebooks
Task-D: Land Change

Import the Datacube Configuration
End of explanation
"""


list_of_products = dc.list_products()
netCDF_products = list_of_products[list_of_products['format'] == 'NetCDF']
netCDF_products

"""
Explanation: Browse the available Data Cubes
End of explanation
"""


# Change the data platform and data cube here

platform = 'LANDSAT_7'
product = 'ls7_usgs_sr_scene'
collection = 'c1'
level = 'l2'

# Get Coordinates
coordinates = api.get_full_dataset_extent(platform = platform, product = product)

"""
Explanation: Pick a product
Use the platform and product names from the previous block to select a Data Cube.
End of explanation
"""


from utils.data_cube_utilities.dc_time import _n64_to_datetime, dt_to_str

extents = api.get_full_dataset_extent(platform = platform, product = product, measurements=[])

latitude_extents = (min(extents['latitude'].values),max(extents['latitude'].values))
longitude_extents = (min(extents['longitude'].values),max(extents['longitude'].values))
time_extents = (min(extents['time'].values),max(extents['time'].values))

print("Latitude Extents:", latitude_extents)
print("Longitude Extents:", longitude_extents)
print("Time Extents:", list(map(dt_to_str, map(_n64_to_datetime, time_extents))))

"""
Explanation: Display Latitude-Longitude and Time Bounds of the Data Cube
End of explanation
"""


## The code below renders a map that can be used to orient yourself with the region.
from utils.data_cube_utilities.dc_display_map import display_map
display_map(latitude = latitude_extents, longitude = longitude_extents)

"""
Explanation: Visualize Data Cube Region
End of explanation
"""


## Vietnam
# Lam Dong Province near reservior 
# lat =  (11.843, 11.922)
# lon =  (107.723, 107.821)
# Field near Gia Nghia
# lat =  (11.90, 11.92)
# lon =  (107.76, 107.78)

## Ghana
# Weija Reservoir
lat = (5.5388, 5.6126)
lon = (-0.4089, -0.3327)
time_extents = ('2001-01-01', '2005-12-31')

display_map(lat, lon)

"""
Explanation: Pick a smaller analysis region and display that region
Try to keep your region to less than 0.02-deg x 0.02-deg for rapid processing. This will give you a region of about 75x75 pixels. You can click on the map above to find the Lat-Lon coordinates of any location. You will want to identify a region with an urban aree or some known vegetation. Pick a time window of 10+ years, as the curvefit algorithm requires a long time series to detect change.
<br><br>
Here is what to expect ... 0.02-deg x 0.02-deg (75x75 pixels) over 10 years will take 10 to 15 minutes to execute. Be patient ...
End of explanation
"""


measurements = ['red', 'green', 'blue', 'nir', 'swir1', 'swir2', 'pixel_qa']
landsat_dataset = \
    dc.load(latitude = lat,
            longitude = lon,
            platform = platform,
             time = time_extents,
             product = product,
             measurements = measurements)

landsat_dataset
#view the dimensions and sample content from the cube

from utils.data_cube_utilities.clean_mask import landsat_clean_mask_full

# obtain the clean mask
clean_mask = landsat_clean_mask_full(dc, landsat_dataset, product=product, platform=platform,
                                     collection=collection, level=level)

landsat_dataset = landsat_dataset.where(clean_mask)

"""
Explanation: Load the dataset and the required spectral bands or other parameters
After loading, you will view the Xarray dataset. Notice the dimensions represent the number of pixels in your latitude and longitude dimension as well as the number of time slices (time) in your time series.
End of explanation
"""


# import utils.data_cube_utilities.dc_ccd as ccd

# %time change_matrix = ccd.process_xarray(landsat_dataset, distributed = True, process = "matrix") #Run process xarray on large dataset

# %time change_volume = (change_matrix.sum(dim='time') - 1).rename('change_volume')

"""
Explanation: PyCCD
The compute time of PyCCD scales sharply with extent sizes. Seemingly small areas can take anywhere from 10-15 minutes to process.
End of explanation
"""


# import matplotlib.pyplot as plt  
# from utils.data_cube_utilities.plotter_utils import figure_ratio

# plt.figure(figsize = figure_ratio(change_volume, fixed_width=8))
# change_volume.plot()
# plt.axes().set_aspect("equal")

"""
Explanation: Plotting change volume
Plot change volume to identify regions/areas of change. The change volume represents the number of land changes for a pixel over the time series.
End of explanation
"""


# %time
# time_map_ccd_product = ccd._nth_occurence_in_ccd_matrix(change_matrix, 1, f = ccd._n64_datetime_to_scalar)

# import datetime
# from matplotlib.ticker import FuncFormatter

# plt.figure(figsize = figure_ratio(time_map_ccd_product, fixed_width=8))
# epochFormatter = FuncFormatter(lambda x, pos: datetime.datetime.utcfromtimestamp(x).strftime('%Y-%m-%d'))
# time_map_ccd_product.plot(cmap = "magma", cbar_kwargs=({'format': epochFormatter}))
# plt.axes().set_aspect("equal")

"""
Explanation: Plot the time of first changes
End of explanation
"""


from utils.data_cube_utilities.dc_rgb import rgb

first = 2 # Acquisition Number ... change here

print( landsat_dataset.time.values[first] )

rgb( landsat_dataset,
     time_index = first,
     bands = ['swir2','nir','green'])

last = 10 # Acquisition Number ... change here

print( landsat_dataset.time.values[last] )

rgb(landsat_dataset,
    time_index = last,
    bands = ['swir2','nir','green'])

"""
Explanation: Validating Change
Use the two images below to review scenes at the beginning of the time series and the end of the time series. 
Review of RGB images
Choose an image from the start and end of the time series. Remember that 0 is the first acquisition and the last acquisition if the number of time steps. You can find that number in the XARRAY report above. Try various combinations of the start and end images to identify the land changes using visual interpretation. You will see some images will have clouds and others will have Landsat-7 "bands". You can also change the RGB image bands to give you an combination you desire.
End of explanation
"""


import utils.data_cube_utilities.trend as trend

import numpy as np

def land_and_water_masking_ls7(dataset):    
    #Create boolean Masks for clear and water pixels
    clear_pixels = dataset.pixel_qa.values == 2 + 64
    water_pixels = dataset.pixel_qa.values == 4 + 64

    a_clean_mask = np.logical_or(clear_pixels, water_pixels)
    return a_clean_mask

def remove_clouds(dataset):  
    return dataset.where(land_and_water_masking_ls7(dataset)).drop('pixel_qa')

def NDVI(dataset):
    return (dataset.nir - dataset.red)/(dataset.nir + dataset.red)

ndvi_trend_product = trend.linear(NDVI(remove_clouds(landsat_dataset)))

(-ndvi_trend_product).plot()

"""
Explanation: Vogelmann NDVI Trend
End of explanation
"""


