Commit 7a240a52 authored by Antoine Berchet's avatar Antoine Berchet
Browse files

Warning when CHIMERE periods are not compatible with simulation window;...

Warning when CHIMERE periods are not compatible with simulation window; implementing becker_ocean fluxes
parent 01615143
#QSUB -s /bin/tcsh #QSUB -s /bin/tcsh
#PBS -q longp #PBS -q xlongp
#PBS -l nodes=1:ppn=8 #PBS -l nodes=1:ppn=10
python -m pycif /home/users/aberchet/CIF/config_files/tuto_chimere/config_chimere_argonaut_n2o_sat_inv_corr.yml python -m pycif /home/users/aberchet/CIF/config_files/RECAPP/config_chimere_fwd_EUROCOM_CO2_AB_6mois_TNO_2008_1.yml
#python -m pycif /homen/users/aberchet/CIF/config_files/tuto_chimere/config_chimere_EUROCOM_satOMI.yml
#QSUB -s /bin/tcsh
#PBS -q xlongp
#PBS -l nodes=1:ppn=10
#PBS -N RECAPP_2008_1
python -m pycif /home/users/aberchet/CIF/config_files/RECAPP/config_chimere_fwd_EUROCOM_CO2_AB_6mois_TNO_2008_1.yml
...@@ -70,7 +70,6 @@ def ini_data(self, **kwargs): ...@@ -70,7 +70,6 @@ def ini_data(self, **kwargs):
workdir = self.workdir workdir = self.workdir
dirchem_ref = "{}/chemical_scheme/{}/".format(workdir, self.schemeid) dirchem_ref = "{}/chemical_scheme/{}/".format(workdir, self.schemeid)
self.dirchem_ref = dirchem_ref self.dirchem_ref = dirchem_ref
shutil.rmtree(dirchem_ref, ignore_errors=True) shutil.rmtree(dirchem_ref, ignore_errors=True)
init_dir(dirchem_ref) init_dir(dirchem_ref)
......
"""
These fluxes were generated for coastal and ocean CO2 sea-air exchanges
"""
from .fetch import fetch
from .get_domain import get_domain
from .read import read
from .write import write
from logging import info
_name = "becker"
_version = "ocean"
_fullname = "Becker coastal fluxes"
input_arguments = {
"file_freq": {
"doc": "Temporal frequency at which files are available",
"default": "1YS",
"accepted": str
}
}
import datetime
import glob
import os
import pandas as pd
import numpy as np
from .....utils import path
from logging import debug
def fetch(ref_dir, ref_file, date_interval, target_dir,
tracer=None, component=None, **kwargs):
"""
Fetch files and dates for Becker coastal fluxes
Args:
ref_dir (str): the path to the input files
ref_file (str): format of the input files
input_dates (list): simulation interval (start and end dates)
target_dir (str): where to copy
tracer: the tracer Plugin, corresponding to the paragraph
:bash:`datavect/components/fluxes/parameters/my_species` in the
configuration yaml; can be needed to fetch extra information
given by the user
component: the component Plugin, same as tracer; corresponds to the paragraph
:bash:`datavect/components/fluxes` in the configuration yaml
Return:
(dict, dict): returns two dictionaries: list_files and list_dates
list_files: for each date that begins a period, a list containing
the names of the files that are available for the dates within this period
list_dates: for each date that begins a period, a list containing
the date intervals (in the form of a list of two dates each)
matching the files listed in list_files
"""
# Reshape input interval to include full years
datei, datef = date_interval
datei = datetime.datetime(year=datei.year, month=1, day=1)
datef = datetime.datetime(year=datef.year + 1, month=1, day=1)
list_period_dates = pd.date_range(datei, datef, freq=tracer.file_freq)
list_dates = {}
list_files = {}
for dd in list_period_dates:
file = dd.strftime("{}/{}".format(ref_dir, ref_file))
file_hours = pd.date_range(
dd, dd + pd.to_timedelta(tracer.file_freq), freq="1H")
list_dates[dd] = [[hh, hh + datetime.timedelta(hours=1)]
for hh in file_hours]
list_files[dd] = (len(file_hours) * [file])
# Fetching
if os.path.isfile(file):
target_file = "{}/{}".format(target_dir, os.path.basename(file))
path.link(file, target_file)
debug(
"Fetched files and dates as follows:\n"
"Dates: {\n" +
"\n".join(["\n".join([" {}:".format(ddi)]
+ [" {}".format(dd)
for dd in list_dates[ddi]])
for ddi in list_dates])
+ "\n}\n\n" +
"Files: {\n" +
"\n".join(["\n".join([" {}:".format(ddi)]
+ [" {}".format(dd)
for dd in list_files[ddi]])
for ddi in list_files])
+ "\n}"
)
return list_files, list_dates
import numpy as np
from logging import debug
from .....utils.classes.setup import Setup
from .....utils.classes.domains import Domain
def get_domain(ref_dir, ref_file, input_dates, target_dir, tracer=None):
"""Read information to define the data horizontal and, if relevant, vertical domain.
There are several possible approaches:
- read a reference file that is necessary in :bash:`ref_dir`
- read a file among the available data files
- read a file specified in the yaml,
by using the corresponding variable name; for instance, tracer.my_file
From the chosen file, obtain the coordinates of the centers and/or the corners
of the grid cells. If corners or centers are not available, deduce them from
the available information.
Warning:
the grid must not be overlapping: e.g for a global grid,
the last grid cell must not be the same as the first
Warning:
Order the centers and corners latitudes and longitudes increasing order
Args:
ref_dir (str): the path to the input files
ref_file (str): format of the input files
input_dates (list): simulation interval (start and end dates)
target_dir (str): where to copy
tracer: the tracer Plugin, corresponding to the paragraph
:bash:`datavect/components/fluxes/parameters/my_species` in the
configuration yaml; can be needed to fetch extra information
given by the user
Return:
Domain: a domain class object, with the definition of the center grid
cells coordinates, as well as corners
"""
# Print some explanation
debug(__doc__)
# For the purpose of demonstration, the domain dimensions are specified by default
# in input_arguments.
# Individual arguments can be modified manually in the yaml
lon_min = tracer.lon_min
lon_max = tracer.lon_max
lat_min = tracer.lat_min
lat_max = tracer.lat_max
nlon = tracer.nlon
nlat = tracer.nlat
debug("if lon and lat are flat vectors in the definition, "
"convert into a grid with:\n"
"zlon, zlat = np.meshgrid(lon, lat)\n")
# Some explanations
debug(
'Here, read the vertical information, from the same file as the horizontal '
'information or from another.\n'
'Get the number of vertical levels.\n'
'Get or deduce the sigma_a/sigma_b coefficients from bottom to top if the '
'vertical '
'extension is in pressure.\n'
'It is possible to specify the vertical extension in m a.g.l. as well by '
'defining the variable heights'
)
nlev = tracer.nlev
punit = "Pa"
sigma_a = np.linspace(0, 1, nlev)
sigma_b = np.linspace(1, 0, nlev)
# Initializes domain
setup = Setup.from_dict(
{
"domain": {
"plugin": {
"name": "dummy",
"version": "std",
"type": "domain",
},
"xmin": lon_min, # minimum longitude for centers
"xmax": lon_max, # maximum longitude for centers
"ymin": lat_min, # minimum latitude for centers
"ymax": lat_max, # maximum latitude for centers
"nlon": nlon, # number of longitudinal cells
"nlat": nlat, # number of latitudinal cells
"nlev": nlev, # number of vertical levels
"sigma_a": sigma_a,
"sigma_b": sigma_b,
"pressure_unit": "Pa" # adapted to sigmas
}
}
)
Setup.load_setup(setup, level=1)
return setup.domain
import numpy as np
import xarray as xr
from logging import debug
def read(
self,
name,
varnames,
dates,
files,
interpol_flx=False,
tracer=None,
model=None,
ddi=None,
**kwargs
):
"""Get fluxes from raw files and load them into a pyCIF
variables.
The list of date intervals and corresponding files is directly provided,
coming from what is returned by the :bash:`fetch` function.
One should loop on dates and files and extract the corresponding temporal slice of data
Warning:
Make sure to optimize the opening of files. There is high chances that
the same file has to be open and closed over and over again to loop on the dates.
If this is the case, make sure not to close it between each date.
Args:
name (str): name of the component
varnames (list[str]): original names of variables to read; use `name`
if `varnames` is empty
dates (list): list of the date intervals to extract
files (list): list of the files matching dates
Return:
xr.DataArray: the actual data with dimension:
time, levels, latitudes, longitudes
"""
# Get domain dimensions for random generation
domain = tracer.domain
nlon = domain.nlon
nlat = domain.nlat
nlev = domain.nlev
# Loop over dates/files and import data
data = []
out_dates = []
for dd, ff in zip(dates, files):
debug(
"Reading the file {} for the date interval {}".format(
ff, dd
)
)
# Generate random values instead of reading
data.append(
np.random.normal(
tracer.average_value, tracer.average_value / 2,
(nlev, nlat, nlon)))
out_dates.append(dd[0])
# if only one level for emissions, create the axis:
xmod = xr.DataArray(
np.array(data),
coords={"time": out_dates},
dims=("time", "lev", "lat", "lon"),
)
return xmod
def write(self, name, flx_file, flx, mode="a"):
raise Exception("Can't write template fluxes so far")
...@@ -2,7 +2,7 @@ import datetime ...@@ -2,7 +2,7 @@ import datetime
import glob import glob
import os import os
import pandas as pd import pandas as pd
import xarray as xr
import numpy as np import numpy as np
import glob import glob
...@@ -48,10 +48,20 @@ def fetch(ref_dir, ref_file, date_interval, target_dir, ...@@ -48,10 +48,20 @@ def fetch(ref_dir, ref_file, date_interval, target_dir,
indout = np.where(deltas == delta_max)[0][0] indout = np.where(deltas == delta_max)[0][0]
else: else:
indout = 0 indout = 0
tmp_dates[dd] = [[dd.to_pydatetime(), file = list_files_avail[indout]
datetime.datetime(dd.year + 1, 1, 1)]] ds = xr.open_dataset(file)
tmp_files[dd] = [list_files_avail[indout]]
if "time" in ds:
dates = ds["time"].values
dt = np.unique(np.diff(dates))[0]
tmp_dates[dd] = [[d - dt, d] for d in dates]
tmp_files[dd] = len(tmp_dates[dd]) * [file]
else:
tmp_dates[dd] = [[dd.to_pydatetime(),
datetime.datetime(dd.year + 1, 1, 1)]]
tmp_files[dd] = [file]
# Fetch to datavect # Fetch to datavect
f = list_files_avail[indout] f = list_files_avail[indout]
......
...@@ -35,7 +35,7 @@ def get_domain(ref_dir, ref_file, input_dates, target_dir, tracer=None): ...@@ -35,7 +35,7 @@ def get_domain(ref_dir, ref_file, input_dates, target_dir, tracer=None):
) )
# Read lon/lat in # Read lon/lat in
nc = xr.open_dataset(domain_file, decode_times=False) nc = xr.open_dataset(domain_file)
lon = nc["lon"] lon = nc["lon"]
lat = nc["lat"] lat = nc["lat"]
......
...@@ -25,8 +25,16 @@ def read( ...@@ -25,8 +25,16 @@ def read(
trcr_flx = [] trcr_flx = []
trcr_dates = [] trcr_dates = []
for dd, file_flx in zip(dates, files): for dd, file_flx in zip(dates, files):
nc = xr.open_dataset(file_flx, decode_times=False) nc = xr.open_dataset(file_flx)
trcr_flx.append(nc[varnames].values)
if "time" in nc:
times = nc["time"].values
ind_data = np.where(times == np.datetime64(dd[1]))[0][0]
trcr_flx.append(nc[varnames][ind_data].values)
else:
trcr_flx.append(nc[varnames].values)
trcr_dates.append(dd[0]) trcr_dates.append(dd[0])
return xr.DataArray( return xr.DataArray(
......
def read( def read(
self, self,
name, name,
tracdir, varnames,
tracfile, dates,
varnames, files,
dates, interpol_flx=False,
interpol_flx=False, comp_type=None,
comp_type=None, ddi=None,
**kwargs **kwargs
): ):
raise Exception( raise Exception(
......
...@@ -85,6 +85,22 @@ def init_components(plugin): ...@@ -85,6 +85,22 @@ def init_components(plugin):
"for the component/tracer {}/{}".format(trcr, comp)) "for the component/tracer {}/{}".format(trcr, comp))
raise e raise e
# Clean empty periods
# Make sure that periods are datetime objects
list_files = {
pd.DatetimeIndex([ddi]).to_pydatetime()[0]:
list_files[ddi] for ddi in list_files
# if list_dates[ddi] != []
}
list_dates = {
pd.DatetimeIndex([ddi]).to_pydatetime()[0]:
[[pd.DatetimeIndex([d]).to_pydatetime()[0]
for d in dd] if type(dd) == list else
pd.DatetimeIndex([dd]).to_pydatetime()[0]
for dd in list_dates[ddi]]
for ddi in list_dates
}
# Crop files and dates depending on global datei and datef # Crop files and dates depending on global datei and datef
try: try:
list_files = { list_files = {
...@@ -104,22 +120,6 @@ def init_components(plugin): ...@@ -104,22 +120,6 @@ def init_components(plugin):
except Exception as e: except Exception as e:
warning(e) warning(e)
pass pass
# Clean empty periods
# Make sure that periods are datetime objects
list_files = {
pd.DatetimeIndex([ddi]).to_pydatetime()[0]:
list_files[ddi] for ddi in list_files
# if list_dates[ddi] != []
}
list_dates = {
pd.DatetimeIndex([ddi]).to_pydatetime()[0]:
[[pd.DatetimeIndex([d]).to_pydatetime()[0]
for d in dd] if type(dd) == list else
pd.DatetimeIndex([dd]).to_pydatetime()[0]
for dd in list_dates[ddi]]
for ddi in list_dates
}
# Attach info to the tracer # Attach info to the tracer
tracer.input_dates = list_dates tracer.input_dates = list_dates
......
...@@ -20,7 +20,19 @@ def ini_periods(self, **kwargs): ...@@ -20,7 +20,19 @@ def ini_periods(self, **kwargs):
# List of sub-simulation windows # List of sub-simulation windows
self.subsimu_dates = date_range(datei, datef, period=self.periods) self.subsimu_dates = date_range(datei, datef, period=self.periods)
# Check that subperiods are all of the same length
if np.unique(np.diff(self.subsimu_dates)).size > 1:
raise Exception("Trying to run CHIMERE on a simulation window not fitting the "
"sub-simulation length. \nPlease check the compatibility between "
"the duration of your sub-simulations ({}) and the simulation "
"window ({} to {}): \n"
"Guessed sub-periods: \n".format(self.periods, datei, datef) +
"\n".join([" - {}: {}".format(d0, d1 - d0)
for d0, d1 in zip(self.subsimu_dates[:-1],
self.subsimu_dates[1:])])
)
# Time steps defined by the user # Time steps defined by the user
nphour_ref = self.nphour_ref nphour_ref = self.nphour_ref
...@@ -44,6 +56,10 @@ def ini_periods(self, **kwargs): ...@@ -44,6 +56,10 @@ def ini_periods(self, **kwargs):
self.list_spec_stop = list_spec_stop self.list_spec_stop = list_spec_stop
self.list_thld_stop = list_thld_stop self.list_thld_stop = list_thld_stop
print(__file__)
import code
code.interact(local=dict(locals(), **globals()))
# Loop over sub simulations # Loop over sub simulations
for dd in self.subsimu_dates[:-1]: for dd in self.subsimu_dates[:-1]:
......
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