Skip to content
GitLab
Commit 36119604 authored by Woude, Auke van der's avatar Woude, Auke van der
Browse files

update for publication

parent 9ba799c0
Expand all Hide whitespace changes
Inline Side-by-side
da/ccffdas/category_info.py
View file @ 36119604
......@@ -49,21 +49,61 @@ categories = {
'CO2.uncertainty': 'n',
'CO': -7.52,
'CO.uncertainty': 'l'},
'cars': {'name': 'cars',
'cars highway': {'name': 'cars highway',
'model': 1,
'spatial': 'Road transport',
'temporal': 't_road',
'fraction_of_total': 1,
'temporal': 't_carhw',
'fraction_of_total': 0.47,
'emission_factors': 36200000,
'CO2' : 1,
'CO2.uncertainty': 'n',
'CO': -4.32,
'CO.uncertainty': 'l'},
'heavy duty': {'name': 'heavy duty highway',
'cars middle road': {'name': 'cars middle road',
'model': 1,
'spatial': 'Road transport',
'temporal': 't_road',
'fraction_of_total': 1.,
'temporal': 't_carmr',
'fraction_of_total': 0.28,
'emission_factors': 36200000,
'CO2' : 1,
'CO2.uncertainty': 'n',
'CO': -4.32,
'CO.uncertainty': 'l'},
'cars urban road': {'name': 'cars urban road',
'model': 1,
'spatial': 'Road transport',
'temporal': 't_carur',
'fraction_of_total': 0.25,
'emission_factors': 36200000,
'CO2' : 1,
'CO2.uncertainty': 'n',
'CO': -4.32,
'CO.uncertainty': 'l'},
'heavy duty highway': {'name': 'heavy duty highway',
'model': 1,
'spatial': 'Road transport',
'temporal': 't_hdvhw',
'fraction_of_total': 0.56,
'emission_factors': 36650000,
'CO2' : 1,
'CO2.uncertainty': 'n',
'CO': -6.11,
'CO.uncertainty': 'l'},
'heavy duty middle road': {'name': 'heavy duty middle road',
'model': 1,
'spatial': 'Road transport',
'temporal': 't_hdvmr',
'fraction_of_total': 0.24,
'emission_factors': 36650000,
'CO2' : 1,
'CO2.uncertainty': 'n',
'CO': -6.11,
'CO.uncertainty': 'l'},
'heavy duty urban': {'name': 'heavy duty urban',
'model': 1,
'spatial': 'Road transport',
'temporal': 't_hdvur',
'fraction_of_total': 0.2,
'emission_factors': 36650000,
'CO2' : 1,
'CO2.uncertainty': 'n',
......
da/ccffdas/emissionmodel.py
View file @ 36119604
This diff is collapsed.
da/ccffdas/energy_use_country.py
View file @ 36119604
energy_use_per_country = {
'AUS': {
'Public power': 161.01,
'Public power ratio gas': 0.78,
'Industry': 158.02,
'Other stationary combustion': 120.04,
'Road transport': 150.8,
'Shipping': 0.14
},
'BEL': {
'Public power': 254.82,
'Public power ratio gas': 0.72,
'Industry': 203.68,
'Other stationary combustion': 383.16,
'Road transport': 172.455,
'Shipping': 5.63
},
'CZE': {
'Public power': 582.9,
'Public power ratio gas': 0.10,
'Industry': 140.13,
'Other stationary combustion': 178.46,
'Road transport': 121.935,
'Shipping': 0.17},
'FRA': {
'Public power': 585.52,
'Public power ratio gas': 0.56,
'Industry': 731.9,
'Other stationary combustion': 1262.1,
'Road transport': 845.07,
'Shipping': 20.0
},
'DEU': {
'Public power': 3515.68,
'Public power ratio gas': 0.19,
'Industry': 1516.34,
'Other stationary combustion': 2069.16,
'Road transport': 1066.98,
'Shipping': 26.11
},
'LUX': {
'Public power': 3.75,
'Public power ratio gas': 0.99,
'Industry': 15.23,
'Other stationary combustion': 25.22,
'Road transport': 37.045,
'Shipping': 0.01
},
'NED': {
'Public power': 868.96,
'Public power ratio gas': 0.47,
'Industry': 427.43,
'Other stationary combustion': 599.52,
'Road transport': 199.74,
'Shipping': 13.93
},
'POL': {
'Public power': 1702.42,
'Public power ratio gas': 0.05,
'Industry': 338.92,
'Other stationary combustion': 618.77,
'Road transport': 362.94,
'Shipping': 0.29
},
'CHE': {
'Public power': 43.41,
'Public power ratio gas': 0.53,
'Industry': 66.5,
'Other stationary combustion': 194.77,
'Road transport': 100.28,
'Shipping': 1.54
},
'GBR': {
'Public power': 1716.77,
'Public power ratio gas': 0.76,
'Industry': 620.58,
'Other stationary combustion': 1485.23,
'Road transport': 787.845,
'Shipping': 72.39}}
'AUS': {'Public power': 161.01,
'Industry': 158.02,
'Other stationary combustion': 120.04,
'Road transport': 150.8,
'Shipping': 0.14},
'BEL': {'Public power': 254.82,
'Industry': 203.68,
'Other stationary combustion': 383.16,
'Road transport': 172.455,
'Shipping': 5.63},
'CZE': {'Public power': 582.9,
'Industry': 140.13,
'Other stationary combustion': 178.46,
'Road transport': 121.935,
'Shipping': 0.17},
'FRA': {'Public power': 585.52,
'Industry': 731.9,
'Other stationary combustion': 1262.1,
'Road transport': 845.07,
'Shipping': 20.0},
'DEU': {'Public power': 3515.68,
'Industry': 1516.34,
'Other stationary combustion': 2069.16,
'Road transport': 1066.98,
'Shipping': 26.11},
'LUX': {'Public power': 3.75,
'Industry': 15.23,
'Other stationary combustion': 25.22,
'Road transport': 37.045,
'Shipping': 0.01},
'NED': {'Public power': 868.96,
'Industry': 427.43,
'Other stationary combustion': 599.52,
'Road transport': 199.74,
'Shipping': 13.93},
'POL': {'Public power': 1702.42,
'Industry': 338.92,
'Other stationary combustion': 618.77,
'Road transport': 362.94,
'Shipping': 0.29},
'CHE': {'Public power': 43.41,
'Industry': 66.5,
'Other stationary combustion': 194.77,
'Road transport': 100.28,
'Shipping': 1.54},
'GBR': {'Public power': 1716.77,
'Industry': 620.58,
'Other stationary combustion': 1485.23,
'Road transport': 787.845,
'Shipping': 72.39}}
da/ccffdas/observationoperator.py
View file @ 36119604
......@@ -94,19 +94,18 @@ recalc_factors = {'CO2': 1, 'CO': 1000}
spname = ['CO2', 'CO']
def run_STILT(dacycle, footprint, datepoint, species, i_member):
def run_STILT(dacycle, footprint, datepoint, species, path, i_member):
"""This function reads in STILT footprints and returns hourly concentration data
Input:
dacycle : dict: Dictionary containing information on the current time
footprint : np.array: the footprint of the station
datepoint : datetime: the datepoint for which the concentrations should be calculated
footprint: np.array: the footprint of the station
datepoint: datepoint: datetime: the datepoint for which the concentrations should be calculated
site : str: name of the location
species : str: the species
path : path: path to the files containing the emissions
i_member : int: index of the ensemblemember
Returns:
float: the concentration increase at the respective location due to the emissions from the respective species"""
# get the emission data:
spatial_emissions = get_spatial_emissions(dacycle, i_member, species, datepoint)
spatial_emissions = get_spatial_emissions(dacycle, i_member, species, path, datepoint)
# multiply footprint with emissions for each hour in the back trajectory. Indices: Time, Category, lat, lon
foot_flux = (footprint[None, :, :, :] * spatial_emissions[:, :, :, :]).sum()
......@@ -114,13 +113,11 @@ def run_STILT(dacycle, footprint, datepoint, species, i_member):
return concentration_increase
@cached
def get_spatial_emissions(dacycle, i_member, species, datepoint):
def get_spatial_emissions(dacycle, i_member, species, path, datepoint):
""" Function that gets the spatial emissions
Input:
dacycle: dict: Dictionary containing information on the current time
i_member: int: the index of the member for which the simulation is run
species: str: the species for which the simulation is run
datepoint: datetime: the date and time of the concentration measurement
Returns:
np.ndarray (3d): the spatial emissions per category, lat and lon"""
#read spatially distributed emissions calculated with the emission model
......@@ -129,7 +126,6 @@ def get_spatial_emissions(dacycle, i_member, species, datepoint):
start_time = (dacycle['time.sample.end'] - backtimes)
indices = get_time_indices(datepoint, start_time)
path = dacycle.dasystem['inputdir']
emisfile = path +'prior_spatial_{0:03d}.nc'.format(i_member) #format: species[SNAP, time, lon, lat]
f = io.ct_read(emisfile,method='read')
emis=f.get_variable(species)
......@@ -207,11 +203,10 @@ class STILTObservationOperator(ObservationOperator):
self.paramdict = rc.read(dacycle.dasystem['paramdict'])
self.pftfile = dacycle.dasystem['file.pft']
with rasterio.Env(logging.getLogger().setLevel(logging.INFO)):
with rasterio.Env(logging.getLogger().setLevel(logging.DEBUG)):
with rasterio.open(self.pftfile) as dataset:
self.pft_shape_orig = dataset.read().shape[-2:]
logging.getLogger().setLevel(logging.DEBUG)
self.lon_start = float(dacycle.dasystem['domain.lon.start'])
self.lon_end = float(dacycle.dasystem['domain.lon.end'])
self.lat_start = float(dacycle.dasystem['domain.lat.start'])
......@@ -354,12 +349,11 @@ class STILTObservationOperator(ObservationOperator):
for i in range(2)])
logging.debug('Shape for paint-by-number: {}'.format(new_shape))
with rasterio.Env(logging.getLogger().setLevel(logging.INFO)):
with rasterio.Env(logging.getLogger().setLevel(logging.DEBUG)):
with rasterio.open(self.pftfile) as dataset:
pftdata = dataset.read(out_shape=new_shape,
resampling=Resampling.mode)
pftdata = pftdata.reshape(pftdata.shape[1:])
logging.getLogger().setLevel(logging.DEBUG)
scaling_factors = np.array(np.array(new_shape) / np.array(self.lu_pref[0].shape), dtype=int)
......@@ -390,8 +384,7 @@ class STILTObservationOperator(ObservationOperator):
for mem, values in enumerate(all_param_values):
param_values = all_param_values[mem]
index_in_state = find_in_state(param_values, 'BIO', str(lutype), return_index=True)
if index_in_state:
param_value = param_values[index_in_state]
if index_in_state: param_value = param_values[index_in_state]
else: param_value = 1
nee_lut_mem = (nee_lut.sum(axis=0) * mask) * param_value
nee_lut_mem_scaled = average2d(nee_lut_mem, domain_shape)
......@@ -407,9 +400,11 @@ class STILTObservationOperator(ObservationOperator):
time_profile_nuc = np.array(worksheet.col_values(19))
self.nuc_time_profile = time_profile_nuc
file = self.model_settings['biosphere_fluxdir']
self.dis_eq_flux = self.get_nc_variable(file, 'TER_dC14', np.float32)[self.indices, :, :]
file = self.model_settings['nuclear_fluxdir']
with nc.Dataset(file) as ds:
self.nuc_flux = ds['E_14CO2_nuc'][:].astype(np.float32)
self.nuc_flux = self.get_nc_variable(file, 'E_14CO2_nuc', np.float32)
def get_foot(self, site, datepoint):
"""Function that gets the footprint for the current time and site.
......@@ -514,22 +509,27 @@ class STILTObservationOperator(ObservationOperator):
# Get the fluxes and multiply them by the footprint and time profile
# Note that the unit is already in Delta notation
nuc_flux = self.nuc_flux # time invariant
nuc_flux = (nuclear_time_profile[:, None, None] * nuc_flux[None, :, :] * self.foot).sum()
nuc_flux = self.nuc_flux
nuc_flux = (nuclear_time_profile[:, None, None] * nuc_flux[None,:] * self.foot).sum()
nuc_flux *= (1./MASS_14CO2) * KMOL2UMOL * PERYR2PERS
# now get the disequilibrium flux
delta_14CO2_ter = (self.dis_eq_flux[indices] * self.foot).sum()
R_14CO2_ter = (delta_14CO2_ter/1000. + 1) * Rstd # [nt,nlat,nlon] 1.639e-12 [1.557e-12 - 1.787e-12]
# Now, calculate to ppm
co2_14_ff = ff_flux * R_14CO2_ff * MASS_C / MASS_C_14
co2_14_bg = background* R_14CO2_bg * MASS_C / MASS_C_14
co2_14_bio= nee * R_14CO2_bg * MASS_C / MASS_C_14
# Calculate the total c14 in the sample (ppm)
co2_14_total = (co2_14_ff + # The C14 due to FF emissions (should be 0)
co2_14_bio + # C14 from the biospehre, ignore disequilibruim flux
nuc_flux + # The C14 from the nuclera power plants
co2_14_bg) # The C14 that was already in the air
(alpha_14CO2_gpp * R_14CO2_bg * gpp) + # The C14 taken up due to GPP
# (note that GPP is negative)
(R_14CO2_ter * ter) + # The C14 from the disequilibrium flux
nuc_flux + # The C14 from the nuclera power plants
co2_14_bg) # The C14 that was already in the air
# Conversion 14CO2 in ppm to DELTA ----
co2_total = ff_flux + nee + nuc_flux + background
co2_total = ff_flux + gpp + ter + nuc_flux + background
Rm = (co2_14_total * MASS_14CO2) / (co2_total * MASS_CO2) # ppm to mass ratio (kg 14CO2 / kgCO2)
A14 = Rm * lambda_14CO2 * 1e3*Nav / MASS_14CO2 # kg14CO2/kgCO2 * 1/s * molecules/kmole / (kg14CO2/kmole14CO2) = molecules 14CO2 / kgCO2 / s
As = A14 * MASS_CO2 / MASS_C # Bq/kgC-CO2
......@@ -537,6 +537,19 @@ class STILTObservationOperator(ObservationOperator):
delta_14CO2 = (Asn * np.exp(lambda_14CO2 * dt_obs) / Aabs -1) * 1000.# per mille
return delta_14CO2
def get_nc_variable(self, file, fluxname, dtype):
"""Helper function that gets the values from an nc file
Input:
file: str: filename of the nc file of which the value should be taken
fluxname: str: name of the variable in the file
Returns:
np.ndarray: the values in the nc file"""
data = nc.Dataset(file)
fluxmap = data[fluxname][:]
fluxmap = fluxmap.astype(dtype)
data.close()
return fluxmap
def run_forecast_model(self,dacycle,do_pseudo,adv):
logging.info('Preparing run...')
self.prepare_run(do_pseudo,adv)
......@@ -628,7 +641,7 @@ class STILTObservationOperator(ObservationOperator):
if not 'C14' in species.upper():
pool = Pool(self.forecast_nmembers)
# Create the function that calculates the conentration
func = partial(run_STILT, self.dacycle, self.foot, datepoint, species)
func = partial(run_STILT, self.dacycle, self.foot, datepoint, species, self.inputdir)
# We need to run over all members
memberlist = list(range(0, self.forecast_nmembers))
......@@ -645,7 +658,7 @@ class STILTObservationOperator(ObservationOperator):
self.c14_dict[site][datepoint] = {}
nee = self.get_biosphere_concentration(self.foot, self.nee_mem, datepoint)
self.c14_dict[site][datepoint]['nee'] = nee[:]
self.c14_dict[site][datepoint]['bio'] = nee[:]
self.c14_dict[site][datepoint]['ff'] = ff_increase[:]
# If not the species is CO2, the NEE is 0.
......
da/ccffdas/stilt-ops_urbanall.rc
View file @ 36119604
!!! Info for the CarbonTracker data assimilation system
basepath : /projects/0/ctdas/awoude/develop/
name :
strategy : CO2C14
strategy : CO2
datadir : /projects/0/ctdas/RINGO/inversions/Data
inputdir : ${basepath}/input/
outputdir : ${basepath}/output/
......@@ -14,20 +14,19 @@ obs.spec.nr : 1
obs.dir : obsfiles${strategy}
do.co : 0
do.c14integrated: 0
do.c14targeted: 1
do.c14targeted: 0
obs.spec.name : CO2
! number of emission categories defined in the emission model
obs.cat.nr : 10
obs.cat.nr : 14
! For Rdam obs
obs.sites.rc : ${datadir}/sites_weights2.rc
! number of parameters
! In the covmatrix and statevector, the ff parameters are first, then the bio parameters!
nffparameters : 24
nbioparameters : 8
nbioparameters : 6
nparameters : ${nffparameters} + ${nbioparameters}
file.pft : ${datadir}/SiBPFTs.nc
file.timeprofs : ${datadir}/CAMS_TEMPO_Tprof_subset.nc
file.pft : /projects/0/ctdas/RINGO/inversions/Data/SiBPFTs.nc
! Settings for the domain:
domain.lon.start : -5.95
......@@ -37,7 +36,7 @@ domain.lat.end : 54.975
domain.lon.num : 260
domain.lat.num : 240
paramdict : ${datadir}/paramdict_sib.rc
paramdict : ${datadir}/paramdict.rc
! 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
......
develop.rc.easy
View file @ 36119604
......@@ -28,7 +28,7 @@
!
! The time for which to start and end the data assimilation experiment in format YYYY-MM-DD HH:MM:SS
time.start : 2016-01-06 00:00:00
time.start : 2016-01-04 00:00:00
time.finish : 2016-01-20 00:00:00
time.fxstart : 2016-01-01 00:00:00
......
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