#!/usr/bin/env python
# coding: utf-8

# # Sandbox Metadata Document
# 
# This document aims to provide an overview of the different data sets available in the Open Data Cube Sandbox. It will provide details on which products are available, what regions are indexed for each product, and the time series available within that extent.
# 
# For information on how to get started with the Datacube, please see the read_me_first.ipynb available [here](https://sandbox.test.frontiersi.io/user/sandbox/notebooks/Read_me_first.ipynb).
# 
# * [Datacube Information](#datacube)  
#     * [ALOS 2](#alos2)  
#     * [Landsat 5](#ls5)  
#     * [Landsat 7](#ls7)  
#     * [Landsat 8](#ls8)  
#     * [Sentinel 1](#s1)
# * [Product coverage illustration](#map)
#     
# * [System Information](#system)

# <a class="anchor" id="datacube"></a>
# ## Datacube Information

# In[1]:


from IPython.core.display import display, HTML
import pandas as pd

import datacube
dc = datacube.Datacube()
dc.list_products()


# <a class="anchor" id="alos2"></a>
# ### alos2_palsar_AMA_ingest
# 
# **Caqueta**  
# ``Time: 2015-01-01 and 2016-01-01 ``   
# ``latitude:  (  0,            3.2)``    
# ``longitude: (-76.26666668, -73.6)``      
# 
# **Vietnam**  
# ``Time: 2015-01-01 and 2016-01-01 ``   
# ``latitude:  (  7.46666667,  12.8)``   
# ``longitude: (103.46666668, 110.4)``     
# 
# *****
# <a class="anchor" id="ls5"></a>
# ### ls5_level1_usgs
# *No Data Indexed*
# *****
# ### ls5_usgs_sr_scene
# 
# ``Time: 1984-05-26 to 2011-11-13``  
# ``y: (-587700.0, 2187300.0)``  
# ``x: (  68100.0,  932400.0)``  
# 
# *****
# <a class="anchor" id="ls7"></a>
# ### ls7_collection1_AMA_ingest
# 
# **Caqueta**  
# ``Time: 1999-08-21 to 2018-03-25``  
# ``latitude:   (  0.000134747292617865, 1.077843593651382)``  
# ``longitude:  (-74.91935994831539,   -73.30266193148462)``  
# ``resolution: ( -0.000269494585236,    0.000269494585236)``  
# 
# **Lake Baringo, Kenya**  
# ``Time: 2005-01-08 to 2016-12-24``  
# ``latitude:   ( 0.4997747685,  0.7495947795)``  
# ``longitude:  (35.9742163305, 36.473586859499996)``  
# ``resolution: (-0.000269493,   0.000269493)``  
# 
# *****
# ### ls7_level1_usgs
# *No Data Indexed*
# 
# *****
# ### ls7_usgs_sr_scene
# 
# ``Time: 1999-07-08 to 2018-11-02``  
# ``y: (-586200.0, 2185200.0)``  
# ``x: (  70800.0,  935100.0)``  
# 
# *****
# <a class="anchor" id="ls8"></a>
# ### ls8_collection1_AMA_ingest
# 
# **Caqueta**  
# ``Time: 2013-04-13 to 2018-03-26``  
# ``latitude:   (  0.000134747292617865, 1.077843593651382)``  
# ``longitude:  (-74.91935994831539,   -73.30266193148462) ``  
# ``resolution: ( -0.000269494585236,    0.000269494585236)``  
# 
# **Vietnam**  
# ``Time: 2014-01-14 to 2016-12-21``  
# ``latitude:   ( 10.513927001104687, 12.611133863411238)``    
# ``longitude:  (106.79005909290998, 108.91906631627438) ``  
# ``resolution: ( -0.000269494585236,  0.000269494585236)``  
# 
# *****
# ### ls8_l1_pc_usgs
# *No Data Indexed*
# 
# *****
# ### ls8_level1_usgs 
# [Global](https://sandbox.test.frontiersi.io/user/Sandbox/view/examples/images/indexed_areas.png)
# 
# 
# *****
# ### ls8_usgs_sr_scene
# ``Time: 2013-03-21 to 2018-11-16``  
# ``y: (-596700.0, 2193300.0)``  
# ``x: (  82500.0,  930300.0)``  
# 
# *****
# <a class="anchor" id="s1"></a>
# ### s1_gamma0_scene
# 
# ``Time: 2015-05-12 to 2018-07-20``
# 
# **Caqueta**   
# ``latitude:  (  1.00018083,  2.)``  
# ``longitude: (-75.,        -74.)``  
# 
# **Samoa**  
# ``latitude:  ( -13.99981917, -13.)``  
# ``longitude: (-171.,        -173.)``    
# 
# 

# <a class="anchor" id="map"></a>
# ### Plotting Indicative Data Extents
# 
# The following code generates an overview of the approximte extents of the product selected on line 1, shown in red. The yellow base image shows the extents of the global ls8_level1_usgs product.

# In[2]:


product = 'ls8_collection1_AMA_ingest'

import matplotlib.pyplot as plt
from pyproj import Proj, transform

import datacube
dc = datacube.Datacube()

outProj = Proj(init='epsg:4326')

x_bounds = []
y_bounds = []
unique_bounds=set()
bounds=[]
times=[]
tiles = dc.find_datasets(product=product)
if max(tiles[0].bounds) > 180: XY = True
else: XY = False

for tile in tiles:
    if XY:
        inProj = Proj(init=tile.crs)
        l_lon,t_lat = transform(inProj,outProj,tile.bounds.left,tile.bounds.top)
        r_lon,b_lat = transform(inProj,outProj,tile.bounds.right,tile.bounds.bottom)
        x_bounds = [l_lon, r_lon, r_lon, l_lon]
        y_bounds = [t_lat, t_lat, b_lat, b_lat]
    
    else:
        x_bounds = [tile.bounds.left, tile.bounds.right, tile.bounds.right, tile.bounds.left]
        y_bounds = [tile.bounds.top, tile.bounds.top, tile.bounds.bottom, tile.bounds.bottom]
        
    times.append(tile.time)
    bounds_key = "".join(str(x) for x in x_bounds) + "".join(str(y) for y in y_bounds)
    if bounds_key not in unique_bounds:
        bounds.append([x_bounds, y_bounds, tile.time])

get_ipython().run_line_magic('matplotlib', 'inline')
fig, ax = plt.subplots(1, figsize = (18,9))
for i in bounds:
    ax.fill(i[0], i[1], c='red')
ax.axis('equal')

plt.title(product)
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
basemap = plt.imread('../images/indexed_areas.png')
ax.imshow(basemap, extent=[-180,180,-90,90])
ax.set_xlim(-180,180)
ax.set_ylim(-90, 90)
plt.show()
print('Time:', min(times[-1]).strftime('%Y-%m-%d'), 'to', max(times[0]).strftime('%Y-%m-%d'))


# <a class="anchor" id="system"></a>
# # System Information

# In[3]:


import sys
print('Python Version', sys.version)


# In[4]:


get_ipython().system('pip3 show datacube xarray numpy pandas gdal matplotlib')


# ## Extent Query Functions

# In[5]:


#Generates text extents of data
product = 's1_gamma0_scene'

import datacube
dc = datacube.Datacube()
import numpy as np

left = np.array([])
bottom = np.array([])
right = np.array([])
top = np.array([])
times = np.array([])

tiles = dc.find_datasets(product=product)
for tile in tiles: #Each tile (unique time for each)
    left = np.append(left, tile.bounds[0])
    right = np.append(right, tile.bounds[2])
    top = np.append(top, tile.bounds[3])
    bottom = np.append(bottom, tile.bounds[1])
    times = np.append(times, tile.time)
    #corners = np.append(corners, [])
print('Max Extent')
print('``Time:', min(times).strftime('%Y-%m-%d'), 'to',  max(times).strftime('%Y-%m-%d')+'``  ')
print('``latitude:  ('+str(min(np.unique(bottom)))+', '+str(max(np.unique(top)))+')``  ')
print('``longitude: ('+str(min(np.unique(left)))+', '+str(max(np.unique(right)))+')``  ')

print('Discrete Extents')
print('``Time:', min(times).strftime('%Y-%m-%d'), 'to',  max(times).strftime('%Y-%m-%d')+'``  ')
print('``latitude:  ('+str(np.unique(bottom))+', '+str(np.unique(top))+')``  ')
print('``longitude: ('+str(np.unique(left))+', '+str(np.unique(right))+')``  ')


# In[6]:


pd.options.display.max_rows = 200
dc.list_measurements()