#!/usr/bin/env python
# obs.py
"""
Author : peters
Revision History:
File created on 28 Jul 2010.
"""
import os
import sys
import logging
import datetime
identifier = 'CarbonTracker CO2'
################### Begin Class MixingRatioSample ###################
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.
"""
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
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]))
################### End Class MixingRatioSample ###################
################### 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 """
from numpy import array, ndarray
def getvalues(self,name,constructor=array):
from numpy import ndarray
result = constructor([getattr(o,name) for o in self])
if isinstance(result,ndarray):
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')