diff --git a/da/stilt/obspack_geocarbon.py b/da/stilt/obspack_geocarbon.py
deleted file mode 100755
index 4dc0fae2cd5a24d72ecb2d702b6cb51ebdd49d27..0000000000000000000000000000000000000000
--- a/da/stilt/obspack_geocarbon.py
+++ /dev/null
@@ -1,638 +0,0 @@
-#!/usr/bin/env python
-
-# obs.py
-
-"""
-Author : peters
-
-Revision History:
-File created on 28 Jul 2010.
-
-"""
-import os
-import sys
-import logging
-
-import datetime as dtm
-from string import strip
-from numpy import array, logical_and
-sys.path.append(os.getcwd())
-sys.path.append('../../')
-
-identifier = 'CarbonTracker CO2 mole fractions'
-version = '0.0'
-
-from da.baseclasses.obs import Observations
-import da.tools.io4 as io
-import da.tools.rc as rc
-################### Begin Class ObsPackObservations ###################
-
-class ObsPackObservations(Observations):
- """ an object that holds data + methods and attributes needed to manipulate mole fraction values """
-
- def setup(self, dacycle):
-
- self.startdate = dacycle['time.sample.start']
- self.enddate = dacycle['time.sample.end']
-
- op_id = dacycle.dasystem['obspack.input.id']
- op_dir = dacycle.dasystem['obspack.input.dir']
-
- if not os.path.exists(op_dir):
- msg = 'Could not find the required ObsPack distribution (%s) ' % op_dir
- logging.error(msg)
- raise IOError, msg
- else:
- self.obspack_dir = op_dir
- self.obspack_id = op_id
-
- self.datalist = []
-
- def add_observations(self):
- """ Returns a MoleFractionList holding individual MoleFractionSample objects for all obs in a file
-
- 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
-
- """
- import string
-
- # Step 1: Read list of available site files in package
-
- infile = os.path.join(self.obspack_dir, 'summary', '%s_dataset_summary.txt' % (self.obspack_id,))
- f = open(infile, 'r')
- lines = f.readlines()
- f.close()
-
- ncfilelist = []
- for line in lines:
- if line.startswith('#'): continue # header
-
- items = line.split()
- #ncfile, lab , start_date, stop_date, data_comparison = items[0:5]
- ncfile, lab , start_date, stop_date, data_comparison= line[:105].split()
-
-
- ncfilelist += [ncfile]
-
- logging.debug("ObsPack dataset info read, proceeding with %d netcdf files" % len(ncfilelist))
-
- for ncfile in ncfilelist:
- logging.info('ncfile %s'%ncfile)
- infile = os.path.join(self.obspack_dir, 'data', 'nc', ncfile + '.nc')
- ncf = io.ct_read(infile, 'read')
- idates = ncf.get_variable('time_components')
- dates = array([dtm.datetime(*d) for d in idates])
-
- subselect = logical_and(dates >= self.startdate , dates <= self.enddate).nonzero()[0]
-
- if len(subselect) == 0:
- ncf.close()
- continue
-
- logging.debug("Trying to add %d observations from file (%s) to the Data list" % (len(subselect), ncfile))
-
- dates = dates.take(subselect, axis=0)
-
- #ccgg_evn = ncf.get_variable('obspack_num').take(subselect) # or should we propagate obs_num which is not unique across datasets??
- ccgg_evn = ncf.get_variable('ccgg_evn').take(subselect,axis=0)
-
- obspackid = ncf.get_variable('obspack_id').take(subselect, axis=0)
- obspackid = [s.tostring().lower() for s in obspackid]
- obspackid = map(strip, obspackid)
-
- datasetname = ncfile # use full name of dataset to propagate for clarity
- lats = ncf.get_variable('latitude').take(subselect, axis=0)
- lons = ncf.get_variable('longitude').take(subselect, axis=0)
- alts = ncf.get_variable('altitude').take(subselect, axis=0)
- obs = ncf.get_variable('value').take(subselect, axis=0)
- species = ncf.get_attribute('dataset_parameter')
- utc2loc = ncf.get_attribute('site_utc2lst')
- flags = ncf.get_variable('obs_flag').take(subselect, axis=0)
- ncf.close()
-
- for n in range(len(dates)):
- used_ids = self.getvalues('id',list)
- if ccgg_evn[n] in used_ids:
- ii = used_ids.index(ccgg_evn[n])
- logging.error("Error when reading from file: %s"%ncfile)
- logging.error("This sample ID (%d) is not unique"%ccgg_evn[n])
- logging.error("Previously used from file: %s"%self.datalist[ii].fromfile)
- logging.error("...skipping")
- raise IOError
- else:
- self.datalist.append(MoleFractionSample(ccgg_evn[n], dates[n], datasetname, obs[n], 0.0, 0.0, 0.0, 0.0, flags[n], alts[n], lats[n], lons[n], obspackid[n], species, utc2loc, 1, 0.0, infile))
-
- logging.info("Observations list now holds %d values" % len(self.datalist))
-
-
- # for k in range(len(dates_foot)):
- # logging.info("k %s" %(k))
- # used_ids = self.getvalues('id',list)
- # if ccgg_evn_foot[k] in used_ids:
- # ii = used_ids.index(ccgg_evn_foot[k])
- # logging.error("Error when reading from file: %s"%ncfile)
- # logging.error("This sample ID (%d) is not unique"%ccgg_evn_foot[k])
- # logging.error("Previously used from file: %s"%self.datalist[ii].fromfile)
- # logging.error("...skipping")
- # raise IOError
- # else:
- # self.datalist.append(MoleFractionSample(ccgg_evn_foot[k], dates_foot[k], datasetname, obs_foot[k], 0.0, 0.0, 0.0, 0.0, flags_foot[k], alts_foot[k], lats_foot[k], lons_foot[k], obspackid_foot[k], species, utc2loc,1, 0.0, infile))
-
- # obs=self.getvalues('obs')
- # logging.info("Obs %s" % obs)
- # #logging.info("datalist %s" % self.datalist[1].obs)
- # logging.info("Observations list now holds %d values" % len(self.datalist))
- #logging.info("Observations list now holds %d values" % len(self.datalist))
-
- def add_simulations(self, filename, silent=False):
- """ 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
- logging.error(msg)
- logging.error("Did the sampling step succeed?")
- logging.error("...exiting")
- raise IOError, msg
-
- ncf = io.ct_read(filename, method='read')
- ids = ncf.get_variable('obs_num')
- simulated = ncf.get_variable('flask')
- ncf.close()
- logging.info("Successfully read data from model sample file (%s)" % filename)
-
- obs_ids = self.getvalues('id',list)
- ids = map(int, ids)
-
- missing_samples = []
-
- for idx, val in zip(ids, simulated):
- if idx in obs_ids:
- index = obs_ids.index(idx)
- self.datalist[index].simulated = val
- else:
- missing_samples.append(idx)
-
- if not silent and missing_samples != []:
- logging.warning('Model samples were found that did not match any ID in the observation list. Skipping them...')
- #msg = '%s'%missing_samples ; logging.warning(msg)
-
- logging.debug("Added %d simulated values to the Data list" % (len(ids) - len(missing_samples)))
-
- def write_sample_coords(self, obsinputfile,proc,nprocs):
- """
- Write the information needed by the observation operator to a file. Return the filename that was written for later use
-
- """
-
- if len(self.datalist) == 0:
- #f.close()
- #return obsinputfile
- logging.debug("No observations found for this time period, nothing written to obs file")
- else:
- frac_even = len(self.datalist)/nprocs
- frac_odd = len(self.datalist) - nprocs*frac_even
- f = io.CT_CDF(obsinputfile, method='create')
- logging.debug('Creating new observations file for ObservationOperator (%s)' % obsinputfile)
-
- if proc < nprocs-1:
- dimid = f.add_dim('obs', len(self.datalist[proc*frac_even:(proc+1)*frac_even]))
- data_id = self.getvalues('id')[proc*frac_even:(proc+1)*frac_even]
- data_xdate = self.getvalues('xdate')[proc*frac_even:(proc+1)*frac_even]
- data_lat = self.getvalues('lat')[proc*frac_even:(proc+1)*frac_even]
- data_lon = self.getvalues('lon')[proc*frac_even:(proc+1)*frac_even]
- data_height = self.getvalues('height')[proc*frac_even:(proc+1)*frac_even]
- data_ss = self.getvalues('samplingstrategy')[proc*frac_even:(proc+1)*frac_even]
- data_evn = self.getvalues('evn')[proc*frac_even:(proc+1)*frac_even]
- data_obs = self.getvalues('obs')[proc*frac_even:(proc+1)*frac_even]
- data_mdm = self.getvalues('mdm')[proc*frac_even:(proc+1)*frac_even]
-
-
- if proc == nprocs-1:
- dimid = f.add_dim('obs', len(self.datalist[proc*frac_even:]))
- data_id = self.getvalues('id')[proc*frac_even:]
- data_xdate = self.getvalues('xdate')[proc*frac_even:]
- data_lat = self.getvalues('lat')[proc*frac_even:]
- data_lon = self.getvalues('lon')[proc*frac_even:]
- data_height = self.getvalues('height')[proc*frac_even:]
- data_ss = self.getvalues('samplingstrategy')[proc*frac_even:]
- data_evn = self.getvalues('evn')[proc*frac_even:]
- data_obs = self.getvalues('obs')[proc*frac_even:]
- data_mdm = self.getvalues('mdm')[proc*frac_even:]
-
-
-
-
- dim200char = f.add_dim('string_of200chars', 200)
- dim10char = f.add_dim('string_of10chars', 10)
- dimcalcomp = f.add_dim('calendar_components', 6)
-
- for key, value in self.site_move.iteritems():
- msg = "Site is moved by %3.2f degrees latitude and %3.2f degrees longitude" % value
- f.add_attribute(key, msg)
-
- data = data_id
-
- savedict = io.std_savedict.copy()
- savedict['name'] = "obs_num"
- savedict['dtype'] = "int"
- savedict['long_name'] = "Unique_Dataset_observation_index_number"
- savedict['units'] = ""
- savedict['dims'] = dimid
- savedict['values'] = data.tolist()
- savedict['comment'] = "Unique index number within this dataset ranging from 0 to UNLIMITED."
- f.add_data(savedict)
-
- data = [[d.year, d.month, d.day, d.hour, d.minute, d.second] for d in data_xdate ]
-
- savedict = io.std_savedict.copy()
- savedict['dtype'] = "int"
- savedict['name'] = "date_components"
- savedict['units'] = "integer components of UTC date/time"
- savedict['dims'] = dimid + dimcalcomp
- savedict['values'] = data
- savedict['missing_value'] = -9
- savedict['comment'] = "Calendar date components as integers. Times and dates are UTC."
- savedict['order'] = "year, month, day, hour, minute, second"
- f.add_data(savedict)
-
- data = data_lat
-
- savedict = io.std_savedict.copy()
- savedict['name'] = "latitude"
- savedict['units'] = "degrees_north"
- savedict['dims'] = dimid
- savedict['values'] = data.tolist()
- savedict['missing_value'] = -999.9
- f.add_data(savedict)
-
- data = data_lon
-
- savedict = io.std_savedict.copy()
- savedict['name'] = "longitude"
- savedict['units'] = "degrees_east"
- savedict['dims'] = dimid
- savedict['values'] = data.tolist()
- savedict['missing_value'] = -999.9
- f.add_data(savedict)
-
- data = data_height
-
- savedict = io.std_savedict.copy()
- savedict['name'] = "altitude"
- savedict['units'] = "meters_above_sea_level"
- savedict['dims'] = dimid
- savedict['values'] = data.tolist()
- savedict['missing_value'] = -999.9
- f.add_data(savedict)
-
- data = data_ss
-
- savedict = io.std_savedict.copy()
- savedict['dtype'] = "int"
- savedict['name'] = "sampling_strategy"
- savedict['units'] = "NA"
- savedict['dims'] = dimid
- savedict['values'] = data.tolist()
- savedict['missing_value'] = -9
- f.add_data(savedict)
-
- data = data_evn
-
- savedict = io.std_savedict.copy()
- savedict['dtype'] = "char"
- savedict['name'] = "obs_id"
- savedict['units'] = "ObsPack datapoint identifier"
- savedict['dims'] = dimid + dim200char
- savedict['values'] = data
- savedict['missing_value'] = '!'
- f.add_data(savedict)
-
- data = data_obs
-
- savedict = io.std_savedict.copy()
- savedict['name'] = "observed"
- savedict['long_name'] = "observedvalues"
- savedict['units'] = "mol mol-1"
- savedict['dims'] = dimid
- savedict['values'] = data.tolist()
- savedict['comment'] = 'Observations used in optimization'
- f.add_data(savedict)
-
- data = data_mdm
-
- savedict = io.std_savedict.copy()
- savedict['name'] = "modeldatamismatch"
- savedict['long_name'] = "modeldatamismatch"
- savedict['units'] = "[mol mol-1]"
- savedict['dims'] = dimid
- savedict['values'] = data.tolist()
- savedict['comment'] = 'Standard deviation of mole fractions resulting from model-data mismatch'
- f.add_data(savedict)
-
- f.close()
-
- logging.debug("Successfully wrote data to obs file")
- logging.info("Sample input file for obs operator now in place [%s]" % obsinputfile)
-
-
-
- def add_model_data_mismatch(self, filename):
- """
- Get the model-data mismatch values for this cycle.
-
- (1) Open a sites_weights file
- (2) Parse the data
- (3) Compare site list against data
- (4) Take care of double sites, etc
-
- """
-
- if not os.path.exists(filename):
- msg = 'Could not find the required sites.rc input file (%s) ' % filename
- logging.error(msg)
- raise IOError, msg
- else:
- self.sites_file = filename
-
- sites_weights = rc.read(self.sites_file)
-
- self.rejection_threshold = int(sites_weights['obs.rejection.threshold'])
- self.global_R_scaling = float(sites_weights['global.R.scaling'])
- self.n_site_categories = int(sites_weights['n.site.categories'])
-
- logging.debug('Model-data mismatch rejection threshold: %d ' % self.rejection_threshold)
- logging.warning('Model-data mismatch scaling factor : %f ' % self.global_R_scaling)
- logging.debug('Model-data mismatch site categories : %d ' % self.n_site_categories)
-
- cats = [k for k in sites_weights.keys() if 'site.category' in k]
-
- site_categories = {}
- for key in cats:
- name, error_summer, error_winter , may_localize, may_reject = sites_weights[key].split(';')
- name = name.strip().lower()
- error_summer = float(error_summer)
- error_winter = float(error_winter)
- may_reject = ("TRUE" in may_reject.upper())
- may_localize = ("TRUE" in may_localize.upper())
- site_categories[name] = {'category': name, 'error_summer': error_summer, 'error_winter': error_winter ,'may_localize': may_localize, 'may_reject': may_reject}
-
- site_info = {}
- site_move = {}
- site_hourly = {} # option added to include only certain hours of the day (for e.g. PAL) IvdL
- site_foot = {} # option added to check for available footprints per observation
- site_incalt = {} # option to increase sampling altitude for sites specified in sites and weights file
- for key, value in sites_weights.iteritems():
- if 'co2_' in key or 'sf6' in key: # to be fixed later, do not yet know how to parse valid keys from rc-files yet.... WP
- sitename, sitecategory = key, value
- sitename = sitename.strip()
- sitecategory = sitecategory.split()[0].strip().lower()
- site_info[sitename] = site_categories[sitecategory]
- if 'site.move' in key:
- identifier, latmove, lonmove = value.split(';')
- site_move[identifier.strip()] = (float(latmove), float(lonmove))
- if 'site.hourly' in key:
- identifier, hourfrom, hourto = value.split(';')
- site_hourly[identifier.strip()] = (int(hourfrom), int(hourto))
- if 'site.foot' in key:
- identifier, foot = value.split(';')
- site_foot[identifier.strip()] = (str(foot))
- if 'site.incalt' in key:
- identifier, incalt = value.split(';')
- site_incalt[identifier.strip()] = (int(incalt))
-
- #for obs in self.datalist:
- do_not_simulate=[]
- do_simulate=[]
- for i,obs in enumerate(self.datalist):
- obs.mdm = 1000.0 # default is very high model-data-mismatch, until explicitly set by script
- obs.flag = 99 # default is do-not-use , until explicitly set by script
- exclude_hourly = False # default is that hourly values are not included
- exclude_footprint = False
-
- identifier = obs.code
- species, site, method, lab, datasetnr = identifier.split('_')
-
- if site_info.has_key(identifier):
- if site_foot.has_key(identifier):
- path_foot = site_foot[identifier]
- dir_foot = os.path.join(path_foot,'%s'%obs.xdate.year,'%02d'%obs.xdate.month)
- files_foot = os.listdir(dir_foot)
- str_id = '%s'%obs.id
- if any(str_id in s for s in files_foot):
- logging.info("id in footprint %s" %str_id)
- else:
- exclude_footprint = True
- logging.info("id not in footprint %s" %str_id)
- if site_hourly.has_key(identifier):
- obs.samplingstrategy = 2
- hourf, hourt = site_hourly[identifier]
- if int(obs.xdate.hour+obs.utc2loc) > hourf and int(obs.xdate.hour+obs.utc2loc) < hourt:
- logging.warning("Observations in hourly dataset INCLUDED, while UTC sampling time %s was between local %s:00-%s:00"%(obs.xdate.time(),hourf,hourt))
- else:
- logging.warning("Observation in hourly dataset EXCLUDED, while UTC sampling time %s was outside local %s:00-%s:00"%(obs.xdate.time(),hourf,hourt))
- exclude_hourly = True
- if site_info[identifier]['category'] == 'do-not-use' or site_info[identifier]['category'] == 'do-not-simulate' or exclude_hourly or exclude_footprint:
- logging.warning("Site found (%s, %d), but data point not used in assimilation !!!" % (identifier, obs.id))
- if int(obs.xdate.month) > 3 and int(obs.xdate.month) < 10:
- obs.mdm = site_info[identifier]['error_summer'] * self.global_R_scaling
- else:
- obs.mdm = site_info[identifier]['error_winter'] * self.global_R_scaling
- obs.may_localize = site_info[identifier]['may_localize']
- obs.may_reject = site_info[identifier]['may_reject']
- obs.flag = 99
- if site_info[identifier]['category'] == 'do-not-simulate':
- do_not_simulate.append(i)
- else:
- logging.debug("Site found (%s, %d)" % (identifier, obs.id))
- if int(obs.xdate.month) > 3 and int(obs.xdate.month) < 10:
- obs.mdm = site_info[identifier]['error_summer'] * self.global_R_scaling
- else:
- obs.mdm = site_info[identifier]['error_winter'] * self.global_R_scaling
- obs.may_localize = site_info[identifier]['may_localize']
- obs.may_reject = site_info[identifier]['may_reject']
- obs.flag = 0
- do_simulate.append(i)
- else:
- logging.warning("Site NOT found (%s, %d), please check sites.rc file (%s) !!!" % (identifier, obs.id, self.sites_file))
-
- if site_move.has_key(identifier):
-
- movelat, movelon = site_move[identifier]
- obs.lat = obs.lat + movelat
- obs.lon = obs.lon + movelon
-
- logging.warning("Observation location for (%s, %d), is moved by %3.2f degrees latitude and %3.2f degrees longitude" % (identifier, obs.id, movelat, movelon))
-
- if site_incalt.has_key(identifier):
-
- incalt = site_incalt[identifier]
- obs.height = obs.height + incalt
-
- logging.warning("Observation location for (%s, %d), is moved by %3.2f meters in altitude" % (identifier, obs.id, incalt))
-
- if len(do_not_simulate) > 0 :
- logging.info("do-not-simulate flags located")
- self.datalist = [self.datalist[k] for k in do_simulate]
- logging.info("After do-not-simulate filter observations list now holds %d values" % len(self.datalist))
- else:
- logging.info("No do-not-simulate flags, continues normally ")
-
- # Add site_info dictionary to the Observations object for future use
-
- self.site_info = site_info
- self.site_move = site_move
- self.site_hourly = site_hourly
- self.site_foot = site_foot
- self.site_incalt = site_incalt
-
- logging.debug("Added Model Data Mismatch to all samples ")
-
- def write_sample_auxiliary(self, auxoutputfile):
- """
- Write selected information contained in the Observations object to a file.
-
- """
-
- f = io.CT_CDF(auxoutputfile, method='create')
- logging.debug('Creating new auxiliary sample output file for postprocessing (%s)' % auxoutputfile)
-
- dimid = f.add_dim('obs', len(self.datalist))
- dim200char = f.add_dim('string_of200chars', 200)
- dim10char = f.add_dim('string_of10chars', 10)
- dimcalcomp = f.add_dim('calendar_components', 6)
-
- if len(self.datalist) == 0:
- f.close()
- #return outfile
-
- for key, value in self.site_move.iteritems():
- msg = "Site is moved by %3.2f degrees latitude and %3.2f degrees longitude" % value
- f.add_attribute(key, msg)
-
- data = self.getvalues('id')
-
- savedict = io.std_savedict.copy()
- savedict['name'] = "obs_num"
- savedict['dtype'] = "int"
- savedict['long_name'] = "Unique_Dataset_observation_index_number"
- savedict['units'] = ""
- savedict['dims'] = dimid
- savedict['values'] = data.tolist()
- savedict['comment'] = "Unique index number within this dataset ranging from 0 to UNLIMITED."
- f.add_data(savedict)
-
- data = [[d.year, d.month, d.day, d.hour, d.minute, d.second] for d in self.getvalues('xdate')]
-
- savedict = io.std_savedict.copy()
- savedict['dtype'] = "int"
- savedict['name'] = "date_components"
- savedict['units'] = "integer components of UTC date/time"
- savedict['dims'] = dimid + dimcalcomp
- savedict['values'] = data
- savedict['missing_value'] = -9
- savedict['comment'] = "Calendar date components as integers. Times and dates are UTC."
- savedict['order'] = "year, month, day, hour, minute, second"
- f.add_data(savedict)
-
- data = self.getvalues('obs')
-
- savedict = io.std_savedict.copy()
- savedict['name'] = "observed"
- savedict['long_name'] = "observedvalues"
- savedict['units'] = "mol mol-1"
- savedict['dims'] = dimid
- savedict['values'] = data.tolist()
- savedict['comment'] = 'Observations used in optimization'
- f.add_data(savedict)
-
- data = self.getvalues('mdm')
-
- savedict = io.std_savedict.copy()
- savedict['name'] = "modeldatamismatch"
- savedict['long_name'] = "modeldatamismatch"
- savedict['units'] = "[mol mol-1]"
- savedict['dims'] = dimid
- savedict['values'] = data.tolist()
- savedict['comment'] = 'Standard deviation of mole fractions resulting from model-data mismatch'
- f.add_data(savedict)
-
- data = self.getvalues('simulated')
- #logging.info("Data (%s)" %data)
- dimmembers = f.add_dim('members', data.shape[1])
-
- savedict = io.std_savedict.copy()
- savedict['name'] = "modelsamples"
- savedict['long_name'] = "modelsamples for all ensemble members"
- savedict['units'] = "mol mol-1"
- savedict['dims'] = dimid + dimmembers
- savedict['values'] = data.tolist()
- savedict['comment'] = 'simulated mole fractions based on optimized state vector'
- f.add_data(savedict)
-
- data = self.getvalues('fromfile')
-
- savedict = io.std_savedict.copy()
- savedict['name'] = "inputfilename"
- savedict['long_name'] = "name of file where original obs data was taken from"
- savedict['dtype'] = "char"
- savedict['dims'] = dimid + dim200char
- savedict['values'] = data
- savedict['missing_value'] = '!'
- f.add_data(savedict)
-
- f.close()
-
- logging.debug("Successfully wrote data to auxiliary sample output file (%s)" % auxoutputfile)
-
- #return outfile
-
-
-
-################### End Class CtObservations ###################
-
-
-
-################### Begin Class MoleFractionSample ###################
-
-class MoleFractionSample(object):
- """
- 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.
-
- """
-
- def __init__(self, idx, xdate, code='XXX', obs=0.0, simulated=0.0, resid=0.0, hphr=0.0, mdm=0.0, flag=0, height=0.0, lat= -999., lon= -999., evn='0000', species='co2',utc2loc=0.0, samplingstrategy=1, sdev=0.0, fromfile='none.nc'):
- self.code = code.strip() # dataset identifier, i.e., co2_lef_tower_insitu_1_99
- self.xdate = xdate # Date of obs
- self.obs = obs # Value observed
- self.simulated = simulated # Value simulated by model
- 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
- self.flag = flag # Flag
- self.height = height # Sample height in masl
- self.lat = lat # Sample lat
- self.lon = lon # Sample lon
- self.id = idx # Obspack ID within distrution (integer), e.g., 82536
- self.evn = evn # Obspack Number within distrution (string), e.g., obspack_co2_1_PROTOTYPE_v0.9.2_2012-07-26_99_82536
- self.sdev = sdev # standard deviation of ensemble
- self.masl = True # Sample is in Meters Above Sea Level
- self.mag = not self.masl # Sample is in Meters Above Ground
- self.species = species.strip()
- self.utc2loc = utc2loc # Hour difference between UTC and local
- self.samplingstrategy = samplingstrategy
- self.fromfile = fromfile # netcdf filename inside ObsPack distribution, to write back later
-
-################### End Class MoleFractionSample ###################
-
-
-if __name__ == "__main__":
- pass
-
-
-
diff --git a/da/stilt/sites_weights_stilt-updatedfoot-ivar.rc b/da/stilt/sites_weights_stilt-updatedfoot-ivar.rc
index 8b05d41693f793f86d98c502a8fc405d9f82f2e4..299327bbd78bbbc6bdd4ba0205fb860216cdf702 100755
--- a/da/stilt/sites_weights_stilt-updatedfoot-ivar.rc
+++ b/da/stilt/sites_weights_stilt-updatedfoot-ivar.rc
@@ -11,7 +11,7 @@ global.R.scaling: 1.e-6
! Number of different site categories:
-n.site.categories: 10
+n.site.categories: 11
! Each site.category.XXX contains four elements:
! 1. (string, len=30) label for the site category.
@@ -29,19 +29,21 @@ site.category.005: pfp_str; 3.0; 3.0
site.category.006: pfp_wbi; 5.0; 3.0; FALSE; TRUE
site.category.007: pfp_wgc; 3.0; 3.0; FALSE; TRUE
site.category.008: pfp_wkt; 3.0; 3.0; FALSE; TRUE
-site.category.009: do-not-use; 1000.00; 1000.00; FALSE; TRUE
-site.category.010: do-not-simulate; 1000.00; 1000.00; FALSE; TRUE
+site.category.009: pfp_air; 1.0; 1.0; FALSE; TRUE
+site.category.010: do-not-use; 1000.00; 1000.00; FALSE; TRUE
+site.category.011: do-not-simulate; 1000.00; 1000.00; FALSE; TRUE
! dataset
co2_amt_surface-pfp_1_allvalid : do-not-simulate noaa 2008-11-23 2012-08-14 no all valid measurements
co2_bao_surface-pfp_1_allvalid : do-not-simulate noaa 2007-08-16 2012-08-20 no all valid measurements
-co2_lef_surface-pfp_1_allvalid : do-not-simulate noaa 2002-08-24 2012-08-22 no all valid measurements
+co2_lef_surface-pfp_1_allvalid : pfp_lef noaa 2002-08-24 2012-08-22 no all valid measurements
co2_sct_surface-pfp_1_allvalid : do-not-simulate noaa 2008-08-14 2012-08-16 no all valid measurements
co2_str_surface-pfp_1_allvalid : do-not-simulate noaa 2007-10-02 2012-08-19 no all valid measurements
-co2_wbi_surface-pfp_1_allvalid : pfp_wbi noaa 2004-09-14 2012-08-12 no all valid measurements
+co2_wbi_surface-pfp_1_allvalid : do-not-simulate noaa 2004-09-14 2012-08-12 no all valid measurements
co2_wgc_surface-pfp_1_allvalid : do-not-simulate noaa 2007-09-21 2012-08-16 no all valid measurements
co2_wkt_surface-pfp_1_allvalid : do-not-simulate noaa 2006-07-07 2012-08-14 no all valid measurements
+co2_lef_aircraft-pfp_1_allvalid : pfp_air noaa 2002-08-24 2012-08-22 no all valid measurements
! Move some stations for better model representivity
@@ -88,7 +90,7 @@ site.move.020: co2_esp_surface-flask_2_representative
! first real is from hour of the day
! second real is to hour of the day
-!n.sites.hourly: 8
+n.sites.hourly: 9
site.hourly.001: co2_amt_surface-pfp_1_allvalid ; 10; 18
site.hourly.002: co2_bao_surface-pfp_1_allvalid ; 10; 18
@@ -98,8 +100,9 @@ site.hourly.005: co2_str_surface-pfp_1_allvalid
site.hourly.006: co2_wbi_surface-pfp_1_allvalid ; 10; 18
site.hourly.007: co2_wgc_surface-pfp_1_allvalid ; 10; 18
site.hourly.008: co2_wkt_surface-pfp_1_allvalid ; 10; 18
+site.hourly.009: co2_lef_aircraft-pfp_1_allvalid ; 10; 18
-n.sites.foot: 8
+n.sites.foot: 9
site.foot.001: co2_amt_surface-pfp_1_allvalid ;/Storage/ctdas-wrfstilt/new_stilt_footprints/flask/
site.foot.002: co2_bao_surface-pfp_1_allvalid ;/Storage/ctdas-wrfstilt/new_stilt_footprints/flask/
@@ -109,3 +112,4 @@ site.foot.005: co2_str_surface-pfp_1_allvalid ;
site.foot.006: co2_wbi_surface-pfp_1_allvalid ;/Storage/ctdas-wrfstilt/new_stilt_footprints/flask/
site.foot.007: co2_wgc_surface-pfp_1_allvalid ;/Storage/ctdas-wrfstilt/new_stilt_footprints/flask/
site.foot.008: co2_wkt_surface-pfp_1_allvalid ;/Storage/ctdas-wrfstilt/new_stilt_footprints/flask/
+site.foot.009: co2_lef_aircraft-pfp_1_allvalid ;/Storage/ctdas-wrfstilt/new_stilt_footprints/flask/
diff --git a/da/stilt/statevector.py b/da/stilt/statevector.py
index b42fdb144bbe124be65c3d17f0ab82acf0043b47..e76b7be88f42c8c944b42f16532bf783c95158b7 100755
--- a/da/stilt/statevector.py
+++ b/da/stilt/statevector.py
@@ -79,57 +79,11 @@ class CO2GriddedStateVector(StateVector):
# Boundary conditions covariance
- cov = np.array([[1.e-18]]) #np.ones((1,1),)
+ cov = np.array([[2]]) #np.ones((1,1),)
covariancematrixlist.append(cov)
- # Test
-
-
-
-
-
- covariancematrixlist = []
- for file in cov_files:
- if not os.path.exists(file):
- msg = "Cannot find the specified file %s" % file
- logging.error(msg)
- raise IOError, msg
- else:
- logging.debug("Using covariance file: %s" % file)
-
- f = io.ct_read(file, 'read')
-
- if 'pco2' in file:
- cov_ocn = f.get_variable('CORMAT')
- cov = cov_ocn*1.e-18
- else:
- cov = f.get_variable('covariance')
- #cov_sf = 10.0/np.sqrt(cov.diagonal().sum()) # this scaling factor makes the total variance close to the value of a single ecoregion
- cov_sf = 20.0 / np.sqrt(cov.diagonal().sum()) # this scaling factor makes the total variance close to the value of a single ecoregion
- cov = cov * cov_sf*1.e-18
-
- f.close()
- covariancematrixlist.append(cov)
-
- # Boundary conditions covariance
-
- cov = np.array([[2.]]) #np.ones((1,1),)
- covariancematrixlist.append(cov)
-
-
-
-
-
-
-
-
-
-
-
-
-
logging.debug("Succesfully closed files after retrieving prior covariance matrices")
@@ -239,6 +193,7 @@ class CO2GriddedStateVector(StateVector):
# Create mean values
new_mean = np.ones(self.nparams, float) # standard value for a new time step is 1.0
+ new_mean[-1] = np.zeros((1),float) # standard value for boundary conditions for new time step is 0.0
# If this is not the start of the filter, average previous two optimized steps into the mix