arrival_time_percentile.py

#!/usr/bin/env python3
#  -*- coding: utf-8 -*-

"""
Compute and plot a distribution of oil spill arrival time at sensistive sites in Qatar based on backtracking simulations with OpenDrift

Author: Thomas Dobbelaere, Earth and Life Institute, UCLouvain, Belgium
Last modified: 19 July 2022
"""

import numpy as np
import os
import matplotlib.pyplot as plt
from map_data import plot_landmark
from postpro import ArrivalTimeMap
from typing import List
import itertools

minlon, maxlon, minlat, maxlat = 50, 53, 24, 27.5
res = 0.01
basedir = "/export/miro/students/tanselain/OpenOil/"
sources = ['ras_laffan','ras_laffan_port', 'abu_fontas', 'umm_al_houl']  
years = ['2017', '2018', '2019', '2020']

def plot_arrival_time(sources: List[str], years: List[str], month: str, outfile: str) -> None:

    atm = ArrivalTimeMap(minlon, maxlon, minlat, maxlat, res)
    for src,yr in itertools.product(sources,years):
        if src == "lusail":
            nc_path = basedir + f'/Output_backward_lusail/{yr}/{month}'
        else:
            nc_path = basedir + f'/Output_backward/{yr}/{src}/{month}'
        for day in range(1,32):
            ncfile = nc_path+f'/out_{month}_{day}.nc'
            if not os.path.isfile(ncfile):  
                continue
            atm.add_arrival_times_from_file(ncfile)

    arrival_time_days = atm.get_arrival_time_days()
    hasdata = ~np.isnan(arrival_time_days)
    arrival_time_days[hasdata] = np.floor(arrival_time_days[hasdata])

    _, ax = plt.subplots(figsize=(12,12))
    ax.pcolormesh(atm.lon, atm.lat, arrival_time_days, cmap="YlOrRd_r")
    cmap = plt.get_cmap('YlOrRd_r')
    colors = cmap(np.arange(5)/5)
    for i,clr in enumerate(colors):
        label = "< 1 day" if i == 0 else f"{i} - {i+1} days"
        ax.add_patch(plt.Rectangle((0,0), 0, 0, color=clr, label=label))
    plot_landmark(ax, sources)
    month_name = month[0].upper()+month[1:]
    ax.text(0, 0, " "+month_name+" ", va="bottom", ha="left", color="k", \
            fontsize=20, transform=ax.transAxes)
    ax.legend( loc="upper right", title="Arrival time\n(5th percentile)", \
                fontsize=14, title_fontsize=14)
    plt.savefig(outfile, dpi=300, bbox_inches="tight")
    plt.close()

if __name__ == "__main__":
    for month in ['april', 'june']:
        outfile = f"fig_arrival_{month}.pdf"
        plot_arrival_time(sources, years, month, outfile)