STOFS3D-Pac scripts: changed climatology data and modified gen_source…

…sink.py to account for year change

src/Utility/Pre-Processing/STOFS-3D-Pac-shadow-VIMS/Source_sink/Clima/Dai_Trenberth/README

@@ -3,7 +3,7 @@ cp /sciclone/schism10/yinglong/NOAA-CSDL/Pacific_project_2021/RUN24e/Source/coas
 cp /sciclone/schism10/yinglong/NOAA-CSDL/Pacific_project_2021/RUN24e/Source/out_river .
 ln -sf /sciclone/schism10/yinglong/NOAA-CSDL/Pacific_project_2021/RUN24e/Source/hgrid.gr3 .
 
-Change "date_forecast" in the gen_clima.py accordingly, then run "python gen_clima.py" to get file:
+Run "python gen_clima.py" to get file:
 Dai_Trenberth_climatology_1990-2018_hourly.csv
 
 

src/Utility/Pre-Processing/STOFS-3D-Pac-shadow-VIMS/Source_sink/Clima/Dai_Trenberth/gen_clima.py

@@ -21,8 +21,8 @@ def read_station_file(station_file_name):
 
 if __name__ == '__main__':
 
-    #year = 2024
-    date_forecast = pd.date_range('2024-01-01', '2025-01-01', freq='MS')
+    #year = 2024 (value is independant with year. date will be changed based on NWM's date, gen_sourcesink.py)
+    date_forecast = pd.date_range('2024-01-01', '2025-02-01', freq='MS')
     river_filename = 'out_river'
     hgrid = 'hgrid.gr3'
 
@@ -52,14 +52,16 @@ def read_station_file(station_file_name):
         if mask.sum() > 0:
             month_mean.loc[mask == 1, i] = 0
 
-    #add January to the end
+    #add January and Feb to the end
     month_mean.loc[13] = month_mean.loc[1]
+    month_mean.loc[14] = month_mean.loc[2]
+
     #climatology monthly
-    df2023 = pd.DataFrame(data=month_mean.values, columns=stations, index=date_forecast)
-    #df_2023_daily = df2023.resample('D').asfreq()
-    #df_2023_daily.interpolate('linear', inplace=True)
-    df_2023_hourly = df2023.resample('H').asfreq()
-    df_2023_hourly.interpolate('linear', inplace=True)
+    df2024 = pd.DataFrame(data=month_mean.values, columns=stations, index=date_forecast)
+    #df_2024_daily = df2024.resample('D').asfreq()
+    #df_2024_daily.interpolate('linear', inplace=True)
+    df_2024_hourly = df2024.resample('H').asfreq()
+    df_2024_hourly.interpolate('linear', inplace=True)
 
     #to daily
     # tmp = month_mean.loc[df.index.month]
@@ -84,7 +86,7 @@ def read_station_file(station_file_name):
     #print(rank)
 
     #select rivers
-    df_clima_hourly = df_2023_hourly.loc[:, rank-1]
+    df_clima_hourly = df_2024_hourly.loc[:, rank-1]
 
     #Find nearest elem
     gd = read_schism_hgrid(hgrid)
@@ -104,8 +106,9 @@ def read_station_file(station_file_name):
     # ie_r = values[15]
     # idxs = np.where(ie_river == ie_r)[0]
     # cols = rank[idxs] - 1
-    # tmp = df_2023_hourly.loc[:, cols[0]]+df_2023_hourly.loc[:, cols[1]]+df_2023_hourly.loc[:, cols[2]]
+    # tmp = df_2024_hourly.loc[:, cols[0]]+df_2024_hourly.loc[:, cols[1]]+df_2024_hourly.loc[:, cols[2]]
     # tmp2 = out_flow_transposed.loc[:, ie_r]
 
 
     out_flow_transposed.to_csv('Dai_Trenberth_climatology_1990-2018_hourly.csv', index_label='date')
+    #out_flow_transposed.to_csv('Dai_Trenberth_climatology_1990-2018_hourly.csv', index_label='date', mode='a', header=False)

src/Utility/Pre-Processing/STOFS-3D-Pac-shadow-VIMS/Source_sink/Clima/Dai_Trenberth/hgrid.gr3

@@ -0,0 +1 @@
+/sciclone/schism10/yinglong/NOAA-CSDL/Pacific_project_2021/RUN24e/Source/hgrid.gr3

src/Utility/Pre-Processing/STOFS-3D-Pac-shadow-VIMS/Source_sink/NWM/Combined/gen_sourcesink.py

@@ -8,6 +8,12 @@
 import pandas as pd
 from netCDF4 import Dataset
 
+def is_leap_year(year):
+    if (year % 4 == 0 and year % 100 != 0) or (year % 400 == 0):
+        return True
+    else:
+        return False
+
 def read_featureID_file(filename):
 
     with open(filename) as f:
@@ -216,8 +222,16 @@ def streamflow_lookup(file, indexes, threshold=-1e-5):
     print(sinks.shape)
 
     #combine with Dai_Trenberth_climatology
-    df = pd.read_csv(clima_filename)
-    df.set_index('date', inplace=True)
+    df = pd.read_csv(clima_filename, parse_dates=['date'], index_col='date')
+
+    if is_leap_year(int(startdate.year)):
+        print(f'Year {startdate.year} is a leap year')
+    else:
+        print('Year {startdate.year} is not a leap year, drop day 02-29')
+        df = df[~(df.index.date == pd.to_datetime(f'2024-02-29').date())]
+    date_forecast = pd.date_range(f'{startdate.year}-01-01', f'{int(startdate.year) + 1}-02-01', freq='H')
+    df.index = date_forecast
+
     mask = (pd.to_datetime(df.index) >= dates[0]) & (pd.to_datetime(df.index) <=dates[-1])
     df_forecast = df[mask]
     df_nwm = pd.DataFrame(data=sources, columns=np.array(eid_sources), index=np.array(dates))