OUTPUT TRACKING ALGORITHM¶

%load_ext autoreload
%autoreload 2

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

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

VAR = 'DBZc'
LEVEL = 5

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

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

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

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

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 

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()

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

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)

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

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

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

gd_t = gd_t.reset_index()

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])

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

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'])

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]]]

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

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

from colormap import rgb2hex

map(int, results)

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])))

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'])

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

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())

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

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()

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"
            )
    ])

width=800, height=400

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
                  ),
                 )

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

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')]

!conda install -c plotly plotly_express -y

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

from shapely import wkt
import numpy as np

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)

geo_frame.tail(2).plot()

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

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

tt.plot()

import geopandas as gpd

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

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

sta.life_cicle()

life_cicle(df)

# 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

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

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

line = 0 #first line 

selected_line = selected_fam.iloc[[line]]
selected_line

## 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!')

selected_line

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)

### 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

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

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');

selected_line

### 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]

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

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)