diff --git a/arrival_time_routes.py b/arrival_time_routes.py index 0bcefe493c9af1062c5768b93c58b4d46ce93b16..dbb3907d5d1a58a9254575e6cd6b7a7c39041dc9 100644 --- a/arrival_time_routes.py +++ b/arrival_time_routes.py @@ -5,53 +5,75 @@ Intersects arrival times with shipping routes Author: Thomas Dobbelaere, Earth and Life Institute, UCLouvain, Belgium -Last modified: 19 July 2022 +Last modified: 9 November 2022 """ import numpy as np +import os import geopandas as gpd +import shapely +from postpro import ArrivalTimeMap import calendar -import geopandas -import common +from typing import List -basedir = "/export/miro/students/tanselain/OpenOil" +# --- 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' -def compute_arrival_routes(sources, routes, year, month, outfile): +# --- USEFUL FUNCTIONS --- # - sources = common.to_iterable(sources) - routes = common.to_iterable(routes) +def compute_arrival_time( + df: gpd.GeoDataFrame, sources: List[str], year: int, month: str, day: int +) -> 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}' + 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)[:] - time_tag = f"{year}{month:02d}" - # read arrival time - shpfile = f"shp/{time_tag}/arrival_time_percentile_{time_tag}.shp" - df_arrival = gpd.read_file(shpfile) - df_arrival = df_arrival[~np.isnan(df_arrival['arrival_ti'])] - # read shipping routes - path = basedir + f'/Indicators/Shipping_selected/Selected_shipping_route_2019' - df_routes = gpd.read_file(path) +def main() -> None: + + # initialize geodataframes + df_routes = gpd.read_file(path_routes) df_routes = df_routes.to_crs("EPSG:4326") - df_arrival.crs = df_routes.crs - # intersection - df_inter = gpd.sjoin(df_arrival, df_routes, op='intersects') - df_inter['dn_total']= df_inter['arrival_ti'] * df_inter['int'] - - res = [] - for route in routes: - df_inter_route = df_inter.loc[df_inter['Route'] == route] - mean_arrival_time = df_inter_route['arrival_ti'].mean() - print(route, mean_arrival_time) - res.append(mean_arrival_time) - return res - -if __name__ == "__main__": - sources = ['ras_laffan', 'abu_fontas', 'umm_al_houl', 'ras_laffan_port'] - routes = ['Barhain', 'Doha', 'Iran', 'Ras_laffan_field', 'Ras_laffan'] - fname = "arrival_routes_201601.png" - a = np.empty((len(routes),12,5)) - for year in range(2016,2021): - for month in range(1,13): - print(calendar.month_name[month], year) - iyr = year-2016 - imn = month-1 - a[:,imn,iyr] = compute_arrival_routes(sources, routes, year, month, fname) - np.save("arrival_by_route.npy", {"routes":routes, "data":a}) + 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("atime_by_routes.csv","w") as csv: + csv.write("year,month,day,route,arrival_time\n") + for yr in range(2016,2021): + for mn in range(1,13): + month = calendar.month_name[mn].lower() + print(f'{month} {yr}') + for d in range(1,32): + compute_arrival_time(df, sources, yr, month, d) + 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},{month},{d},{route},{v}\n") + +if __name__ == '__main__': + main() \ No newline at end of file