update script intersecting arrival time with routes

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