RG filter and making netcdf and tiff with same rounding results

src/lisfloodutilities/gridding/configuration/1arcmin/config_rg.txt

@@ -12,6 +12,8 @@ DATA_TYPE_PACKED = i2
 STANDARD_NAME = integral_wrt_time_of_surface_downwelling_shortwave_flux_in_air
 LONG_NAME = Daily Calculated Radiation
 
+KIWIS_FILTER_PLUGIN_CLASSES = {'SolarRadiationLimitsKiwisFilter': {'EXCLUDE_BELLOW_LATITUDE': 66.0, 'EXCLUDE_BELLOW_VALUE': 0.0}}
+
 [VAR_TIME]
 
 UNIT = days since 1990-01-02 00:00:00.0

src/lisfloodutilities/gridding/lib/filters.py

@@ -286,6 +286,37 @@ def set_filtered_stations(self, df: pd.DataFrame):
             self.filtered_station_ids[station_id] = ''
 
 
+class SolarRadiationLimitsKiwisFilter(KiwisFilter):
+    """
+    Class to filter Solar Radiation Kiwis files whose data coordinates are both less equal a defined latitude
+    and values less equal a defined threshold.
+    This was developed to avoid wrong values of zero daily solar radiation in EFAS domain bellow 66 degrees Latitude (empirical).
+    
+    Expects to have in filter_args a dictionary containing the definition of the limits using the
+    keys EXCLUDE_BELLOW_LATITUDE and EXCLUDE_BELLOW_VALUE
+    """
+
+    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)
+        # Calculating the radius in decimal degrees
+        self.threshold_max_latitude = 72.0
+        if 'EXCLUDE_BELLOW_LATITUDE' in self.args:
+            self.threshold_max_latitude = self.args['EXCLUDE_BELLOW_LATITUDE']
+        self.threshold_min_value = 0.0
+        if 'EXCLUDE_BELLOW_VALUE' in self.args:
+            self.threshold_min_value = self.args['EXCLUDE_BELLOW_VALUE']
+
+    def apply_filter(self, df: pd.DataFrame) -> pd.DataFrame:
+        df = super().apply_filter(df)
+        # Convert to float so it can be compared to the thresholds
+        df[self.COL_LAT] = df[self.COL_LAT].astype(float)
+        df[self.COL_VALUE] = df[self.COL_VALUE].astype(float)
+        # Filter values
+        df = df[~((df[self.COL_LAT] <= self.threshold_max_latitude) & (df[self.COL_VALUE] <= self.threshold_min_value))]
+        self.print_statistics(df)
+        return df
+
+
 class DowgradedDailyTo6HourlyObservationsKiwisFilter(ObservationsKiwisFilter):
     """
     Class to filter Kiwis files metadata for stations whose daily data was down graded to 6hourly data

src/lisfloodutilities/gridding/lib/utils.py

@@ -355,8 +355,8 @@ def do_height_correction(self) -> bool:
         return self.height_correction_factor != 0.0
 
     @property
-    def height_correction_factor(self) -> float:
-        return float(self.get_config_field('PROPERTIES', 'HEIGHT_CORRECTION_FACTOR'))
+    def height_correction_factor(self) -> np.float32:
+        return np.float32(self.get_config_field('PROPERTIES', 'HEIGHT_CORRECTION_FACTOR'))
 
     @property
     def truncate_negative_values(self) -> bool:
@@ -402,8 +402,7 @@ def __init__(self, conf: Config, quiet_mode: bool = False, use_broadcasting: boo
         super().__init__(quiet_mode)
         self.conf = conf
         self.use_broadcasting = use_broadcasting
-        self.unit_conversion = float(self.conf.get_config_field('PROPERTIES', 'UNIT_CONVERSION'))
-        # self.compressed_NaN = (self.conf.VALUE_NAN - self.conf.add_offset) / self.conf.scale_factor
+        self.unit_conversion = np.float32(self.conf.get_config_field('PROPERTIES', 'UNIT_CONVERSION'))
 
     def correct_height(self, df: pd.DataFrame) -> pd.DataFrame:
         if self.conf.do_height_correction:
@@ -432,16 +431,17 @@ def reset_height(self, result: ma.MaskedArray) -> np.ndarray:
         return result
 
     def prepare_grid(self, result: np.ndarray, grid_shape: np.ndarray) -> np.ndarray:
-        # Compress data
+        # Pack data
         if self.unit_conversion != 1.0:
             result = result * self.unit_conversion
         result[np.where(result == self.conf.VALUE_NAN)] = np.nan
         result[np.where(result < self.conf.value_min)] = np.nan
         result[np.where(result > self.conf.value_max)] = np.nan
-        result = np.round(result.astype(float), 1)
-        result = (result - self.conf.add_offset) / self.conf.scale_factor
+        result = result.astype(np.float32)
+        result = np.round(result, 1)
+        result[~np.isnan(result)] -= self.conf.add_offset
+        result[~np.isnan(result)] /= self.conf.scale_factor
         result[np.isnan(result)] = self.conf.VALUE_NAN
-        # Reshape grid
         grid_data = result.reshape(grid_shape)
         return grid_data
 
@@ -455,15 +455,16 @@ def check_grid_nan(self, filename: Path, result: np.ndarray):
         current_grid[np.where(current_grid > self.conf.value_max)] = np.nan
         count_grid_nan = np.count_nonzero(np.isnan(current_grid))
         if count_dem_nan != count_grid_nan:
-            print(f"WARNING: The grid interpolated from file {filename.name} contains different NaN values ({count_grid_nan}) than the DEM ({count_dem_nan}). diff: {count_grid_nan - count_dem_nan}")
+            self.print_msg(f"WARNING: The grid interpolated from file {filename.name} contains different NaN values ({count_grid_nan}) than the DEM ({count_dem_nan}). diff: {count_grid_nan - count_dem_nan}")
 
     def read_tiff(self, tiff_filepath: Path) -> np.ndarray:
         # This method could be implemented on a GDALReader class
         src = rasterio.open(tiff_filepath)
         data = src.read(1)
         scale = src.scales[0]
         offset = src.offsets[0]
-        data = data * scale + offset
+        # data = data * scale + offset
+        data = data.astype(np.float32) * np.float32(scale) + np.float32(offset)
         src.close()
         return self.prepare_grid(data, self.conf.dem.lons.shape)
 

src/lisfloodutilities/gridding/lib/writers.py

@@ -66,6 +66,8 @@ def setNaN(self, value, defaultNaN=np.nan):
             value[value==31082] = defaultNaN
         except Exception as e:
             self.print_msg(f"value==31082 : {str(e)}")
+        value[value < self.conf.value_min_packed] = defaultNaN
+        value[value > self.conf.value_max_packed] = defaultNaN
         return value
 
     def open(self, out_filename: Path):
@@ -136,14 +138,13 @@ def open(self, out_filename: Path):
 
     def setup_grid(self, grid: np.ndarray, print_stats: bool = True) -> np.ndarray:
         values = self.setNaN(copy.deepcopy(grid))
-        values[values < self.conf.value_min_packed] = np.nan
-        values[values > self.conf.value_max_packed] = np.nan
-        values[values != self.conf.VALUE_NAN] *= self.conf.scale_factor
-        values[values != self.conf.VALUE_NAN] += self.conf.add_offset
+        values = np.trunc(values)
+        values[~np.isnan(values)] *= self.conf.scale_factor
+        values[~np.isnan(values)] += self.conf.add_offset
         if print_stats:
             self.print_grid_statistics(values)
         values[np.isnan(values)] = self.conf.VALUE_NAN * self.conf.scale_factor + self.conf.add_offset
-        return values   
+        return values
 
     def write(self, grid: np.ndarray, timestamp: datetime = None, print_stats: bool = True):
         timestep = -1