Merge branch 'main' of https://github.com/rodekruis/river-flood-data-…

…analysis

GloFAS/GloFAS_analysis/probability_calculator.py

@@ -11,7 +11,7 @@
 import GloFAS.GloFAS_prep.configuration as cfg
 
 class FloodProbabilityProcessor:
-    def __init__(self, leadtime, adminLevel, forecastType, measure='max', start_date=None, end_date=None, nrCores=4):
+    def __init__(self, leadtime, adminLevel, forecastType, measure='max', start_date=None, end_date=None, nrCores=4, comparisonShape='polygon'):
         # start and end_date are only necessary if your forecastType is not 'reforecast' but 'forecast'
         self.leadtime = leadtime
         self.adminLevel = adminLevel
@@ -22,6 +22,7 @@ def __init__(self, leadtime, adminLevel, forecastType, measure='max', start_date
         self.area = cfg.MaliArea
         self.lakesPath = cfg.lakesPath
         self.forecastType = forecastType
+        self.comparisonShape = comparisonShape
         if forecastType == 'reforecast':
             self.MONTHSDAYS = get_monthsdays()
         elif forecastType ==  'forecast':
@@ -32,11 +33,18 @@ def __init__(self, leadtime, adminLevel, forecastType, measure='max', start_date
         self.threshold_da = openThreshold(self.DataDir, self.crs, self.RPyr, self.area, self.reference_Q_da)
         self.threshold_gdf = aggregation(self.threshold_da, adminPath, 'polygon', measure=cfg.measure)
         self.nrCores = nrCores
-
+        
     def exceedance(self, Q_da, threshold_gdf, nrEnsemble):
-        '''Check exceedance of threshold values for a single ensemble member'''
-        GloFAS_gdf = aggregation(Q_da, threshold_gdf, 'polygon', nrEnsemble, measure=self.measure)
-        exceedance = GloFAS_gdf[f'{self.measure}'] > threshold_gdf[f'{self.measure}']
+        '''Check exceedance of threshold values for a single ensemble member
+        threshold_gdf should be in same shape as eventual outcme, so if comparisonShape=polygon, threshold_gdf should be in polygon'''
+        if self.comparisonShape=='polygon':
+            GloFAS_gdf = aggregation(Q_da, threshold_gdf, 'polygon', nrEnsemble, measure=self.measure)
+            exceedance = GloFAS_gdf[f'{self.measure}'] > threshold_gdf[f'{self.measure}']
+        elif self.comparisonShape=='point': # assuming this then refers to aggregation to a point
+            GloFAS_gdf = aggregation (Q_da, threshold_gdf, 'point', nrEnsemble)
+            exceedance = GloFAS_gdf ['rastervalue'] > threshold_gdf['rastervalue']
+        else: 
+            raise ValueError (f"Not a valid shape : {self.comparisonShape}, pick between 'polygon' and 'point'")
         return exceedance
 
     def probability(self, Q_da, threshold_gdf, nrCores):
@@ -92,6 +100,7 @@ def concatenate_prob(self, Q_da, nrCores):
         floodProb_gdf_concat = gpd.GeoDataFrame(floodProb_gdf_concat, geometry='geometry')
 
         # Save to file
-        output_path = f"{self.DataDir}/floodedRP{self.RPyr}yr_leadtime{self.leadtime}_ADM{self.adminLevel}.gpkg"
+        # make the save to file so that 
+        output_path = f"{self.DataDir}/floodedRP{self.RPyr}yr_leadtime{self.leadtime}.gpkg"
         floodProb_gdf_concat.to_file(output_path, driver="GPKG")
         return floodProb_gdf_concat

GloFAS/GloFAS_prep/Q_timeseries_station.py

@@ -0,0 +1,53 @@
+# creates a timeseries of the discharge (or no timeseries but just static value for discharge threshold), 
+# for the location corresponding to a certain station (DNH station)
+import numpy as np
+import pandas as pd
+import geopandas as gpd
+from datetime import datetime, timedelta
+from GloFAS.GloFAS_prep.aggregation import aggregation
+from GloFAS.GloFAS_data_extractor.reforecast_dataFetch import get_monthsdays
+from GloFAS.GloFAS_data_extractor.forecast_dataFetch import compute_dates_range
+from GloFAS.GloFAS_data_extractor.threshold_opener import openThreshold
+from GloFAS.GloFAS_data_extractor.unzip_open import openGloFAS, unzipGloFAS
+import GloFAS.GloFAS_prep.configuration as cfg
+
+def Q_over_time(Q_da, forecastType):
+    '''Concatenate exceedance probability calculations over multiple time steps'''
+    all_years_data = []
+    if forecastType=='reforecast':
+        for month, day in self.MONTHSDAYS:
+            print(f"Starting for {month}, {day}")
+            Q_da_rasterpath = unzipGloFAS (self.DataDir, self.leadtime, month, day)
+            Q_da = openGloFAS(Q_da_rasterpath, self.leadtime)  # Open data for the specific month/day
+            for timestamp in Q_da.time:  # Iterate over years
+                startLoop = datetime.now()
+                valid_timestamp = pd.to_datetime(str(timestamp.values)) + timedelta(hours=self.leadtime)
+                # Calculate probability
+                flooded_df = self.probability(Q_da, self.threshold_gdf, nrCores)
+                flooded_df['ValidTime'] = valid_timestamp  # Add timestamp
+                all_years_data.append(flooded_df)
+                print(f"Time step {valid_timestamp} done, took: {datetime.now() - startLoop}")
+
+    if self.forecastType=='forecast': 
+        for date in self.MONTHSDAYSYEARS:
+            startLoop = datetime.now()
+            year  = date.strftime('%Y')
+            month = date.strftime('%m')
+            day   = date.strftime('%d')
+            Q_da_rasterpath = unzipGloFAS(self.DataDir, self.leadtime, month, day, year)
+            Q_da = openGloFAS(Q_da_rasterpath, lakesPath=self.lakesPath, crs=self.crs)
+            valid_timestamp = pd.to_datetime(str(timestamp.values)) + timedelta(hours=self.leadtime)
+            # Calculate probability
+            flooded_df = self.probability(Q_da, self.threshold_gdf, nrCores)
+            flooded_df['ValidTime'] = valid_timestamp  # Add timestamp
+            all_years_data.append(flooded_df)
+            print(f"Time step {valid_timestamp} done, took: {datetime.now() - startLoop}")
+    # Concatenate results for all time steps
+    floodProb_gdf_concat = pd.concat(all_years_data, ignore_index=True)
+    floodProb_gdf_concat = gpd.GeoDataFrame(floodProb_gdf_concat, geometry='geometry')
+
+    # Save to file
+    # make the save to file so that 
+    output_path = f"{self.DataDir}/floodedRP{self.RPyr}yr_leadtime{self.leadtime}.gpkg"
+    floodProb_gdf_concat.to_file(output_path, driver="GPKG")
+    return floodProb_gdf_concat

GloFAS/GloFAS_prep/station_to_glofas.py

@@ -25,7 +25,7 @@
         ups_area_point_m=row['Catchment area (km2)'] * 1e6,  # Convert from km² to m²
         stationName=row['Station names'],
         glofas_ups_area_xr=glofas_ups_area_da,
-        radius_m=50000
+        radius_m=100000
     )
     stations_upsarea_df.loc[i, 'Glofas_Point_X'] = glofas_x
     stations_upsarea_df.loc[i, 'Glofas_Point_Y'] = glofas_y

comparison/pointMatching.py

@@ -164,7 +164,7 @@ def find_corresponding_point_within_box(station_lon, station_lat, ups_area_point
         linewidth=1.5,  # Border thickness
         label=f"GloFAS location corresponding to minimal upstream area difference "
         )
-    plt.plot ([],[], label= f'found within {radius_m/1e3:.0f} km radius (difference={best_area_diff/1e6:.1f} km²)')
+    plt.plot ([],[], ' ', label= f'found within {radius_m/1e3:.0f} km radius (difference={best_area_diff/1e6:.1f} km²)')
     # Add the colorbar
     # sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
     # sm.set_array([])  # Required for ScalarMappable