emissionmodel.py 19.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#!/usr/bin/env python
# stilt_tools.py

"""
Author : I. Super

Revision History:
Newly developed code, September 2017

This module holds an emission model that prepares emission files used by the observation operator and
to create pseudo-data

"""

import shutil
import os
import logging
import datetime as dtm
import numpy as np
from numpy import array, logical_and
import da.tools.io4 as io
import math
23
import pytz
24
25
26
27

import da.tools.rc as rc
from da.tools.general import create_dirs, to_datetime
import netCDF4 as nc
28
29
30
31

from multiprocessing import Pool
from functools import partial

32
33
34
identifier = 'EmissionModel ensemble '
version = '1.0'

35
#from da.ccffdas.observationoperator import get_time_indices
36
# Improve: should not be a python file? Defenitely not be stored in this folder!
37
from da.ccffdas import energy_use_country
38
39
energy_use_per_country = energy_use_country.energy_use_per_country 

40
from da.ccffdas import category_info
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
categories = category_info.categories
################### Begin Class Emission model ###################

class EmisModel(object):

    def __init__(self, dacycle=None):
        if dacycle != None:
            self.dacycle = dacycle
        else:
            self.dacycle = {}

    def setup(self, dacycle):
        self.dacycle = dacycle
        self.startdate = self.dacycle['time.fxstart']
        self.enddate = self.dacycle['time.finish']
        self.emisdir = dacycle.dasystem['datadir']
        self.inputdir = dacycle.dasystem['inputdir']
        self.proxyfile = dacycle.dasystem['emis.input.spatial']
        self.nrspc = int(dacycle.dasystem['obs.spec.nr'])
        self.species = dacycle.dasystem['obs.spec.name'].split(',')
        self.nrcat = int(dacycle.dasystem['obs.cat.nr'])
        self.nmembers = int(dacycle['da.optimizer.nmembers'])
        self.paramfile = dacycle.dasystem['emis.paramfiles']
        self.countries = [country.strip() for country in dacycle.dasystem['countries'].split(';')]
        areafile = dacycle.dasystem['area.file']
        self.area = nc.Dataset(areafile)['Area'][:]
67
        self.time_prof_file = dacycle.dasystem['file.timeprofs']
68
69
70

        self.energy_use_per_country = energy_use_per_country
        self.categories = categories
71
        self.paramdict = rc.read(dacycle.dasystem['paramdict'])
72
        self.country_mask_file = dacycle.dasystem['country.mask']
73
74


75
76
77
78
79
80
81
82
83
84
85
86
87
        self.countries_dict = { # From ISO3166-1 alpha-3 to alpha-2
                 'AUS': 'at',
                 'BEL': 'be',
                 'CZE': 'cz',
                 'FRA': 'fr',
                 'DEU': 'de',
                 'LUX': 'lu',
                 'NED': 'nl',
                 'POL': 'pl',
                 'CHE': 'cz',
                 'GBR': 'gb'}

        logging.debug('Emismodel has been set-up')
88

89
    def find_in_state(self, params, station, cat, name=None, return_index=False):
90
91
92
93
94
        """Function that finds the index in the state vector"""
        if not name:
            key = station + '.' + cat
        else:
            key = station + '.' + cat + '.' + name
Woude, Auke van der's avatar
Woude, Auke van der committed
95
        #print('Looking for {}'.format(key))
96
97
98
99
100

        if key in  self.paramdict:
            i_in_state = int(self.paramdict[key])
            if return_index: return i_in_state
            else:
101
                value = params[i_in_state]
102
103
                return value
        elif return_index: return False
Woude, Auke van der's avatar
Woude, Auke van der committed
104
        else: return 1
105

106
    def get_emis(self, dacycle, samples, indices, do_pseudo):
107
108
        """set up emission information for pseudo-obs (do_pseudo=1) and ensemble runs (do_pseudo=0)"""
        
109
110
111
        self.timestartkey = self.dacycle['time.sample.start']
        self.timefinishkey = self.dacycle['time.sample.end']
        time_profiles = self.make_time_profiles(indices=indices)
112
        self.time_profiles=  time_profiles
113
114
115
116
117
118
119
120
121
122

        pool = Pool(self.nmembers)
        # Create the function that calculates the conentration
        func = partial(self.get_emissions, dacycle, time_profiles)
        # We need to run over all members
        memberlist = list(range(0, self.nmembers))
        _ = pool.map(func, memberlist)
        pool.close()
        pool.join()
        logging.debug('Succesfully wrote prior emission files')
123
        
124
    def get_yearly_emissions(self, params):
Woude, Auke van der's avatar
Woude, Auke van der committed
125
        yremis = np.zeros((len(self.categories), len(self.countries), len(self.species)), dtype=np.float32)
126
127
128
129
        for i_country, country in enumerate(self.countries):
            for i_cat, (cat, values) in enumerate(self.categories.items()):
                emission_factor = values['emission_factors']

130
                emission_factor *= self.find_in_state(params, country, cat, 'emission_factors')
131
132
133
134
135
136
                if cat == 'Public power gas':
                    fraction_of_total = energy_use_per_country[country]['Public power ratio gas']
                    
                elif cat == 'Public power coal':
                    fraction_of_total = 1 - energy_use_per_country[country]['Public power ratio gas']
                else: fraction_of_total = values['fraction_of_total']
137
                fraction_of_total *= self.find_in_state(params, country, cat, 'fraction_of_total')
138
139
140
141
               
                e_use = self.energy_use_per_country[country][values['spatial']] 
                for i_species, specie in enumerate(self.species):
                    emission_ratio = values[specie]
142
                    emission_ratio *= self.find_in_state(params, country, cat, specie)
143
144
145
146
147
                    uncertainty = values[specie+'.uncertainty']
                    if uncertainty == 'l':
                        emission_ratio = np.exp(emission_ratio)
                    emission = e_use * fraction_of_total * emission_factor * emission_ratio
                    yremis[i_cat, i_country, i_species] = emission
148
        return yremis
149

150
    def get_emissions(self, dacycle, time_profiles, member):
151
        """read in proxy data used for spatial distribution of the gridded emissions, disaggregate yearly totals for the area"""
152
153
154
155
156
        prmfile=os.path.join(dacycle['dir.input'],'parameters.%03d.nc'%member)
        f = io.ct_read(prmfile, 'read')
        params = f.get_variable('parametervalues')
        f.close()
        yremis = self.get_yearly_emissions(params)
157
        # Create a recalculation factor from kg/km2/yr to umol/m2/sec
158
159
        M_mass = np.array([44e-9, 28e-9][:self.nrspc])
        sec_year = 24*3600.*self.ndays #seconds in a year (leapyear)
160
161
162
163
164
165
166
167
        kgperkmperyr2umolperm2pers = np.array(M_mass)[:, None, None] * sec_year * self.area[None, :, :] 

        #read in proxy data for spatial disaggregation
        infile = os.path.join(self.emisdir, self.proxyfile)
        proxy = io.ct_read(infile, method='read')
        proxy_category_names = proxy['emis_category_names'][:]
        proxy_category_names = [b''.join(category_name).strip().decode() for category_name in proxy_category_names]
         
Woude, Auke van der's avatar
Woude, Auke van der committed
168
        spatial_distributions = np.zeros((self.nrcat, len(self.countries), *self.area.shape), dtype=np.float32)
169
170
171
172
173
        for i, category in enumerate(self.categories):
            spatial_name = self.categories[category]['spatial']
            cat_index = proxy_category_names.index(spatial_name)
            for country in self.countries:
                country_index = self.countries.index(country)
174
175
176
177
178
                # Get the emission distribution for the category
                category_distribution_country = proxy['proxy_maps'][cat_index, country_index, :, :]
                spatial_distributions[i, country_index, :, :] = category_distribution_country
                 
        # Multiply spatial distributions with the yearly emissions in the country to get spatially distributed emissions
179
        spatial_emissions = spatial_distributions[:, :, None, :, :] * yremis[:, :, :, None, None] # cat, country, species, lat, lon
180
181
        # Sum over the countries to overlay them.
        spatial_emissions = spatial_emissions.sum(axis=1) # Indices: category, species, lat, lon
182
183
184
185
186
187

        emissions = []
        for i, category in enumerate(self.categories):
            spatial_name = self.categories[category]['spatial']
            temporal_name = self.categories[category]['temporal']
            temporal_profile = time_profiles[temporal_name]
188
            emissions.append(spatial_emissions[i, :, :, :] * temporal_profile[None, :, :, :])
189

Woude, Auke van der's avatar
Woude, Auke van der committed
190
191
        self.emissions = emissions
        emissions = np.array(emissions)
192
193
        emissions = np.swapaxes(emissions, 0,1) # Indices: species, category, time, lat, lon
        
194
        # Recalculate spatial emissions to umol/sec/m2
195
        emissions = emissions / kgperkmperyr2umolperm2pers[:, None, None, :, :]
196
197
198
199
200
        ## create output file
        prior_file = os.path.join(self.inputdir, 'prior_spatial_{0:03d}.nc'.format(member))
        f = io.CT_CDF(prior_file, method='create')
        dimid = f.add_dim('ncat', self.nrcat)
        dimid2 = f.add_dim('ops',2 )
201
        dimtime= f.add_dim('time', emissions.shape[2])
202
203
204
205
206
207
208
209
210
        dimlat = f.add_dim('lat', self.area.shape[0])
        dimlon = f.add_dim('lon', self.area.shape[1])
    
        #loop over all tracers
        for i, species in enumerate(self.species): 
            savedict = io.std_savedict.copy() 
            savedict['name'] = species
            savedict['long_name'] = "Spatially distributed emissions"
            savedict['units'] = "micromole/m2/s"
211
            savedict['dims'] = dimid + dimtime + dimlat + dimlon
212
213
            savedict['values'] = emissions[i, :, :, :]
            savedict['dtype'] = 'float'
214
215
            f.add_data(savedict)
        f.close()
216
217
218
219
220
221
222
223
224
225
            
    def make_time_profiles(self, indices):
        """Function that calculates the time profiles based on pre-specified
        monthly, daily and hourly profiles. Temperature and radiation affect 
        household heating and powerplant time profiles.
        Input: 
            year: int: The year for which the profiles should be calculated
        Returns: 
            dict of np.arrays:
                The temporal profiles (one for each hour) for each gridcel and timestep """
Woude, Auke van der's avatar
Woude, Auke van der committed
226
227
        # --- Settings
        year = self.startdate.year
228
229
230
231
232
        ndays_year = (dtm.datetime(year, 12, 31) - dtm.datetime(year - 1, 12, 31)).days
        self.ndays = ndays_year
        times = np.array([dtm.datetime(year, 1, 1, 0, 0, 0) + dtm.timedelta(hours = i) for i in range(ndays_year*24)])
        times_add1 = times[indices.start: indices.stop + 1]
        times = times[indices]
Woude, Auke van der's avatar
Woude, Auke van der committed
233
        self.times = times
234
        numdays = (times[-1] - times[0]).days + 1
Woude, Auke van der's avatar
Woude, Auke van der committed
235
236
237
        datapath   = '/projects/0/ctdas/RINGO/EmissionInventories/DynEmMod_TimeProfiles'
        infilename = '{}/RINGO_ECMWF_DailyMeteo{}.nc'.format(datapath, year)
        with nc.Dataset(infilename) as infile:
238
239
240
241
242
243
            T2myr    = infile.variables['T2m'     ][:ndays_year] - 273.15         # T2m           K --> oC
            U10      = infile.variables['U10m'    ][:ndays_year]                  # u             m/s
            Rinyr    = infile.variables['Rsin'    ][:ndays_year] / (24. * 1.0e4)  # Radiation (?) J/m2/day --> J/cm2/hr
            T2myr_av = infile.variables['T2m_avg' ][:ndays_year] - 273.15         #
            U10_av   = infile.variables['U10m_avg'][:ndays_year]
            Rinyr_av = infile.variables['Rsin_avg'][:ndays_year]
244
    
Woude, Auke van der's avatar
Woude, Auke van der committed
245
        ndays, nlat, nlon = T2myr.shape
246

247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
        # --- 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
        fr_coal = 0.8  #    constant emission  for coal-fired power plants
        fr_gas  = 0.1  #    constant emission  for gas-fired power plants
    
        T0_cons = 18.  # temperature threshold for consumer heating
        T0_gls  = 15.  # temperature threshold for glasshouse heating
        T0_coal = 25.  # temperature threshold for coal-fired power plants
        U0_gas  = 10.  # wind speed  threshold for gas-fired power plants
        R0_gas  = 50.  # radiation   threshold for gas-fired power plants
    
        HDC_cons    = np.empty((ndays, nlat, nlon)) # Heating Demand Category 1 (Household heating)
        HDC_gls     = np.empty((ndays, nlat, nlon)) # Heating Demand Category 2 (Glasshouse heating)
        HDC_coal    = np.empty((ndays, nlat, nlon)) # Heating Demand Category 3 (Power plants)
        HDC_gas     = np.empty((ndays, nlat, nlon)) # Heating Demand Category 4 (Renewable activity)
    
264
265
266
267
268
269
270
        for day in range(ndays): #
            HDC_cons[day] = np.fmax(T0_cons - T2myr[day, :, :], 0) # Daily demand for consumers / household heating
            HDC_gls [day] = np.fmax(T0_gls  - T2myr[day, :, :], 0) # Daily demand for glasshouse heating
            HDC_coal[day] = np.fmax(T0_coal - T2myr[day, :, :], 0) # Daily demand for coal-fired powerplant productvity
            wind_engy   = np.fmax(U0_gas  -   U10[day, :, :], 0) # Wind energy
            solar_engy  = np.fmax(R0_gas  - Rinyr[day, :, :], 0) # Solar energy
            HDC_gas [day] = wind_engy * solar_engy
Woude, Auke van der's avatar
Woude, Auke van der committed
271

272
273
274
275
276
277
        HC_cons         = HDC_cons.mean(axis=0)
        HC_gls          = HDC_gls.mean(axis=0)
        HC_coal         = HDC_coal.mean(axis=0)
        HC_gas          = HDC_gas.mean(axis=0)
    
        t_consd  = (HDC_cons + fr_cons * HC_cons) / ((1 + fr_cons) * HC_cons) # daily time profile for consumer/household heating
278
        t_glsd   = (HDC_gls  + fr_gls  * HC_gls)  / ((1 + fr_gls ) * HC_gls)  #                        glasshouse
279
        t_coald  = (HDC_coal + fr_coal * HC_coal) / ((1 + fr_coal) * HC_coal) #                        coal-fired powerplant
280
281
282
283
284
285
286
        t_gasd   = (HDC_gas  + fr_gas  * HC_gas)  / ((1 + fr_gas ) * HC_gas)  #                        gas-fired  powerplant

        #### Get the time profiles for all sectors:
        with nc.Dataset(self.time_prof_file) as ds:
            public_power_monthly= ds['FM_A'][:]
            public_power_weekly = ds['FW_A'][:]
            public_power_hourly = ds['FH_A'][:]
287
    
288
289
290
291
292
293
294
295
296
297
298
299
            industry_monthly = ds['FM_B'][:]
            industry_weekly  = np.array([1.02] * 6 + [0.88]) # TNO https://atmosphere.copernicus.eu/sites/default/files/2019-07/MACC_TNO_del_1_3_v2.pdf
            industry_hourly  = 1
        
            other_stationary_monthly = ds['FM_C'][:]
            other_stationary_hourly  = ds['FH_C'][:]
            
            road_transport_monthly    = ds['FM_F'][:]
            road_transport_weekly     = ds['FW_F'][:]
            road_transport_hourly_wkd = ds['FHwd_F'][:]
            road_transport_hourly_sat = ds['FHst_F'][:]
            road_transport_hourly_sun = ds['FHsn_F'][:]
300
    
301
302
303
            shipping_monthly = np.array([0.88, 0.92, 0.98, 1.03, 1.05, 1.06, 1.01, 1.02, 1.06, 1.05, 1.01, 0.93]) # TNO
            shipping_weekly = 1
            shipping_daily  = 1
304
    
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
        t_public_power_coal = np.zeros((len(times_add1), *self.area.shape))
        t_public_power_gas  = np.zeros_like(t_public_power_coal)
        t_industry     = np.zeros_like(t_public_power_coal)
        t_other_stat_cons   = np.zeros_like(t_public_power_coal)
        t_other_stat_gls    = np.zeros_like(t_public_power_coal)
        t_road         = np.zeros_like(t_public_power_coal)
        t_ship         = np.zeros_like(t_public_power_coal)

        for i, t in enumerate(times_add1):
            month = t.month -1 # Make index
            day = (t.day -1) % 7 # Make weekly index
            day_ind = t.day - 1  # Make index, starting at time 0
            hour = t.hour      # Hours start at 0
            weekday = t.weekday() < 6
            saturday = t.weekday() == 6
            sunday = t.weekday() == 7
            
            # Monthly time profiles            
            public_power_month_mul = public_power_monthly[month, :, :]
            industry_month_mul     = industry_monthly[month, :, :]
            other_stationary_month_mul = other_stationary_monthly[month, :, :]
            road_transport_month_mul   = road_transport_monthly[month, :, :]
            shipping_month_mul = shipping_monthly[month]
            
            # Weekly ('daily') profiles:
            public_power_day_mul_coal = t_coald[day_ind, :, :]
            public_power_day_mul_gas  = t_gasd[day_ind, :, :]
            industry_day_mul          = industry_weekly[day]
            other_stat_day_mul_cons   = t_consd[t.day, :, :] # Index should start at startday, check!
            other_stat_day_mul_gls    = t_glsd[t.day, :, :]
            road_transport_day_mul = road_transport_weekly[day, :, :]
            shipping_day_mul = shipping_weekly
            
            # Hourly profiles: 
            public_power_hour_mul = public_power_hourly[hour, :, :]
            industry_hour_mul = industry_hourly
            other_stationary_hour_mul = other_stationary_hourly[hour, :, :]
            if weekday: road_hour_mul = road_transport_hourly_wkd[hour, :, :]
            elif saturday: road_hour_mul = road_transport_hourly_sat[hour, :, :]
            elif sunday: road_hour_mul = road_transport_hourly_sun[hour, :, :]
            shipping_hour_mul = shipping_daily
            
            public_power_coal_tprof = public_power_month_mul * public_power_day_mul_coal * public_power_hour_mul
            public_power_gas_tprof  = public_power_month_mul * public_power_day_mul_gas  * public_power_hour_mul
            industry_tprof = industry_month_mul * industry_day_mul * industry_hour_mul
            other_stat_cons_tprof = other_stationary_month_mul * other_stat_day_mul_cons * other_stationary_hour_mul 
            other_stat_gls_tprof  = other_stationary_month_mul * other_stat_day_mul_gls  * other_stationary_hour_mul
            road_tprof = road_transport_month_mul * road_transport_day_mul * road_hour_mul
            shipping_tprof = shipping_month_mul * shipping_day_mul * shipping_hour_mul

            t_public_power_coal[i, :, :] = public_power_coal_tprof
            t_public_power_gas[i, :, :]  = public_power_gas_tprof
            t_industry[i, :, :] = industry_tprof 
            t_other_stat_cons[i, :, :] = other_stat_cons_tprof
            t_other_stat_gls[i, :, :]  = other_stat_gls_tprof
            t_road[i, :, :] = road_tprof
            t_ship[i, :, :] = shipping_tprof

        time_profiles = {
            't_gas': t_public_power_gas,
            't_coal': t_public_power_coal,
            't_ind': t_industry,
            't_cons': t_other_stat_cons,
            't_gls': t_other_stat_gls,
            't_road': t_road,
            't_ship': t_ship}

        # Roll the time profiles to be consisten with time zones
        with nc.Dataset(self.country_mask_file) as ds:
            masks = ds['country_mask'][:]
            country_names = [b''.join(c).decode() for c in  ds['country_names'][:]]

        for sect, profile in time_profiles.items():
            new_profile = np.zeros((len(times), *profile.shape[1:]))
            for country in self.countries:
                mask = masks[country_names.index(country)]
                country_times = mask * profile 
                
                timezone = pytz.timezone(pytz.country_timezones[self.countries_dict[country]][0])
                offset = timezone.utcoffset(times[0]).seconds//3600
        
                rolled = np.roll(country_times, -offset, 0)[:len(times)]
                new_profile += rolled

            new_profile[new_profile == 0] = 1    
            time_profiles[sect] = new_profile.astype(np.float32) 
391
392
    
        logging.debug('Time profiles created')
393

394
395
        return time_profiles

396
397
398
399
400
401
################### End Class Emission model ###################


if __name__ == "__main__":
    pass