Commit e707cf68 authored by Peters, Wouter's avatar Peters, Wouter
Browse files

modified the baseclasses to have simpler init functions, and new observation object

parent 9c523b10
......@@ -13,48 +13,56 @@ import sys
import logging
import datetime
identifier = 'CarbonTracker CO2'
identifier = 'Observation baseclass'
version = '0.0'
################### Begin Class MixingRatioSample ###################
################### Begin Class Observation ###################
class MixingRatioSample(object):
"""
Holds the data that defines a Mixing Ratio 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.
class Observation(object):
""" an object that holds data + methods and attributes needed to manipulate observations for a DA cycle """
"""
def __init__(self):
self.Identifier = identifier
self.Version = version
self.Data = ObservationList([]) # Initialize with an empty list of obs
def __init__(self,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',sdev=0.0):
self.code = code.strip() # Site code
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.mdm = mdm # Model data mismatch
self.flag = flag # Flag
self.height = height # Sample height
self.lat = lat # Sample lat
self.lon = lon # Sample lon
self.id = evn # Event number
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
msg = '%s object initialized'%self.Identifier ; logging.debug(msg)
msg = '%s version: %s'%(self.Identifier,self.Version) ; logging.info(msg)
def __str__(self):
day=self.xdate.strftime('%Y-%m-%d %H:%M:%S')
return ' '.join(map(str,[self.code,day,self.obs,self.flag,self.id]))
""" Prints a list of Observation objects"""
return "This is a %s object, version %s"%(self.Identifier,self.Version)
def Initialize(self):
""" Perform all steps needed to start working with observational data, this can include moving data, concatenating files,
selecting datasets, etc.
"""
def Validate(self):
""" Make sure that data needed for the ObservationOperator (such as observation input lists, or parameter files)
are present.
"""
################### End Class MixingRatioSample ###################
def AddObs(self):
""" Add the observation data to the Observation object. This is in a form of a list that goes under the name Data. The
list has as only requirement that it can return the observed+simulated values through a method "getvalues"
"""
def AddSimulations(self):
""" Add the observation data to the Observation object. This is in a form of a list that goes under the name Data. The
list has as only requirement that it can return the observed+simulated values through a method "getvalues"
"""
################### Begin Class MixingRatioList ###################
class MixingRatioList(list):
""" This is a special type of list that holds MixingRatioSample objects. It has methods to extract data from such a list easily """
################### End Class Observations ###################
################### Begin Class ObservationList ###################
class ObservationList(list):
""" This is a special type of list that holds observed sample objects. It has methods to extract data from such a list easily """
from numpy import array, ndarray
def getvalues(self,name,constructor=array):
......@@ -64,189 +72,6 @@ class MixingRatioList(list):
return result.squeeze()
else:
return result
def unflagged(self):
return MixingRatioSample([o for o in self if o.flag == 0])
def flagged(self):
return MixingRatioSample([o for o in self if o.flag != 0])
def selectsite(self,site='mlo'):
l=[o for o in self if o.code == site]
return MixingRatioSample(l)
def selecthours(self,hours=range(1,24)):
l=[o for o in self if o.xdate.hour in hours]
return MixingRatioSample(l)
################### End Class MixingRatioList ###################
################### Begin Class CtObservations ###################
class CtObservations(object):
""" an object that holds data + methods and attributes needed to manipulate mixing ratio values """
def __init__(self,DaSystem,sampledate):
sfname = DaSystem.da_settings['obs.input.fname']+'.%s'%(sampledate.strftime('%Y%m%d'),)+'.nc'
#msg = 'This filename is hardcoded but needs replacement' ; logging.warning(msg)
filename = os.path.join(DaSystem.da_settings['obs.input.dir'],sfname)
if not os.path.exists(filename):
msg = 'Could not find the required observation input file (%s) ' % filename ; logging.error(msg)
raise IOError,msg
else:
self.ObsFilename = filename
self.MrData = MixingRatioList([])
def __str__(self):
""" Prints a list of MixingRatioSample objects"""
return self.ObsFilename
def AddObs(self):
""" Returns a MixingRatioList holding individual MixingRatioSample objects for all obs in a file
Note that this routine was not yet ported to Nio because it cannot read string arrays !!!
"""
import pyhdf.SD as HDF
import datetime as dtm
from string import strip, join
ncf = HDF.SD(self.ObsFilename)
idates = ncf.select('date_components').get()
ids = ncf.select('id').get()
sites = ncf.select('site').get()
sites = [s.tostring().lower() for s in sites]
sites = map(strip,sites)
lats = ncf.select('lat').get()
lons = ncf.select('lon').get()
alts = ncf.select('alt').get()
obs = ncf.select('obs').get()
species = ncf.select('species').get()
date = ncf.select('date').get()
window = ncf.select('sampling_strategy').get()
flags = ncf.select('NOAA_QC_flags').get()
flags = [s.tostring().lower() for s in flags]
flags = map(strip,flags)
dummy = ncf.end()
msg = "Successfully read data from obs file (%s)"%self.ObsFilename ; logging.debug(msg)
dates = [dtm.datetime(*d) for d in idates]
if self.MrData == None:
self.MrData = MixingRatioList([])
for n in range(len(dates)):
self.MrData.append( MixingRatioSample(dates[n],sites[n],obs[n],0.0,0.0,0.0,0.0,flags[n],alts[n],lats[n],lons[n],ids[n],0.0) )
msg = "Added %d observations to the MrData list" % len(dates) ; logging.debug(msg)
def AddSimulations(self, filename,silent=True):
""" Adds model simulated values to the mixing ratio objects """
import pyhdf.SD as HDF
import Nio
if not os.path.exists(filename):
msg = "Could not find required model sample file (%s)"%filename ; logging.debug(msg)
raise IOError,msg
ncf = Nio.open_file(filename,format='hdf')
ids = ncf.variables['id'].get_value()
simulated = ncf.variables['sampled_mole_frac'].get_value()*1e6
dummy = ncf.close()
msg = "Successfully read data from model sample file (%s)"%filename ; logging.info(msg)
obs_ids = self.MrData.getvalues('id')
obs_ids = obs_ids.tolist()
ids = map(int,ids)
missing_samples = []
for id,val in zip(ids,simulated):
if id in obs_ids:
index = obs_ids.index(id)
dummy = self.MrData[index].simulated = val
else:
missing_samples.append(id)
if not silent and missing_samples != []:
msg = 'Model samples were found that did not match any ID in the observation list. Skipping them...' ; logging.warning(msg)
msg = '%s'%missing_samples ; logging.warning(msg)
msg = "Added %d simulated values to the MrData list" % (len(ids)-len(missing_samples) ) ; logging.debug(msg)
################### End Class CtObservations ###################
################# Stand alone methods that manipulate the objects above ##############
def PrepareObs(CycleInfo,lag):
""" PrepareObs """
import shutil
from da.tools.general import AdvanceTime
# First, we make a CarbonTracker Observation object which holds all the information related to the mixing ratios going in
# and coming out of the model. This object also holds the observations themselves, and the simulated value + statistics
sampledate = CycleInfo.da_settings['time.start']
sampledate = AdvanceTime(sampledate,lag*int(CycleInfo.da_settings['cyclelength']))
samples = CtObservations(CycleInfo.DaSystem, sampledate)
dummy = samples.AddObs()
#
# Next should follow a section where the observations that are read are further processes, and written to a NetCDF file for
# the transport model to use for x,y,z,t smapling.
#
# For carbontracker the observations already come in netcdf format thanks to Andy and Ken. We will therefore simply copy
# the observations.nc file to the pre-cooked obs files
#
filename= samples.ObsFilename
saveas = os.path.join(CycleInfo.da_settings['dir.input'],'observations.nc')
dummy = shutil.copy2(filename,saveas)
msg = 'Successfully copied file (%s) to (%s) ' % (filename,saveas) ; logging.info(msg)
return samples
if __name__ == "__main__":
sys.path.append('../../')
import os
import sys
from da.tools.general import StartLogger
from da.tools.initexit import CycleControl
import numpy as np
import da.tools.rc as rc
opts = ['-v']
args = {'rc':'../../da.rc','logfile':'da_initexit.log','jobrcfilename':'test.rc'}
StartLogger()
DaCycle = CycleControl(opts,args)
DaCycle.Initialize()
print DaCycle
samples = CtObservations(DaCycle.DaSystem,datetime.datetime(2005,3,5))
dummy = samples.AddObs()
dummy = samples.AddSimulations('/Users/peters/tmp/test_da/output/samples.000.nc')
print samples.MrData.getvalues('code')
print samples.MrData.getvalues('obs')
mlo=samples.MrData.selectsite('mlo_01d0')
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