arrival_time_routes.py

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

"""
Intersects arrival times with shipping routes

Author: Thomas Dobbelaere, Earth and Life Institute, UCLouvain, Belgium
Last modified: 9 November 2022
"""

import numpy as np
import os
import geopandas as gpd
import shapely
from postpro import ArrivalTimeMap
import calendar
from typing import List
import itertools

# --- USEFUL VARIABLES --- #
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']  
routes = ['Barhain', 'Doha', 'Iran', 'Ras_laffan_field', 'Ras_laffan']
path_routes = basedir + f'/Indicators/Shipping_selected/Selected_shipping_route_2019'
route_names = {
    "Barhain" : "Bahrain",
    "Iran" : "Qatar-Iran",
    "Ras_laffan_field": 'Ras Laffan to North Field',
    "Ras_laffan": "Ras Laffan LNG export",
    "Doha" : 'Doha-Mesaieed'
}


# --- USEFUL FUNCTIONS --- #

def compute_arrival_time(
   df: gpd.GeoDataFrame, sources: List[str], year: int, month: str
) -> None:
    atm = ArrivalTimeMap(minlon, maxlon, minlat, maxlat, res)
    for src in sources:
        if src == "lusail":
            nc_path = basedir + f'Output_backward_lusail/{year}/{month}'
        else:
            nc_path = basedir + f'Output_backward/{year}/{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)
    a = atm.get_arrival_time_days()
    df.loc[:,'atime_days'] = a.reshape(-1)[:]

def main() -> None:

    # initialize geodataframes
    df_routes = gpd.read_file(path_routes)
    df_routes = df_routes.to_crs("EPSG:4326")
    data = []
    lons = np.arange(minlon, maxlon+res, res)
    lats = np.arange(minlat, maxlat+res, res)
    for i in range(lats.size-1):
        for j in range(lons.size-1):
            px = lons[[j,   j, j+1, j+1]]
            py = lats[[i, i+1, i+1,   i]]
            poly = shapely.geometry.Polygon(np.column_stack([px,py]))
            data.append((0.0, poly))
    df = gpd.GeoDataFrame(data, columns=['atime_days','geometry'], crs=df_routes.crs)
    with open("arrival_time_route.csv","w") as csv:
        csv.write("year,month,route,arrival_time\n")
        years = np.arange(2016,2021)
        months = np.arange(1,13)
        for yr,mn in itertools.product(years,months):
            month = calendar.month_name[mn].lower()
            print(f"{month} {yr}")
            compute_arrival_time(df, sources, yr, month)
            df_arrival = df[~np.isnan(df['atime_days'])]
            if df_arrival.shape[0] == 0:
                continue
            df_inter = gpd.sjoin(df_arrival, df_routes, op='intersects')
            df_inter['dn_total']= df_inter['atime_days'] * df_inter['int']
            for route in routes:
                df_inter_route = df_inter.loc[df_inter['Route'] == route]
                v = df_inter_route['atime_days'].mean()
                csv.write(f"{yr},{calendar.month_abbr[mn]},{route_names[route]},{v}\n")
                
if __name__ == '__main__':
    main()