Commit df419572 authored by Isabelle Pison's avatar Isabelle Pison
Browse files

netcdf CAMS provisory fetch

parent 671605f8
......@@ -13,7 +13,7 @@ def fetch(ref_dir, ref_file, input_dates, target_dir, tracer=None, **kwargs):
for datei in input_dates:
tmp_files = []
tmp_dates = []
#print('DDDDDDDDDDD',datei)
print('ddddddddddddd',datei)
for dd in input_dates[datei]:
dir_dd = dd.strftime(ref_dir)
dir_dd_next = (dd+datetime.timedelta(days=1)).strftime(ref_dir)
......@@ -34,5 +34,9 @@ def fetch(ref_dir, ref_file, input_dates, target_dir, tracer=None, **kwargs):
unique_dates, unique_index = np.unique(tmp_dates, return_index=True)
list_files[datei] = np.array(tmp_files)[unique_index]
list_dates[datei] = unique_dates
# print('lllllllll',len(list_files))
# print('LLLLLLLLLL',len(list_dates))
return list_files, list_dates
......@@ -24,10 +24,10 @@ def get_domain(ref_dir, ref_file, input_dates, target_dir, tracer=None):
pv = grib_file_reader(files_3d[0], [], attribute="pv")
hybrid = grib_file_reader(files_3d[0], ["hybrid"])[0]
lon_min = lon.min() - (lon[1] - lon[0]) / 2
lon_max = lon.max() + (lon[-1] - lon[-2]) / 2
lat_min = lat.min() - (lat[1] - lat[0]) / 2
lat_max = lat.max() + (lat[-1] - lat[-2]) / 2
lon_min = lon.min() #- (lon[1] - lon[0]) / 2
lon_max = lon.max() #+ (lon[-1] - lon[-2]) / 2
lat_min = lat.min() #- (lat[1] - lat[0]) / 2
lat_max = lat.max() #+ (lat[-1] - lat[-2]) / 2
nlon = lon.size
nlat = lat.size
......
......@@ -30,12 +30,14 @@ def read(
list_files = tracfile[:]
xout = []
print('dates',dates)
for dd, dd_file in zip(dates, list_files):
dir_dd = dd.strftime(tracdir)
jscan = grib_file_reader("{}/{}".format(dir_dd, dd_file), [], 'jScansPositively')
spec = grib_file_reader("{}/{}".format(dir_dd, dd_file), varnames)
xout.append(spec)
# Flipping upside down to follow convention from other models
if( jscan == 0) :
xmod = xr.DataArray(
......
......@@ -9,28 +9,31 @@ def fetch(ref_dir, ref_file, input_dates, target_dir, tracer=None, **kwargs):
list_files = {}
list_dates = {}
#print('fffffffffffff',input_dates)
for datei in input_dates:
print('ddddddddd',datei)
tmp_files = []
tmp_dates = []
for dd in input_dates[datei]:
#print('ppppppppppp',dd)
# WARNING: does not work if change in the year -> see grib ECMWF for adaptation
files_cams, dates_cams = find_valid_file(ref_dir, ref_file, dd)
print(files_cams,dates_cams)
#print(files_cams,dates_cams)
tmp_files.extend(files_cams)
tmp_dates.extend(dates_cams)
# print('FFFFFFFFFFFF',tmp_dates,tmp_files)
# Fetching
local_files = []
#print('RRRRRRRRRR',ref_file)
for f, dd in zip(tmp_files, tmp_dates):
target_file = "{}/{}".format(target_dir, dd.strftime(ref_file))
path.link(ref_dir+'/'+f, target_file)
print('ttt',target_file)
path.link(f, target_file)
local_files.append(target_file)
unique_dates, unique_index = np.unique(tmp_dates, return_index=True)
list_files[datei] = np.array(tmp_files)[unique_index]
list_dates[datei] = unique_dates
#print('oooooooooooooooo',list_dates)
print('lllllllll',len(list_files))
print('LLLLLLLLLL',len(list_dates))
return list_files, list_dates
......@@ -33,13 +33,15 @@ def get_domain(ref_dir, ref_file, input_dates, target_dir, tracer=None):
"CAMS domain could not be initialized as no file was found"
)
print('Domain file for CAMS BCs:',domain_file)
nc = xr.open_dataset(domain_file, decode_times=False)
# Read lon/lat in domain_file
lon = nc["longitude"]
lat = nc["latitude"]
nlon = lon.size
nlat = lat.size
nlat = lat.size
#print(nlon,nlat)
#print(lon[0])
lon_min = lon.min()
......@@ -48,13 +50,16 @@ def get_domain(ref_dir, ref_file, input_dates, target_dir, tracer=None):
lat_max = lat.max()
print(lon_min)
# Compute corners #WARNING: only valid for regular grid in lat/lon -> to generalize
dx = (lon[1] - lon[0]) / 2
dy = (lat[1] - lat[0]) / 2
lonc = np.linspace(lon_min - dx/2., lon_max + dx/2., nlon + 1)
latc = np.linspace(lat_min - dy/2., lat_max + dy/2., nlat + 1)
# Read vertical information in domain_file
sigma_a = nc["ap"].values
sigma_b = nc["bp"].values
#print('SSSSSSSSSSSSSSSS',sigma_a)
#print('ssssssssssssssss',sigma_b)
nlevs = sigma_a.size - 1
#print('LLLLLLLLLLLLLLLLLLLLLL',nlevs)
# Initializes domain
......@@ -81,5 +86,13 @@ def get_domain(ref_dir, ref_file, input_dates, target_dir, tracer=None):
)
Setup.load_setup(setup, level=1)
zlon, zlat = np.meshgrid(lon, lat)
zlonc, zlatc = np.meshgrid(lonc, latc)
setup.domain.zlon = zlon
setup.domain.zlat = zlat
setup.domain.zlonc = zlonc
setup.domain.zlatc = zlatc
return setup.domain
......@@ -29,7 +29,7 @@ def read(
comp_type: type of boundary conditions to generate
"""
print('ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ')
if type(tracfile) == str:
tracfile = [dd.strftime(tracfile) for dd in dates]
else:
......@@ -42,45 +42,36 @@ def read(
)
list_files = tracfile[:]
print('LLLLLLLLLLLLLL',list_files)
# Reading fields for periods within the simulation window
trcr_conc = []
times = []
xout = []
print('dates',dates) #devrai etre liste des dates toutes les 3 heures
mdmdfs
for dd, dd_file in zip(dates, list_files):
file_conc = dd.strftime(dd_file)
dir_conc = dd.strftime(tracdir)
if not os.path.isfile("{}/{}".format(dir_conc, file_conc)) and getattr(
self, "closest_year", False
):
info(
"Warning: could not find correct year for CAMS; "
"using closest available one"
)
list_dates = [
datetime.datetime.strptime(os.path.basename(f), tracfile)
for f in glob.glob("{}/z_cams_l_nilu_*nc".format(dir_conc))
]
delta_dates = np.abs(dd - np.array(list_dates))
file_conc = list_dates[np.argmin(delta_dates)].strftime(tracfile)
nc = xr.open_dataset(
"{}/{}".format(dir_conc, file_conc), decode_times=False
)
# Convert number of hours since first day of the month to datetime
times.extend(
[datetime.datetime(dd.year, dd.month, 1) + datetime.timedelta(hours=t)
for t in nc['time'].values]
)
trcr_conc.append(nc[varnames]['time' == dates].values)
print('ffff',file_conc)
print('ddddd',dd)
conc = readnc(file_conc, [varnames])
#print(conc.shape)
time = readnc(file_conc, ['time'])
#print(time)
date_file = (file_conc.split('nilu_')[1]).split('_v')[0]
print(date_file)
date_ref = datetime.datetime(year = int(date_file[0:4]), month = int(date_file[4:6]), day =1 )
print(date_ref)
#print(date_ref)
dates_file = [ date_ref + datetime.timedelta(hours = t ) for t in time ]
print(dates_file)
d = dates_file.index(dd)
#print(d)
xout.append(conc[ d, :, :, :])
#conc = xr.open_dataset(file_conc, decode_times = True)
#spec = conc.where( datetime.datetime(conc.time).year == dd.year ) #.dropna('time')
#print('ssss',spec[varnames])
#xout.append(spec[varnames])
# Putting in DataArray for the right period = combien d'heures??? Juste dates???
# XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
xmod = xr.DataArray(
np.array(trcr_conc),
np.array(xout)[0,:,:,:],
coords={"time": dates},
dims=("time", "lev", "lat", "lon"),
)
......
import datetime
import glob
import os
import calendar
import numpy as np
def find_valid_file(ref_dir, file_format, dd):
# Get all files and dates matching the file and format
list_files_orig = os.listdir(ref_dir)
list_dates = []
list_files = []
list_files_avail = [] #list available files
# if dd first hour of the month, needs file of the month before
for f in list_files_orig:
try:
list_dates.append(datetime.datetime.strptime(f, file_format))
list_files.append(f)
except:
try:
check = datetime.datetime.strptime(f, file_format)
list_files_avail.append(f)
except:
continue
print(list_files_avail)
# actually in each file, information every 3 hour
delta_t = 3
list_files = []
list_dates = []
for f in list_files_avail:
date_beg = datetime.datetime.strptime(f, file_format)
date_end = date_beg + datetime.timedelta( days = calendar.mdays[date_beg.month] + calendar.isleap(date_beg.year))
nb_deltas = int((date_end - date_beg).days * 24 / delta_t)
print(nb_deltas)
list_dates_covered = [ date_beg + ( k + 1 ) * datetime.timedelta( hours = delta_t ) for k in range(nb_deltas) ]
print(list_dates_covered)
for d in list_dates_covered:
list_files.append(f)
list_dates.append(d)
list_files = np.array(list_files)
list_dates = np.array(list_dates)
......@@ -29,17 +46,24 @@ def find_valid_file(ref_dir, file_format, dd):
raise Exception("Did not find any valid CAMS files in {} "
"with format {}. Please check your yml file"
.format(ref_dir, file_format))
#print('yyyyyyyyyyyyyyy',dd,list_dates,list_files)
# si date = premiere heure du mois -> fichier = celui du mois d'avant
# sinon, fichier = celui du mois
if dd.hour == 0 and dd.day == 1:
file_ref1 = list_files[0]
dd1 = list_dates[0]
else:
file_ref1 = list_files[1]
dd1 = list_dates[1]
file_ref2 = list_files[1]
dd2 = list_dates[1]
# Convert ref date
ref_date = datetime.datetime.strptime(dd.strftime(file_format), file_format)
print('ref_date',ref_date)
# Compute deltas
mask = (list_dates - ref_date) <= datetime.timedelta(0)
file_ref1 = ref_dir + list_files[mask][np.argmax(list_dates[mask])]
date_ref1 = list_dates[mask].max()
print(file_ref1, date_ref1)
date_ref2 = date_ref1 + datetime.timedelta(hours = delta_t)
print(list_files[np.where(list_dates == date_ref2)][0])
file_ref2 = ref_dir + list_files[np.where(list_dates == date_ref2)][0]
print(file_ref2, date_ref2)
# Reconvert to original date
dd1 = dd + (date_ref1 - ref_date)
dd2 = dd + (date_ref2 - ref_date)
return [file_ref1, file_ref2], [dd1, dd2]
......@@ -84,6 +84,7 @@ def make_obs(self, ddi, datastore, runsubdir, mode, tracer, do_simu=True):
# extract first level
mask = np.isnan(dataout["level"]) | (dataout["level"] <= 0)
dataout.loc[mask, "level"] = 1
#print(verfi dataout = level du monitor.nc)
# Add column with index of current hour
dataout.loc[:, "hour"] = \
......@@ -112,4 +113,4 @@ def make_obs(self, ddi, datastore, runsubdir, mode, tracer, do_simu=True):
data_prior.to_csv(f, sep=" ", index=False, header=False)
self.iniobs = True
\ No newline at end of file
......@@ -34,6 +34,7 @@ def init_reindex(
# Temporal re-indexing if any
tinterp = getattr(param, "time_interpolation", None)
#if interp == 'grut': returnXXXXXXXXXXXXx
yml_dict = {
"plugin": {
"name": "time_interpolation",
......
......@@ -29,7 +29,7 @@ def calc_indexes(mapper, trid_out):
in_dates_start = np.array(in_dates[ddi][:-1])
in_dates_end = np.array(in_dates[ddi][1:])
in_dates_all = np.array(in_dates[ddi])
#print(les in_dates, les out_dates: Faire un systeme de correspondance entre 2005 et 2011 modifiant de force les in_dates)
if len(in_dates_all) <= 1:
continue
......
......@@ -31,6 +31,7 @@ def adjoint(
# Stop here if do not need to compute the full adjoint
sparse_data = mapper["outputs"][trid].get("sparse_data", False)
if sparse_data:
#print(xmod[trid] -> level identique a celui du monitor.nc d'entree??)
inout_datastore["inputs"][trid] = xmod[trid]
continue
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment