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

Paint by number

parent 54d116f9
......@@ -68,21 +68,7 @@ class EmisModel(object):
self.categories = categories
self.paramdict = rc.read(dacycle.dasystem['paramdict'])
# --- Settings
year = self.startdate.year
ndays = (dtm.datetime(year, 12, 31) - dtm.datetime(year - 1, 12, 31)).days
self.ndays = ndays
times = np.array([dtm.datetime(year, 1, 1, 0, 0, 0) + dtm.timedelta(hours = i) for i in range(ndays*24)])
self.times = times
datapath = '/projects/0/ctdas/RINGO/EmissionInventories/DynEmMod_TimeProfiles'
infilename = '{}/RINGO_ECMWF_DailyMeteo{}.nc'.format(datapath, year)
with nc.Dataset(infilename) as infile:
self.T2myr = infile.variables['T2m' ][:ndays] - 273.15 # T2m K --> oC
self.U10 = infile.variables['U10m' ][:ndays] # u m/s
self.Rinyr = infile.variables['Rsin' ][:ndays] / (24. * 1.0e4) # Radiation (?) J/m2/day --> J/cm2/hr
self.T2myr_av = infile.variables['T2m_avg' ][:ndays] - 273.15 #
self.U10_av = infile.variables['U10m_avg'][:ndays]
self.Rinyr_av = infile.variables['Rsin_avg'][:ndays]
logging.debug('Emismodel has been set-up')
def find_in_state(self, params, station, cat, name=None, return_index=False):
......@@ -91,6 +77,7 @@ class EmisModel(object):
key = station + '.' + cat
else:
key = station + '.' + cat + '.' + name
#print('Looking for {}'.format(key))
if key in self.paramdict:
i_in_state = int(self.paramdict[key])
......@@ -99,7 +86,7 @@ class EmisModel(object):
value = params[i_in_state]
return value
elif return_index: return False
return 1
else: return 1
def get_emis(self, dacycle, samples, indices, do_pseudo):
"""set up emission information for pseudo-obs (do_pseudo=1) and ensemble runs (do_pseudo=0)"""
......@@ -119,7 +106,7 @@ class EmisModel(object):
logging.debug('Succesfully wrote prior emission files')
def get_yearly_emissions(self, params):
yremis = np.zeros((len(self.categories), len(self.countries), len(self.species)))
yremis = np.zeros((len(self.categories), len(self.countries), len(self.species)), dtype=np.float32)
for i_country, country in enumerate(self.countries):
for i_cat, (cat, values) in enumerate(self.categories.items()):
emission_factor = values['emission_factors']
......@@ -157,10 +144,7 @@ class EmisModel(object):
proxy_category_names = proxy['emis_category_names'][:]
proxy_category_names = [b''.join(category_name).strip().decode() for category_name in proxy_category_names]
proxy_country_names = proxy['country_names'][:]
proxy_country_names = [b''.join(country_name).strip().decode() for country_name in proxy_country_names]
spatial_distributions = np.zeros((self.nrcat, len(self.countries), self.area.shape[0], self.area.shape[1]))
spatial_distributions = np.zeros((self.nrcat, len(self.countries), *self.area.shape), dtype=np.float32)
for i, category in enumerate(self.categories):
spatial_name = self.categories[category]['spatial']
cat_index = proxy_category_names.index(spatial_name)
......@@ -174,7 +158,6 @@ class EmisModel(object):
spatial_emissions = spatial_distributions[:, :, None, :, :] * yremis[:, :, :, None, None] # cat, country, species, lat, lon
# Sum over the countries to overlay them.
spatial_emissions = spatial_emissions.sum(axis=1) # Indices: category, species, lat, lon
self.spatial_emissions = spatial_emissions
emissions = []
for i, category in enumerate(self.categories):
......@@ -184,7 +167,8 @@ class EmisModel(object):
temporal_profile = time_profiles[temporal_name]
emissions.append(spatial_emissions[cat_index, :, :, :] * temporal_profile[None, :, :, :])
emissions = np.asarray(emissions)
self.emissions = emissions
emissions = np.array(emissions)
emissions = np.swapaxes(emissions, 0,1) # Indices: species, category, time, lat, lon
# Recalculate spatial emissions to umol/sec/m2
......@@ -219,8 +203,23 @@ class EmisModel(object):
Returns:
dict of np.arrays:
The temporal profiles (one for each hour) for each gridcel and timestep """
# --- Settings
year = self.startdate.year
ndays = (dtm.datetime(year, 12, 31) - dtm.datetime(year - 1, 12, 31)).days
self.ndays = ndays
times = np.array([dtm.datetime(year, 1, 1, 0, 0, 0) + dtm.timedelta(hours = i) for i in range(ndays*24)])
self.times = times
datapath = '/projects/0/ctdas/RINGO/EmissionInventories/DynEmMod_TimeProfiles'
infilename = '{}/RINGO_ECMWF_DailyMeteo{}.nc'.format(datapath, year)
with nc.Dataset(infilename) as infile:
T2myr = infile.variables['T2m' ][:ndays] - 273.15 # T2m K --> oC
U10 = infile.variables['U10m' ][:ndays] # u m/s
Rinyr = infile.variables['Rsin' ][:ndays] / (24. * 1.0e4) # Radiation (?) J/m2/day --> J/cm2/hr
T2myr_av = infile.variables['T2m_avg' ][:ndays] - 273.15 #
U10_av = infile.variables['U10m_avg'][:ndays]
Rinyr_av = infile.variables['Rsin_avg'][:ndays]
ndays, nlat, nlon = self.T2myr.shape
ndays, nlat, nlon = T2myr.shape
# --- calculate degree day sum and corresponding time profiles
fr_cons = 0.2 # constant emission for consumers (cooking, warm water)
fr_gls = 0. # no constant emission for glasshouses
......@@ -239,11 +238,11 @@ class EmisModel(object):
HDC_gas = np.empty((ndays, nlat, nlon)) # Heating Demand Category 4 (Renewable activity)
for i in range(ndays):
HDC_cons[i] = np.fmax(T0_cons - self.T2myr[i, :, :], 0) # Daily demand for consumers / household heating
HDC_gls [i] = np.fmax(T0_gls - self.T2myr[i, :, :], 0) # Daily demand for glasshouse heating
HDC_coal[i] = np.fmax(T0_coal - self.T2myr[i, :, :], 0) # Daily demand for coal-fired powerplant productvity
dum1 = np.fmax(U0_gas - self.U10[i, :, :], 0) # Wind energy
dum2 = np.fmax(R0_gas - self.Rinyr[i, :, :], 0) # Solar energy
HDC_cons[i] = np.fmax(T0_cons - T2myr[i, :, :], 0) # Daily demand for consumers / household heating
HDC_gls [i] = np.fmax(T0_gls - T2myr[i, :, :], 0) # Daily demand for glasshouse heating
HDC_coal[i] = np.fmax(T0_coal - T2myr[i, :, :], 0) # Daily demand for coal-fired powerplant productvity
dum1 = np.fmax(U0_gas - U10[i, :, :], 0) # Wind energy
dum2 = np.fmax(R0_gas - Rinyr[i, :, :], 0) # Solar energy
HDC_gas [i] = dum1 * dum2
HC_cons = HDC_cons.mean(axis=0)
......@@ -381,7 +380,7 @@ class EmisModel(object):
datepoint = self.startdate
enddate = self.enddate
t_ind = np.ones((nt, nlat, nlon))
t_ind = np.ones((indices.stop, nlat, nlon), dtype=np.float32)
# Make container (dict) containing all time profiles
time_profiles = {'t_gas': t_gas[indices].astype(np.float32),
't_coal': t_coal[indices].astype(np.float32),
......
This diff is collapsed.
......@@ -10,7 +10,7 @@ restartdir : ${basepath}/restart/
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 : 2
obs.spec.nr : 1
obs.dir : obsfiles${strategy}
do.co : 0
do.c14integrated: 0
......@@ -26,6 +26,16 @@ nffparameters : 24
nbioparameters : 6
nparameters : ${nffparameters} + ${nbioparameters}
file.pft : /projects/0/ctdas/RINGO/inversions/Data/SiBPFTs.nc
! Settings for the domain:
domain.lon.start : -5.95
domain.lon.end : 19.95
domain.lat.start : 43.025
domain.lat.end : 54.975
domain.lon.num : 260
domain.lat.num : 240
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
......@@ -56,7 +66,7 @@ run.obsflag : 0
! back trajectory time of STILT footprints, also applied to OPS (in hours)
run.backtime : 24
biosphere_fluxdir : /projects/0/ctdas/RINGO/EmissionInventories/True/RINGO_ORCHIDEE_GPP_TER_dC14_old.nc
biosphere_fluxdir : /projects/0/ctdas/RINGO/inversions/Data/SiBfile.nc
files_startdate : 2016-01-01 00:00:00
! choose propagation scheme:
......
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