Skip to content
Snippets Groups Projects
Commit 73bd5c60 authored by Christoph.Fischer's avatar Christoph.Fischer
Browse files

w2w cluster aew vort comp

parent 6f455d39
No related tags found
No related merge requests found
......@@ -10,16 +10,26 @@ import os
# lonW = -100
# lonE = 45
in_data = "/home/ws/he7273/phd_all/data/coll_oper/jjaso2021.nc"
in_wts = "/home/ws/he7273/phd_all/data/coll_oper/jjaso2021.json"
in_data = "/project/meteo/w2w/C3/fischer/data/coll_oper/2021/jjaso2021.nc"
in_wts = "/project/meteo/w2w/C3/fischer/belanger/out/jjaso2021.json"
circles_file = "/home/ws/he7273/phd_all/data/aew/circles.nc"
fig_dir = '/home/ws/he7273/phd_all/data/vortices/figs'
# in_data = "/home/ws/he7273/phd_all/data/coll_oper/jjaso2021.nc"
# in_wts = "/home/ws/he7273/phd_all/data/coll_oper/jjaso2021.json"
circles_file = "/project/meteo/w2w/C3/fischer/data/circles.nc"
fig_dir = '/project/meteo/w2w/C3/fischer/data/vortices/figs'
# circles_file = "/home/ws/he7273/phd_all/data/aew/circles.nc"
# fig_dir = '/home/ws/he7273/phd_all/data/vortices/figs'
# kw_rain_cm_file = '/home/he7273/phd_all/data/tracked/colorpalette_dyamond_prec_rate.txt'
kw_rain_cm_file = '/project/meteo/w2w/C3/fischer/belanger/colorpalette_dyamond_prec_rate.txt'
# plot_dir = '/home/ws/he7273/phd_all/data/coll_oper/' # '/project/meteo/w2w/C3/fischer/belanger/plots/' # join('/home/ws/he7273/phd_all/data/aew/plots/') # '/project/meteo/w2w/C3/fischer/belanger/plots/'
field = 'rv'
bpf = False
bpf = True
fig_dir = fig_dir + '_' + ('bpf' if bpf else 'no_bpf') + '_' + field + '/'
try:
......@@ -36,8 +46,9 @@ base_threshold = {'cv': 3e-5, 'cva': 1e-5, 'rv': 3e-5}
# time of interest, if None all
# june-oct is AEW season
start_date = '2021-08-27T00:00' # '2022-08-01T00:00' # None # '2008-08-01T00:00' # # '2008-08-01T00:00'
end_date = '2021-08-28T10:00' # '2022-08-15T00:00' # None # '2008-08-15T00:00' # None # '2008-08-03T00:00'
## TODO CHANGE DATE AND BPF 2-6D
start_date = '2021-09-01T00:00' # '2022-08-01T00:00' # None # '2008-08-01T00:00' # # '2008-08-01T00:00'
end_date = '2021-09-15T00:00' # '2022-08-15T00:00' # None # '2008-08-15T00:00' # None # '2008-08-03T00:00'
radius = 500000 # 500 km
prominence_radius = 400000 # 300 km
......
......@@ -120,7 +120,8 @@ class AEWVortexIdentification(IdentificationStrategy):
dataset['foci_c'] = xr.zeros_like(dataset[self.field], dtype=int)
l2 = xr.zeros_like(dataset[u_name], dtype=float)
# TODO compute CPs in parallel in identify()
"""
for time in dataset.time.values:
for lev in dataset.level.values:
print(time)
......@@ -150,7 +151,8 @@ class AEWVortexIdentification(IdentificationStrategy):
dt['saddle'] = cp
dt['foci_c'] = cp_c
dataset.loc[dict(level=lev, time=time)] = dt
"""
dataset = create_wt_troughs_and_areas(dataset, self.json_desc, self.circles, False, d=self.config.radius)
"""
......@@ -190,6 +192,33 @@ class AEWVortexIdentification(IdentificationStrategy):
def identify(self, data_chunk: xr.Dataset, **kwargs):
objs = []
cp = data_chunk.saddle
cp_c = data_chunk.foci_c
u_name = self.config.u_dim
v_name = self.config.v_dim
sx = ScalarField()
sx.import_from_arrays(data_chunk.latitude.values.tolist(), data_chunk.longitude.values.tolist(), data_chunk[u_name].data)
sy = ScalarField()
sy.import_from_arrays(data_chunk.latitude.values.tolist(), data_chunk.longitude.values.tolist(), data_chunk[v_name].data)
v = VectorField()
v.import_from_sfs(sx, sy)
cpp = IMTreatment.vortex_detection.vortex_detection.get_critical_points(v)
foci = cpp.sadd[0].xy
foci_c = cpp.foc_c[0].xy
for i in range(foci.shape[0]):
cp.loc[dict(longitude=round_off_rating(foci[i][1]), latitude=round_off_rating(foci[i][0]))] = 1
for i in range(foci_c.shape[0]):
cp_c.loc[dict(longitude=round_off_rating(foci_c[i][1]), latitude=round_off_rating(foci_c[i][0]))] = 1
data_chunk['saddle'] = cp
data_chunk['foci_c'] = cp_c
# v = data_chunk.vortices
cp = data_chunk.saddle
......
......@@ -10,12 +10,10 @@ from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from PIL import Image
from datetime import datetime
def get_kitweather_rain_cm():
def get_kitweather_rain_cm(rain_cm_file):
rgb_colors = []
pathtotxtfile = '/home/he7273/phd_all/data/tracked/'
filename_colorpalette = 'colorpalette_dyamond_prec_rate.txt'
with open(pathtotxtfile + filename_colorpalette, 'r') as f:
with open(rain_cm_file, 'r') as f:
lines = f.readlines()
for i, line in enumerate(lines):
rgb_colors.append([float(line[0:3])/255, float(line[4:7])/255, float(line[8:11])/255, 1])
......@@ -79,8 +77,8 @@ def plot_kw_style(dataset, dataset_desc, config, lev=700):
subplot_kw=dict(projection=ccrs.PlateCarree()))
extent = [-75, 45, -10, 40]
levels_rain, rain_cm, norm = get_kitweather_rain_cm()
levels_rain, rain_cm, norm = get_kitweather_rain_cm(config.kw_rain_cm_file)
distance_plot_to_cbar = 0.010
axins = ax.inset_axes([1 + distance_plot_to_cbar, 0.05, 0.015, 0.93],
transform=ax.transAxes)
......@@ -131,9 +129,9 @@ def plot_kw_style(dataset, dataset_desc, config, lev=700):
# plot vortices
# ds_t.sel(level=700).vortices.plot.contourf(levels=[-0.5,0.5,99], colors=('#00000000', 'blue'), subplot_kws={'transform_first': True}, add_colorbar=False)
ds_t.sel(level=lv).saddle.plot.contourf(levels=[-0.5,0.5,99], colors=('#00000000', 'green'), subplot_kws={'transform_first': True}, add_colorbar=False)
ds_t.sel(level=lv).foci_c.plot.contourf(levels=[-0.5, 0.5, 99], colors=('#00000000', 'grey'),
subplot_kws={'transform_first': True}, add_colorbar=False)
ds_t.sel(level=lv).saddle.plot.contourf(levels=[-0.5,0.5,99], colors=('#00000000', 'orange'), subplot_kws={'transform_first': True}, add_colorbar=False)
ds_t.sel(level=lv).foci_c.plot.contourf(levels=[-0.5, 0.5, 99], colors=('#00000000', 'orange'),
subplot_kws={'transform_first': True}, add_colorbar=False) #### TODO all points seem to end up here!!
# ds_t.prec_rate_rea.plot.contourf(levels=levels_rain, extend='max', subplot_kws={'transform_first': True},
# cmap=rain_cm, norm=norm, add_colorbar=False)
......
......@@ -255,8 +255,8 @@ def bpf(dataset, u_str, v_str, min_days=2, max_days=6): # 2-6day bpf as default
time_res_days = time_res_hrs / 24.0 # 6/24 = 1/4
print("Execute bpf for " + str(min_days) + " to " + str(max_days) + " days.")
u_filtered = butter_bandpass_filter(dataset[u_str].data, time_res_days / 6, time_res_days / 2, 1, axis=0) # time is dim0
v_filtered = butter_bandpass_filter(dataset[v_str].data, time_res_days / 6, time_res_days / 2, 1, axis=0) # time is dim0
u_filtered = butter_bandpass_filter(dataset[u_str].data, time_res_days / max_days, time_res_days / min_days, 1, axis=0) # time is dim0
v_filtered = butter_bandpass_filter(dataset[v_str].data, time_res_days / max_days, time_res_days / min_days, 1, axis=0) # time is dim0
u_bpf = xr.zeros_like(dataset[u_str])
v_bpf = xr.zeros_like(dataset[v_str])
......
......@@ -13,11 +13,13 @@ data_lat = (3, 35)
data_lon = (-100, 45)
aew_clim_dir = '/project/meteo/w2w/C3/fischer/belanger/aew_clim/cv_clim_era5.nc' # '/home/ws/he7273/phd_all/data/aew/clim/cv_clim_era5.nc' # # '/lsdf/MOD/Gruppe_Transregio/Gruppe_Knippertz/kitweather/data/era5/cv_clim_era5.nc' 'C:\\Users\\Christoph\\phd\\data\\enstools-feature\\cv_clim_era5.nc' # '/home/christoph/phd/data/aew/clim/cv_clim_era5.nc' # '/home/christoph/phd/data/framework_example_ds/aew/' # '/project/meteo/w2w/C3/fischer/belanger/aew_clim/' #
in_files = '/project/meteo/w2w/C3/fischer/data/coll_oper/2021/jjaso2021.nc' # '/home/ws/he7273/phd_all/data/coll_oper/jja2021/jja2021.nc' # 'C:\\Users\\Christoph\\phd\\data\\enstools-feature\\2008_sum_uv.nc' # '/home/christoph/phd/data/framework_example_ds/aew/cv_aug_08.nc'
out_dir = '/project/meteo/w2w/C3/fischer/belanger/out/' # join('/home/ws/he7273/phd_all/data/aew/out/') # '/project/meteo/w2w/C3/fischer/belanger/out/'
in_files = '/project/meteo/w2w/C3/fischer/data/aew_new/2009/jjaso2009_uv_05.nc' #
'/home/ws/he7273/phd_all/data/coll_oper/jja2021/jja2021.nc' # 'C:\\Users\\Christoph\\phd\\data\\enstools-feature\\2008_sum_uv.nc' # '/home/christoph/phd/data/framework_example_ds/aew/cv_aug_08.nc'
out_dir = '/project/meteo/w2w/C3/fischer/data/aew_new/2009/' #
join('/home/ws/he7273/phd_all/data/aew/out/') # '/project/meteo/w2w/C3/fischer/belanger/out/'
out_json_path = out_dir + 'jjaso2021.json'
out_data_path = out_dir + 'jjaso2021_wts.nc'
out_json_path = out_dir + 'jjaso2009.json'
out_data_path = out_dir + 'jjaso2009_wts.nc'
generate_output = True
plot_dir = '/home/ws/he7273/phd_all/data/coll_oper/' # '/project/meteo/w2w/C3/fischer/belanger/plots/' # join('/home/ws/he7273/phd_all/data/aew/plots/') # '/project/meteo/w2w/C3/fischer/belanger/plots/'
......
......@@ -176,6 +176,12 @@ if enable_out:
### add WTs to orig DS. field tracks and field WTs
ds = add_wts_to_ds(ds, ob)
print(ds)
if 'u' in ds.data_vars:
ds = ds.drop_vars(['u'])
if 'v' in ds.data_vars:
ds = ds.drop_vars(['v'])
if 'cv' in ds.data_vars:
ds = ds.drop_vars(['cv'])
pipeline.set_data(ds)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment