"""
PQPF (Hong Kong)
========================================================
This example demonstrates how to perform
operational PQPF up to three hours from
raingauge and radar data, using data from Hong Kong.
"""

###########################################################
# Definitions
# --------------------------------------------------------
#

import os
import numpy as np
import pandas as pd
from pyresample import utils
import xarray
import cartopy.feature as cfeature
import cartopy.crs as ccrs
import matplotlib
import imageio
import matplotlib.pyplot as plt
from matplotlib.colors import BoundaryNorm, LinearSegmentedColormap

from swirlspy.rad.iris import read_iris_grid
from swirlspy.qpe.utils import dbz2rr, rr2rf, locate_file, timestamps_ending
from swirlspy.qpe.utils import temporal_interp, multiple_acc
from swirlspy.obs.rain import Rain
from swirlspy.qpf.rover import rover
from swirlspy.qpf.sla import sla
from swirlspy.qpf.rover import rover
from swirlspy.core.resample import grid_resample

plt.switch_backend('agg')
THIS_DIR = os.getcwd()
os.chdir(THIS_DIR)

start_time = pd.Timestamp.now()

#############################################################
# Initialising
# ---------------------------------------------------
#
# This section demonstrates extracting
# radar reflectivity data.
#
# Step 1: Define a basetime.
#

# Supply basetime
basetime = pd.Timestamp('201902190800')

##############################################################################
# Step 2: Using basetime, generate timestamps of desired radar files
# timestamps_ending() and locate files using locate_file().

# Obtain radar files
dir = THIS_DIR + '/../tests/samples/iris/ppi'
located_files = []
radar_ts = timestamps_ending(
    basetime,
    duration=pd.Timedelta(minutes=60)
)

for timestamp in radar_ts:
    located_files.append(locate_file(dir, timestamp))

##############################################################################
# Step 3: Read data from radar files into xarray.DataArray
# using read_iris_grid().

reflectivity_list = []  # stores reflec from read_iris_grid()
for filename in located_files:
    reflec = read_iris_grid(
        filename
    )
    reflectivity_list.append(reflec)

########################################################################
# Step 4: Define the target grid as a pyresample AreaDefinition.

# Defining target grid
area_id = "hk1980_250km"
description = ("A 250 m resolution rectangular grid "
               "centred at HKO and extending to 250 km "
               "in each direction in HK1980 easting/northing coordinates")
proj_id = 'hk1980'
projection = ('+proj=tmerc +lat_0=22.31213333333334 '
              '+lon_0=114.1785555555556 +k=1 +x_0=836694.05 '
              '+y_0=819069.8 +ellps=intl +towgs84=-162.619,-276.959,'
              '-161.764,0.067753,-2.24365,-1.15883,-1.09425 +units=m '
              '+no_defs')
x_size = 500
y_size = 500
area_extent = (587000, 569000, 1087000, 1069000)
area_def_tgt = utils.get_area_def(
    area_id, description, proj_id, projection, x_size, y_size, area_extent
)

##############################################################################
# Step 5: Reproject the radar data from read_iris_grid() from Centered
# Azimuthal (source) projection to HK 1980 (target) projection.
#

# Extracting the AreaDefinition of the source projection
area_def_src = reflectivity_list[0].attrs['area_def']

# Reprojecting
reproj_reflectivity_list = []
for reflec in reflectivity_list:
    reproj_reflec = grid_resample(
        reflec, area_def_src, area_def_tgt,
        coord_label=['easting', 'northing']
    )
    reproj_reflectivity_list.append(reproj_reflec)

############################################################################
# Step 6: Assigning reflectivity xarrays at the last three timestamps to
# variables for use during ROVER QPF.

xarray1 = reproj_reflectivity_list[-3]
xarray2 = reproj_reflectivity_list[-2]
xarray3 = reproj_reflectivity_list[-1]

initialising_time = pd.Timestamp.now()

###########################################################
# Running ROVER and Semi-Lagrangian Advection
# -------------------------------------------
#
# 1. Concatenate required reflectivity xarrays along time dimension.
# 2. Run ROVER on all members, with the concatenated xarray as the input.
# 3. Store motion field xarrays as a list.
# 4. Perform Semi-Lagrangian Advection on all members using the motion fields
#    from rover.
# 5. Store forecasted reflectivities as list.


# Combining two reflectivity DataArrays
# the order of the coordinate keys is now ['y', 'x', 'time']
# as opposed to ['time', 'x', 'y']
reflec_concat_6min = xarray.concat([xarray2, xarray3], dim='time')
reflec_concat_12min = xarray.concat([xarray1, xarray3], dim='time')

# Running rover on 4 members
# Mem-1
u1, v1 = rover(
    reflec_concat_6min,
    start_level=1, max_level=7, sigma=2.5
)
# Mem-2
u2, v2 = rover(
    reflec_concat_12min,
    start_level=2, max_level=7, sigma=2.5
)
# Rover-A
u3, v3 = rover(
    reflec_concat_6min
)
# Mem-4
u4, v4 = rover(
    reflec_concat_12min
)

# Storing motion fields for quiver plot
motion_list = [[u1, v1], [u2, v2], [u3, v3], [u4, v4]]

rover_time = pd.Timestamp.now()

# Running SLA on all members

z1 = sla(
    reflec_concat_6min,
    u1, v1, steps=30
)

z2 = sla(
    reflec_concat_12min,
    u2, v2, steps=15
)

z3 = sla(
    reflec_concat_6min,
    u3, v3, steps=30
)

z4 = sla(
    reflec_concat_12min,
    u4, v4, steps=15
)

# appending all reflectivities to list

z_sla_list = [z1, z2, z3, z4]

sla_time = pd.Timestamp.now()

##############################################################
# Concatenating observed and forecasted reflectivities
# --------------------------------------------------------
#
# 1. Add forecasted reflectivity to reproj_reflectivity_list.
# 2. Concatenate observed and forecasted reflectivity
#    xarray.DataArrays along the time dimension.
# 3. Concatenate reflectivities of different members along
#    a fourth dimension.

z_cat_list = []
for reflectivity in z_sla_list:
    z_all = reproj_reflectivity_list + [reflectivity[1:, ...]]
    z_cat = xarray.concat(z_all, dim='time')
    z_cat_list.append(z_cat)

# Concatenating different members
z_ens_6min = xarray.concat(
    [z_cat_list[0],
     z_cat_list[2]],
    xarray.IndexVariable('member', ['Mem-1', 'Mem-3'])
)

z_ens_12min = xarray.concat(
    [z_cat_list[1],
     z_cat_list[3]],
    xarray.IndexVariable('member', ['Mem-2', 'Mem-4'])
)
concat_time = pd.Timestamp.now()


##########################################################################
# Accumulating hourly rainfall for 3-hour forecast
# ------------------------------------------------
#
# Hourly accumulated rainfall is calculated every 30 minutes, the first
# endtime is the basetime i.e. T+0min.
#
# | Step 1: Convert reflectivity in dBZ to rainrates in mm/h with dbz2rr().
# | Step 2: Replacing time coordinates with end times instead of start times.
# | Step 3: Interpolate rainrates temporally to ensure that rainrates are
# |         equally spaced along the time axis.
# | Step 4: Convert rainrates to rainfalls in 6 and 12 minutes with rr2rf().

# 6 minute ensembles
rainrates_ens_6min = dbz2rr(z_ens_6min, a=58.53, b=1.56)
rainrates_ens_6min_endtime = rainrates_ens_6min.copy()
rainrates_ens_6min_endtime.coords['time'] = \
    [
        pd.Timestamp(t) + pd.Timedelta(minutes=6)
        for t in rainrates_ens_6min_endtime.coords['time'].values
    ]
rainfalls_ens_6min = rr2rf(rainrates_ens_6min_endtime)

# 12 minute ensembles
rainrates_ens_12min = dbz2rr(z_ens_12min, a=58.53, b=1.56)
rainrates_ens_12min_endtime = rainrates_ens_12min.copy()
rainrates_ens_12min_endtime.coords['time'] = \
    [
        pd.Timestamp(t) + pd.Timedelta(minutes=12)
        for t in rainrates_ens_12min_endtime.coords['time'].values
    ]
rainrates_ens_12min_interp = temporal_interp(
    rainrates_ens_12min,
    pd.Timestamp('201902190700'),
    pd.Timestamp('201902191100')
)
rainfalls_ens_12min = rr2rf(rainrates_ens_12min_interp)

######################################################################
# Step 3: Compute hourly accumulated rainfall every 30 minutes.

acc_rf_6min = multiple_acc(
    rainfalls_ens_6min, basetime, basetime+pd.Timedelta(hours=3)
)

acc_rf_12min = multiple_acc(
    rainfalls_ens_12min, basetime, basetime+pd.Timedelta(hours=3)
)

acc_time = pd.Timestamp.now()

######################################################################
# Calculating Probability Exceeding Threshold
# ----------------------------------------------
#
# In this example, thresholds are 0.5mm, 5mm,
# 10mm, 30mm, 50mm and 70 mm.
#
# Probabilities are 0%, 25%, 50%, 75% and 100%, as there are
# four members.
#
# Steps:
#
# 1. Define threshold.
# 2. Concatenate rainfalls from 6 minute and 12 minute ensembles
#    along the member dimension.
# 3. Use a loop to calculate the probability of exceeding threshold
#    for each gridcell and store results in list.
# 4. Concatenate the contents of the list. Result is an xarray with
#    dimensions (threshold, time, y, x).
# 5. Plot the results using xarray.plot(). In this example,
#    only the probability exceeding 0.5mm rainfall
#    every hour will be plotted.

# Define threshold
threshold = [0.5, 5., 10., 30., 50., 70.]

# Concatenating rainfalls_ens_6min and rainfalls_ens_12min
acc_rf = xarray.concat(
    [acc_rf_6min, acc_rf_12min], dim='member'
    )

# Define list to store probabilities of exceeding rainfall thresholds
# List corresponds to threshold
prob_list = []

# Calculate probability
for th in threshold:
    bool_forecast = acc_rf >= th
    count = bool_forecast.sum(dim='member')
    prob = count/len(bool_forecast.coords['member'].values) * 100
    prob_list.append(prob)

# Generate coordinate xarray for threshold
th_xarray = xarray.IndexVariable(
    'threshold', threshold
)

# concatenate
prob_rainfall = xarray.concat(
    prob_list,
    dim=th_xarray
)

prob_rainfall.attrs['long_name'] = "Probability of Exceeding Threshold"
prob_rainfall.attrs['units'] = "%"

# Plot the results

# Extracting the crs
crs = area_def_tgt.to_cartopy_crs()

# Defining cmap
cmap = LinearSegmentedColormap.from_list(
    'custom blue', ['#FFFFFF', '#000099']
)

# Obtaining xarray slices to be plotted
threshold_list = [0.5]
timelist = [basetime + pd.Timedelta(hours=i) for i in range(4)]
da_plot = prob_rainfall.sel(
    threshold=threshold_list, time=timelist
)

# Defining coastlines
hires = cfeature.GSHHSFeature(
   levels=[1],
   scale='h',
   edgecolor='k'
)

for th in threshold_list:
    # Plotting
    p = da_plot.sel(threshold=th).plot(
        col='time', col_wrap=2,
        subplot_kws={'projection': crs},
        cbar_kwargs={
            'ticks': [0, 25, 50, 75, 100],
            'format': '%.3g'
        },
        cmap=cmap
        )

    for idx, ax in enumerate(p.axes.flat):
        ax.add_feature(hires)  # coastlines
        ax.gridlines()  # gridlines
        ax.set_title(
            f"% Exceeding {th}mm rainfall\n"
            f"Based @ {basetime.strftime('%H:%MH')}",
            loc='left',
            fontsize=8
        )
        ax.set_title(
            ''
        )
        ax.set_title(
            f"{basetime.strftime('%Y-%m-%d')} \n"
            f"Valid @ {timelist[idx].strftime('%H:%MH')} ",
            loc='right',
            fontsize=8
        )
    plt.savefig(
        THIS_DIR +
        f"/../tests/outputs/p_{th}.png",
        dpi=300
        )

prob_time = pd.Timestamp.now()

###############################################################################
# Rainfall percentiles
# ----------------------------
#
# 1. Using xarray.DataArray.mean(),
#    calculate the mean rainfall of all gridpoints.
# 2. Using xarray.DataArray.min(),
#    find the minimum rainfall of all gridpoints.
# 3. Using xarray.DataArray.max(),
#    find the maximum rainfall of all gridpoints.
# 4. Using xarray.DataArray.quantile()
#    find the 25th, 50th and 75th percentile rainfall
#    of all gridpoints.
# 5. Concatenate rainfall along percentile dimension.
# 6. Plot results using xarray.plot(). In this example
#    only the maximum percentile rainfall
#    every hour will be plotted.

# mean
mean_rainfall = acc_rf.mean(dim='member', keep_attrs=True)
# max
max_rainfall = acc_rf.max(dim='member', keep_attrs=True)
# min
min_rainfall = acc_rf.min(dim='member', keep_attrs=True)
# quartiles
q_rainfall = acc_rf.quantile(
    [.25, .5, .75], dim='member',
    interpolation='nearest', keep_attrs=True)

# generate index
percentile = xarray.IndexVariable(
    'percentile',
    ['Minimum',
     '25th Percentile',
     '50th Percentile',
     'Mean',
     '75th Percentile',
     'Maximum']
    )

# concatenating

p_rainfall = xarray.concat(
    [min_rainfall,
     q_rainfall.sel(quantile=.25).squeeze().drop('quantile'),
     q_rainfall.sel(quantile=.50).squeeze().drop('quantile'),
     mean_rainfall,
     q_rainfall.sel(quantile=.75).squeeze().drop('quantile'),
     max_rainfall],
    dim=percentile
)

p_rainfall.attrs['long_name'] = "Hourly Accumulated Rainfall"
p_rainfall.attrs['units'] = "mm"

# Defining levels

# Defining the colour scheme
levels = [
    0, 0.5, 2, 5, 10, 20,
    30, 40, 50, 70, 100, 150,
    200, 300, 400, 500, 600, 700
]

cmap = matplotlib.colors.ListedColormap([
    '#ffffff', '#9bf7f7', '#00ffff', '#00d5cc', '#00bd3d', '#2fd646',
    '#9de843', '#ffdd41', '#ffac33', '#ff621e', '#d23211', '#9d0063',
    '#e300ae', '#ff00ce', '#ff57da', '#ff8de6', '#ffe4fd'
])

norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

# Obtaining xarray slices to be plotted
plot_positions = ['Maximum']
timelist = [basetime + pd.Timedelta(hours=i) for i in range(4)]
da_plot = p_rainfall.sel(percentile=plot_positions, time=timelist)

for pos in plot_positions:
    # Plotting
    p = da_plot.sel(percentile=pos).plot(
        col='time', col_wrap=2,
        subplot_kws={'projection': crs},
        cbar_kwargs={
            'ticks': levels,
            'format': '%.3g'
        },
        cmap=cmap,
        norm=norm
        )

    for idx, ax in enumerate(p.axes.flat):
        ax.add_feature(hires)  # using GSHHS coastlines defined previously
        ax.gridlines()
        ax.set_title(
            f"{pos} Radar-Based Rainfall\n"
            f"Based @ {basetime.strftime('%H:%MH')}",
            loc='left',
            fontsize=8
        )
        ax.set_title(
            ''
        )
        ax.set_title(
            f"{basetime.strftime('%Y-%m-%d')} \n"
            f"Valid @ {timelist[idx].strftime('%H:%MH')} ",
            loc='right',
            fontsize=8
        )
        position = pos.split(" ")[0]
    plt.savefig(
        THIS_DIR +
        f"/../tests/outputs/rainfall_{position}.png",
        dpi=300
    )
ptime = pd.Timestamp.now()

##############################################################################
# Extract the rainfall values at a specified location
# ------------------------------------------------------------------
#
# In this example, the rainfall values at the location is assumed
# to be the same as the nearest gridpoint.
#
# 1. Read information regarding the radar stations into a pandas.DataFrame.
# 2. Extract the rainfall values at the nearest gridpoint to location
#    for given timesteps (in this example, 30 minute intervals).
# 3. Store rainfall values over time in an xarray.DataArray.
# 4. Plot the time series of rainfall with boxplots at desired station.
#    In this case, the 15th percentile member is plotted.

# Getting radar station coordinates
df = pd.read_csv(
    os.path.join(THIS_DIR, "../tests/samples/hk_raingauge.csv"),
    usecols=[0, 1, 2, 3, 4]
)

# Extract rainfall values at gridpoint closest to the
# location specified for given timesteps and storing it
# in xarray.DataArray.
station_rf_list = []
station_name = []
for index, row in df.iterrows():
    station_rf_list.append(p_rainfall.sel(
        northing=row[1], easting=row[2],
        method='nearest'
    ).drop('northing').drop('easting'))
    station_name.append(row[0])

station_name_index = xarray.IndexVariable(
    'ID', station_name
    )
station_rf = xarray.concat(
    station_rf_list,
    dim=station_name_index
)

# Extracting the 15th ranked station
xr_15_percentile = station_rf.quantile(
    .15, dim='ID', interpolation='nearest').drop('quantile')

# Plotting
_, tax = plt.subplots(figsize=(20, 15))
plt.plot(
    np.arange(1, len(xr_15_percentile.coords['time'].values) + 1),
    xr_15_percentile.loc['Mean'].values,
    'ko-'
)  # plot line

# Storing percentiles as dictionary to call ax.bxp for boxplot
stats_list = []
for i in range(len(xr_15_percentile.coords['time'].values)):
    stats = {
        'med': xr_15_percentile.loc['50th Percentile'].values[i],
        'q1': xr_15_percentile.loc['25th Percentile'].values[i],
        'q3': xr_15_percentile.loc['75th Percentile'].values[i],
        'whislo': xr_15_percentile.loc['Maximum'].values[i],
        'whishi': xr_15_percentile.loc['Minimum'].values[i]
    }
    stats_list.append(stats)

# Plot boxplot
tax.bxp(
    stats_list, showfliers=False
)

# Labels
xcoords = xr_15_percentile.coords['time'].values
xticklabels = [pd.to_datetime(str(t)).strftime("%-H:%M") for t in xcoords]
tax.set_xticklabels(xticklabels)
tax.xaxis.set_tick_params(labelsize=20)
tax.yaxis.set_tick_params(labelsize=20)
plt.title('Time Series of Hourly Accumulated Rainfall', fontsize=25)
plt.ylabel("Hourly Accumulated Rainfall [mm]", fontsize=22)
plt.xlabel("Time", fontsize=18)
plt.savefig(THIS_DIR+"/../tests/outputs/pqpf_time_series.png")

extract_time = pd.Timestamp.now()

######################################################################
# Checking run time of each component
# --------------------------------------------------------------------
#

print(f"Start time: {start_time}")
print(f"Initialising time: {initialising_time}")
print(f"Rover time: {rover_time}")
print(f"SLA time: {sla_time}")
print(f"Concatenating time: {concat_time}")
print(f"Accumulating rainfall time: {acc_time}")
print(
    "Calculate and plot probability exceeding threshold: "
    f"{prob_time}"
    )
print(
    f"Plotting rainfall maps: {ptime}"
)
print(f"Extracting and plotting time series time: {extract_time}")

print(f"Time to initialise: {initialising_time-start_time}")
print(f"Time to run rover: {rover_time-initialising_time}")
print(f"Time to perform SLA: {sla_time-rover_time}")
print(f"Time to concatenate xarrays: {concat_time - sla_time}")
print(f"Time to accumulate rainfall: {acc_time - concat_time}")
print(
    "Time to calculate and plot probability exceeding threshold: "
    f"{prob_time-acc_time}"
    )
print(f"Time to plot rainfall maps: {ptime-prob_time}")
print(
    f"Time to extract station data and plot time series: "
    f"{extract_time-ptime}"
)