"""
Quantitative Precipitation Estimate (QPE)
====================================================

This example demonstrates how to perform QPE,
using raingauge data from Manila.
"""

import os
import pandas
import numpy as np
from swirlspy.qpe.rfmap import RainGauge
from swirlspy.qpe.rfmap import to_map_coordinates
from swirlspy.qpe.rfmap import epsilon
from swirlspy.qpe.rfmap import rbf_interpolation
from swirlspy.qpe.rfmap import save_fig


###############################################################################
# Read files and declare raingauge object
# ---------------------------------------
#
# Read files in consecutive timesteps to obtain accumulated rainfall,
# return as a raingauge object.

THIS_DIR = '../tests/qpe'
os.chdir(THIS_DIR)
    
filelist = ["rf60m_201808102340_qced", "rf60m_201808102350_qced", 
            "rf60m_201808110000_qced", "rf60m_201808110010_qced", 
            "rf60m_201808110020_qced", "rf60m_201808110030_qced",
            "rf60m_201808110040_qced"]
table = pandas.read_csv(filepath_or_buffer=filelist[0], header=None, 
						delim_whitespace=True, skiprows=[0], 
						usecols=[1,2,3], dtype=float, na_values=3276.7, 
						low_memory=False) #override nan values
nrow, ncol = table.shape
accu_rf = np.zeros(nrow)
#accumulate rainfall
for i in range (0,len(filelist)):
    file = open(filelist[i])
    if i == 0:
        start_time = file.readline()
    elif i == len(filelist)-1:
        end_time = file.readline()
    table = pandas.read_csv(filepath_or_buffer=filelist[i], header=None,  
							delim_whitespace=True,  skiprows=[0],
							usecols=[1,2,3], dtype=float, na_values=3276.7, 
							low_memory=False) #override nan values
    table = table.dropna()
    gauge_lat = table[1].values
    gauge_lon = table[2].values
    accu_rf = accu_rf + table[3].values
    
#Declare raingauge object
rg_object = RainGauge(start_time, end_time, gauge_lat, 
					    gauge_lon, accu_rf)   


###############################################################################
# Interpolation
# -------------
#
# Perform RBF interpolation.
lllon = 119.8
urlon = 122.0
lllat = 13.6
urlat = 16.5
resolution = 1000
mode= "rbf"
m_rg_lon, m_rg_lat, m_lllon, m_lllat, m_urlon, m_urlat, \
rbf_gx, rbf_gy = to_map_coordinates(lllon, urlon, lllat, urlat, 
                                    rg_object.gauge_lon, 
                                    rg_object.gauge_lat,
                                    resolution) 
epsilon_value = epsilon(urlon, lllon, resolution, rg_object.rg_number)
interpolated_rf = rbf_interpolation(m_rg_lon, m_rg_lat, rbf_gx, 
                                    rbf_gy, rg_object.rainfall, 
                                    epsilon_value) 
###############################################################################
# Generate rainfall map
# ---------------------
#
# Here, we create a rainfall map with Cartopy basemap and
# showing raingauge values. Please refer to 
# `swirlspy.qpe.rfmap.save_fig` for customization.

#Options here: showing raingauge values, Cartopy basemap
save_fig(mode, rg_object.start_time, rg_object.end_time, lllon, urlon, lllat, 
        urlat, m_rg_lon, m_rg_lat, rg_object.rainfall, 
        rbf_gx, rbf_gy, interpolated_rf, True, True)