diff --git a/mlair/data_handler/data_handler_wrf_chem.py b/mlair/data_handler/data_handler_wrf_chem.py
index 5820fa95a768fc8279774f8496436cda1de3e9bb..112f8fba1e37527d83356c4dd6cbab9886c67c49 100644
--- a/mlair/data_handler/data_handler_wrf_chem.py
+++ b/mlair/data_handler/data_handler_wrf_chem.py
@@ -21,6 +21,11 @@ from typing import Tuple, Union, List, Dict
import logging
import cartopy.crs as ccrs
+# ToDo
+# test imports for circle plots
+from cartopy.geodesic import Geodesic
+import shapely
+
float_np_xr = Union[float, np.ndarray, xr.DataArray, xr.Dataset]
@@ -244,7 +249,7 @@ class SingleGridColumnWrfChemDataLoader(BaseWrfChemDataLoader):
return {k: list(v.values) for k, v in self._nearest_coords.items()}
-class DataHandlerSingleGridColumn2(DataHandlerSingleStation):
+class DataHandlerSingleGridColumn(DataHandlerSingleStation):
_requirements = remove_items(inspect.getfullargspec(DataHandlerSingleStation).args, ["self", "station"])
def __init__(self, *args, external_coords_file, **kwargs):
@@ -252,10 +257,13 @@ class DataHandlerSingleGridColumn2(DataHandlerSingleStation):
super().__init__(*args, **kwargs)
@staticmethod
- def str2coords(str_coords: str, sep='-') -> Tuple[float_np_xr, float_np_xr]:
+ def coord_str2coords(str_coords: str, sep='__', dec_marker='_') -> Tuple[float_np_xr, float_np_xr]:
if isinstance(str_coords, list) and len(str_coords) == 1:
str_coords = str_coords[0]
- lat, lon = str_coords.split(sep=sep)
+ _, lat, lon = str_coords.split(sep=sep)
+ lat = lat.replace(dec_marker, '.')
+ lon = lon.replace(dec_marker, '.')
+
return np.array(float(lat)), np.array(float(lon))
def setup_data_path(self, data_path: str, sampling: str):
@@ -264,7 +272,7 @@ class DataHandlerSingleGridColumn2(DataHandlerSingleStation):
def load_data(self, path, station, statistics_per_var, sampling, station_type=None, network=None,
store_data_locally=False, data_origin: Dict = None, start=None, end=None):
- lat, lon = self.str2coords(station)
+ lat, lon = self.coord_str2coords(station)
sgc_loader = SingleGridColumnWrfChemDataLoader((lat, lon),
data_path=path,
external_coords_file=self.external_coords_file,
@@ -289,6 +297,15 @@ class DataHandlerSingleGridColumn2(DataHandlerSingleStation):
# raise NotImplementedError
return data.chunk({self.time_dim:-1}), meta
+ def get_X(self, upsampling=False, as_numpy=False):
+ if as_numpy is True:
+ return None
+ elif as_numpy is False:
+ return self.get_transposed_history()
+
+ # def get_Y(self, upsampling=False, as_numpy=False):
+ # raise NotImplementedError
+
# def set_inputs_and_targets(self):
# # inputs = self._data.sel({self.target_dim: helpers.to_list(self.variables)})
# # targets = self._data.sel(
@@ -323,8 +340,8 @@ class DataHandlerSingleGridColumn2(DataHandlerSingleStation):
class DataHandlerWRF(DefaultDataHandler):
"""Data handler using CDC."""
- data_handler = DataHandlerSingleGridColumn2
- data_handler_transformation = DataHandlerSingleGridColumn2
+ data_handler = DataHandlerSingleGridColumn
+ data_handler_transformation = DataHandlerSingleGridColumn
_requirements = data_handler.requirements()
@@ -332,6 +349,7 @@ class DataHandlerWRF(DefaultDataHandler):
if __name__ == '__main__':
def plot_map_proj(data, xlim=None, ylim=None, filename=None, point=None):
+
crs_latlon = ccrs.PlateCarree()
if ylim is None:
ylim = [-90, 90]
@@ -339,6 +357,8 @@ if __name__ == '__main__':
xlim = [-180, 180]
if filename is None:
filename = 'test_fig.pdf'
+ if point is not None:
+ lat, lon = point
# plt.figure(figsize=(14, 6))
plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree())
@@ -346,9 +366,23 @@ if __name__ == '__main__':
data.squeeze().plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree(), x='XLONG', y='XLAT',
cmap=plt.cm.Reds, )
ax.coastlines(resolution='50m')
+
+ ##### Idea to draw circle from https://stackoverflow.com/a/58735566
+ # circle_points = cartopy.geodesic.Geodesic().circle(lon=lon, lat=lat, radius=radius_in_meters,
+ # n_samples=n_points, endpoint=False)
+ # geom = shapely.geometry.Polygon(circle_points)
+ # ax.add_geometries((geom,), crs=cartopy.crs.PlateCarree(), facecolor='red', edgecolor='none', linewidth=0)
+ #### own try:
+ circle_points = Geodesic().circle(lon=lon, lat=lat,
+ radius=100000, # radius_in_meters,
+ )
+ geom = shapely.geometry.Polygon(circle_points)
+ ax.add_geometries((geom,), crs=ccrs.PlateCarree(), facecolor='none', edgecolor='blue', linewidth=1)
+ #
+ #####
ax.set_ylim(ylim)
ax.set_xlim(xlim)
- plt.plot(point[1], point[0], 'bo', markersize=7, transform=crs_latlon)
+ plt.plot(point[1], point[0], '+k', markersize=7, transform=crs_latlon)
plt.tight_layout()
plt.savefig(filename)
plt.close('all')
@@ -364,13 +398,13 @@ 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/',
common_file_starter='wrfout_d0',
- time_dim_name='Time',
+ time_dim_name='XTIME',
logical_x_coord_name='west_east',
logical_y_coord_name='south_north',
logical_z_coord_name='bottom_top',
- rechunk_values={'Time': 1, 'bottom_top': 2},
+ rechunk_values={'XTIME': 1, 'bottom_top': 2},
external_coords_file='/home/felix/Data/WRF-Chem/upload_aura_2021-02-24/coords.nc',
)
@@ -396,9 +430,9 @@ if __name__ == '__main__':
dist_xr.attrs.update(dict(units='km'))
dist_xr_set = xr.Dataset({'dist': dist_xr})
- for i, (data, xlim, ylim) in enumerate(((wrf_new._data.T2.isel({'Time': 0}), [-42, 66], [23, 80]),
+ for i, (data, xlim, ylim) in enumerate(((wrf_new._data.T2.isel({'XTIME': 0}), [-42, 66], [23, 80]),
(dist_xr_set.dist, [-42, 66], [23, 80]),
- (wrf_new._data.T2.isel({'Time': 0}).where(
+ (wrf_new._data.T2.isel({'XTIME': 0}).where(
dist_xr.sel({'points': 0}).drop('points') <= 100), [2, 15],
[45, 58]),
(dist_xr_set.dist.where(
@@ -410,9 +444,9 @@ if __name__ == '__main__':
point=[lat_np, lon_np], filename=f'Example_dist{i}.pdf')
for i, (data, xlim, ylim) in enumerate(
- ((wrf_new._data.o3.isel({'Time': 0, 'bottom_top': 0}), [-42, 66], [23, 80]),
+ ((wrf_new._data.o3.isel({'XTIME': 0, 'bottom_top': 0}), [-42, 66], [23, 80]),
(dist_xr_set.dist, [-42, 66], [23, 80]),
- (wrf_new._data.o3.isel({'Time': 0, 'bottom_top': 0}).where(
+ (wrf_new._data.o3.isel({'XTIME': 0, 'bottom_top': 0}).where(
dist_xr.sel({'points': 0}).drop('points') <= 100), [2, 15], [45, 58]),
(dist_xr_set.dist.where(dist_xr.sel({'points': 0}).drop('points') <= 100), [2, 15], [45, 58]),
)):
@@ -423,13 +457,13 @@ if __name__ == '__main__':
######################### # Larger 4D data
use_second_dummy_dataset = False
if use_second_dummy_dataset:
- wrf_dh_4d = BaseWrfChemDataLoader(data_path='/home/felix/Data/WRF-Chem/upload_aura/2009/2009',
+ wrf_dh_4d = BaseWrfChemDataLoader(data_path='/home/felix/Data/WRF-Chem/upload_aura/2009/',
common_file_starter='wrfout_d01_2009', time_dim_name='Time')
wrf_dh_4d.open_data()
- wrf_dh_4d.rechunk_data({"Time": 1, "bottom_top": 34, "south_north": 36, "west_east": 40})
+ wrf_dh_4d.rechunk_data({"XTIME": 1, "bottom_top": 34, "south_north": 36, "west_east": 40})
lat_np = np.array([50.73333])
lon_np = np.array([7.1])
- wrf_dh_4d._data = wrf_dh_4d._data.assign_coords(wrf_dh._data.coords)
+ wrf_dh_4d._data = wrf_dh_4d._data.assign_coords(wrf_dh_4d._data.coords)
icoords = dask.compute(wrf_dh_4d.compute_nearest_icoordinates(lat_np, lon_np))[0]
dist_xr = wrf_dh_4d.get_distances(lat_xr, lon_xr)