Correct radius calculation based in latitude

src/lisfloodutilities/gridding/lib/filters.py

@@ -5,7 +5,8 @@
 from datetime import datetime as dt
 from typing import List, Tuple
 from scipy.spatial import cKDTree
-import geopy.distance
+from geopy.distance import geodesic
+import math
 
 
 class KiwisFilter:
@@ -195,6 +196,8 @@ class ObservationsKiwisFilter(KiwisFilter):
     """
 
     CLUSTER_COLLAPSE_RADIUS = np.float32(0.011582073434000193) # decimal degrees (1287 m)
+    EARTH_RADIUS_IN_METERS = 6371000
+    EARTH_RADIUS_IN_KILOMETERS = 6371
 
     def __init__(self, filter_columns: dict = {}, filter_args: dict = {}, var_code: str = '', quiet_mode: bool = False):
         super().__init__(filter_columns, filter_args, var_code, quiet_mode)
@@ -203,14 +206,41 @@ def __init__(self, filter_columns: dict = {}, filter_args: dict = {}, var_code:
         self.provider_radius = {}
         for provider_id in self.args:
             radius_km = self.args[provider_id]
-            radius = self.kilometers2degrees(radius_km)
-            self.provider_radius[provider_id] = radius
+            # radius = self.kilometers2degrees(radius_km)
+            self.provider_radius[provider_id] = radius_km
 
     @staticmethod
-    def kilometers2degrees(km: np.float32) -> np.float32:
+    def kilometers2degrees_approximated(radius_km: np.float32) -> np.float32:
         # Convert km to degrees of latitude
-        delta_lat = km * np.float32(0.00899928005)
+        delta_lat = radius_km * np.float32(0.00899928005)
         return delta_lat
+
+    @staticmethod
+    def kilometers2degrees(radius_km: np.float32) -> np.float32:
+        # Convert the radius from km to radians (using the Haversine formula)
+        radius_rad = geodesic(kilometers=radius_km).meters / ObservationsKiwisFilter.EARTH_RADIUS_IN_METERS
+        return np.float32(radius_rad)
+
+    @staticmethod
+    def calculate_radius(lat: np.float32, lon: np.float32, radius_km: np.float32) -> Tuple[np.float32, np.float32]:
+        """
+        Calculate the radius around a point in decimal degrees.
+        
+        Args:
+            lat (float): Latitude of the point in decimal degrees.
+            lon (float): Longitude of the point in decimal degrees.
+            km (float): Radius in kilometers.
+        
+        Returns:
+            tuple: Radius in decimal degrees for latitude and longitude.
+        """
+        # Convert kilometers to radians
+        km_to_rad = radius_km / ObservationsKiwisFilter.EARTH_RADIUS_IN_KILOMETERS
+        # Calculate the radius in decimal degrees for latitude
+        lat_degrees = km_to_rad * (180 / math.pi)
+        # Calculate the radius in decimal degrees for longitude at the given latitude
+        lon_degrees = km_to_rad * (180 / math.pi) / math.cos(math.radians(lat))
+        return abs(lat_degrees), abs(lon_degrees)
 
     def apply_filter(self, df: pd.DataFrame) -> pd.DataFrame:
         df = super().apply_filter(df)
@@ -230,7 +260,10 @@ def has_neighbor_within_radius_from_other_providers(self, row: pd.Series, tree:
         cur_provider_id = row[self.COL_PROVIDER_ID]
         if cur_provider_id == provider_id:
             location = (row[self.COL_LON], row[self.COL_LAT])
-            nearby_points = tree.query_ball_point(location, radius)
+            radius_lat_degrees, radius_lon_degrees = self.calculate_radius(lat=np.float32(row[self.COL_LAT]),
+                                                                           lon=np.float32(row[self.COL_LON]),
+                                                                           radius_km=radius)
+            nearby_points = tree.query_ball_point(location, radius_lon_degrees)
             return len(nearby_points) > 0
         return False
 
@@ -424,7 +457,10 @@ def get_decumulated_24h_value_for_missing_6h_values(self, row: pd.Series, tree:
         decumulated_value = cur_24h_value
         if cur_provider_id == provider_id:
             location = (row[self.COL_LON], row[self.COL_LAT])
-            nearby_points = tree.query_ball_point(location, radius)
+            radius_lat_degrees, radius_lon_degrees = self.calculate_radius(lat=np.float32(row[self.COL_LAT]),
+                                                                           lon=np.float32(row[self.COL_LON]),
+                                                                           radius_km=radius)
+            nearby_points = tree.query_ball_point(location, radius_lon_degrees)
             stations_6h = stations_6h_df.loc[stations_6h_df.index.isin(nearby_points)][['sum_6h_values', 'count_6h_slots']]
             if stations_6h.empty:
                 return decumulated_value

src/lisfloodutilities/gridding/tools/get_stats_by_maskarea_and_pixelarea.py

@@ -0,0 +1,149 @@
+import netCDF4 as nc
+import numpy as np
+from pathlib import Path
+
+OUTPUT_PATH = '/mnt/nahaUsers/gomesgo/ERA5_vs_EFAS/var_statistics2'
+
+data = {
+    # 'EFAS': ('/mnt/nahaUsers/grimast/meteo_template/meteo_from_francesca/EMO-1arcmin-{v}_{yyyy}_monsum.nc', ['pr6', 'e0', 'et'], ['2019', '2020', '2022'], 0, 12),
+    # 'EFAS_FRANCESCA_ATOS_2023': ('/mnt/nahaUsers/grimast/meteo_template/meteo_from_francesca/{v}_{yyyy}_from_francesca_and_ATOS_monsum.nc', ['pr', 'e0', 'et'], ['2023'], 0, 7),
+    # 'EFAS_FRANCESCA_ATOS_2024': ('/mnt/nahaUsers/grimast/meteo_template/meteo_from_francesca/{v}_{yyyy}_from_francesca_and_ATOS_monsum.nc', ['pr', 'e0', 'et'], ['2024'], 0, 3),
+    # 'EFAS_KISTERS_REFRESH_2023': ('/mnt/nahaUsers/grimast/meteo_template/meteo_from_francesca/{v}_EFASv5_Kisters_refresh_202307_202403_mm_per_day_monsum.nc', ['pr6'], ['2023'], 0, 6),
+    # 'EFAS_KISTERS_REFRESH_2024': ('/mnt/nahaUsers/grimast/meteo_template/meteo_from_francesca/{v}_EFASv5_Kisters_refresh_202307_202403_mm_per_day_monsum.nc', ['pr6'], ['2024'], 6, 9),
+    'EFAS_KISTERS_REFRESH_202404': ('/mnt/nahaUsers/grimast/meteo_template/meteo_from_francesca/{v}_EFASv5_Kisters_refresh_202404_mm_per_day_monsum.nc', ['pr6'], ['2024'], 0, 1),
+    }
+
+pixarea_filename = '/mnt/nahaUsers/gomesgo/ERA5_vs_EFAS/pixarea.nc'
+flip_pixarea_updown = False
+# pixarea_filename = '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/PixAREA_01degrees.tif.nc'
+# flip_pixarea_updown = True
+
+pixarea_var = 'Band1'
+
+
+mask_var = 'Band1'
+
+
+
+mask_filenames = [
+    '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/mask_Elbe_1arcmin.tif.nc',
+    '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/mask_Vistula_1arcmin.tif.nc',
+    '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/mask_Europe_commonareas_1arcmin.tif.nc',
+    '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/Po_mask_1arcmin_LARGE.tif.nc',
+    '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/mask_Oder_1arcmin.tif.nc',
+    '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/mask_EmiliaRomagna_1arcmin.tif.nc',
+    '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/mask_Piemonte_1arcmin.tif.nc',
+    '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/mask_Lombardia_1arcmin.tif.nc',
+    ]
+
+
+# mask_filenames = [
+#     '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/mask_Elbe_01degrees.tif.nc',
+#     '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/mask_Vistula_01degrees.tif.nc',
+#     '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/mask_Europe_commonareas_01degrees.tif.nc',
+#     '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/mask_Oder_01degrees.tif.nc',
+#     '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/test_maskmap_PoValley_01degrees.tif.nc',
+#     '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/mask_EmiliaRomagna_01degrees.tif.nc',
+#     '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/mask_Piemonte_01degrees.tif.nc',
+#     '/mnt/nahaUsers/grimast/meteo_template/Goncalo_monthly_statistics/masks/mask_Lombardia_01degrees.tif.nc',
+#     ]
+
+
+for mask_filename in mask_filenames:
+    f=Path(mask_filename)
+    output_basename=f.stem
+
+    # Open the mask NetCDF file
+    mask_nc = nc.Dataset(mask_filename, 'r')
+    # Read the mask data
+    mask_data = mask_nc.variables[mask_var][:]
+    mask_data = np.flipud(mask_data)
+
+    # Open the Pixel Area NetCDF file
+    pixarea_nc = nc.Dataset(pixarea_filename, 'r')
+    # Read the pixel area data
+    pixarea_data = pixarea_nc.variables[pixarea_var][:]
+    if flip_pixarea_updown:
+        pixarea_data = np.flipud(pixarea_data)
+    masked_pixarea = np.ma.masked_where(mask_data !=1, pixarea_data)
+
+    for dataset_name in data:
+        data_filename_pattern, vars_array, years, start_step, end_step = data[dataset_name]
+
+        for v in vars_array:
+            output_file_path = Path(OUTPUT_PATH, f'{v}_{output_basename}.csv')
+            with open(output_file_path, 'a') as out_file:
+                for yyyy in years:
+                    data_filename = data_filename_pattern.format(yyyy=yyyy, v=v)
+
+                    # Open the data NetCDF file
+                    data_nc = nc.Dataset(data_filename, 'r')
+
+                    # Assuming the variable containing the monthly data is named 'monthly_data'
+                    # and the mask variable is named 'mask'
+                    monthly_data_var = v
+
+                    # data_lat = data_nc.variables['lat'][:]
+                    # data_lon = data_nc.variables['lon'][:]
+                    # mask_lat = mask_nc.variables['lat'][:]
+                    # mask_lon = mask_nc.variables['lon'][:]
+            #        print('data_lat', data_lat)
+            #        print('data_lon', data_lon)
+            #        print('mask_lat', mask_lat)
+            #        print('mask_lon', mask_lon)
+                    # pixarea_lat = pixarea_nc.variables['lat'][:]
+                    # pixarea_lon = pixarea_nc.variables['lon'][:]
+            #        print('pixarea_lat', pixarea_lat)
+            #        print('pixarea_lon', pixarea_lon)
+
+
+                    # Initialize an empty list to store the yearly totals
+                    yearly_totals = []
+
+                    # Loop through each timeslice (assuming there are 12 for each year)
+                    for i in range(start_step, end_step):
+                        # Read the data for the current timeslice
+                        monthly_data = data_nc.variables[monthly_data_var][i, :, :]
+
+                        # Apply the mask to the data subset
+                        # masked_monthly_data = np.where(mask_data, data_subset, 0)
+                        masked_monthly_data = np.ma.masked_where(mask_data !=1, monthly_data)
+
+                        monthly_sum_m = masked_monthly_data / 1000 / 4
+
+                        monthly_sum_m3 = monthly_sum_m * masked_pixarea
+
+                        monthly_totals_m3 = np.sum(monthly_sum_m3)
+                        # monthly_totals_m3 = np.mean(masked_monthly_data)
+
+                        # print(f"Total values of {v} yyyy-mm {yyyy}-{i}: {monthly_totals_m3}")
+                        # out_file.write(f"{dataset_name};{yyyy};{i+1};{monthly_totals_m3}\n")
+
+                        out_file.write(f"EFAS;{yyyy};{i+1+3};{monthly_totals_m3}\n")
+                        # if yyyy == '2023':
+                        #     out_file.write(f"EFAS;{yyyy};{i+1+6};{monthly_totals_m3}\n")
+                        # elif yyyy == '2024':
+                        #     out_file.write(f"EFAS;{yyyy};{i+1-6};{monthly_totals_m3}\n")
+
+                        # If it's the first month, initialize the yearly sum with the masked data
+                        if i == 0:
+                            yearly_sum = monthly_totals_m3
+                        else:
+                            # Otherwise, add the masked data to the running yearly sum
+                            yearly_sum += monthly_totals_m3
+
+                    # Add the yearly sum to the list of yearly totals
+                    yearly_totals.append(yearly_sum)
+
+                    # Convert the list to a numpy array if necessary
+                    yearly_totals = np.array(yearly_totals)
+
+                    total_values = np.sum(yearly_totals)
+
+                    # Close the NetCDF files
+                    data_nc.close()
+
+                    # print(f"Total values of {v} year {yyyy}: {total_values}")
+
+    mask_nc.close()
+    pixarea_nc.close()