-
Notifications
You must be signed in to change notification settings - Fork 0
/
combine_aeolus_rayleigh.py
116 lines (94 loc) · 5.42 KB
/
combine_aeolus_rayleigh.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
import xarray
import datetime
import os
import shutil
import re
import logging
import matplotlib.pyplot as plt
logging.basicConfig(level=logging.DEBUG)
mylogger = logging.getLogger(__name__)
mylogger.debug('This message should go to the log file')
mylogger.info('So should this')
#logging.warning('And this, too')
mylogger.error('And non-ASCII stuff, too, like Øresund and Malmö')
mylogger.info ("is this working")
class NetCDFInterFaceClass():
def __init__(self, logger=False, analysis=None, path=None):
self.logger=logger
self.AnalysisTime= analysis
self.project_dynamic_data_partition = os.path.join("/home/stephen/","BlackForestProjectData/", "dynamic")
this_path=os.path.join("/home/stephen/","BlackForestProjectData/AUG19/Cycle/CTRL/20190817/0600/obsproc/production")
self.path=dict(obsproc=dict(production=this_path))
def _read_single_file_(self, single_netcdf_file):
self.logger.debug("opening {}".format(single_netcdf_file))
self.prof_data = xarray.open_dataset(single_netcdf_file, group='rayleigh_profile_data')
self.wind_data = xarray.open_dataset(single_netcdf_file, group='rayleigh_wind_data')
# note the date is since datetime.datetime(2000, 1, 1)
epoch = datetime.datetime(2000, 1, 1)
Tstart = []
Tcog = []
for key in zip(self.wind_data.rayleigh_wind_result_start_time,
self.wind_data.rayleigh_wind_result_COG_time):
Tstart.append(epoch + datetime.timedelta(seconds=key[0].item(0)))
Tcog.append(epoch + datetime.timedelta(seconds=key[1].item(0)))
# make a figure
figure_name = self.AnalysisTime.strftime("Start_time_%Y-%m-%d.png")
fig, ax = plt.subplots(1,1,figsize=(8,8))
ax.plot(Tstart, label="{}".format('start'))
ax.plot(Tcog, label="{}".format('COG'))
plt.legend()
plt.savefig(figure_name)
plt.show()
a=1
def _combine_daily_aeolus_(self):
"""
the Aeolus data is available as a daily dataset
we need to combine this into a dataset for the specific assimilation window
obsproc
"""
FILE_MASK = 'L2B_rayleigh_wind_(?P<start>\d{4}-\d{2}-\d{2})-'
FILE_MASK += '(?P<end>\d{4}-\d{2}-\d{2})\.nc'
aeolus_repository = os.path.join(self.project_dynamic_data_partition, "satellite/aeolus/netcdf/daily")
selected_aeolus_files = []
Analysis_day = self.AnalysisTime.strftime("%Y-%m-%d")
Analysis_next_day = (self.AnalysisTime + datetime.timedelta(days=1)).strftime("%Y-%m-%d")
aeolus_file = "L2B_rayleigh_wind_"+Analysis_day+"-"+Analysis_next_day+".nc"
#for aeolus_file in daily_aeolus_files:
# file_matches_mask = re.match(FILE_MASK, aeolus_file)
# if file_matches_mask:
# data_granula_start = datetime.datetime.strptime(file_matches_mask.group("start"), "%Y-%m-%d")
# data_granula_end = datetime.datetime.strptime(file_matches_mask.group("end"), "%Y-%m-%d")
# # we pick all files within a window of 24 hours around analysis time
# if abs((data_granula_start - self.AnalysisTime).total_seconds()) < 86400 or abs(
# (data_granula_end - self.AnalysisTime).total_seconds()) < 86400:
# selected_aeolus_files.append(os.path.join(aeolus_repository, 'daily', aeolus_file))
self.aeolus_data_file = "_".join(['aeolus_data_file', self.AnalysisTime.strftime("%Y%m%d%H%M")]) + ".nc"
self.full_path_aeolus_data_file = os.path.join(self.path["obsproc"]["production"], self.aeolus_data_file)
if os.path.exists(self.full_path_aeolus_data_file):
os.remove(self.full_path_aeolus_data_file)
src = os.path.join(aeolus_repository, aeolus_file)
shutil.copyfile(src, self.full_path_aeolus_data_file)
#if (len(selected_aeolus_files) >= 0):
# self.logger.debug("for {} we found {} aeolus files".format(self.AnalysisTime, len(selected_aeolus_files)))
# self.logger.debug("{}".format(selected_aeolus_files))
# ds_rayleigh = {'rayleigh_profile_data': [], "rayleigh_wind_data": []}
# for groupno, group in enumerate(ds_rayleigh):
# for aeolus_file_key in selected_aeolus_files:
# ds_rayleigh[group].append(xarray.open_dataset(aeolus_file_key, group=group))
# ds = xarray.concat(ds_rayleigh[group], dim=group)
# if groupno == 0:
# ds.to_netcdf(self.full_path_aeolus_data_file, 'w', group=group)
# else:
# ds.to_netcdf(self.full_path_aeolus_data_file, 'a', group=group)
if __name__ == '__main__':
"""
Given the assimilation date we combine the Aeolus rayleigh observations around this date into one file
"""
analysis = datetime.datetime.strptime("201908170600","%Y%m%d%H%M")
aeolus_data = NetCDFInterFaceClass(analysis=analysis, logger=mylogger)
aeolus_data_file = "L2B_rayleigh_wind_2019-08-17-2019-08-18.nc"
# aeolus_data._combine_daily_aeolus_()
single_netcdf_file = os.path.join(aeolus_data.project_dynamic_data_partition, "satellite/aeolus/netcdf/daily",
aeolus_data_file)
single_netcdf_file = os.path.join("/home/stephen/BlackForestProjectData/AUG19/Cycle/CTRL/20190817/0600/obsproc/production/aeolus_data_file_201908170600.nc")
aeolus_data._read_single_file_(single_netcdf_file)