Rainfall Model#

import warnings
warnings.filterwarnings('ignore')

import os
import os.path as op
import sys

import matplotlib.pyplot as plt
import xarray as xr
import pandas as pd
import numpy as np
import geopandas as gpd
import utm
sys.path.insert(0, os.path.join(os.path.abspath(''), '..','..', '..' ))

from bluemath.rainfall_forecast.rainfall_forecast import *
from bluemath.rainfall_forecast.graffiti_rain import *

Warning: cannot import cf-json, install "metocean" dependencies for full functionality

# database
p_data = r'/media/administrador/HD2/SamoaTonga/data'
site, site_, nm = 'Upolu', 'upolu', 'up'
p_site = op.join(p_data, site)

p_deliv = op.join(p_site, 'd09_rainfall_forecast')
p_figs = op.join(p_deliv, 'figs')

#Basins path
p_basins = op.join(p_deliv, 'basins', nm + '_db' +  '_fixed.shp')

#Som fitting
p_som_fit = os.path.join(p_deliv, 'som_fit_30x30_4d.nc')

#Parameters from classification
p_params_fit = os.path.join(p_deliv, 'params_fit_30x30_4d_clust.nc')

#Historical track
n = 166
p_hist_track = os.path.join(p_deliv,'sp_{0}.nc'.format(n))
times_sel = [150, 385] #Times to cut that specific track to our area

#output files
p_rain_out = os.path.join(p_deliv,'rain_sp_{0}.nc'.format(n))

# make data folders
for f in [p_deliv, p_figs]:
    if not op.isdir(f):
        os.makedirs(f)

#Parameters

area_interest = [[187.848, 188.63],[-14.078, -13.756]] #upolu
discretization = 0.05 #local grid discretization

bas_raw = gpd.read_file(p_basins)

Regional Scale#

A model to estimate rainfall caused by the action of tropical cyclones (TCs) according to atmospheric, geographical, and oceanic variables has been developed. For this purpose, three databases have been used: precipitation measured by satellite GPM mission with half-hourly temporal resolution and 0.1° of spatial resolution, TCs historical tracks from IBTrACS and sea surface temperature (SST) from the NOAA OI SST V2 High Resolution Dataset with daily data at 0.25° spatial resolution.

The model assumes a circular shape of the TC influence area and defines the precipitation clustered in 900 groups following a Self-Organizing Maps (SOM) algorithm as a function of:

TC track minimum pressure
TC track sea surface temperature
TC track latitude
\(\Delta\)SST : difference in sea surface temperature between the current day and the previous day averaged in an circular area of radius equal to 200 km from the TC track
For each of the 900 clusters, the mean precipitacion over the radius dimension is calculated, and the curves are shown below.

from IPython.display import Image
Image('./resources/Rainfall_SOM.png', width=1000)    

Load needed files#

To run the model, we need to load the data corresponding to the SOM clustering (som_fit dataset) and the parameters associated to each of the 900 clusters (params_fit dataset).

som_fit = xr.open_dataset(p_som_fit)
params_fit = xr.open_dataset(p_params_fit)

Rainfall footprint#

For a historical TC track:

tco = xr.open_dataset(p_hist_track)
tc = process_tc_file(tco)

TC: AMOS, starting day 2016-04-17

#We can select the dates we are interested in
tc = tc.isel(date_time=slice(times_sel[0], times_sel[1]))

Plot TC track

The track is coloured according to the minimum pressure value

fig_tc = Plot_track(tc, figsize=[20,10], cmap='jet_r', m_size=10)

Obtain rainfall fields

interp=1 #interp: 1:activate / 0:deactivate 
interp_factor=1 #multiply by interp_factor the number of points along the track, in this case we already have the data each half hour in the TCs database
min_points=5 #minimum points to consider TC track and process it

tc_r = Generate_Rain_TC_4D(tc, som_fit, params_fit, interp, interp_factor, min_points)

time step 209

Once the model is constructed, from any TC track with its corresponding values of SST, latitude, minimum pressure and \(\Delta\)SST, a precipitation field from the closest cluster can be assigned to each position of the TC track.

Plot rain field

rmax = 20
figsize = [25,12]
fig_rf = plot_rain_fields(tc,tc_r,figsize,rmax)

Total storm rainfall#

Here we grid the results into a regular grid with the resolution given below (discretization), to obtain the resultant rainfall grid in mm/hr associated to each temporal moment.
Then we aggregate the results to obtain the resultant rainfall footprint for all the storm.

Obtain the grid for each temporal moment

discretization = 0.5
l_ds = tc_grid(tc_r,discretization)

time step 209

#Dataset containing for each temporal moment the resultant rainfall grid
dst = xr.concat(l_ds,dim='time')

#Dataset with the resultant rainfall footprint for all the TC duration
dss = dst.sum(dim='time')*0.5

dss.to_netcdf(p_rain_out)

By aggregating the rainfall at each point of the region, the rainfall SWATH for the duration of the storm can be obtained

fig_total_storm_tc = tc_total_rain(tc, dss, figsize=[25,12])

Local Scale#

Define area of interest and gris parameters

area_interest_l = [[area_interest[0][0]-discretization*3, area_interest[0][1]+discretization*3],[area_interest[1][0]-discretization*3, area_interest[1][1]+discretization*3]] #upolu

l_ds_area = tc_grid(tc_r,discretization, area=area_interest_l)
dst_area = xr.concat(l_ds_area,dim='time')
dss_area = dst_area.sum(dim='time')*0.5

time step 209

Basins discretization#

The discretization of the island in its watersheds have been carried out by means of QGIS, an open-source cross-platform desktop geographic information system (GIS) application that supports viewing, editing, and analysis of geospatial data. From the 5m resolution digital elevation model, a number of basins have been obtained, as in the following figure:

basins = Obtain_Basins(bas_raw, site)
Plot_Basins(basins)

Hyetograph for each basin#

The rainfall over time has been then obtained for each of the following basins, as it will be the input of the flooding model.

D_5 = Plot_Rain_Hyetograph_Basins(dst_area, dss_area, basins, INFO, site, thr_plot=1, thr=5)
D_5.to_netcdf(os.path.join(p_deliv, 'Rain_Real_Basins_' + site + '_Tc_' + str(tc.name.values) + '_sp' + str(n) + '.nc'))

Rainfall Model

Contents