Skip to content
Snippets Groups Projects
Commit b7fe5224 authored by Felix Kleinert's avatar Felix Kleinert
Browse files

first version which runs (in serial mode)

parent 49498099
No related branches found
No related tags found
1 merge request!259Draft: Resolve "WRF-Datahandler should inherit from SingleStationDatahandler"
Show whitespace changes
Inline Side-by-side
mlair/data_handler/data_handler_wrf_chem.py
+ 51
− 17
View file @ b7fe5224
......@@ -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)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment