diff --git a/mlair/data_handler/data_handler_wrf_chem.py b/mlair/data_handler/data_handler_wrf_chem.py index 5869d4a7711bbf6267e9f695779267e72c5b056b..fd8b8a4df42a132515598d6e95a1e14510bf3c3a 100644 --- a/mlair/data_handler/data_handler_wrf_chem.py +++ b/mlair/data_handler/data_handler_wrf_chem.py @@ -754,26 +754,14 @@ class DataHandlerSectorGrid(DataHandlerSingleGridColumn): wind_dir_of_interest = self.compute_wind_dir_of_interest() sector_allocation = self.windsector.get_sect_of_value(value=wind_dir_of_interest, external_edges=ws_edges) existing_sectors = np.unique(sector_allocation.data) - sector_history = xr.ones_like(self.history) - sector_history_var_names = [f"{var}_sect" for var in sector_history.coords[self.target_dim].values] + sector_history, sector_history_var_names = self.setup_history_like_xr_and_var_names() with self.loader as loader, TimeTracking(name="loader in modify history"): # setup sector history grid_data = self.preselect_and_transform_neighbouring_data_based_on_radius(loader) for i, sect in enumerate(existing_sectors, start=1): with TimeTracking(name=f"processing sector {sect} ({i}/{len(existing_sectors)}) for {self.station}"): - # select data in wind sector - sec_data = self.get_section_data_from_circle(grid_data, loader, sect) - sec_data = self.apply_aggregation_method_on_sector_data(sec_data, loader) - sec_data_history = self._make_and_return_history_window(dim_name_of_shift=self.time_dim, - window=self.window_history_size, - data=sec_data.to_array(self.target_dim) - ) - sec_data_history = sec_data_history.broadcast_like(self.history) - sec_data_history = sec_data_history.transpose(*self.history.dims) - time_index_of_sector = sector_history[self.time_dim].where((sector_allocation.squeeze() == sect), - drop=True).values - sector_history.loc[{self.time_dim: time_index_of_sector}] = sec_data_history.sel( - {self.time_dim: time_index_of_sector}).values + sector_history = self.get_wind_sectordata_in_history_format(grid_data, loader, sect, + sector_allocation, sector_history) sector_history = sector_history.assign_coords({self.target_dim: sector_history_var_names}) combined_history = xr.concat([self.history, sector_history], dim=self.target_dim) @@ -782,6 +770,48 @@ class DataHandlerSectorGrid(DataHandlerSingleGridColumn): else: return self.history + def setup_history_like_xr_and_var_names(self, var_name_suffix="sect"): + """ + Returns ones_like xarray from self.history and list of variable names which can be modified by passing a + var_name_suffix + + :param var_name_suffix: + :type var_name_suffix: + :return: + :rtype: + """ + sector_history = xr.ones_like(self.history) + sector_history_var_names = [f"{var}_{var_name_suffix}" for var in sector_history.coords[self.target_dim].values] + return sector_history, sector_history_var_names + + def get_wind_sectordata_in_history_format(self, grid_data, loader, sect, sector_allocation, sector_history, + external_selector: dict = None): + # select data in wind sector + sec_data = self.get_section_data_from_circle(grid_data, loader, sect) + sec_data = self.apply_aggregation_method_on_sector_data(sec_data, loader) + sec_data_history = self._make_and_return_history_window(dim_name_of_shift=self.time_dim, + window=self.window_history_size, + data=sec_data.to_array(self.target_dim) + ) + sec_data_history = sec_data_history.broadcast_like(self.history) + sec_data_history = sec_data_history.transpose(*self.history.dims) + if external_selector is None: + time_index_of_sector = self.extract_time_index_of_sector(sect, sector_allocation, sector_history) + selector = {self.time_dim: time_index_of_sector} + else: + selector = external_selector + + sector_history.loc[selector] = sec_data_history.sel(selector).values + # sector_history.loc[{self.time_dim: time_index_of_sector}] = sec_data_history.sel( + # {self.time_dim: time_index_of_sector}).values + + return sector_history + + def extract_time_index_of_sector(self, sect, sector_allocation, sector_history): + time_index_of_sector = sector_history[self.time_dim].where((sector_allocation.squeeze() == sect), + drop=True).values + return time_index_of_sector + def get_section_data_from_circle(self, grid_data, loader, sect, compute=True): sec_data = grid_data.where( self.windsector.is_in_sector(sect, loader.geo_infos.bearing.drop('points').squeeze())) @@ -841,6 +871,69 @@ class DataHandlerMainSectWRF(DefaultDataHandler): _requirements = data_handler.requirements() +class DataHandler3SectorGrid(DataHandlerSectorGrid): + _requirements = DataHandlerSectorGrid.requirements() + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + @TimeTrackingWrapper + def modify_history(self): + if self.transformation_is_applied: + ws_edges = self.get_applied_transdormation_on_wind_sector_edges() + wind_dir_of_interest = self.compute_wind_dir_of_interest() + sector_allocation = self.windsector.get_sect_of_value(value=wind_dir_of_interest, external_edges=ws_edges) + existing_sectors = np.unique(sector_allocation.data) + sector_history, sector_history_var_names = self.setup_history_like_xr_and_var_names() + sector_history_left, sector_history_var_names_left = self.setup_history_like_xr_and_var_names(var_name_suffix="sect_left") + sector_history_right, sector_history_var_names_right = self.setup_history_like_xr_and_var_names(var_name_suffix="sect_right") + with self.loader as loader, TimeTracking(name="loader in modify history"): + # setup sector history + grid_data = self.preselect_and_transform_neighbouring_data_based_on_radius(loader) + for i, sect in enumerate(existing_sectors, start=1): + sect_left, sect_right = self.get_left_and_right_sector_of_main(sect) + with TimeTracking( + name=f"processing sector {sect} ({i}/{len(existing_sectors)}) for {self.station}"): + sector_history = self.get_wind_sectordata_in_history_format(grid_data, loader, sect, + sector_allocation, sector_history) + time_index_of_sector = self.extract_time_index_of_sector(sect, sector_allocation, sector_history) + + minor_selector = {self.time_dim: time_index_of_sector} + sector_history_left = self.get_wind_sectordata_in_history_format(grid_data, loader, sect_left, + sector_allocation, + sector_history_left, + external_selector=minor_selector) + sector_history_right = self.get_wind_sectordata_in_history_format(grid_data, loader, sect_right, + sector_allocation, + sector_history_right, + external_selector=minor_selector + ) + + sector_history = sector_history.assign_coords({self.target_dim: sector_history_var_names}) + sector_history_left = sector_history.assign_coords({self.target_dim: sector_history_var_names_left}) + sector_history_right = sector_history.assign_coords({self.target_dim: sector_history_var_names_right}) + combined_history = xr.concat([self.history, sector_history, + sector_history_left, sector_history_right], + dim=self.target_dim) + + return combined_history + else: + return self.history + + def get_left_and_right_sector_of_main(self, main_sect): + left_sect_idx = self.wind_sectors.index(main_sect) - 1 + # cyclic index: start at idx 0 if idx + 1 is larger then list of wind sectors + right_sect_idx = (self.wind_sectors.index(main_sect) + 1) % len(self.wind_sectors) + return self.wind_sectors[left_sect_idx], self.wind_sectors[right_sect_idx] + + +class DataHandlerMainMinorSectWRF(DefaultDataHandler): + """Data handler using DataHandlerSectorGrid.""" + data_handler = DataHandler3SectorGrid + data_handler_transformation = DataHandler3SectorGrid + _requirements = data_handler.requirements() + + if __name__ == '__main__': import cartopy.crs as ccrs