diff --git a/mlair/data_handler/data_handler_wrf_chem.py b/mlair/data_handler/data_handler_wrf_chem.py index 5536f8d6dacb747128270474d392cef2b6d55697..5869d4a7711bbf6267e9f695779267e72c5b056b 100644 --- a/mlair/data_handler/data_handler_wrf_chem.py +++ b/mlair/data_handler/data_handler_wrf_chem.py @@ -324,7 +324,8 @@ class SingleGridColumnWrfChemDataLoader(BaseWrfChemDataLoader): target_dim=DEFAULT_TARGET_DIM, target_var=DEFAULT_TARGET_VAR, window_history_size=DEFAULT_WINDOW_HISTORY_SIZE, window_lead_time=DEFAULT_WINDOW_LEAD_TIME, - external_coords_file: str = None, **kwargs): + external_coords_file: str = None, + wind_sectors=None, **kwargs): super().__init__(**kwargs) self._set_coords(coords) self.target_dim = target_dim @@ -334,6 +335,7 @@ class SingleGridColumnWrfChemDataLoader(BaseWrfChemDataLoader): self._nearest_icoords = None self._nearest_coords = None self.external_coords_file = external_coords_file + self.wind_sectors = wind_sectors logging.debug("SingleGridColumnWrfChemDataLoader Initialised") @@ -378,7 +380,7 @@ class SingleGridColumnWrfChemDataLoader(BaseWrfChemDataLoader): # set wind_sectors # initialize WindSector() providing helper methods for wind sector calculations - ws = WindSector() + ws = WindSector(wind_sectors=self.wind_sectors) # expand dataset by new dataarray self._geo_infos['wind_sectors'] = xr.full_like(bearing, fill_value=np.nan) # xr.DataArray(coords=bearing.coords, dims=bearing.dims) @@ -414,8 +416,7 @@ class SingleGridColumnWrfChemDataLoader(BaseWrfChemDataLoader): else: raise TypeError(f"`coords' must be a tuple of floats or a dict, but is of type: {type(coords)}") - def \ - get_coordinates(self, as_arrays=False) -> Union[Tuple[np.ndarray, np.ndarray], dict]: + def get_coordinates(self, as_arrays=False) -> Union[Tuple[np.ndarray, np.ndarray], dict]: if as_arrays: return np.array(self.__coords['lat']), np.array(self.__coords['lon']) else: @@ -560,7 +561,9 @@ class DataHandlerSingleGridColumn(DataHandlerSingleStation): data = data.transpose(self.iter_dim, self.time_dim, self.target_dim, ...) data = self._slice_prep(data, start=start, end=end) # ToDo add metadata - meta = None + _meta = {'station_lon': self.loader.get_coordinates()['lon'].tolist(), + 'station_lat': self.loader.get_coordinates()['lat'].tolist()} + meta = pd.DataFrame(_meta, index=station).T return data, meta @@ -612,10 +615,11 @@ class DataHandlerSingleGridColumn(DataHandlerSingleStation): {self._logical_z_coord_name: self.targetvar_logical_z_coord_selector}) self.remove_nan(self.time_dim) - self.history = self.modify_history() + self.history = self.modify_history() self.label = self.modify_label() self.observation = self.modify_observation() + self.remove_nan(self.time_dim) # self.input_data = self.input_data.compute() # self.label = self.label.compute() @@ -623,7 +627,7 @@ class DataHandlerSingleGridColumn(DataHandlerSingleStation): # self.target_data = self.target_data.compute() # self._data.close() - def modify_history(self, **kwargs): + def modify_history(self): """ Place holder for more user spec. processing of history. :param **kwargs: @@ -633,7 +637,7 @@ class DataHandlerSingleGridColumn(DataHandlerSingleStation): """ return self.history - def modify_label(self, **kwargs): + def modify_label(self): """ Place holder for more user spec. processing of label. :param **kwargs: @@ -744,7 +748,7 @@ class DataHandlerSectorGrid(DataHandlerSingleGridColumn): return trafo @TimeTrackingWrapper - def modify_history(self, **kwargs): + 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() @@ -863,6 +867,8 @@ if __name__ == '__main__': filename = 'test_fig.pdf' if point is not None: lat, lon = point + else: + lat, lon = None, None # plt.figure(figsize=(14, 6)) plt.figure() ax = plt.axes(projection=ccrs.PlateCarree()) @@ -875,7 +881,7 @@ if __name__ == '__main__': ##### Idea to draw circle from https://stackoverflow.com/a/58735566 #### own try: - if circle: + if circle and ((lon is not None) and (lat is not None)): circle_points = Geodesic().circle(lon=lon, lat=lat, radius=radius*1000., # radius_in_meters, ) @@ -886,7 +892,8 @@ if __name__ == '__main__': ax.set_ylim(ylim) ax.set_xlim(xlim) - plt.plot(lon, lat, '+k', markersize=7, transform=crs_latlon) + if (lon is not None) and (lat is not None): + plt.plot(lon, lat, '+k', markersize=7, transform=crs_latlon) if set_background: ax.stock_img() # plt.tight_layout() @@ -931,9 +938,9 @@ if __name__ == '__main__': # lat_np = np.array([50.73333, 47.73333]) # lon_np = np.array([7.1, 8.6]) - - lat_np = np.array([50.73333]) - lon_np = np.array([7.1]) + # 48_8479__10_0963 + lat_np = np.array([48.8479]) + lon_np = np.array([10.0963]) lat_xr = xr.DataArray(lat_np, dims=["points"], coords={'points': range(len(lat_np))}) lon_xr = xr.DataArray(lon_np, dims=["points"], coords={'points': range(len(lon_np))}) @@ -942,31 +949,33 @@ if __name__ == '__main__': if use_first_dummy_dataset: wrf_new = SingleGridColumnWrfChemDataLoader((lat_xr, lon_xr), - data_path='/home/felix/Data/WRF-Chem/upload_aura_2021-02-24/2009/', + # data_path='/home/felix/Data/WRF-Chem/upload_aura_2021-02-24/2009/', + data_path="/media/felix/INTENSO/WRF_CHEM/JFM_2009/", common_file_starter='wrfout_d01_', time_dim_name='XTIME', logical_x_coord_name='west_east', logical_y_coord_name='south_north', logical_z_coord_name='bottom_top', rechunk_values={'XTIME': 24, 'bottom_top': 10}, - external_coords_file='/home/felix/Data/WRF-Chem/upload_aura_2021-02-24/coords.nc', + # external_coords_file='/home/felix/Data/WRF-Chem/upload_aura_2021-02-24/coords.nc', + external_coords_file="/media/felix/INTENSO/WRF_CHEM/monthly/coords.nc", + wind_sectors=['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW'], ) + windsect = WindSector(wind_sectors=['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW']) + radius_from_point = 100. # in km - # wrf_gridcol = SingleGridColumnWrfChemDataLoader((lat_xr, lon_xr), - # data_path='/home/felix/Data/WRF-Chem/', - # common_file_starter='wrfout_d01_2010-', - # time_dim_name='Time', - # - # ) - # wrf_gridcol.get_nearest_coords() - # wrf_dh = WrfChemDataHandler(data_path='/home/felix/Data/WRF-Chem/', - # common_file_starter='wrfout_d01_2010-', - # time_dim_name='Time', - # ) - # wrf_dh.open_data() - # wrf_dh.rechunk_data({"XTIME": 1, "y": 36, "x": 40}) - # T2 = wrf_dh._data.T2 + x_de_extract = [2, 15] + y_de_extract = [45, 60] + + x_eu_extract = [-4, 20] + y_eu_extract = [40, 60] + + with wrf_new as loader: + o3 = loader.data.o3 + plot_map_proj(data=o3.isel({'XTIME': 0, 'bottom_top': 0}), xlim=x_de_extract, ylim=y_de_extract, + filename="ozone_map_eur.pdf", cbar_kwargs={'orientation': 'horizontal', 'pad': 0.01}, + circle=False) icoords = dask.compute(wrf_new.compute_nearest_icoordinates(lat_np, lon_np))[0] dist_np = wrf_new.get_distances(lat_np, lon_np) @@ -978,8 +987,7 @@ if __name__ == '__main__': bearing_xr.attrs.update(dict(units='°')) geo_info = xr.Dataset({'dist': dist_xr, 'bearing': bearing_xr}) - windsect = WindSector() - radius_from_point = 150. # in km + ###test for rotation vecrot = VectorRotateLambertConformal2latlon(xlat=wrf_new.data.XLAT.data, xlong=wrf_new.data.XLONG.data) @@ -994,58 +1002,80 @@ if __name__ == '__main__': ustg_dim='west_east_stag', vstg_dim='south_north_stag') ### end test rotation - - for i, (data, xlim, ylim, kwargs) in enumerate(((wrf_new._data.T2.isel({'XTIME': 0}), [-42, 66], [23, 80], - {'circle': True, 'cbar_kwargs': {'orientation': 'horizontal', 'pad': 0.01}}), + point = [wrf_new.get_nearest_coords()['lat'][0], wrf_new.get_nearest_coords()['lon'][0]] + for i, (data, xlim, ylim, point, kwargs) in enumerate(((wrf_new._data.o3.isel({'XTIME': 0,'bottom_top': 0}), + x_eu_extract, y_eu_extract, point, + {'circle': True, + 'cbar_kwargs': {'orientation': 'horizontal', + 'pad': 0.01, 'shrink': .7}}), (wrf_new._data.o3.isel({'XTIME': 0, 'bottom_top': 0}), - [-42, 66], [23, 80], {'circle': True, 'cbar_kwargs': {'orientation': 'horizontal', 'pad': 0.01}}), - (wrf_new.geo_infos.dist, [-42, 66], [23, 80], - {'cmap': plt.cm.get_cmap('Greys'), 'circle': True, 'cbar_kwargs': {'orientation': 'horizontal', 'pad': 0.01}}), - (wrf_new.geo_infos.bearing, [-42, 66], [23, 80], - {'cmap': plt.cm.get_cmap('twilight'), 'circle': True, 'cbar_kwargs': {'orientation': 'horizontal', 'pad': 0.01}}), - (wrf_new.geo_infos.wind_sectors, [-42, 66], [23, 80], { - 'cmap': plt.cm.get_cmap('tab20',16), - 'circle': True, + x_eu_extract, y_eu_extract, point, + {'circle': False, + 'cbar_kwargs': {'orientation': 'horizontal', + 'pad': 0.01, 'shrink': .7}}), + (wrf_new.geo_infos.dist, x_eu_extract, y_eu_extract, point, + {'cmap': plt.cm.get_cmap('Greys'), 'circle': False, + 'cbar_kwargs': {'orientation': 'horizontal', + 'pad': 0.01, 'shrink': .7}}), + (wrf_new.geo_infos.bearing, x_eu_extract, y_eu_extract, point, + {'cmap': plt.cm.get_cmap('twilight'), 'circle': False, + 'cbar_kwargs': {'orientation': 'horizontal', + 'pad': 0.01, 'shrink': .7}}), + (wrf_new.geo_infos.wind_sectors, x_eu_extract, y_eu_extract, point, + { + 'cmap': plt.cm.get_cmap('tab20', len(wrf_new.wind_sectors)), + 'circle': False, 'cbar_kwargs': {'orientation': 'horizontal', 'pad': 0.01} #'set_background': True, # 'levels': 16, }), - (wrf_new._data.T2.isel({'XTIME': 0}).where( + (wrf_new._data.o3.isel({'XTIME': 0,'bottom_top': 0}).where( windsect.is_in_sector('N', wrf_new.geo_infos['bearing'])), - [-42, 66], [23, 80], {'circle': False, 'cbar_kwargs': {'orientation': 'horizontal', 'pad': 0.01}}), - (wrf_new._data.T2.isel({'XTIME': 0}).where( + x_eu_extract, y_eu_extract, point, + {'circle': False, + 'cbar_kwargs': {'orientation': 'horizontal', 'pad': 0.01}}), + (wrf_new._data.o3.isel({'XTIME': 0,'bottom_top': 0}).where( wrf_new.geo_infos.dist.sel({'points': 0}).drop('points') <= radius_from_point), - [2, 15], - [45, 58], {'circle': False, 'cbar_kwargs': {'orientation': 'horizontal', 'pad': 0.01}}), + x_de_extract, y_de_extract, None, + {'circle': True, + 'cbar_kwargs': {'orientation': 'horizontal', + 'pad': 0.01}}), (geo_info.dist.where( wrf_new.geo_infos.dist.sel({'points': 0}).drop('points') <= radius_from_point), - [2, 15], - [45, 58], {'circle': False, 'cbar_kwargs': {'orientation': 'horizontal', 'pad': 0.01}}), + x_de_extract, y_de_extract, None, + {'cmap': plt.cm.get_cmap('Greys'), 'circle': False, + 'cbar_kwargs': {'orientation': 'horizontal', 'pad': 0.01, + 'shrink': .6}}), (wrf_new.geo_infos.bearing.where( dist_xr.sel({'points': 0}).drop('points') <= radius_from_point), - [2, 15], - [45, 58], {'circle': False, 'cbar_kwargs': {'orientation': 'horizontal', 'pad': 0.01}}), + x_de_extract, y_de_extract, None, + {'cmap': plt.cm.get_cmap('twilight'), 'circle': False, + 'cbar_kwargs': {'orientation': 'horizontal', + 'pad': 0.01, 'shrink': .65}}), (wrf_new.geo_infos.wind_sectors.where( dist_xr.sel({'points': 0}).drop('points') <= radius_from_point), - [2, 15], - [45, 58], {'cmap': plt.cm.get_cmap('tab20', 16), + x_de_extract,y_de_extract, None, + {'cmap': plt.cm.get_cmap('tab20', len(wrf_new.wind_sectors)), 'circle': False, 'cbar_kwargs': {'orientation': 'horizontal', - 'shrink': .6, 'pad': 0.01} + 'shrink': .6, 'pad': 0.01, + } #'set_background': True, }), - (wrf_new.geo_infos.dist.where( + (wrf_new.data.o3.isel({'XTIME': 0,'bottom_top': 0}).where( wrf_new.geo_infos.dist.sel({'points': 0}).drop('points') <= radius_from_point).where( windsect.is_in_sector('N', wrf_new.geo_infos['bearing'])), - [2, 15], - [45, 58], {'circle': False, 'cbar_kwargs': {'orientation': 'horizontal', - 'fraction':.1, 'pad': 0.01}}), + x_de_extract, y_de_extract, None, + {'circle': False, 'cbar_kwargs': {'orientation': 'horizontal', + 'fraction':.1, 'pad': 0.01, + }}), )): plot_map_proj(data, xlim=xlim, ylim=ylim, # point=[lat_np, lon_np], - point=[wrf_new.get_nearest_coords()['lat'][0], wrf_new.get_nearest_coords()['lon'][0]], - filename=f'EEExample_dist{i}.png', + # point=[wrf_new.get_nearest_coords()['lat'][0], wrf_new.get_nearest_coords()['lon'][0]], + point=point, + filename=f'Example_sector{i}.pdf', radius=radius_from_point, # cbar_kwargs={'orientation': 'horizontal'}, **kwargs, @@ -1055,8 +1085,8 @@ if __name__ == '__main__': ((wrf_new._data.o3.isel({'XTIME': 0, 'bottom_top': 0}), [-42, 66], [23, 80]), (geo_info.dist, [-42, 66], [23, 80]), (wrf_new._data.o3.isel({'XTIME': 0, 'bottom_top': 0}).where( - dist_xr.sel({'points': 0}).drop('points') <= radius_from_point), [2, 15], [45, 58]), - (geo_info.dist.where(dist_xr.sel({'points': 0}).drop('points') <= radius_from_point), [2, 15], [45, 58]), + dist_xr.sel({'points': 0}).drop('points') <= radius_from_point), x_de_extract, y_de_extract), + (geo_info.dist.where(dist_xr.sel({'points': 0}).drop('points') <= radius_from_point), x_de_extract, y_de_extract), )): plot_map_proj(data, xlim=xlim, ylim=ylim, @@ -1083,11 +1113,11 @@ if __name__ == '__main__': for i, (data, xlim, ylim) in enumerate(((wrf_dh_4d._data.T2, [-42, 66], [23, 80]), (geo_info.dist, [-42, 66], [23, 80]), (wrf_dh_4d._data.T2.where( - dist_xr.sel({'points': 0}).drop('points') <= 100), [2, 15], - [45, 58]), + dist_xr.sel({'points': 0}).drop('points') <= 100), x_de_extract, + y_de_extract), (geo_info.dist.where( - dist_xr.sel({'points': 0}).drop('points') <= 100), [2, 15], - [45, 58]), + dist_xr.sel({'points': 0}).drop('points') <= 100), x_de_extract, + y_de_extract), )): plot_map_proj(data, xlim=xlim, ylim=ylim, diff --git a/run_wrf_dh.py b/run_wrf_dh.py index d278269ef876745108b0a2da025dfcc2ec44ab5e..717ce02b1732daf92c1059c18d1f3236f1b0872a 100644 --- a/run_wrf_dh.py +++ b/run_wrf_dh.py @@ -7,6 +7,8 @@ from mlair.data_handler.data_handler_wrf_chem import DataHandlerWRF, DataHandler from mlair.workflows import DefaultWorkflow from mlair.helpers import remove_items from mlair.configuration.defaults import DEFAULT_PLOT_LIST + +from mlair.model_modules.model_class import MyPaperModel, MyLSTMModel, MyCNNModel import os @@ -33,20 +35,36 @@ def main(parser_args): train_min_length=1, val_min_length=1, test_min_length=1, # data_handler=DataHandlerSingleStation, # data_handler=DataHandlerSingleGridColumn, + upsampeling=True, epochs=20, window_lead_time=2, window_history_size=6, - stations=["coords__50_7333__7_1000", "coords__50_0000__0_0000", - # "coords__48_7444__7_6000", "coords__50_0000__1_0000", + # ('Germany', (5.98865807458, 47.3024876979, 15.0169958839, 54.983104153)) + stations=["coords__48_8479__10_0963", "coords__51_8376__14_1417", + "coords__50_7536__7_0827", "coords__51_4070__6_9656", + "coords__49_8421__7_8662", "coords__49_7410__7_1935", + "coords__51_1566__11_8182", "coords__51_4065__6_9660", + # "coords__50_7536__7_0827", "coords__54_0000__10_000", + + # "coords__50_7333__7_1000", "coords__50_0000__8_0000", + # "coords__48_7444__7_6000", "coords__51_0000__11_0000", # "coords__52_7555__8_1000", "coords__50_0000__2_0000", # "coords__51_7666__8_6000", "coords__50_0000__3_0000", # "coords__45_7777__9_1000", "coords__50_0000__4_0000", ], - # data_handler=DataHandlerWRF, - data_handler=DataHandlerMainSectWRF, #, - data_path="/p/project/deepacf/inbound_data/IASS_upload/20??", + + data_handler=DataHandlerWRF, + # data_handler=DataHandlerMainSectWRF, #, + # data_path='/home/felix/Data/WRF-Chem/upload_aura_2021-02-24/2009/', + # data_path='/home/felix/Data/WRF-Chem/test_cut_nc/', + # data_path='/home/felix/Data/WRF-Chem/test_cut_nc_joint', + data_path="/home/felix/Data/WRF-Chem/test_cut_nc_joint/short_test", + # data_path="/media/felix/INTENSO/WRF_CHEM/JFM_2009", + date_format_of_nc_file="%Y-%m", + time_dim='XTIME', - external_coords_file='/p/project/deepacf/inbound_data/IASS_upload/coords.nc', + external_coords_file='/home/felix/Data/WRF-Chem/test_cut_nc/coords.nc', + transformation={ "T2": {"method": "standardise"}, "Q2": {"method": "standardise"}, @@ -62,7 +80,7 @@ def main(parser_args): 'CLDFRA': {"method": "min_max"}, }, variables=['T2', 'o3', 'wdir10ll', 'wspd10ll', 'no', 'no2', 'co', 'PSFC', 'PBLH', 'CLDFRA'], - target_var=['T2'], + target_var='o3', statistics_per_var={'T2': None, 'o3': None, 'wdir10ll': None, 'wspd10ll': None, 'no': None, 'no2': None, 'co': None, 'PSFC': None, 'PBLH': None, 'CLDFRA': None, }, # variables=['T2', 'Q2', 'PBLH', 'U10ll', 'V10ll', 'wdir10ll', 'wspd10ll'], @@ -73,7 +91,6 @@ def main(parser_args): var_logical_z_coord_selector=0, targetvar_logical_z_coord_selector=0, aggregation_dim='bottom_top', - radius=100, # km start='2009-01-01', @@ -89,7 +106,11 @@ def main(parser_args): test_end='2009-01-04', sampling='hourly', - use_multiprocessing=False, + + interpolation_limit=0, + # as_image_like_data_format=False, + # model=MyLSTMModel, + model=MyCNNModel, **parser_args.__dict__) workflow.run()