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

# # OUTPUT TRACKING ALGORITHM
# 
# 
# ---
# Author: **Helvecio B. Leal Neto** & **Alan J. P. Calheiros**\
# **National Institute for Space Research - Brazil - (2021)**
# 
# 
# 
# ## About
# 
# This notebook is designed for viewing the tracking results of the storm/precipitation tracking algorithm beta version. The results presented here refer to the tracking of clusters via radar data provided by the GoAmazon project, for the following periods:
# 
# **Start**: 2014-09-07 00:00:00
# 
# **End**: 2014-09-09 00:00:00
# 
# The tracking threshold is:
# 
# * **20** dBZ
# * inner 1 - ***35*** dBZ
# * inner 2 - ***40*** dBZ
# 
# Minimum size threshold per cluster:
# 
# * **30** pixels
# * inner 1 - ***15*** pixels
# * inner 2 - ***10*** pixels

# ## Dependencies libraries

# In[2]:


get_ipython().run_line_magic('load_ext', 'autoreload')
get_ipython().run_line_magic('autoreload', '2')

import sys
sys.path.append("../")
import stanalyzer as sta


# In[2]:


import pandas as pd
import geopandas as gpd
import numpy as np
from shapely import wkt
import matplotlib.pyplot as plt # standard graphics library
from matplotlib.colors import LinearSegmentedColormap
import gzip
import netCDF4
import cartopy.io.img_tiles as cimgt
import cartopy.crs as ccrs # cartographic coordinate reference system
import io
from urllib.request import urlopen, Request
from PIL import Image
from math import floor
from matplotlib import patheffects


# In[103]:


VAR = 'DBZc'
LEVEL = 5

CPT = sta.loadCPT('../stanalyzer/cpt/humidity.cpt')
CMAP = LinearSegmentedColormap('cpt', CPT)
VMIN = 0
VMAX = 65


# In[4]:


pd.read_csv


# In[5]:


pd.read_pickle


# In[3]:


track_frame = sta.read_file('../output/tracking_compressed.pkl')


# In[7]:


track_frame.nc_file


# In[8]:


## This function returns the duration of events
lifes = sta.life_cicle(track_frame,sort=True)
lifes


# In[9]:


### Filter by time
TIME_MIN = 0
TIME_MAX = 15
UNIT = 'h'


# In[10]:


## Apply filter by time
df_filter1 = sta.time_filter(track_frame,TIME_MIN,TIME_MAX,UNIT)


# In[11]:


sta.life_cicle(df_filter1,sort=True)


# In[12]:


fam_test = df_filter1.query('uid == 118')


# In[13]:


fam_test.head()


# In[104]:


def get_cmap(n, name='hsv'):
    '''Returns a function that maps each index in 0, 1, ..., n-1 to a distinct 
    RGB color; the keyword argument name must be a standard mpl colormap name.'''
    return plt.cm.get_cmap(name, n)

def open_ncdf(nc_file):
    ## OPEN RADAR DATA COMPRESS
    with gzip.open(nc_file) as gz:
        with netCDF4.Dataset('dummy', mode='r', memory=gz.read()) as nc:
            data = nc.variables[VAR][0][LEVEL][:].filled()
            extent = [nc.variables['lon0'][:].min(),nc.variables['lon0'][:].max(),nc.variables['lat0'][:].min(),nc.variables['lat0'][:].max()]
    return data,extent

## Image
def image_spoof(self, tile): # this function pretends not to be a Python script
    url = self._image_url(tile) # get the url of the street map API
    req = Request(url) # start request
    req.add_header('User-agent','Anaconda 3') # add user agent to request
    fh = urlopen(req) 
    im_data = io.BytesIO(fh.read()) # get image
    fh.close() # close url
    img = Image.open(im_data) # open image with PIL
    img = img.convert(self.desired_tile_form) # set image format
    return img, self.tileextent(tile), 'lower' # reformat for cartopy

## Scale bar
def scale_bar(ax, proj, length, location=(0.5, 0.05), linewidth=3,
              units='km', m_per_unit=1000):
    # find lat/lon center to find best UTM zone
    x0, x1, y0, y1 = ax.get_extent(proj.as_geodetic())
    # Projection in metres
    utm = ccrs.UTM(utm_from_lon((x0+x1)/2))
    # Get the extent of the plotted area in coordinates in metres
    x0, x1, y0, y1 = ax.get_extent(utm)
    # Turn the specified scalebar location into coordinates in metres
    sbcx, sbcy = x0 + (x1 - x0) * location[0], y0 + (y1 - y0) * location[1]
    # Generate the x coordinate for the ends of the scalebar
    bar_xs = [sbcx - length * m_per_unit/2, sbcx + length * m_per_unit/2]
    # buffer for scalebar
    buffer = [patheffects.withStroke(linewidth=5, foreground="w")]
    # Plot the scalebar with buffer
    ax.plot(bar_xs, [sbcy, sbcy], transform=utm, color='k',
        linewidth=linewidth, path_effects=buffer)
    # buffer for text
    buffer = [patheffects.withStroke(linewidth=3, foreground="w")]
    # Plot the scalebar label
    t0 = ax.text(sbcx, sbcy, str(length) + ' ' + units, transform=utm,
        horizontalalignment='center', verticalalignment='bottom',
        path_effects=buffer, fontsize=14, zorder=4)
    left = x0+(x1-x0)*(0.9)
    # Plot the N arrow
    t1 = ax.text(left, sbcy, u'\u25B2\nN', transform=utm,
        horizontalalignment='center', verticalalignment='bottom',
        path_effects=buffer, fontsize=14, zorder=4)
    # Plot the scalebar without buffer, in case covered by text buffer
    ax.plot(bar_xs, [sbcy, sbcy], transform=utm, color='k',
        linewidth=linewidth, zorder=6)

def utm_from_lon(lon):
    return floor( ( lon + 180 ) / 6) + 1 


# In[105]:


def plot(dframe,t = 0):
    
    ## GET TIMESTAMP LIST
    time_list = [np.unique(sorted(dframe.time))[t]]
    locked_by_time = dframe.query('time <= @time_list')

    geo_cols = []
    list_cols = locked_by_time.columns.to_list()
    tresholds = []
    for c in list_cols:
        if 'geom_' in c or 'line' in c or 'traj' in c or 'vect' in c:
            geo_cols.append(c)
        if 'geom_' in c and 'geom_i' not in c:
            tresholds.append(int(c[-2:]))
            
    geo_cols = sorted(geo_cols)
    
    no_geo_cols = [item for item in list_cols if item not in geo_cols]
    ### CREATE GEODATAFRAME
    geo_frame = gpd.GeoDataFrame(locked_by_time[no_geo_cols])

    ### OPEN MATRIX
    nc_file = locked_by_time.loc[locked_by_time['time'] == time_list[-1]]['nc_file'].unique()[0]
    data,extent = open_ncdf(nc_file)
    ## FILTER
    data[data <= 0 ] = np.nan
    
    ## FIGURE
    cimgt.Stamen.get_image = image_spoof # reformat web request for street map spoofing
    osm_img = cimgt.Stamen('terrain-background'); # spoofed, downloaded street map
    
    plt.figure(figsize=(10,10),dpi=100)
    ax = plt.axes(projection = osm_img.crs)
    ax.set_extent(extent)
    
    ### BACKGROUND
    zoom = 1
    scale = np.ceil(-np.sqrt(2)*np.log(np.divide(zoom,350.0))) # empirical solve for scale based on zoom
    scale = (scale<20) and scale or 19 # scale cannot be larger than 19
    ax.add_image(osm_img, int(scale)) # add OSM with zoom specification
    
    color_range = ['red','C1','C2','C3','C4','C5']
    
    ## Geometries
    for g in range(len(geo_cols)):
        ## Plot Boundary
        if 'geom_' in geo_cols[g] and 'geom_i' not in geo_cols[g]:
            last_time = time_list[-1]
            last_series = locked_by_time.query('time == @last_time')[geo_cols[g]].apply(wkt.loads)
            geo_frame['geometry'] = last_series
            geo_frame.exterior.plot(ax=ax,color=color_range[g],linewidth=1,label=geo_cols[g],transform = ccrs.PlateCarree())
          
        ## Plot Trajectory
        if 'traj' in geo_cols[g]:
            lock_uids = locked_by_time.query('time == @last_time').uid.values.tolist()           
            traject = locked_by_time.query('uid == @lock_uids')[geo_cols[g]].apply(wkt.loads)
            geo_frame['geometry'] = traject
            geo_frame.plot(ax=ax,color='black',label=geo_cols[g],transform = ccrs.PlateCarree())
            
#         # Plot vectores
#         if 'vect' in geo_cols[g]:
#             lock_uids = locked_by_time.query('time == @last_time').uid.values.tolist()           
#             traject = locked_by_time.query('uid == @lock_uids')[geo_cols[g]].apply(wkt.loads)
#             geo_frame['geometry'] = traject
#             geo_frame.plot(ax=ax,color='gray',label=geo_cols[g],transform = ccrs.PlateCarree())

    ### Plot INFO
    buffer = [patheffects.withStroke(linewidth=3, foreground="w")]
    transform = ccrs.PlateCarree()._as_mpl_transform(ax)
    
    for i,row in locked_by_time.query('time == @last_time').iterrows():
        print(row.p0,row.p1,row.uid)
        ax.annotate(s=' UID \n'+str(row.uid),
                        xy=(float(row.lon),float(row.lat)),
                        xycoords=transform, weight='bold', size=10, path_effects=buffer)
        
    
    #### PLOT DATA
    cs = ax.contourf(data, cmap = CMAP, vmin = VMIN, vmax = VMAX, extent = extent,alpha=0.9, transform = ccrs.PlateCarree())
    
    ### LINES LABELS
    gl = ax.gridlines( draw_labels=True,
              linewidth=1, color='gray', alpha=0.5, linestyle='--') # label all axes
    gl.top_labels = False # turn off top label
    gl.right_labels = False # turn off right label
    
    ## PLOT SCALE BAR
    scale_bar(ax, ccrs.Mercator(), 50)
    ## COLOR BAR
    sm = plt.cm.ScalarMappable(cmap=CMAP,norm=plt.Normalize(VMIN,VMAX))
    sm._A = []
    cbar = plt.colorbar(sm,ax=ax,extend='both',pad=0.05,shrink=0.6,orientation='horizontal')
    cbar.set_label('Refletividade (dBZ)')
    
    ## LEGEND
    ax.legend()


# In[106]:


df_filter1.columns


# In[107]:


df_filter1.query('time == 5')


# In[108]:


plot(df_filter1,5)


# In[109]:


import plotly.express as px
import plotly.graph_objects as go


# In[110]:


geo_cols = []
for c in df_filter1.columns:
    if 'geom_' in c or 'line' in c or 'traj' in c or 'vect' in c:
        geo_cols.append(c)


# In[111]:


no_geo_cols = [item for item in df_filter1.columns if item not in geo_cols]


# In[112]:


geo_df = gpd.GeoDataFrame(df_filter1)
geo_df['geometry'] = df_filter1['geom_20'].apply(wkt.loads)


# In[113]:


geo_df.query('time <= 5')


# In[114]:


gd_t = geo_df.query('time == 5')


# In[115]:


gd_t = gd_t.reset_index()


# In[ ]:





# In[116]:


def ncdf(nc_file):
    ## OPEN RADAR DATA COMPRESS
    with gzip.open(nc_file) as gz:
        with netCDF4.Dataset('dummy', mode='r', memory=gz.read()) as nc:
            data = nc.variables[VAR][0][LEVEL][:].filled()
            lon = nc.variables['lon0'][:].filled()
            lat = nc.variables['lat0'][:].filled()
    return data,lon,lat

data_,lon_,lat_ = ncdf(gd_t.nc_file.iloc[0])


# In[117]:


data_[data_ == -9999] = np.nan


# In[118]:


img_pts,lons,lats = [],[],[]

for i in range(len(data_)):
    for j in range(len(data_)):
        if data_[i][j] != np.nan:
            img_pts.append(data_[i][j])
            lons.append(lon_[i][j])
            lats.append(lat_[i][j])
            
img_df = pd.DataFrame(list(zip(lons,lats,img_pts)),
               columns =['lon','lat','metric'])


# In[119]:


import datashader.transfer_functions as tf
from colorcet import rainbow
import datashader as ds
import json

cvs = ds.Canvas(plot_width=data_.shape[0], plot_height=data_.shape[1])
agg = cvs.points(img_df, x='lon', y='lat')
coords_lat, coords_lon = agg.coords['lat'].values, agg.coords['lon'].values

# Corners of the image, which need to be passed to mapbox
coordinates = [[coords_lon[0], coords_lat[0]],
               [coords_lon[-1], coords_lat[0]],
               [coords_lon[-1], coords_lat[-1]],
               [coords_lon[0], coords_lat[-1]]]


# In[306]:


np.nanmin(agg.values)


# In[301]:


# data_[data_ <= 20 ] = np.nan

agg.values = data_
img = tf.shade(agg, cmap=complet_rbb['HEX'].tolist(),alpha=255)[::-1].to_pil()


# In[151]:


from colormap import rgb2hex


# In[293]:


np.nanmin(data_)


# In[ ]:


map(int, results)


# In[280]:


tab = pd.read_table('../stanalyzer/cpt/NEXRAD.cpt',sep=';',names=['RGB','DBZ'])
tab['HEX'] = tab['RGB'].apply(lambda x: rgb2hex(int(x.split(",")[0]),int(x.split(",")[1]),int(x.split(",")[2])))


# In[281]:


tab2 = pd.DataFrame(columns=['RGB','DBZ','HEX'])

tab2['DBZ'] = np.arange(-186,-30)
tab2['HEX'] = ['#000000'] * len(tab2['DBZ'])
tab2['RGB'] = ['0, 0, 0'] * len(tab2['DBZ'])


# In[283]:


complet_rbb = pd.concat([tab2,tab])


# In[308]:


complet_rbb


# In[307]:


CMAP


# In[303]:


complet_rbb['HEX'].tolist()


# In[302]:


rainbow


# In[88]:


geo1 = json.loads(gpd.GeoSeries(gd_t['geom_20'].apply(wkt.loads)).to_json())
geo2 = json.loads(gpd.GeoSeries(gd_t['geom_35'].apply(wkt.loads)).to_json())
geo3 = json.loads(gpd.GeoSeries(gd_t['geom_40'].apply(wkt.loads)).to_json())
traj = json.loads(gpd.GeoSeries(gd_t['trajectory'].apply(wkt.loads)).to_json())


# In[300]:


img


# In[102]:


img


# In[69]:


img


# In[1]:


layers1 = [              
         {'source': geo1,
           'type': "line",
           'color': "black",
           'line':dict(width=1.5)
         },
    
         {'source': geo2,
           'type': "line",
           'color': "royalblue",
           'name':'ds',
           'line':dict(width=0.8)
          },
    
        {'source': geo3,
           'type': "line",
           'color': "red",
           'name':'ds',
           'line':dict(width=0.8)
          },
    
        {'source': traj,
           'type': "line",
           'color': "red",
           'name':'ds',
           'line':dict(width=1.8)
          },
    
        {"sourcetype": "image",
          'below': "traces",
          "source": img,
          "coordinates": coordinates,
        }
        ]


# In[90]:


colorbar=dict(thickness=20,
                           ticklen=3, tickcolor='orange',
                           tickfont=dict(size=14, color='orange'))

showscale=True,
colorscale=[[0, 'green'],[1, 'red']],
cmin=0,
cmax=10


# In[317]:


mapbox_access_token = 'pk.eyJ1Ijoic3BsaW50ZXJzdG0iLCJhIjoiY2tydWxjMzdlMTRxYTJwcGZlYmc0aWJyYSJ9.3-nO2w18a1fbjmAXrH7kEA'

hover_info = dict(bgcolor="white",
                font_size=14,
                font_family="Rockwell")

hover_temp = "<b>%{customdata[0]} </b><br><br>" + \
             "longitude: %{lon}<br>" + \
             "latitude: %{lat}<br>" + \
             "UID: %{customdata[1]}<extra></extra>"

### FIGURE
data = go.Scattermapbox(customdata=gd_t[['timestamp','uid','time']],
                        lat=list(gd_t["lat"]), lon=list(gd_t["lon"]),
                        textfont=dict(size=16, color='black'),
                        hoverlabel=hover_info,
                        mode='markers',
                        marker=dict(size=10, color='black',
                                showscale=True,colorscale=[[0, 'green'],[1, 'red']],cmin=0,cmax=50,
                                colorbar=dict(thickness=20,
                                   ticklen=3, tickcolor='orange',
                                   tickfont=dict(size=14, color='orange'))),
                        hovertemplate = hover_temp,
                      )

layout = dict(margin=dict(l=0, t=0, r=0, b=0, pad=0),
              mapbox=dict(layers=layers1,
                          accesstoken=mapbox_access_token,
                          center=dict(lat=-3.14, lon=-60),
                          style='light',
                          zoom=6),
              
             )

fig = go.Figure(data=data, layout=layout)

fig.update_mapboxes()

fig.show()


# In[ ]:


fig = go.Figure()

fig = px.scatter_mapbox(gd_t, lat="lat", lon="lon", hover_name="uid", hover_data=["time", "uid","max_ref_20"],
                        color_discrete_sequence=["fuchsia"],width=900, height=700)

fig.update_layout(mapbox=dict(layers=layers,
                  style="stamen-terrain",
                  center =  {'lon': -60, 'lat': -3.14},
                  zoom = 6))

# buttons for updatemenu
buttons = [dict(method='restyle',
                label='linear',
                visible=True,
                args=[{'label': 'linear',
                       'visible':[True, False],}]),
           
           dict(method='restyle',
                label='sqrt',
                visible=True,
                args=[{'label': 'linear',
                       'visible':[False, True],}])]

# specify updatemenu        
um = [{'buttons':buttons, 'direction': 'down'}]

fig.update_layout(updatemenus=um)
fig.show()


# In[ ]:


fig.update_layout(
    updatemenus=[
            dict(
            buttons=list([
                dict(label="High",
                     method="update",               
                     args=[{"visible": [False]},
                          {"shapes": layers1[3]}])
            ]),
                
            type = "buttons",
            direction="right",
            pad={"r": 0, "t": 0},
            showactive=True,
            x=0.38,
            xanchor="left",
            y=1,
            yanchor="top"
            )
    ])


# In[ ]:


width=800, height=400


# In[ ]:


fig = px.density_mapbox(img_df,lat="lat", lon="lon",z='metric',radius=5)

fig.update_layout(mapbox=dict(layers=layers,
                  style="stamen-terrain",
                  center =  {'lon': -60, 'lat': -3.14},
                  zoom = 6
                  ),
                 )

fig.update_layout(mapbox=dict(layers=layers,
                  style="stamen-terrain",
                  center =  {'lon': -60, 'lat': -3.14},
                  zoom = 6
                  ),
                 )


# In[ ]:





# In[ ]:





# In[ ]:


x_states = []
y_states = []

for i,row in gd_t.iterrows():
    x = row.geometry.exterior.xy[0].tolist()
    y = row.geometry.exterior.xy[0].tolist()
    for segment in range(len(x)):
        x_states = x_states + x
        y_states = y_states + y


# In[ ]:


state_trace = [dict(
            type = 'scatter',
            legendgroup = "States",
            line = dict(color='white', width=1),
            x = x_states,
            y = y_states,
            hoverinfo='none',
            showlegend=False,
            mode='lines')]


# In[52]:


dir(go.Choroplethmapbox)


# In[39]:


fig = go.Figure(go.Choroplethmapbox(geojson=gd_t.geometry,
                                    locations=gd_t.index,
                                    z=gd_t.uid,
                                    colorscale='Viridis',
                                    zauto=True,
                                    marker_opacity=0.5,
                                    marker_line_width=0.5)
                )
# fig.update_layout(mapbox_style='white-bg',
#                   #mapbox_accesstoken=mapbox_token,
#                   mapbox_zoom=12,
#                   mapbox_center={'lat': 45.41117, 'lon': -75.69812})
# fig.update_layout(margin={'r':0, 't':0, 'l':0, 'b':0})
# pio.renderers.default = 'browser'
fig.show()


# In[41]:


df


# In[43]:


import plotly.express as px

df = px.data.gapminder()
fig = px.choropleth(df, locations="iso_alpha",
                    color="lifeExp", # lifeExp is a column of gapminder
                    hover_name="country", # column to add to hover information
                    color_continuous_scale=px.colors.sequential.Plasma,
                    animation_frame='year')

# fig2 = px.scatter_geo(df, locations="iso_alpha",
#                     size="lifeExp", # lifeExp is a column of gapminder
#                     hover_name="country", # column to add to hover information
#                     color_continuous_scale=px.colors.sequential.Plasma,
#                     animation_frame='year')
fig.add_trace(fig2.data[0])

# for i, frame in enumerate(fig.frames):
#     fig.frames[i].data += (fig2.frames[i].data[0],)
fig.show()


# In[ ]:


get_ipython().system('conda install -c plotly plotly_express -y')


# In[ ]:





# In[ ]:


open_ncdf(df_filter1['nc_file'].iloc[0])


# In[ ]:





# In[ ]:





# In[ ]:


from shapely import wkt
import numpy as np


# In[ ]:


fam_test = df_filter1.query('uid == 118')
fam_test['geom_35'].fillna(value='GEOMETRYCOLLECTION EMPTY', inplace=True)


geo_frame = gpd.GeoDataFrame()
geo_frame['geometry'] = fam_test['geom_35'].apply(wkt.loads)


# In[ ]:


geo_frame.tail(2).plot()


# In[ ]:


fam_test['geom_35'].apply(wkt.loads)


# In[ ]:


tt = gpd.GeoDataFrame()
tt['geometry'] = fam_test['geom_intersect'].apply(wkt.loads)


# In[ ]:


tt.plot()


# In[ ]:


import geopandas as gpd


# In[ ]:


### PATH
PATH_FILE = '../tracks/S201409070000_E201409100000_VDBZc_T20_L5.pkl'


# In[ ]:


### Read tracking file
df = sta.read_file(PATH_FILE)
df.head()


# In[ ]:


sta.life_cicle()


# In[ ]:


life_cicle(df)


# In[ ]:


# Dataframe library
import pandas as pd
# Numerical Python library
import numpy as np
# netCDF4 library
import netCDF4
# Import gzip to open netCDF
import gzip
# Visualization library
import matplotlib.pyplot as plt


# ## Variables
# 
# **Fam_Nº**-> Refers to the number of the Tracked Family.
# <br>
# 
# **timestamp** ->A digital record of the time of occurrence of a particular event.
# <br>
# **time** -> Refers to the tracking time in the algorithm.
# <br>
# **uid** -> Unique IDentifier, it is used to generate the families.
# <br>
# **id_t** -> Referring cluster identifier at the time of tracking occurrence. From the DBSCAN clustering algorithm.
# <br>
# **lat** -> Refers latitude centroid, taken from the reference matrix of the original nc files.
# <br>
# **lon** -> Refers longitude centroid, taken from the reference matrix of the original nc files.
# <br>
# **p0** -> The first coordinate point of centroid in matrix (clusters or nc_file): (p0,p1)=(x,y)=(lon,lat).
# <br>
# **p1** -> The second coordinate point of centroid in matrix (clusters or nc_file): (p0,p1)=(x,y)=(lon,lat).
# <br>
# **size_%THRESHOLD** -> Total number of Pixels in the main cluster. Each point depends on the sensor's spatial resolution (pixel size): RADAR 2x2km.
# <br>
# **mean_ref_%THRESHOLD** -> Averaged reflectivity of the cluster. Value in dBZ.
# <br>
# **max_ref_%THRESHOLD** -> Max reflectivity of the cluster. Value in dBZ.
# <br>
# **angle_%THRESHOLD_orig** -> Original displacement angle of the cluster at the current time.
# <br>
# **angle_%THRESHOLD_cor** -> Corrected displacement angle of the cluster at the current time. 
# <br>
# **vel_%THRESHOLD_orig** -> Original displacement velocity of the cluster at the current time in kilometers per hour (km/h).
# <br>
# **vel_%THRESHOLD_cor** -> Corrected displacement velocity of the cluster at the current time in kilometers per hour (km/h).
# <br>
# **mean_total_ref_%THRESHOLD** -> Average reflectivity of the inner clusters by threshold (Value in dBZ).
# <br>
# **total_size_%THRESHOLD** -> Total size of inner clusters by threshold (number of pixels).
# <br>
# **n_cluster_%THRESHOLD** -> Total number of inner clusters by Threshold.
# <br>
# **avg_angle_%THRESHOLD** -> Averaged angle for the inner cluster by threshold (Value in degree).
# <br>
# **avg_vel_%THRESHOLD** -> Averaged velocity for inner clusters by threshold (Value in km/h).
# <br>
# **status** -> Status of occurrence, type: NEW-> New cluster; CONT-> Continous cluster; SPLT -> Splitted cluster; MERG -> Merged Cluster.
# <br>
# **delta_t** -> Time interval for cluster life cycle.
# <br>
# **nc_file** -> Path of netCDF file.
# <br>
# **cluster_file** -> Path of cluster file (From DBSCAN).
# <br>
# **dsize_%THRESHOLD** -> Difference between the sizes of two consecutive clusters (in Pixel).
# <br>
# **dmean_ref_%THRESHOLD** -> Difference between the mean reflectivities of two consecutive clusters for main threshold (in dBZ).
# <br>
# **dmean_total_ref_%THRESHOLD** -> Difference between the mean reflectivities of all clusters between two consecutive times for an inner threshold (in dBZ).
# <br>
# **dtotal_size_%THRESHOLD** -> Difference between the total size (in pixel) of all clusters between two consecutive times for an inner threshold (values in pixel).

# ## Read tracking file
# 
# Tracking DataFrame.

# In[ ]:


fam_tracking = pd.read_pickle("./S201409070000_E201409100000_VDBZc_T20_L5.pkl")
fam_tracking


# ### Example how to select a FAM by uid

# In[ ]:


uid = 97
selected_fam = fam_tracking.query('uid == @uid')
selected_fam


# ## Example how to select a cluster in the family

# In[ ]:


line = 0 #first line 

selected_line = selected_fam.iloc[[line]]
selected_line


# ## Example how to open the cluster file and the original data to extract reflectivity values

# In[ ]:


## OPEN CLUSTERS
def open_cluster(path):
    try:
        cluster = np.load(path['cluster_file'].values[0])['arr_0']
        cluster[cluster == 0] = np.NAN
        return cluster
    except:
        print('File not found!')


# In[ ]:


selected_line


# In[ ]:


cluster_matrix_all = open_cluster(selected_line)
print('Original dimensions of cluster->',cluster_matrix_all.shape)

THRESHOLD_LEVEL = 0 #to select the main threshold (ex: 0-20dBZ,1-35dBZ,2-40dBZ)
cluster_matrix = cluster_matrix_all[:,:,THRESHOLD_LEVEL]
print('Selected dimensions of cluster->',cluster_matrix.shape)


# In[ ]:


### OPEN NETCDF
def open_file(file_path):
    VAR_NAME = 'DBZc'
    LEVEL = 5 #2.5km height
    THRESHOLDS = [20,35,40] #dBZ
    with gzip.open(file_path['nc_file'].values[0]) as gz:
        with netCDF4.Dataset('dummy', mode='r', memory=gz.read()) as nc:
            data = nc.variables[VAR_NAME][0][LEVEL][:].filled()
            data[data == -9999.] = np.NAN
    data[data < THRESHOLDS[0]] = np.nan
    return data


# In[ ]:


nc_matrix = open_file(selected_line)
print('NetCDF Max/Min values (thresholded):\n',np.nanmax(nc_matrix),np.nanmin(nc_matrix))


# In[ ]:


fig, (ax,ax1) = plt.subplots(1,2, figsize=(15,6))

ax.imshow(nc_matrix)
ax1.imshow(cluster_matrix);
ax.set_title('Original file')
ax1.set_title('Cluster file');


# ## Extracting reflectivities from the selected cluster
# 
# To extract the reflectivity values of an individual cluster, you will need to choose the tracking 'id_t', this should be done as follows:
# 
# Visualization of individual line.

# In[ ]:


selected_line


# In[ ]:


### This line shows that id_t is equal to 20. 
selected_id_t = selected_line.id_t.values[0]

### Get XY coordinates from cluster matrix
x,y = np.where(cluster_matrix == selected_id_t)

### Get reflectivities values from nc_file cluster
dbz_list = nc_matrix[x,y]


# In[ ]:


print('List with reflectivity values of an individual cluster.\n',dbz_list)


# ### Cluster location view

# In[ ]:


fig, (ax,ax1) = plt.subplots(1,2, figsize=(15,6))
ax.imshow(nc_matrix)
ax1.imshow(cluster_matrix);
ax.set_title('Original file')
ax1.set_title('Cluster file');

ax.scatter(selected_line.p0,selected_line.p1,marker='x',color='r',s=100)
ax1.scatter(selected_line.p0,selected_line.p1,marker='x',color='r',s=100)


# In[ ]: