Commit f363a0e5 authored by Woude, Auke van der's avatar Woude, Auke van der
Browse files

QoL fixes

parent 107cec70
......@@ -89,7 +89,7 @@ alpha_14CO2_gpp = 1 + epsilon_14CO2_gpp/1000.
model_settings = rc.read('/projects/0/ctdas/RINGO/inversions/ffdas/exec/da/rc/ffdas/stilt.rc')
recalc_factors = [1, 1000]
spname = ['CO2', 'CO']
def run_STILT(footprint, datepoint, site, i_species, i_member):
def run_STILT(footprint, datepoint, site, i_species, path, i_member):
"""This function reads in STILT footprints and returns hourly concentration data
Input:
footprint: np.array: the footprint of the station
......@@ -100,8 +100,8 @@ def run_STILT(footprint, datepoint, site, i_species, i_member):
Returns:
float: the concentration increase at the respective location due to the emissions from the respective species"""
# get the date:
temp_profiles = get_temporal_profiles(datepoint, site, i_member)
spatial_emissions = get_spatial_emissions(i_member, i_species)
temp_profiles = get_temporal_profiles(datepoint, site, path, i_member)
spatial_emissions = get_spatial_emissions(i_member, i_species, path)
# multiply footprint with emissions for each hour in the back trajectory. Indices: Time, Category, lat, lon
foot_flux = (footprint[:, None, :, :] * spatial_emissions[None, :, :, :] * temp_profiles[:, :, None, None]).sum()
......@@ -110,7 +110,7 @@ def run_STILT(footprint, datepoint, site, i_species, i_member):
return concentration_increase
@cached
def get_temporal_profiles(datepoint, station_name, i_member):
def get_temporal_profiles(datepoint, station_name, path, i_member):
""" Function that reads in the temporal profiles for the current timestep
Input:
datepoint: datepoint: datetime for which the temporal profile should be found
......@@ -120,7 +120,7 @@ def get_temporal_profiles(datepoint, station_name, i_member):
np.array (2): the time profiles for all categories. Indices: time, category"""
#read temporal profiles for the times within the footprint
time_indices = get_time_indices(datepoint)
temporal_prior = os.path.join(model_settings['inputdir'],'prior_temporal_{0}_{1:03}.nc'.format(station_name, i_member))
temporal_prior = path + 'prior_temporal_{0}_{1:03}.nc'.format(station_name, i_member)
temporal_prior = io.ct_read(temporal_prior, 'read')
temp_profiles = []
for category, values in categories.items():
......@@ -132,7 +132,7 @@ def get_temporal_profiles(datepoint, station_name, i_member):
return temp_profiles.T # Transpose to be similar to spatial data in dimensions
@cached
def get_spatial_emissions(i_member, i_species):
def get_spatial_emissions(i_member, i_species, path):
""" Function that gets the spatial emissions
Input:
i_member: int: the index of the member for which the simulation is run
......@@ -140,7 +140,7 @@ def get_spatial_emissions(i_member, i_species):
Returns:
np.ndarray (3d): the spatial emissions per category, lat and lon"""
#read spatially distributed emissions calculated with the emission model
emisfile = os.path.join(model_settings['inputdir'],'prior_spatial_{0:03d}.nc'.format(i_member)) #format: species[SNAP,lon,lat]
emisfile = path +'prior_spatial_{0:03d}.nc'.format(i_member) #format: species[SNAP,lon,lat]
f = io.ct_read(emisfile,method='read')
emis=f.get_variable(spname[i_species])
f.close()
......@@ -174,28 +174,25 @@ def get_time_indices(datepoint, startdate=None):
def get_biosphere_concentration(foot, gpp_mem, ter_mem):
"""Function that calculates the atmospheric increase due to the exchange of fluxes over the footprint
Input:
i_member: int: member number for which the biosphere concentration should be calculated
site: str: location for which the concentration should be calculated
datepoint: datetime.datetime: the datepoint for which the concentration should be calculated
total: bool, optional: Whether to returned summed values or gpp and ter seperately as tuple
foot: np.array(times, lat, lon):
the footprint of the station
gpp_mem: np.array(nmmebers, time, lat, lon)
map of the gpp for each member
ter_mem: np.array(nmembers, time, lat, lon)
map of the ter for each member
Returns:
if total == True:
float: The total biosphere flux in umol/s
else:
tuple of 2 floats: GPP and TER in umol/s """
tuple of 2 floats: GPP and TER in umol/s """
# First, get the footprint
# First, recalculate to good units
gpp = gpp_mem[:, :, :, :] * (1./MASS_C) * MOL2UMOL * PERHR2PERS
ter = ter_mem[:, :, :, :] * (1./MASS_C) * MOL2UMOL * PERHR2PERS
# Multiply with the footprint
gpp_increase = - (gpp * foot[None, :, :, :]).sum()
ter_increase = (ter * foot[None, :, :, :]).sum()
# logging.debug('GGP flux = {0:.3}; TER flux = {1:.3}'.format(gpp_increase, ter_increase))
return gpp_increase, ter_increase
DO_RINGO = 0
class STILTObservationOperator(ObservationOperator):
def __init__(self, filename, dacycle=None): #only the filename used to specify the location of the stavector file for wrf-stilt runs
""" The instance of an STILTObservationOperator is application dependent """
......@@ -207,8 +204,7 @@ class STILTObservationOperator(ObservationOperator):
logging.info('Observation Operator initialized: %s (%s)' % (self.ID, self.version))
self.load_rc(filename) # load the specified rc-file
self.filename = filename
self.model_settings = rc.read(filename)
if dacycle != None:
self.dacycle = dacycle
else:
......@@ -216,9 +212,6 @@ class STILTObservationOperator(ObservationOperator):
self.startdate=None
if DO_RINGO: logging.warning('DO_RINGO is used, so the time indices +=1')
def setup(self, dacycle):
""" Execute all steps needed to prepare the ObsOperator for use inside CTDAS,
only done at the very first cycle normally """
......@@ -234,6 +227,7 @@ class STILTObservationOperator(ObservationOperator):
self.bgfile = dacycle.dasystem['obs.background']
self.btime = int(dacycle.dasystem['run.backtime'])
self.categories = categories
self.inputdir = dacycle.dasystem['inputdir']
self.paramdict = rc.read('/projects/0/ctdas/RINGO/inversions/Data/paramdict.rc')
biosphere_fluxes = nc.Dataset(dacycle.dasystem['biosphere_fluxdir'])
......@@ -244,27 +238,11 @@ class STILTObservationOperator(ObservationOperator):
with nc.Dataset(self.dacycle.dasystem['country.mask']) as countries:
self.masks = countries['country_mask'][:]
self.country_names = countries['country_names'][:]
self.nffparams = int( self.dacycle.dasystem['nffparameters'])
self.nbioparams =int( self.dacycle.dasystem['nbioparameters'])
self.noise = {'CO2': 2.2, 'CO': 8}
self.nffparams = int(self.dacycle.dasystem['nffparameters'])
self.nbioparams = int(self.dacycle.dasystem['nbioparameters'])
self.noise = {'CO2': 2.2, 'CO': 8, 'C14': 2, 'C14_PSEUDO': 0, 'C14_INTEGRATED': 2, 'C14targeted': 2}
logging.info('Noise is hardcoded to be: {}'.format(self.noise))
def load_rc(self, name):
"""This method loads a STILT rc-file with settings for this simulation
Args:
name: str: name of the .rc file that will be loaded
Should include: - num_backtimes
- files_startdate
- footdir
- datadir
Returns
None
"""
self.rcfile = rc.RcFile(name)
self.model_settings = rc.read(name)
logging.debug('rc-file loaded successfully')
def get_initial_data(self, samples):
"""Function that loads the initial data to the observation operator"""
self.allsamples = samples
......@@ -297,8 +275,7 @@ class STILTObservationOperator(ObservationOperator):
# Set the number of ensemble members
self.forecast_nmembers = int(self.dacycle['da.optimizer.nmembers'])
# Set the inputdir
self.param_inputdir = self.dacycle['dir.input']
self.param_inputdir = self.dacycle.dasystem['inputdir']
#list all observation locations and species
self.lat = []
self.lon = []
......@@ -712,13 +689,16 @@ class STILTObservationOperator(ObservationOperator):
# Get the background concentration
background = self.get_background_orig(species.upper())
# First, add noise (this represents errors in the transport model
noise = np.random.normal(0, self.noise[sample.species.upper()])
# Some different cases for different species.
# First case: Species is not c14:
# Then calculate the increase in conc due to ff emissions
if not 'C14' in species.upper():
pool = Pool(self.forecast_nmembers)
# Create the function that calculates the conentration
func = partial(run_STILT, self.foot, datepoint, site, i_species)
func = partial(run_STILT, self.foot, datepoint, site, i_species, self.inputdir)
# We need to run over all members
memberlist = list(range(0, self.forecast_nmembers))
......@@ -801,12 +781,10 @@ class STILTObservationOperator(ObservationOperator):
ter = self.c14_dict[site][datepoint]['bio'][1][mem]
ff_flux = self.c14_dict[site][datepoint]['ff'][mem]
c14concentrations[mem] = self.get_c14_concentration(datepoint, gpp, ter, ff_flux, self.get_background_orig('CO2'))
return c14concentrations
return c14concentrations + noise
# If the species is NOT c14, return the atmospheric concentration.
# First, add noise (this represents errors in the transport model
noise = np.random.normal(0, self.noise[sample.species.upper()])
new_conc = nee + ff_increase + background + noise
self.new_conc = new_conc
return new_conc
......
!!! Info for the CarbonTracker data assimilation system
datadir : /projects/0/ctdas/RINGO/inversions/Data
inputdir : /projects/0/ctdas/RINGO/inversions/ffdas/input/
basepath : /projects/0/ctdas/RINGO/inversions/
name :
strategy : CO2C14
datadir : ${basepath}/Data
inputdir : ${basepath}/${name}${strategy}/input/
outputdir : ${basepath}/${name}${strategy}/output/
restartdir : ${basepath}/${name}${strategy}/restart/
! list of all observation sites
obs.input.id : obsfiles.csv
! number of observation sites included; number of species included and to be used in inversion
obs.input.nr : 100
obs.spec.nr : 1
obs.dir : obsfiles1b
obs.spec.nr : 2
obs.dir : obsfiles${strategy}
do.co : 0
do.c14integrated: 0
do.c14targeted: 0
......@@ -18,9 +22,9 @@ obs.cat.nr : 14
obs.sites.rc : ${datadir}/sites_weights.rc
! number of parameters
! In the covmatrix and statevector, the ff parameters are first, then the bio parameters!
nffparameters : 32
nffparameters : 16
nbioparameters : 4
nparameters : 36
nparameters : 20
! set fixed seed for random number generator, or use 0 if you want to use any random seed
random.seed : 4385
!file with prior estimate of scaling factors (statevector) and covariances
......
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