Example on 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/test_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)

/opt/conda/envs/swirlspy/lib/python3.6/site-packages/scipy/interpolate/rbf.py:257: LinAlgWarning: Ill-conditioned matrix (rcond=2.85228e-19): result may not be accurate.
  self.nodes = linalg.solve(self.A, self.di)

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)