Commit 41f34822 authored by Ingrid Luijkx's avatar Ingrid Luijkx
Browse files

Added option to transition from 25 to 34 layers for OD meteo after December...

Added option to transition from 25 to 34 layers for OD meteo after December 2005. Added option to read initial random seed from file. Localization method is now read in from rc file. Changed terminology of mixing ratio to mole fraction.
parent f3dbd3b9
#!/usr/bin/env python
# expand_fluxes.py
import sys
import os
import getopt
import shutil
import logging
import netCDF4
import numpy as np
from string import join
from datetime import datetime, timedelta
from pylab import date2num, num2date
sys.path.append('../../')
from da.tools.general import create_dirs
import da.tools.io4 as io
"""
Author: Wouter Peters (Wouter.Peters@wur.nl)
Revision History:
File created on 11 May 2012.
"""
def write_mixing_ratios(dacycle):
"""
Write Sample information to NetCDF files. These files are organized by site and
have an unlimited time axis to which data is appended each cycle.
The needed information is obtained from the sample_auxiliary.nc files and the original input data files from ObsPack.
The steps are:
(1) Create a directory to hold timeseries output files
(2) Read the sample_auxiliary.nc file for this cycle and get a list of original files they were obtained from
(3) For each file, copy the original data file from ObsPack (if not yet present)
(4) Open the copied file, find the index of each observation, fill in the simulated data
"""
dirname = create_dirs(os.path.join(dacycle['dir.analysis'], 'data_molefractions'))
#
# Some help variables
#
dectime0 = date2num(datetime(2000, 1, 1))
dt = dacycle['cyclelength']
startdate = dacycle['time.start']
enddate = dacycle['time.end']
logging.debug("DA Cycle start date is %s" % startdate.strftime('%Y-%m-%d %H:%M'))
logging.debug("DA Cycle end date is %s" % enddate.strftime('%Y-%m-%d %H:%M'))
dacycle['time.sample.stamp'] = "%s_%s" % (startdate.strftime("%Y%m%d%H"), enddate.strftime("%Y%m%d%H"),)
# Step (1): Get the posterior sample output data file for this cycle
infile = os.path.join(dacycle['dir.output'], 'sample_auxiliary_%s.nc' % dacycle['time.sample.stamp'])
ncf_in = io.ct_read(infile, 'read')
obs_num = ncf_in.get_variable('obs_num')
obs_val = ncf_in.get_variable('observed')
simulated = ncf_in.get_variable('modelsamples')
infilename = ncf_in.get_variable('inputfilename')
infiles = netCDF4.chartostring(infilename).tolist()
#infiles = [join(s.compressed(),'') for s in infilename]
ncf_in.close()
# Step (2): Get the prior sample output data file for this cycle
infile = os.path.join(dacycle['dir.output'], 'optimizer.%s.nc' % startdate.strftime('%Y%m%d'))
if os.path.exists(infile):
optimized_present = True
else:
optimized_present = False
if optimized_present:
ncf_fc_in = io.ct_read(infile, 'read')
fc_obs_num = ncf_fc_in.get_variable('obspack_num')
fc_obs_val = ncf_fc_in.get_variable('observed')
fc_simulated = ncf_fc_in.get_variable('modelsamplesmean_prior')
fc_simulated_ens = ncf_fc_in.get_variable('modelsamplesdeviations_prior')
fc_flag = ncf_fc_in.get_variable('flag')
if not dacycle.dasystem.has_key('opt.algorithm'):
fc_r = ncf_fc_in.get_variable('modeldatamismatchvariance')
fc_hphtr = ncf_fc_in.get_variable('totalmolefractionvariance')
elif dacycle.dasystem['opt.algorithm'] == 'serial':
fc_r = ncf_fc_in.get_variable('modeldatamismatchvariance')
fc_hphtr = ncf_fc_in.get_variable('totalmolefractionvariance')
elif dacycle.dasystem['opt.algorithm'] == 'bulk':
fc_r = ncf_fc_in.get_variable('modeldatamismatchvariance').diagonal()
fc_hphtr = ncf_fc_in.get_variable('totalmolefractionvariance').diagonal()
filesitecode = ncf_fc_in.get_variable('sitecode')
fc_sitecodes = netCDF4.chartostring(filesitecode).tolist()
#fc_sitecodes = [join(s.compressed(),'') for s in filesitecode]
ncf_fc_in.close()
# Expand the list of input files with those available from the forecast list
infiles_rootdir = os.path.split(infiles[0])[0]
infiles.extend(os.path.join(infiles_rootdir, f + '.nc') for f in fc_sitecodes)
#Step (2): For each observation timeseries we now have data for, open it and fill with data
for orig_file in set(infiles):
if not os.path.exists(orig_file):
logging.error("The original input file (%s) could not be found, continuing to next file..." % orig_file)
continue
copy_file = os.path.join(dirname, os.path.split(orig_file)[-1])
if not os.path.exists(copy_file):
shutil.copy(orig_file, copy_file)
logging.debug("Copied a new original file (%s) to the analysis directory" % orig_file)
ncf_out = io.CT_CDF(copy_file, 'write')
# Modify the attributes of the file to reflect added data from CTDAS properly
ncf_out.Caution = '==================================================================================='
try:
ncf_out.History += '\nOriginal observation file modified by user %s on %s\n' % (os.environ['USER'], datetime.today().strftime('%F'),)
except:
ncf_out.History = '\nOriginal observation file modified by user %s on %s\n' % (os.environ['USER'], datetime.today().strftime('%F'),)
ncf_out.CTDAS_info = 'Simulated values added from a CTDAS run by %s on %s\n' % (os.environ['USER'], datetime.today().strftime('%F'),)\
+ '\nCTDAS was run on platform %s' % (os.environ['HOST'],)\
+ '\nCTDAS job directory was %s' % (dacycle['dir.da_run'],)\
+ '\nCTDAS Da System was %s' % (dacycle['da.system'],)\
+ '\nCTDAS Da ObsOperator was %s' % (dacycle['da.obsoperator'],)
ncf_out.CTDAS_startdate = dacycle['time.start'].strftime('%F')
ncf_out.CTDAS_enddate = dacycle['time.finish'].strftime("%F")
ncf_out.original_file = orig_file
# get nobs dimension
if ncf_out.dimensions.has_key('id'):
dimidob = ncf_out.dimensions['id']
dimid = ('id',)
elif ncf_out.dimensions.has_key('obs'):
dimidob = ncf_out.dimensions['obs']
dimid = ('obs',)
if dimidob.isunlimited:
nobs = ncf_out.inq_unlimlen()
else:
nobs = len(dimid)
# add nmembers dimension
dimmembersob = ncf_out.createDimension('nmembers', size=simulated.shape[1])
dimmembers = ('nmembers',)
nmembers = len(dimmembers)
# Create empty arrays for posterior samples, as well as for forecast sample statistics
savedict = io.std_savedict.copy()
savedict['name'] = "flag_forecast"
savedict['long_name'] = "flag_for_obs_model in forecast"
savedict['units'] = "None"
savedict['dims'] = dimid
savedict['comment'] = 'Flag (0/1/2/99) for observation value, 0 means okay, 1 means QC error, 2 means rejected, 99 means not sampled'
ncf_out.add_variable(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modeldatamismatch"
savedict['long_name'] = "modeldatamismatch"
savedict['units'] = "[mol mol-1]^2"
savedict['dims'] = dimid
savedict['comment'] = 'Variance of mole fractions resulting from model-data mismatch'
ncf_out.add_variable(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "totalmolefractionvariance_forecast"
savedict['long_name'] = "totalmolefractionvariance of forecast"
savedict['units'] = "[mol mol-1]^2"
savedict['dims'] = dimid
savedict['comment'] = 'Variance of mole fractions resulting from prior state and model-data mismatch'
ncf_out.add_variable(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modelsamplesmean"
savedict['long_name'] = "mean modelsamples"
savedict['units'] = "mol mol-1"
savedict['dims'] = dimid
savedict['comment'] = 'simulated mixing ratios based on optimized state vector'
ncf_out.add_variable(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modelsamplesmean_forecast"
savedict['long_name'] = "mean modelsamples from forecast"
savedict['units'] = "mol mol-1"
savedict['dims'] = dimid
savedict['comment'] = 'simulated mixing ratios based on prior state vector'
ncf_out.add_variable(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modelsamplesstandarddeviation"
savedict['long_name'] = "standard deviaton of modelsamples over all ensemble members"
savedict['units'] = "mol mol-1"
savedict['dims'] = dimid
savedict['comment'] = 'std dev of simulated mixing ratios based on optimized state vector'
ncf_out.add_variable(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modelsamplesstandarddeviation_forecast"
savedict['long_name'] = "standard deviaton of modelsamples from forecast over all ensemble members"
savedict['units'] = "mol mol-1"
savedict['dims'] = dimid
savedict['comment'] = 'std dev of simulated mixing ratios based on prior state vector'
ncf_out.add_variable(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modelsamplesensemble"
savedict['long_name'] = "modelsamples over all ensemble members"
savedict['units'] = "mol mol-1"
savedict['dims'] = dimid + dimmembers
savedict['comment'] = 'ensemble of simulated mixing ratios based on optimized state vector'
ncf_out.add_variable(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modelsamplesensemble_forecast"
savedict['long_name'] = "modelsamples from forecast over all ensemble members"
savedict['units'] = "mol mol-1"
savedict['dims'] = dimid + dimmembers
savedict['comment'] = 'ensemble of simulated mixing ratios based on prior state vector'
ncf_out.add_variable(savedict)
else:
logging.debug("Modifying existing file (%s) in the analysis directory" % copy_file)
ncf_out = io.CT_CDF(copy_file, 'write')
# Get existing file obs_nums to determine match to local obs_nums
if ncf_out.variables.has_key('id'):
file_obs_nums = ncf_out.get_variable('id')
elif ncf_out.variables.has_key('obspack_num'):
file_obs_nums = ncf_out.get_variable('obspack_num')
# Get all obs_nums related to this file, determine their indices in the local arrays
selected_obs_nums = [num for infile, num in zip(infiles, obs_num) if infile == orig_file]
# Optimized data 1st: For each index, get the data and add to the file in the proper file index location
for num in selected_obs_nums:
model_index = obs_num.tolist().index(num)
file_index = file_obs_nums.tolist().index(num)
#var = ncf_out.variables['modeldatamismatch'] # Take from optimizer.yyyymmdd.nc file instead
#var[file_index] = mdm[model_index]
var = ncf_out.variables['modelsamplesmean']
var[file_index] = simulated[model_index, 0]
var = ncf_out.variables['modelsamplesstandarddeviation']
var[file_index] = simulated[model_index, 1:].std()
var = ncf_out.variables['modelsamplesensemble']
var[file_index] = simulated[model_index, :]
# Now forecast data too: For each index, get the data and add to the file in the proper file index location
if optimized_present:
selected_fc_obs_nums = [num for sitecode, num in zip(fc_sitecodes, fc_obs_num) if sitecode in orig_file]
for num in selected_fc_obs_nums:
model_index = fc_obs_num.tolist().index(num)
file_index = file_obs_nums.tolist().index(num)
var = ncf_out.variables['modeldatamismatch']
var[file_index] = np.sqrt(fc_r[model_index])
var = ncf_out.variables['modelsamplesmean_forecast']
var[file_index] = fc_simulated[model_index]
var = ncf_out.variables['modelsamplesstandarddeviation_forecast']
var[file_index] = fc_simulated_ens[model_index, 1:].std()
var = ncf_out.variables['modelsamplesensemble_forecast']
var[file_index] = fc_simulated_ens[model_index, :]
var = ncf_out.variables['totalmolefractionvariance_forecast']
var[file_index] = fc_hphtr[model_index]
var = ncf_out.variables['flag_forecast']
var[file_index] = fc_flag[model_index]
# close the file
status = ncf_out.close()
return None
if __name__ == "__main__":
import logging
from da.tools.initexit import CycleControl
from da.carbondioxide.dasystem import CO2DaSystem
sys.path.append('../../')
logging.root.setLevel(logging.DEBUG)
dacycle = CycleControl(args={'rc':'../../ctdas-od-gfed2-glb6x4-obspack-full.rc'})
dacycle.setup()
dacycle.parse_times()
dasystem = CO2DaSystem('../rc/carbontracker_ct09_opfnew.rc')
dacycle.dasystem = dasystem
while dacycle['time.start'] < dacycle['time.finish']:
outfile = write_mixing_ratios(dacycle)
dacycle.advance_cycle_times()
sys.exit(0)
......@@ -183,7 +183,7 @@ def timehistograms_new(fig, infile):
pydates = np.array([dt.datetime(1970, 1, 1) + dt.timedelta(seconds=int(d)) for d in date])
sampled = (np.ma.getmask(simulated) == False)
sampled = (np.ma.getmaskarray(simulated) == False)
if len(sampled.nonzero()[0]) < 2:
logging.warning("Too few simulated values found, continuing...")
......@@ -376,7 +376,7 @@ def timevssite_new(fig, infile):
f.close()
pydates = np.array([dt.datetime(1970, 1, 1) + dt.timedelta(seconds=int(d)) for d in date])
sampled = (np.ma.getmask(simulated) == False)
sampled = (np.ma.getmaskarray(simulated) == False)
if len(sampled.nonzero()[0]) < 2:
logging.warning("Too few simulated values found, continuing...")
......@@ -589,7 +589,7 @@ def residuals_new(fig, infile):
f.close()
pydates = np.array([dt.datetime(1970, 1, 1) + dt.timedelta(seconds=int(d)) for d in date])
sampled = (np.ma.getmask(simulated) == False)
sampled = (np.ma.getmaskarray(simulated) == False)
if len(sampled.nonzero()[0]) < 2:
logging.warning("Too few simulated values found, continuing...")
......
......@@ -205,7 +205,7 @@ class Optimizer(object):
savedict['units'] = "mol mol-1"
savedict['dims'] = dimobs
savedict['values'] = self.Hx.tolist()
savedict['comment'] = '%s mean mixing ratios based on %s state vector' % (type, type)
savedict['comment'] = '%s mean mole fractions based on %s state vector' % (type, type)
f.add_data(savedict)
savedict = io.std_savedict.copy()
......@@ -214,7 +214,7 @@ class Optimizer(object):
savedict['units'] = "mol mol-1"
savedict['dims'] = dimobs + dimmembers
savedict['values'] = self.HX_prime.tolist()
savedict['comment'] = '%s mixing ratio deviations based on %s state vector' % (type, type)
savedict['comment'] = '%s mole fraction deviations based on %s state vector' % (type, type)
f.add_data(savedict)
# Continue with prior only data
......@@ -380,7 +380,7 @@ class Optimizer(object):
self.X_prime[:, :] = np.dot(I, self.X_prime) - np.dot(Kw, self.HX_prime) # HX' = I - K~ * HX'
# Now do the adjustments of the modeled mixing ratios using the linearized ensemble. These are not strictly needed but can be used
# Now do the adjustments of the modeled mole fractions using the linearized ensemble. These are not strictly needed but can be used
# for diagnosis.
part3 = np.dot(HPH, np.transpose(la.inv(sHPHR))) # HPH(sqrt(HPH+R))^-1
......
......@@ -19,7 +19,7 @@ from numpy import array, logical_and
sys.path.append(os.getcwd())
sys.path.append('../../')
identifier = 'CarbonTracker CO2 mixing ratios'
identifier = 'CarbonTracker CO2 mole fractions'
version = '0.0'
from da.baseclasses.obs import Observations
......@@ -29,7 +29,7 @@ import da.tools.rc as rc
################### Begin Class CO2Observations ###################
class CO2Observations(Observations):
""" an object that holds data + methods and attributes needed to manipulate mixing ratio values """
""" an object that holds data + methods and attributes needed to manipulate mole fraction values """
def setup(self, dacycle):
self.startdate = dacycle['time.sample.start']
......@@ -50,9 +50,9 @@ class CO2Observations(Observations):
self.datalist = []
def add_observations(self):
""" Returns a MixingRatioList holding individual MixingRatioSample objects for all obs in a file
""" Returns a MoleFractionList holding individual MoleFractionSample objects for all obs in a file
The CarbonTracker mixing ratio files are provided as one long list of obs for all possible dates. So we can
The CarbonTracker mole fraction files are provided as one long list of obs for all possible dates. So we can
either:
(1) read all, and the subselect the data we will use in the rest of this cycle
......@@ -94,11 +94,11 @@ class CO2Observations(Observations):
logging.debug("Successfully read data from obs file (%s)" % self.obs_filename)
for n in range(len(dates)):
self.datalist.append(MixingRatioSample(ids[n], dates[n], sites[n], obs[n], 0.0, 0.0, 0.0, 0.0, flags[n], alts[n], lats[n], lons[n], evn[n], species[n], strategy[n], 0.0, self.obs_filename))
self.datalist.append(MoleFractionSample(ids[n], dates[n], sites[n], obs[n], 0.0, 0.0, 0.0, 0.0, flags[n], alts[n], lats[n], lons[n], evn[n], species[n], strategy[n], 0.0, self.obs_filename))
logging.debug("Added %d observations to the Data list" % len(dates))
def add_simulations(self, filename, silent=True):
""" Adds model simulated values to the mixing ratio objects """
""" Adds model simulated values to the mole fraction objects """
if not os.path.exists(filename):
......@@ -383,7 +383,7 @@ class CO2Observations(Observations):
savedict['units'] = "mol mol-1"
savedict['dims'] = dimid + dimmembers
savedict['values'] = data.tolist()
savedict['comment'] = 'simulated mixing ratios based on optimized state vector'
savedict['comment'] = 'simulated mole fractions based on optimized state vector'
f.add_data(savedict)
data = self.getvalues('fromfile')
......@@ -419,11 +419,11 @@ class CO2Observations(Observations):
################### Begin Class MixingRatioSample ###################
################### Begin Class MoleFractionSample ###################
class MixingRatioSample(object):
class MoleFractionSample(object):
"""
Holds the data that defines a Mixing Ratio Sample in the data assimilation framework. Sor far, this includes all
Holds the data that defines a mole fraction Sample in the data assimilation framework. Sor far, this includes all
attributes listed below in the __init__ method. One can additionally make more types of data, or make new
objects for specific projects.
......@@ -434,8 +434,8 @@ class MixingRatioSample(object):
self.xdate = xdate # Date of obs
self.obs = obs # Value observed
self.simulated = simulated # Value simulated by model
self.resid = resid # Mixing ratio residuals
self.hphr = hphr # Mixing ratio prior uncertainty from fluxes and (HPH) and model data mismatch (R)
self.resid = resid # Mole fraction residuals
self.hphr = hphr # Mole fraction prior uncertainty from fluxes and (HPH) and model data mismatch (R)
self.mdm = mdm # Model data mismatch
self.may_localize = True # Whether sample may be localized in optimizer
self.may_reject = True # Whether sample may be rejected if outside threshold
......@@ -452,7 +452,7 @@ class MixingRatioSample(object):
self.samplingstrategy = samplingstrategy
self.fromfile = fromfile # netcdf filename inside observation distribution, to write back later
################### End Class MixingRatioSample ###################
################### End Class MoleFractionSample ###################
if __name__ == "__main__":
......
......@@ -23,7 +23,7 @@ import da.tools.rc as rc
from da.baseclasses.obs import Observations
from da.tools.general import name_convert
identifier = 'CarbonTracker CO2 mixing ratios'
identifier = 'CarbonTracker CO2 mole fractions'
version = '0.0'
......@@ -31,7 +31,7 @@ version = '0.0'
################### Begin Class ObsPackObservations ###################
class ObsPackObservations(Observations):
""" an object that holds data + methods and attributes needed to manipulate mixing ratio values """
""" an object that holds data + methods and attributes needed to manipulate mole fraction values """
def setup(self, dacycle):
self.startdate = dacycle['time.sample.start']
......@@ -51,9 +51,9 @@ class ObsPackObservations(Observations):
self.datalist = []
def add_observations(self):
""" Returns a MixingRatioList holding individual MixingRatioSample objects for all obs in a file
""" Returns a MoleFractionList holding individual MoleFractionSample objects for all obs in a file
The ObsPack mixing ratio files are provided as time series per site with all dates in sequence.
The ObsPack mole fraction files are provided as time series per site with all dates in sequence.
We will loop over all site files in the ObsPackage, and subset each to our needs
"""
......@@ -106,14 +106,14 @@ class ObsPackObservations(Observations):
ncf.close()
for n in range(len(dates)):
self.datalist.append(MixingRatioSample(ids[n], dates[n], site, obs[n], 0.0, 0.0, 0.0, 0.0, flags[n], alts[n], lats[n], lons[n], evn[n], species, strategy, 0.0, infile))
self.datalist.append(MoleFractionSample(ids[n], dates[n], site, obs[n], 0.0, 0.0, 0.0, 0.0, flags[n], alts[n], lats[n], lons[n], evn[n], species, strategy, 0.0, infile))
logging.debug("Added %d observations from file (%s) to the Data list" % (len(dates), ncfile))
logging.info("Observations list now holds %d values" % len(self.datalist))
def add_simulations(self, filename, silent=False):
""" Adds model simulated values to the mixing ratio objects """
""" Adds model simulated values to the mole fraction objects """
if not os.path.exists(filename):
msg = "Sample output filename for observations could not be found : %s" % filename
......@@ -405,7 +405,7 @@ class ObsPackObservations(Observations):
savedict['units'] = "mol mol-1"
savedict['dims'] = dimid + dimmembers
savedict['values'] = data.tolist()
savedict['comment'] = 'simulated mixing ratios based on optimized state vector'
savedict['comment'] = 'simulated mole fractions based on optimized state vector'
f.add_data(savedict)
data = self.getvalues('fromfile')
......@@ -430,11 +430,11 @@ class ObsPackObservations(Observations):
################### Begin Class MixingRatioSample ###################
################### Begin Class MoleFractionSample ###################
class MixingRatioSample(object):
class MoleFractionSample(object):
"""
Holds the data that defines a Mixing Ratio Sample in the data assimilation framework. Sor far, this includes all
Holds the data that defines a mole fraction Sample in the data assimilation framework. Sor far, this includes all
attributes listed below in the __init__ method. One can additionally make more types of data, or make new
objects for specific projects.
......@@ -445,8 +445,8 @@ class MixingRatioSample(object):
self.xdate = xdate # Date of obs
self.obs = obs # Value observed
self.simulated = simulated # Value simulated by model
self.resid = resid # Mixing ratio residuals
self.hphr = hphr # Mixing ratio prior uncertainty from fluxes and (HPH) and model data mismatch (R)
self.resid = resid # Mole fraction residuals
self.hphr = hphr # Mole fraction prior uncertainty from fluxes and (HPH) and model data mismatch (R)
self.mdm = mdm # Model data mismatch
self.may_localize = True # Whether sample may be localized in optimizer
self.may_reject = True # Whether sample may be rejected if outside threshold
......@@ -463,7 +463,7 @@ class MixingRatioSample(object):
self.samplingstrategy = samplingstrategy
self.fromfile = fromfile # netcdf filename inside ObsPack distribution, to write back later
################### End Class MixingRatioSample ###################
################### End Class MoleFractionSample ###################
if __name__ == "__main__":
......
......@@ -18,7 +18,7 @@ from numpy import array, logical_and
sys.path.append(os.getcwd())
sys.path.append('../../')
identifier = 'CarbonTracker CO2 mixing ratios'
identifier = 'CarbonTracker CO2 mole fractions'
version = '0.0'
from da.baseclasses.obs import Observations
......@@ -27,7 +27,7 @@ import da.tools.rc as rc
################### Begin Class ObsPackObservations ###################
class ObsPackObservations(Observations):
""" an object that holds data + methods and attributes needed to manipulate mixing ratio values """
""" an object that holds data + methods and attributes needed to manipulate mole fraction values """
def setup(self, dacycle):
......@@ -48,9 +48,9 @@ class ObsPackObservations(Observations):
self.datalist = []
def add_observations(self):
""" Returns a MixingRatioList holding individual MixingRatioSample objects for all obs in a file
""" Returns a MoleFractionList holding individual MoleFractionSample objects for all obs in a file
The ObsPack mixing ratio files are provided as time series per site with all dates in sequence.
The ObsPack mole fraction files are provided as time series per site with all dates in sequence.
We will loop over all site files in the ObsPackage, and subset each to our needs
"""
......@@ -100,14 +100,14 @@ class ObsPackObservations(Observations):
ncf.close()
for n in range(len(dates)):
self.datalist.append(MixingRatioSample(obspacknum[n], dates[n], datasetname, obs[n], 0.0, 0.0, 0.0