From a93fad0fb4ec24df3783ef1a6cf2f10bc4322cd9 Mon Sep 17 00:00:00 2001 From: "Ochsner, David" Date: Thu, 20 Jun 2019 10:20:49 +0200 Subject: [PATCH] Format using black --- README.md | 8 +- config.py | 62 +-- config_CHE.py | 53 +-- config_Carbosense_TNO.py | 58 ++- config_Carbosense_carbocount.py | 92 ++-- config_Carbosense_maiolica.py | 47 +- config_Carbosense_meteotest.py | 21 +- config_EDGAR.py | 67 ++- config_VPRM.py | 31 +- config_VPRM_5km.py | 31 +- config_d1_ch.py | 56 ++- config_d1_eleni.py | 59 +-- config_d1_tno.py | 83 ++-- config_d1_zane.py | 59 +-- config_d2_ch.py | 56 ++- config_d2_tno.py | 83 ++-- config_hes-main_TNO.py | 62 ++- config_hes-main_carbocount.py | 95 ++-- config_hes-nest_TNO.py | 62 ++- config_hes-nest_carbocount.py | 95 ++-- config_tno.py | 47 +- countries.py | 29 +- country_code.py | 388 ++++++++-------- grid.py | 61 ++- hourly_emissions/produce_hourly_emi.py | 610 ++++++++++++++----------- main_ch.py | 135 +++--- main_edgar.py | 71 +-- main_tno.py | 240 +++++----- main_tno_example.py | 136 +++--- main_vprm.py | 123 ++--- make_online_emissions.py | 501 +++++++++++--------- merge_inventories.py | 80 ++-- profiles/temporal_profiles.py | 312 +++++++------ profiles/vertical_profiles.py | 115 ++--- 34 files changed, 2283 insertions(+), 1745 deletions(-) diff --git a/README.md b/README.md index d0f062e..ce270f3 100644 --- a/README.md +++ b/README.md @@ -20,6 +20,12 @@ https://github.com/C2SM-RCM/cosmo-emission-processing. It deliberately is not a fork of any of the existing repositories since that would introduce an unwanted dependency. However, the commit history will be kept intact as much as possible. +## Formatting + +This code is formatted using [black](https://black.readthedocs.io/en/stable/). +Run `find . -name '*.py' -exec black -l 80 --target-version py36 {} +` before commiting +to format your changes before commiting. + ## License -This work is licensed under the Creative Commons Attribution 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA. \ No newline at end of file +This work is licensed under the Creative Commons Attribution 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA. diff --git a/config.py b/config.py index 24d8b03..642884c 100644 --- a/config.py +++ b/config.py @@ -1,39 +1,40 @@ # "constant" paths and values for TNO tnoCamsPath = "/input/TNOMACC/CO2/TNO_CAMS_CO2_emissions_" tnoMACCIIIPath = "/input/TNOMACC/MACCIII/TNO_MACC_III_emissions_" -tno_xmin = -30. -tno_xmax = 60. -tno_ymin = 30. -tno_ymax = 72. -tno_dx = 1/8. -tno_dy = 1/16. -#tno_lons = np.arange(tno_startlon,tno_endlon+tno_dlon,tno_dlon) -#tno_lats = np.arange(tno_startlat,tno_endlat+tno_dlat,tno_dlat) -#tno_nx = round((tno_endlon-tno_startlon)/tno_dlon) + 1. -#tno_ny = round((tno_endlat-tno_startlat)/tno_dlat) + 1. - - - -#case specific parameters -species = ['CO2']#, 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx'] #among 'CO2', 'PM2.5', 'CO', 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx' - -cat_kind="SNAP" -snap = [1,2,34,5,6,70,8,9,10] #70 corresponds to all 7* -#maccversion = 'III' # use this version for TNO/MACC data -tno_snap = [ 1,2,34,5,6,71,72,73,8,9,10] +tno_xmin = -30.0 +tno_xmax = 60.0 +tno_ymin = 30.0 +tno_ymax = 72.0 +tno_dx = 1 / 8.0 +tno_dy = 1 / 16.0 +# tno_lons = np.arange(tno_startlon,tno_endlon+tno_dlon,tno_dlon) +# tno_lats = np.arange(tno_startlat,tno_endlat+tno_dlat,tno_dlat) +# tno_nx = round((tno_endlon-tno_startlon)/tno_dlon) + 1. +# tno_ny = round((tno_endlat-tno_startlat)/tno_dlat) + 1. + + +# case specific parameters +species = [ + "CO2" +] # , 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx'] #among 'CO2', 'PM2.5', 'CO', 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx' + +cat_kind = "SNAP" +snap = [1, 2, 34, 5, 6, 70, 8, 9, 10] # 70 corresponds to all 7* +# maccversion = 'III' # use this version for TNO/MACC data +tno_snap = [1, 2, 34, 5, 6, 71, 72, 73, 8, 9, 10] year = 2015 -gridname = 'Berlin-coarse_redo' +gridname = "Berlin-coarse_redo" -is60= True +is60 = True if is60: - output_path ="./testdata/temp_to_del_60/" + output_path = "./testdata/temp_to_del_60/" else: - output_path ="./testdata/temp_to_del_64/" -#invs = ['CH4_TNO','CO2_TNO','CO_TNO','NOx_TNO','Berlin'] + output_path = "./testdata/temp_to_del_64/" +# invs = ['CH4_TNO','CO2_TNO','CO_TNO','NOx_TNO','Berlin'] # Domain -#Berlin-coarse +# Berlin-coarse dx = 0.1 dy = 0.1 if is60: @@ -42,10 +43,9 @@ nx = 70 ny = 60 else: - xmin = -1.4-2*dx - ymin = 2.5-2*dy - nx = 70+4 - ny = 60+4 + xmin = -1.4 - 2 * dx + ymin = 2.5 - 2 * dy + nx = 70 + 4 + ny = 60 + 4 pollon = -170.0 pollat = 43.0 - diff --git a/config_CHE.py b/config_CHE.py index c3dcd91..9390878 100644 --- a/config_CHE.py +++ b/config_CHE.py @@ -1,46 +1,43 @@ # "constant" paths and values for TNO -#tnoCamsPath = "/project/hjm/CHE/TNO_Anthropogenic_emissions/v1_1_2018_12/TNO_6x6_GHGco_v1_1/TNO_GHGco_v1_1_year2015.nc" +# tnoCamsPath = "/project/hjm/CHE/TNO_Anthropogenic_emissions/v1_1_2018_12/TNO_6x6_GHGco_v1_1/TNO_GHGco_v1_1_year2015.nc" tnoCamsPath = "/input/TNOMACC/TNO_GHGco/Future_years_emissions/TNO_GHGco_v1_1_CIRCE_BAU_year2030.nc" tnoMACCIIIPath = tnoCamsPath -tno_xmin = -30. -tno_xmax = 60. -tno_ymin = 30. -tno_ymax = 72. -tno_dx = 1/10. -tno_dy = 1/20. +tno_xmin = -30.0 +tno_xmax = 60.0 +tno_ymin = 30.0 +tno_ymax = 72.0 +tno_dx = 1 / 10.0 +tno_dy = 1 / 20.0 -#case specific parameters -species = ['CO2']#,'CH4','CO','NOX','NMVOC'] +# case specific parameters +species = ["CO2"] # ,'CH4','CO','NOX','NMVOC'] -cat_kind="NFR_BAU" -snap = [ "A", "B", "C", "D", "E", "F", - "G", "H", "I", "J", "K", "L" ] +cat_kind = "NFR_BAU" +snap = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L"] # snap = [ "J", "L" ] -#tno_snap = [ "A", "B", "C", "D", "E", "F1","F2","F3", +# tno_snap = [ "A", "B", "C", "D", "E", "F1","F2","F3", # "G", "H", "I", "J", "K", "L" ] -tno_snap = [ "A", "B", "C", "D", "E", "F1","F2","F3", - "G", "H", "I", "J", "L" ] +tno_snap = ["A", "B", "C", "D", "E", "F1", "F2", "F3", "G", "H", "I", "J", "L"] year = 2030 -gridname = 'Europe_BAU' -output_path ="./testdata/CHE_TNO_v1_1_2018_12/CHE_TNO_online/" +gridname = "Europe_BAU" +output_path = "./testdata/CHE_TNO_v1_1_2018_12/CHE_TNO_online/" -offline=False +offline = False # Domain -#CHE_Europe domain +# CHE_Europe domain dx = 0.05 dy = 0.05 pollon = -170.0 pollat = 43.0 if not offline: - xmin = -17#-2*dx - ymin = -11#-2*dy - nx = 760#+4 - ny = 610#+4 + xmin = -17 # -2*dx + ymin = -11 # -2*dy + nx = 760 # +4 + ny = 610 # +4 else: - xmin = -17-2*dx - ymin = -11-2*dy - nx = 760+4 - ny = 610+4 - + xmin = -17 - 2 * dx + ymin = -11 - 2 * dy + nx = 760 + 4 + ny = 610 + 4 diff --git a/config_Carbosense_TNO.py b/config_Carbosense_TNO.py index 2733bee..f5377d6 100644 --- a/config_Carbosense_TNO.py +++ b/config_Carbosense_TNO.py @@ -1,29 +1,44 @@ # "constant" paths and values for TNO -tnoCamsPath = "/input/TNOMACC/TNO_GHGco/TNO_6x6_GHGco_v1_1/TNO_GHGco_v1_1_year2015.nc" +tnoCamsPath = ( + "/input/TNOMACC/TNO_GHGco/TNO_6x6_GHGco_v1_1/TNO_GHGco_v1_1_year2015.nc" +) tnoMACCIIIPath = tnoCamsPath -tno_xmin = -30. -tno_xmax = 60. -tno_ymin = 30. -tno_ymax = 72. -tno_dx = 1/10. -tno_dy = 1/20. - -#case specific parameters -species = ['CO2', 'CH4', 'CO'] - -cat_kind="NFR" - -snap = [ "A", "B", "C", "D", "E", "F", - "G", "H", "I", "J", "K", "L" ] - -tno_snap = [ "A", "B", "C", "D", "E", "F1", "F2", "F3", - "G", "H", "I", "J", "K", "L" ] +tno_xmin = -30.0 +tno_xmax = 60.0 +tno_ymin = 30.0 +tno_ymax = 72.0 +tno_dx = 1 / 10.0 +tno_dy = 1 / 20.0 + +# case specific parameters +species = ["CO2", "CH4", "CO"] + +cat_kind = "NFR" + +snap = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L"] + +tno_snap = [ + "A", + "B", + "C", + "D", + "E", + "F1", + "F2", + "F3", + "G", + "H", + "I", + "J", + "K", + "L", +] year = 2018 -gridname = 'tno_CO2_CO_CH4_1km_CH_only' -output_path ="./testdata" +gridname = "tno_CO2_CO_CH4_1km_CH_only" +output_path = "./testdata" -offline=False +offline = False # Carbosense COSMO-1 Domain dx = 0.01 @@ -40,4 +55,3 @@ ymin -= 2 * dy nx += 4 ny += 4 - diff --git a/config_Carbosense_carbocount.py b/config_Carbosense_carbocount.py index 0c6f7be..b6ab3af 100644 --- a/config_Carbosense_carbocount.py +++ b/config_Carbosense_carbocount.py @@ -1,30 +1,48 @@ # For CarboCount Swiss inventory, unit m, x is easterly, y is northly input_path = "/input/CH_EMISSIONS/CarboCountCO2/einzelgrids/" -ch_xn = 760 +ch_xn = 760 ch_yn = 500 -ch_xll = 470000 +ch_xll = 470000 ch_yll = 60000 ch_cell = 500 nodata_value = -9999 -origin = 'carbocount' -gridname = origin + '_CO2_1km' - -species = ['CO2'] - -ch_cat = [ 'bm', 'cf', 'df', 'hf', 'hk', 'if', 'ka', 'ki', 'ks', 'kv', - 'la', 'lf', 'lw', 'mi', 'mt', 'nf', 'pf', 'pq', 'rf', 'vk', - 'zp', 'zv' ] - +origin = "carbocount" +gridname = origin + "_CO2_1km" +species = ["CO2"] +ch_cat = [ + "bm", + "cf", + "df", + "hf", + "hk", + "if", + "ka", + "ki", + "ks", + "kv", + "la", + "lf", + "lw", + "mi", + "mt", + "nf", + "pf", + "pq", + "rf", + "vk", + "zp", + "zv", +] year = 2018 -output_path ="./testdata" +output_path = "./testdata" -offline=False +offline = False # Carbosense COSMO-1 Domain dx = 0.01 @@ -43,27 +61,27 @@ ny += 4 -mapping = { 'bm': 'B', - 'cf': 'B', - 'df': 'C', - 'hf': 'C', - 'hk': 'C', - 'if': 'B', - 'ka': 'J', - 'ki': 'J', - 'ks': 'J', - 'kv': 'J', - 'la': 'J', - 'lf': 'L', - 'lw': 'L', - 'mi': 'B', - 'mt': 'C', - 'nf': 'B', - 'pf': 'B', - 'pq': 'B', - 'rf': 'B', - 'vk': 'F', - 'zp': 'B', - 'zv': 'F', - } - +mapping = { + "bm": "B", + "cf": "B", + "df": "C", + "hf": "C", + "hk": "C", + "if": "B", + "ka": "J", + "ki": "J", + "ks": "J", + "kv": "J", + "la": "J", + "lf": "L", + "lw": "L", + "mi": "B", + "mt": "C", + "nf": "B", + "pf": "B", + "pq": "B", + "rf": "B", + "vk": "F", + "zp": "B", + "zv": "F", +} diff --git a/config_Carbosense_maiolica.py b/config_Carbosense_maiolica.py index 1f7aa5c..6f58d05 100644 --- a/config_Carbosense_maiolica.py +++ b/config_Carbosense_maiolica.py @@ -1,30 +1,34 @@ # for Maiolica Swiss inventory, unit m, x is easterly, y is northly input_path = "/project/jae/MaiolicaSynthesisCH4" -ch_xn = 704 +ch_xn = 704 ch_yn = 442 ch_xll = 484000 ch_yll = 75000 ch_cell = 500 nodata_value = -9999 -origin = 'maiolica' -gridname = origin + '_CH4_1km' - -species = ['CH4'] - -ch_cat = ['ara', 'dep_scaled', 'forest', 'gas_scaled', 'lakes', 'lwi_scaled', - 'wetlands', 'wildanimals'] - - +origin = "maiolica" +gridname = origin + "_CH4_1km" +species = ["CH4"] +ch_cat = [ + "ara", + "dep_scaled", + "forest", + "gas_scaled", + "lakes", + "lwi_scaled", + "wetlands", + "wildanimals", +] year = 2018 -output_path ="./testdata" +output_path = "./testdata" -offline=False +offline = False # Carbosense COSMO-1 Domain dx = 0.01 @@ -43,12 +47,13 @@ ny += 4 # Don't scale for CH4 (i.e., use time and height constant category 'I') -mapping = { 'ara': 'I', - 'dep_scaled': 'I', - 'forest': 'I', - 'gas_scaled': 'I', - 'lakes': 'I', - 'lwi_scaled': 'I', - 'wetlands': 'I', - 'wildanimals': 'I', - } +mapping = { + "ara": "I", + "dep_scaled": "I", + "forest": "I", + "gas_scaled": "I", + "lakes": "I", + "lwi_scaled": "I", + "wetlands": "I", + "wildanimals": "I", +} diff --git a/config_Carbosense_meteotest.py b/config_Carbosense_meteotest.py index 8c3c8ca..a69bcb6 100644 --- a/config_Carbosense_meteotest.py +++ b/config_Carbosense_meteotest.py @@ -1,29 +1,25 @@ # for MeteoTest Swiss inventory, unit m, x is easterly, y is northly input_path = "/input/CH_EMISSIONS/emiskat15/" -ch_xll = 480000 +ch_xll = 480000 ch_yn = 1200 ch_yll = 60000 -ch_xn = 1800 +ch_xn = 1800 ch_cell = 200 nodata_value = -9999 -origin = 'meteotest' -gridname = origin + '_CO_1km' - -species = ['CO'] - -ch_cat = [ "A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M" ] - - +origin = "meteotest" +gridname = origin + "_CO_1km" +species = ["CO"] +ch_cat = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M"] year = 2018 -output_path ="./testdata" +output_path = "./testdata" -offline=False +offline = False # Carbosense COSMO-1 Domain dx = 0.01 @@ -40,4 +36,3 @@ ymin -= 2 * dy nx += 4 ny += 4 - diff --git a/config_EDGAR.py b/config_EDGAR.py index b8b2a63..51dadf6 100644 --- a/config_EDGAR.py +++ b/config_EDGAR.py @@ -1,41 +1,60 @@ # "constant" paths and values for TNO input_path = "/input/EDGAR/v432_FT_CHE/" -edgar_xmin = -30. -edgar_xmax = 60. -edgar_ymin = 30. -edgar_ymax = 69. +edgar_xmin = -30.0 +edgar_xmax = 60.0 +edgar_ymin = 30.0 +edgar_ymax = 69.0 edgar_dx = 0.1 edgar_dy = 0.1 -edgar_nx = (edgar_xmax - edgar_xmin)/edgar_dx -edgar_ny = (edgar_ymax - edgar_ymin)/edgar_dy +edgar_nx = (edgar_xmax - edgar_xmin) / edgar_dx +edgar_ny = (edgar_ymax - edgar_ymin) / edgar_dy -#case specific parameters -species = 'CO2' +# case specific parameters +species = "CO2" -gnfr = [ "A", "B", "C", "D", "E", "F", - "G", "H", "I", "J", "K", "L" ] -edgar_cat = [ "AGS", "CHE", "ENE", "FFF", "IND", "IRO","NEU","NFE","NMM", "PRO", "PRU_SOL", "RCO", "REF_TRF", "SWD_INC", "TNR_Aviation_CDS", "TNR_Aviation_CRS", "TNR_Aviation_LTO","TNR_Other","TNR_Ship","TRO"] +gnfr = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L"] +edgar_cat = [ + "AGS", + "CHE", + "ENE", + "FFF", + "IND", + "IRO", + "NEU", + "NFE", + "NMM", + "PRO", + "PRU_SOL", + "RCO", + "REF_TRF", + "SWD_INC", + "TNR_Aviation_CDS", + "TNR_Aviation_CRS", + "TNR_Aviation_LTO", + "TNR_Other", + "TNR_Ship", + "TRO", +] year = 2015 -gridname = 'Europe' -output_path ="./testdata/EDGAR/" +gridname = "Europe" +output_path = "./testdata/EDGAR/" -offline=True +offline = True # Domain -#CHE_Europe domain +# CHE_Europe domain dx = 0.05 dy = 0.05 pollon = -170.0 pollat = 43.0 if not offline: - xmin = -17#-2*dx - ymin = -11#-2*dy - nx = 760#+4 - ny = 610#+4 + xmin = -17 # -2*dx + ymin = -11 # -2*dy + nx = 760 # +4 + ny = 610 # +4 else: - xmin = -17-2*dx - ymin = -11-2*dy - nx = 760+4 - ny = 610+4 - + xmin = -17 - 2 * dx + ymin = -11 - 2 * dy + nx = 760 + 4 + ny = 610 + 4 diff --git a/config_VPRM.py b/config_VPRM.py index 882a1db..c37d835 100644 --- a/config_VPRM.py +++ b/config_VPRM.py @@ -1,31 +1,30 @@ # "constant" paths and values for TNO vprm_path = "/scratch/snx3000/haussaij/VPRM/single_files/vprm_fluxes_EU%s_GPP_2015010110.nc" -tno_xmin = -30. -tno_xmax = 60. -tno_ymin = 30. -tno_ymax = 72. +tno_xmin = -30.0 +tno_xmax = 60.0 +tno_ymin = 30.0 +tno_ymax = 72.0 tno_dx = 1000 tno_dy = 1000 -output_path ="./testdata/VPRM/output/" +output_path = "./testdata/VPRM/output/" -offline=True +offline = True # Domain -#Berlin-coarse +# Berlin-coarse dx = 0.05 dy = 0.05 pollon = -170.0 pollat = 43.0 if not offline: - xmin = -17#-2*dx - ymin = -11#-2*dy - nx = 760#+4 - ny = 610#+4 + xmin = -17 # -2*dx + ymin = -11 # -2*dy + nx = 760 # +4 + ny = 610 # +4 else: - xmin = -17-2*dx - ymin = -11-2*dy - nx = 760+4 - ny = 610+4 - + xmin = -17 - 2 * dx + ymin = -11 - 2 * dy + nx = 760 + 4 + ny = 610 + 4 diff --git a/config_VPRM_5km.py b/config_VPRM_5km.py index 9ac395f..5fb431e 100644 --- a/config_VPRM_5km.py +++ b/config_VPRM_5km.py @@ -1,31 +1,30 @@ # "constant" paths and values for TNO vprm_path = "/scratch/snx3000/haussaij/VPRM/input/20150101/gpp_2015010110.nc" -tno_xmin = -30. -tno_xmax = 60. -tno_ymin = 30. -tno_ymax = 72. +tno_xmin = -30.0 +tno_xmax = 60.0 +tno_ymin = 30.0 +tno_ymax = 72.0 tno_dx = 5000 tno_dy = 5000 -output_path ="../output/" +output_path = "../output/" -offline=True +offline = True # Domain -#Berlin-coarse +# Berlin-coarse dx = 0.05 dy = 0.05 pollon = -170.0 pollat = 43.0 if not offline: - xmin = -17#-2*dx - ymin = -11#-2*dy - nx = 760#+4 - ny = 610#+4 + xmin = -17 # -2*dx + ymin = -11 # -2*dy + nx = 760 # +4 + ny = 610 # +4 else: - xmin = -17-2*dx - ymin = -11-2*dy - nx = 760+4 - ny = 610+4 - + xmin = -17 - 2 * dx + ymin = -11 - 2 * dy + nx = 760 + 4 + ny = 610 + 4 diff --git a/config_d1_ch.py b/config_d1_ch.py index 9c4ea9d..6ae3a4b 100644 --- a/config_d1_ch.py +++ b/config_d1_ch.py @@ -1,44 +1,50 @@ # for MeteoTest Swiss inventory, unit m, x is easterly, y is northly input_path = "/input/CH_EMISSIONS/emiskat15/" -ch_xll = 480000 +ch_xll = 480000 ch_yn = 1200 ch_yll = 60000 -ch_xn = 1800 +ch_xn = 1800 ch_cell = 200 nodata_value = -9999 -species = ['BC','CO','NH3','NMVOC','NOX','PM10','PM25','SO2']#, 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx'] #among 'CO2', 'PM2.5', 'CO', 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx' +species = [ + "BC", + "CO", + "NH3", + "NMVOC", + "NOX", + "PM10", + "PM25", + "SO2", +] # , 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx'] #among 'CO2', 'PM2.5', 'CO', 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx' -cat_kind="NFR" -#output cat -gnfr = [ "A", "B", "C", "D", "E", "F", - "G", "H", "I", "J", "K", "L", "M" ] -#input cat -ch_cat = [ "A", "B", "C", "D", "E", "F", - "G", "H", "I", "J", "K", "L", "M" ] +cat_kind = "NFR" +# output cat +gnfr = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M"] +# input cat +ch_cat = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M"] year = 2015 -gridname = 'd1_ch' -#output_path ="./testdata/d1_offline/ch/" -output_path ="./testdata/d1_online/ch/" +gridname = "d1_ch" +# output_path ="./testdata/d1_offline/ch/" +output_path = "./testdata/d1_online/ch/" -#offline=True -offline=False +# offline=True +offline = False # Domain -#Europe domain, rotated pole coordinate +# Europe domain, rotated pole coordinate dx = 0.12 dy = 0.12 pollon = -170.0 pollat = 43.0 if not offline: - xmin = -16.08#-2*dx - ymin = -9.54#-2*dy - nx = 192#+4 - ny = 164#+4 + xmin = -16.08 # -2*dx + ymin = -9.54 # -2*dy + nx = 192 # +4 + ny = 164 # +4 else: - xmin = -16.08-2*dx - ymin = -9.54-2*dy - nx = 192+4 - ny = 164+4 - + xmin = -16.08 - 2 * dx + ymin = -9.54 - 2 * dy + nx = 192 + 4 + ny = 164 + 4 diff --git a/config_d1_eleni.py b/config_d1_eleni.py index a6ef7c0..083a09f 100644 --- a/config_d1_eleni.py +++ b/config_d1_eleni.py @@ -1,41 +1,48 @@ tnoCamsPath = "/input/TNOMACC/MACCIII/TNO_MACC_III_emissions_2011.nc" -tno_xmin = -30. -tno_xmax = 60. -tno_ymin = 30. -tno_ymax = 72. -tno_dx = 1/8. -tno_dy = 1/16. +tno_xmin = -30.0 +tno_xmax = 60.0 +tno_ymin = 30.0 +tno_ymax = 72.0 +tno_dx = 1 / 8.0 +tno_dy = 1 / 16.0 -species = ['CO','NOX','NMVOC','SO2','NH3','PM10','PM25']#, 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx'] #among 'CO2', 'PM2.5', 'CO', 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx' +species = [ + "CO", + "NOX", + "NMVOC", + "SO2", + "NH3", + "PM10", + "PM25", +] # , 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx'] #among 'CO2', 'PM2.5', 'CO', 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx' -cat_kind="SNAP" -#output cat -snap = [ 1, 2, 34, 5, 6, 71, 72, 73, 74, 75, 8, 9, 10] -#input cat +cat_kind = "SNAP" +# output cat +snap = [1, 2, 34, 5, 6, 71, 72, 73, 74, 75, 8, 9, 10] +# input cat tno_snap = [1, 2, 34, 5, 6, 71, 72, 73, 74, 75, 8, 9, 10] year = 2013 -gridname = 'd1_eleni' -output_path ="./testdata/eleni/d1_online/" -#output_path ="./testdata/d1_online/tno/" +gridname = "d1_eleni" +output_path = "./testdata/eleni/d1_online/" +# output_path ="./testdata/d1_online/tno/" -#offline=True -offline=False +# offline=True +offline = False # Domain -#Europe domain, rotated pole coordinate +# Europe domain, rotated pole coordinate dx = 0.025 dy = 0.025 pollon = -156.5 pollat = 52.3 if not offline: - xmin = -4.8#-2*dx - ymin = -4.0#-2*dy - nx = 384#+4 - ny = 330#+4 + xmin = -4.8 # -2*dx + ymin = -4.0 # -2*dy + nx = 384 # +4 + ny = 330 # +4 else: - xmin = -4.8-2*dx - ymin = -4.0-2*dy - nx = 384+4 - ny = 330+4 - + xmin = -4.8 - 2 * dx + ymin = -4.0 - 2 * dy + nx = 384 + 4 + ny = 330 + 4 diff --git a/config_d1_tno.py b/config_d1_tno.py index 105ddac..6c8d6e0 100644 --- a/config_d1_tno.py +++ b/config_d1_tno.py @@ -1,45 +1,68 @@ # "constant" paths and values for TNO, regular lat/lon # for MeteoTest Swiss inventory, use calculated regular domain in the code -tnoCamsPath = "/input/TNOMACC/CAMS-REG-AP_v2_2/CAMS-REG-AP_v2_2_1_emissions_year2015.nc" -tno_xmin = -30. -tno_xmax = 60. -tno_ymin = 30. -tno_ymax = 72. -tno_dx = 1/10. -tno_dy = 1/20. +tnoCamsPath = ( + "/input/TNOMACC/CAMS-REG-AP_v2_2/CAMS-REG-AP_v2_2_1_emissions_year2015.nc" +) +tno_xmin = -30.0 +tno_xmax = 60.0 +tno_ymin = 30.0 +tno_ymax = 72.0 +tno_dx = 1 / 10.0 +tno_dy = 1 / 20.0 -species = ['CO','NOX','NMVOC','SO2','NH3','PM10','PM25']#, 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx'] #among 'CO2', 'PM2.5', 'CO', 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx' +species = [ + "CO", + "NOX", + "NMVOC", + "SO2", + "NH3", + "PM10", + "PM25", +] # , 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx'] #among 'CO2', 'PM2.5', 'CO', 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx' -cat_kind="NFR" -#output cat -snap = [ "A", "B", "C", "D", "E", "F", - "G", "H", "I", "J", "K", "L" ] -#input cat -tno_snap = [ "A", "B", "C", "D", "E", "F1","F2","F3","F4", - "G", "H", "I", "J", "K", "L" ] +cat_kind = "NFR" +# output cat +snap = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L"] +# input cat +tno_snap = [ + "A", + "B", + "C", + "D", + "E", + "F1", + "F2", + "F3", + "F4", + "G", + "H", + "I", + "J", + "K", + "L", +] year = 2015 -gridname = 'd1_tno' -#output_path ="./testdata/d1_offline/tno/" -output_path ="./testdata/d1_online/tno/" +gridname = "d1_tno" +# output_path ="./testdata/d1_offline/tno/" +output_path = "./testdata/d1_online/tno/" -#offline=True -offline=False +# offline=True +offline = False # Domain -#Europe domain, rotated pole coordinate +# Europe domain, rotated pole coordinate dx = 0.12 dy = 0.12 pollon = -170.0 pollat = 43.0 if not offline: - xmin = -16.08#-2*dx - ymin = -9.54#-2*dy - nx = 192#+4 - ny = 164#+4 + xmin = -16.08 # -2*dx + ymin = -9.54 # -2*dy + nx = 192 # +4 + ny = 164 # +4 else: - xmin = -16.08-2*dx - ymin = -9.54-2*dy - nx = 192+4 - ny = 164+4 - + xmin = -16.08 - 2 * dx + ymin = -9.54 - 2 * dy + nx = 192 + 4 + ny = 164 + 4 diff --git a/config_d1_zane.py b/config_d1_zane.py index 9692170..d786660 100644 --- a/config_d1_zane.py +++ b/config_d1_zane.py @@ -1,41 +1,48 @@ tnoCamsPath = "/input/TNOMACC/MACCIII/TNO_MACC_III_emissions_2011.nc" -tno_xmin = -30. -tno_xmax = 60. -tno_ymin = 30. -tno_ymax = 72. -tno_dx = 1/8. -tno_dy = 1/16. +tno_xmin = -30.0 +tno_xmax = 60.0 +tno_ymin = 30.0 +tno_ymax = 72.0 +tno_dx = 1 / 8.0 +tno_dy = 1 / 16.0 -species = ['CO','NOX','NMVOC','SO2','NH3','PM10','PM25']#, 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx'] #among 'CO2', 'PM2.5', 'CO', 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx' +species = [ + "CO", + "NOX", + "NMVOC", + "SO2", + "NH3", + "PM10", + "PM25", +] # , 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx'] #among 'CO2', 'PM2.5', 'CO', 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx' -cat_kind="SNAP" -#output cat -snap = [ 1, 2, 34, 5, 6, 71, 72, 73, 74, 75, 8, 9, 10] -#input cat +cat_kind = "SNAP" +# output cat +snap = [1, 2, 34, 5, 6, 71, 72, 73, 74, 75, 8, 9, 10] +# input cat tno_snap = [1, 2, 34, 5, 6, 71, 72, 73, 74, 75, 8, 9, 10] year = 2011 -gridname = 'd1_zane' -output_path ="./testdata/d1_offline/zane/" -#output_path ="./testdata/d1_online/tno/" +gridname = "d1_zane" +output_path = "./testdata/d1_offline/zane/" +# output_path ="./testdata/d1_online/tno/" -offline=True -#offline=False +offline = True +# offline=False # Domain -#Europe domain, rotated pole coordinate +# Europe domain, rotated pole coordinate dx = 0.0625 dy = 0.0625 pollon = -171.0 pollat = 42.5 if not offline: - xmin = -24.96875#-2*dx - ymin = -21.84375#-2*dy - nx = 800#+4 - ny = 700#+4 + xmin = -24.96875 # -2*dx + ymin = -21.84375 # -2*dy + nx = 800 # +4 + ny = 700 # +4 else: - xmin = -24.96875-2*dx - ymin = -21.84375-2*dy - nx = 800+4 - ny = 700+4 - + xmin = -24.96875 - 2 * dx + ymin = -21.84375 - 2 * dy + nx = 800 + 4 + ny = 700 + 4 diff --git a/config_d2_ch.py b/config_d2_ch.py index 5946905..b67fb85 100644 --- a/config_d2_ch.py +++ b/config_d2_ch.py @@ -1,44 +1,50 @@ # for MeteoTest Swiss inventory, unit m, x is easterly, y is northly input_path = "/input/CH_EMISSIONS/emiskat15/" -ch_xll = 480000 +ch_xll = 480000 ch_yn = 1200 ch_yll = 60000 -ch_xn = 1800 +ch_xn = 1800 ch_cell = 200 nodata_value = -9999 -species = ['BC','CO','NH3','NMVOC','NOX','PM10','PM25','SO2']#, 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx'] #among 'CO2', 'PM2.5', 'CO', 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx' +species = [ + "BC", + "CO", + "NH3", + "NMVOC", + "NOX", + "PM10", + "PM25", + "SO2", +] # , 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx'] #among 'CO2', 'PM2.5', 'CO', 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx' -cat_kind="NFR" -#output cat -gnfr = [ "A", "B", "C", "D", "E", "F", - "G", "H", "I", "J", "K", "L", "M" ] -#input cat -ch_cat = [ "A", "B", "C", "D", "E", "F", - "G", "H", "I", "J", "K", "L", "M" ] +cat_kind = "NFR" +# output cat +gnfr = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M"] +# input cat +ch_cat = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M"] year = 2015 -gridname = 'd2_ch' -#output_path ="./testdata/d2_offline/ch/" -output_path ="./testdata/d2_online/ch/" +gridname = "d2_ch" +# output_path ="./testdata/d2_offline/ch/" +output_path = "./testdata/d2_online/ch/" -#offline=True -offline=False +# offline=True +offline = False # Domain -#Europe domain, rotated pole coordinate +# Europe domain, rotated pole coordinate dx = 0.02 dy = 0.02 pollon = -170.0 pollat = 43.0 if not offline: - xmin = -3.48#-2*dx - ymin = -1.74#-2*dy - nx = 220#+4 - ny = 142#+4 + xmin = -3.48 # -2*dx + ymin = -1.74 # -2*dy + nx = 220 # +4 + ny = 142 # +4 else: - xmin = -3.48-2*dx - ymin = -1.74-2*dy - nx = 220+4 - ny = 142+4 - + xmin = -3.48 - 2 * dx + ymin = -1.74 - 2 * dy + nx = 220 + 4 + ny = 142 + 4 diff --git a/config_d2_tno.py b/config_d2_tno.py index 3846ebb..2b143dc 100644 --- a/config_d2_tno.py +++ b/config_d2_tno.py @@ -1,46 +1,69 @@ # "constant" paths and values for TNO, regular lat/lon # for MeteoTest Swiss inventory, use calculated regular domain in the code -tnoCamsPath = "/input/TNOMACC/CAMS-REG-AP_v2_2/CAMS-REG-AP_v2_2_1_emissions_year2015.nc" +tnoCamsPath = ( + "/input/TNOMACC/CAMS-REG-AP_v2_2/CAMS-REG-AP_v2_2_1_emissions_year2015.nc" +) tnoMACCIIIPath = tnoCamsPath -tno_xmin = -30. -tno_xmax = 60. -tno_ymin = 30. -tno_ymax = 72. -tno_dx = 1/10. -tno_dy = 1/20. +tno_xmin = -30.0 +tno_xmax = 60.0 +tno_ymin = 30.0 +tno_ymax = 72.0 +tno_dx = 1 / 10.0 +tno_dy = 1 / 20.0 -species = ['CO','NOX','NMVOC','SO2','NH3','PM10','PM25']#, 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx'] #among 'CO2', 'PM2.5', 'CO', 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx' +species = [ + "CO", + "NOX", + "NMVOC", + "SO2", + "NH3", + "PM10", + "PM25", +] # , 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx'] #among 'CO2', 'PM2.5', 'CO', 'PM10', 'CH4', 'SO2', 'NMVOC', 'NH3', 'NOx' -cat_kind="NFR" -#output cat -snap = [ "A", "B", "C", "D", "E", "F", - "G", "H", "I", "J", "K", "L" ] -#input cat -tno_snap = [ "A", "B", "C", "D", "E", "F1","F2","F3","F4", - "G", "H", "I", "J", "K", "L" ] +cat_kind = "NFR" +# output cat +snap = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L"] +# input cat +tno_snap = [ + "A", + "B", + "C", + "D", + "E", + "F1", + "F2", + "F3", + "F4", + "G", + "H", + "I", + "J", + "K", + "L", +] year = 2015 -gridname = 'd2_tno' -output_path ="./testdata/d2_offline/tno/" -#output_path ="./testdata/d2_online/tno/" +gridname = "d2_tno" +output_path = "./testdata/d2_offline/tno/" +# output_path ="./testdata/d2_online/tno/" -offline=True -#offline=False +offline = True +# offline=False # Domain -#Europe domain, rotated pole coordinate +# Europe domain, rotated pole coordinate dx = 0.02 dy = 0.02 pollon = -170.0 pollat = 43.0 if not offline: - xmin = -3.48#-2*dx - ymin = -1.74#-2*dy - nx = 220#+4 - ny = 142#+4 + xmin = -3.48 # -2*dx + ymin = -1.74 # -2*dy + nx = 220 # +4 + ny = 142 # +4 else: - xmin = -3.48-2*dx - ymin = -1.74-2*dy - nx = 220+4 - ny = 142+4 - + xmin = -3.48 - 2 * dx + ymin = -1.74 - 2 * dy + nx = 220 + 4 + ny = 142 + 4 diff --git a/config_hes-main_TNO.py b/config_hes-main_TNO.py index 13542e0..5b3f8a5 100644 --- a/config_hes-main_TNO.py +++ b/config_hes-main_TNO.py @@ -1,35 +1,50 @@ # "constant" paths and values for TNO -tnoCamsPath = "/input/TNOMACC/TNO_GHGco/TNO_6x6_GHGco_v1_1/TNO_GHGco_v1_1_year2015.nc" +tnoCamsPath = ( + "/input/TNOMACC/TNO_GHGco/TNO_6x6_GHGco_v1_1/TNO_GHGco_v1_1_year2015.nc" +) tnoMACCIIIPath = tnoCamsPath -tno_xmin = -30. -tno_xmax = 60. -tno_ymin = 30. -tno_ymax = 72. -tno_dx = 1/10. -tno_dy = 1/20. - -#case specific parameters -species = ['CO2'] - -cat_kind="NFR" - -snap = [ "A", "B", "C", "D", "E", "F", - "G", "H", "I", "J", "K", "L" ] - -tno_snap = [ "A", "B", "C", "D", "E", "F1", "F2", "F3", - "G", "H", "I", "J", "K", "L" ] +tno_xmin = -30.0 +tno_xmax = 60.0 +tno_ymin = 30.0 +tno_ymax = 72.0 +tno_dx = 1 / 10.0 +tno_dy = 1 / 20.0 + +# case specific parameters +species = ["CO2"] + +cat_kind = "NFR" + +snap = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L"] + +tno_snap = [ + "A", + "B", + "C", + "D", + "E", + "F1", + "F2", + "F3", + "G", + "H", + "I", + "J", + "K", + "L", +] year = 2018 -gridname = 'tno_CO2_FLEXPART_main' -output_path ="./testdata" +gridname = "tno_CO2_FLEXPART_main" +output_path = "./testdata" -offline=False +offline = False # FLEXPART Main dx = 0.16 dy = 0.12 -pollon = 0.0 # non-rotated grid -pollat = 90.0 # non-rotated grid +pollon = 0.0 # non-rotated grid +pollat = 90.0 # non-rotated grid xmin = -11.92 ymin = 36.06 nx = 207 @@ -40,4 +55,3 @@ ymin -= 2 * dy nx += 4 ny += 4 - diff --git a/config_hes-main_carbocount.py b/config_hes-main_carbocount.py index d6fc75a..2e71391 100644 --- a/config_hes-main_carbocount.py +++ b/config_hes-main_carbocount.py @@ -1,35 +1,54 @@ # For CarboCount Swiss inventory, unit m, x is easterly, y is northly input_path = "/input/CH_EMISSIONS/CarboCountCO2/einzelgrids/" -ch_xn = 760 +ch_xn = 760 ch_yn = 500 -ch_xll = 470000 +ch_xll = 470000 ch_yll = 60000 ch_cell = 500 nodata_value = -9999 -origin = 'carbocount' -gridname = origin + '_CO2_FLEXPART_main' - -species = ['CO2'] - -ch_cat = [ 'bm', 'cf', 'df', 'hf', 'hk', 'if', 'ka', 'ki', 'ks', 'kv', - 'la', 'lf', 'lw', 'mi', 'mt', 'nf', 'pf', 'pq', 'rf', 'vk', - 'zp', 'zv' ] +origin = "carbocount" +gridname = origin + "_CO2_FLEXPART_main" +species = ["CO2"] +ch_cat = [ + "bm", + "cf", + "df", + "hf", + "hk", + "if", + "ka", + "ki", + "ks", + "kv", + "la", + "lf", + "lw", + "mi", + "mt", + "nf", + "pf", + "pq", + "rf", + "vk", + "zp", + "zv", +] year = 2018 -output_path ="./testdata" +output_path = "./testdata" -offline=False +offline = False # FLEXPART Main dx = 0.16 dy = 0.12 -pollon = 180.0 # non-rotated grid -pollat = 90.0 # non-rotated grid +pollon = 180.0 # non-rotated grid +pollat = 90.0 # non-rotated grid xmin = -11.92 ymin = 36.06 nx = 207 @@ -42,27 +61,27 @@ ny += 4 -mapping = { 'bm': 'B', - 'cf': 'B', - 'df': 'C', - 'hf': 'C', - 'hk': 'C', - 'if': 'B', - 'ka': 'J', - 'ki': 'J', - 'ks': 'J', - 'kv': 'J', - 'la': 'J', - 'lf': 'L', - 'lw': 'L', - 'mi': 'B', - 'mt': 'C', - 'nf': 'B', - 'pf': 'B', - 'pq': 'B', - 'rf': 'B', - 'vk': 'F', - 'zp': 'B', - 'zv': 'F', - } - +mapping = { + "bm": "B", + "cf": "B", + "df": "C", + "hf": "C", + "hk": "C", + "if": "B", + "ka": "J", + "ki": "J", + "ks": "J", + "kv": "J", + "la": "J", + "lf": "L", + "lw": "L", + "mi": "B", + "mt": "C", + "nf": "B", + "pf": "B", + "pq": "B", + "rf": "B", + "vk": "F", + "zp": "B", + "zv": "F", +} diff --git a/config_hes-nest_TNO.py b/config_hes-nest_TNO.py index 5f1fe16..a51e4dc 100644 --- a/config_hes-nest_TNO.py +++ b/config_hes-nest_TNO.py @@ -1,35 +1,50 @@ # "constant" paths and values for TNO -tnoCamsPath = "/input/TNOMACC/TNO_GHGco/TNO_6x6_GHGco_v1_1/TNO_GHGco_v1_1_year2015.nc" +tnoCamsPath = ( + "/input/TNOMACC/TNO_GHGco/TNO_6x6_GHGco_v1_1/TNO_GHGco_v1_1_year2015.nc" +) tnoMACCIIIPath = tnoCamsPath -tno_xmin = -30. -tno_xmax = 60. -tno_ymin = 30. -tno_ymax = 72. -tno_dx = 1/10. -tno_dy = 1/20. - -#case specific parameters -species = ['CO2'] - -cat_kind="NFR" - -snap = [ "A", "B", "C", "D", "E", "F", - "G", "H", "I", "J", "K", "L" ] - -tno_snap = [ "A", "B", "C", "D", "E", "F1", "F2", "F3", - "G", "H", "I", "J", "K", "L" ] +tno_xmin = -30.0 +tno_xmax = 60.0 +tno_ymin = 30.0 +tno_ymax = 72.0 +tno_dx = 1 / 10.0 +tno_dy = 1 / 20.0 + +# case specific parameters +species = ["CO2"] + +cat_kind = "NFR" + +snap = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L"] + +tno_snap = [ + "A", + "B", + "C", + "D", + "E", + "F1", + "F2", + "F3", + "G", + "H", + "I", + "J", + "K", + "L", +] year = 2018 -gridname = 'tno_CO2_FLEXPART_Nest' -output_path ="./testdata" +gridname = "tno_CO2_FLEXPART_Nest" +output_path = "./testdata" -offline=False +offline = False # FLEXPART Nest dx = 0.02 dy = 0.015 -pollon = 180.0 # non-rotated grid -pollat = 90.0 # non-rotated grid +pollon = 180.0 # non-rotated grid +pollat = 90.0 # non-rotated grid xmin = 4.97 ymin = 45.4875 nx = 305 @@ -40,4 +55,3 @@ ymin -= 2 * dy nx += 4 ny += 4 - diff --git a/config_hes-nest_carbocount.py b/config_hes-nest_carbocount.py index 4939861..2633dfa 100644 --- a/config_hes-nest_carbocount.py +++ b/config_hes-nest_carbocount.py @@ -1,35 +1,54 @@ # For CarboCount Swiss inventory, unit m, x is easterly, y is northly input_path = "/input/CH_EMISSIONS/CarboCountCO2/einzelgrids/" -ch_xn = 760 +ch_xn = 760 ch_yn = 500 -ch_xll = 470000 +ch_xll = 470000 ch_yll = 60000 ch_cell = 500 nodata_value = -9999 -origin = 'carbocount' -gridname = origin + '_CO2_FLEXPART_Nest' - -species = ['CO2'] - -ch_cat = [ 'bm', 'cf', 'df', 'hf', 'hk', 'if', 'ka', 'ki', 'ks', 'kv', - 'la', 'lf', 'lw', 'mi', 'mt', 'nf', 'pf', 'pq', 'rf', 'vk', - 'zp', 'zv' ] +origin = "carbocount" +gridname = origin + "_CO2_FLEXPART_Nest" +species = ["CO2"] +ch_cat = [ + "bm", + "cf", + "df", + "hf", + "hk", + "if", + "ka", + "ki", + "ks", + "kv", + "la", + "lf", + "lw", + "mi", + "mt", + "nf", + "pf", + "pq", + "rf", + "vk", + "zp", + "zv", +] year = 2018 -output_path ="./testdata" +output_path = "./testdata" -offline=False +offline = False # FLEXPART Nest dx = 0.02 dy = 0.015 -pollon = 180.0 # non-rotated grid -pollat = 90.0 # non-rotated grid +pollon = 180.0 # non-rotated grid +pollat = 90.0 # non-rotated grid xmin = 4.97 ymin = 45.4875 nx = 305 @@ -42,27 +61,27 @@ ny += 4 -mapping = { 'bm': 'B', - 'cf': 'B', - 'df': 'C', - 'hf': 'C', - 'hk': 'C', - 'if': 'B', - 'ka': 'J', - 'ki': 'J', - 'ks': 'J', - 'kv': 'J', - 'la': 'J', - 'lf': 'L', - 'lw': 'L', - 'mi': 'B', - 'mt': 'C', - 'nf': 'B', - 'pf': 'B', - 'pq': 'B', - 'rf': 'B', - 'vk': 'F', - 'zp': 'B', - 'zv': 'F', - } - +mapping = { + "bm": "B", + "cf": "B", + "df": "C", + "hf": "C", + "hk": "C", + "if": "B", + "ka": "J", + "ki": "J", + "ks": "J", + "kv": "J", + "la": "J", + "lf": "L", + "lw": "L", + "mi": "B", + "mt": "C", + "nf": "B", + "pf": "B", + "pq": "B", + "rf": "B", + "vk": "F", + "zp": "B", + "zv": "F", +} diff --git a/config_tno.py b/config_tno.py index 435c4a5..1948251 100644 --- a/config_tno.py +++ b/config_tno.py @@ -1,17 +1,31 @@ # TNO inventory tnofile = "/input/TNOMACC/TNO_GHGco/Future_years_emissions/TNO_GHGco_v1_1_CIRCE_BAU_year2030.nc" -tno_xmin = -30. -tno_xmax = 60. -tno_ymin = 30. -tno_ymax = 72. -tno_dx = 1/10. -tno_dy = 1/20. - -cat_kind="NFR" - -tno_cat = [ "A", "B", "C", "D", "E", "F1","F2","F3", - "G", "H", "I", "J", "K","L" ] +tno_xmin = -30.0 +tno_xmax = 60.0 +tno_ymin = 30.0 +tno_ymax = 72.0 +tno_dx = 1 / 10.0 +tno_dy = 1 / 20.0 + +cat_kind = "NFR" + +tno_cat = [ + "A", + "B", + "C", + "D", + "E", + "F1", + "F2", + "F3", + "G", + "H", + "I", + "J", + "K", + "L", +] # COSMO domain xmin = -17 @@ -23,17 +37,16 @@ pollon = -170.0 pollat = 43.0 -offline=False +offline = False if offline: - xmin -= 2*dx - ymin -= 2*dy + xmin -= 2 * dx + ymin -= 2 * dy nx += 4 ny += 4 -species = ['co2_ff'] -output_cat = [ "A", "B"] +species = ["co2_ff"] +output_cat = ["A", "B"] output_path = "./testdata/oae_paper/" output_name = "tno.nc" - diff --git a/countries.py b/countries.py index 0380a15..2f25543 100644 --- a/countries.py +++ b/countries.py @@ -3,50 +3,57 @@ from osgeo import ogr + class Point(object): """ Wrapper for ogr point """ + def __init__(self, lat, lng): """ Coordinates are in degrees """ self.point = ogr.Geometry(ogr.wkbPoint) self.point.AddPoint(lng, lat) - + def getOgr(self): return self.point + ogr = property(getOgr) + class Country(object): """ Wrapper for ogr country shape. Not meant to be instantiated directly. """ + def __init__(self, shape): self.shape = shape - + def getIso(self): - return self.shape.GetField('ISO2') + return self.shape.GetField("ISO2") + iso = property(getIso) - + def __str__(self): - return self.shape.GetField('NAME') - + return self.shape.GetField("NAME") + def contains(self, point): return self.shape.geometry().Contains(point.ogr) + class CountryChecker(object): """ Loads a country shape file, checks coordinates for country location. """ - + def __init__(self, country_file): - driver = ogr.GetDriverByName('ESRI Shapefile') + driver = ogr.GetDriverByName("ESRI Shapefile") self.countryFile = driver.Open(country_file) self.layer = self.countryFile.GetLayer() - + def getCountry(self, point): """ Checks given gps-incoming coordinates for country. Output is either country shape index or None """ - + for i in range(self.layer.GetFeatureCount()): country = self.layer.GetFeature(i) if country.geometry().Contains(point.ogr): return Country(country) - + # nothing found return None diff --git a/country_code.py b/country_code.py index a376976..b112314 100644 --- a/country_code.py +++ b/country_code.py @@ -1,197 +1,197 @@ -#constants +# constants country_codes = { - 'ALB': 1, - 'AUT': 2, - 'BEL': 3, - 'BGR': 4, - 'DNK': 6, - 'FIN': 7, - 'FRA': 8, - 'FGD': 9, - 'FFR': 10, - 'GRC': 11, - 'HUN': 12, - 'IRL': 14, - 'ITA': 15, - 'LUX': 16, - 'NLD': 17, - 'NOR': 18, - 'POL': 19, - 'PRT': 20, - 'ROM': 21, - 'ROU': 21, - 'ESP': 22, - 'SWE': 23, - 'CHE': 24, - 'TUR': 25, - 'GBR': 27, - 'BLR': 39, - 'UKR': 40, - 'MKD': 41, - 'MDA': 42, - 'EST': 43, - 'LVA': 44, - 'LTU': 45, - 'CZE': 46, - 'SVK': 47, - 'SVN': 48, - 'HRV': 49, - 'BIH': 50, - 'YUG': 51, - 'GEO': 54, - 'MLT': 57, - 'DEU': 60, - 'RUS': 61, - 'ARM': 56, - 'AZE': 58, - 'CYP': 55, - 'ISL': 18, # added for AQMEII, assumed to be same as Norway + "ALB": 1, + "AUT": 2, + "BEL": 3, + "BGR": 4, + "DNK": 6, + "FIN": 7, + "FRA": 8, + "FGD": 9, + "FFR": 10, + "GRC": 11, + "HUN": 12, + "IRL": 14, + "ITA": 15, + "LUX": 16, + "NLD": 17, + "NOR": 18, + "POL": 19, + "PRT": 20, + "ROM": 21, + "ROU": 21, + "ESP": 22, + "SWE": 23, + "CHE": 24, + "TUR": 25, + "GBR": 27, + "BLR": 39, + "UKR": 40, + "MKD": 41, + "MDA": 42, + "EST": 43, + "LVA": 44, + "LTU": 45, + "CZE": 46, + "SVK": 47, + "SVN": 48, + "HRV": 49, + "BIH": 50, + "YUG": 51, + "GEO": 54, + "MLT": 57, + "DEU": 60, + "RUS": 61, + "ARM": 56, + "AZE": 58, + "CYP": 55, + "ISL": 18, # added for AQMEII, assumed to be same as Norway # ISO 2 (from EMEP site) - 'AL': 1, # Albania - 'AT': 2, # Austria - 'BE': 3, # Belgium - 'BG': 4, # Bulgaria - 'FCS': 5, # Former Czechoslovakia - 'DK': 6, # Denmark - 'FI': 7, # Finland - 'FR': 8, # France - 'FGD': 9, # Former German Democratic Republic - 'FFR': 10, # Former Federal Republic of Germany - 'GR': 11, # Greece - 'HU': 12, # Hungary - 'IS': 13, # Iceland - 'IE': 14, # Ireland - 'IT': 15, # Italy - 'LU': 16, # Luxembourg - 'NL': 17, # Netherlands - 'NO': 18, # Norway - 'PL': 19, # Poland - 'PT': 20, # Portugal - 'RO': 21, # Romania - 'ES': 22, # Spain - 'AD': 22, # Andorra (assigned to Spain) - 'SE': 23, # Sweden - 'CH': 24, # Switzerland - 'TR': 25, # Turkey - 'FSU': 26, # Former USSR - 'GB': 27, # United Kingdom - 'VOL': 28, # Volcanic emissions - 'REM': 29, # Remaining land Areas - 'BAS': 30, # Baltic Sea - 'NOS': 31, # North Sea - 'ATL': 32, # Remaining North-East Atlantic Ocean - 'MED': 33, # Mediterranean Sea - 'BLS': 34, # Black Sea - 'NAT': 35, # Natural marine emissions - 'RUO': 36, # Kola & Karelia - 'RUP': 37, # St.Petersburg & Novgorod-Pskov - 'RUA': 38, # Kaliningrad - 'BY': 39, # Belarus - 'UA': 40, # Ukraine - 'MD': 41, # Republic of Moldova - 'RUR': 42, # Rest of the Russian Federation - 'EE': 43, # Estonia - 'LV': 44, # Latvia - 'LT': 45, # Lithuania - 'CZ': 46, # Czech Republic - 'SK': 47, # Slovakia - 'SI': 48, # Slovenia - 'HR': 49, # Croatia - 'BA': 50, # Bosnia and Herzegovina - 'CS': 51, # Serbia and Montenegro - 'MK': 52, # The former Yugoslav Republic of Macedonia - 'KZ': 53, # Kazakhstan in the former official EMEP domain - 'GE': 54, # Georgia - 'CY': 55, # Cyprus - 'CYN': 55, # Cyprus (HJM) - 'AM': 56, # Armenia - 'MT': 57, # Malta - 'ASI': 58, # Remaining Asian areas - 'LI': 59, # Liechtenstein - 'DE': 60, # Germany - 'RU': 61, # Russian Federation in the former official EMEP domain - 'MC': 62, # Monaco - 'NOA': 63, # North Africa - 'MAR': 63, # Maroko (HJM) - 'TUN': 63, # Tunisia - 'DZA': 63, # Algeria (HJM) - 'EU': 64, # European Community - 'US': 65, # United States - 'CA': 66, # Canada - 'BIC': 67, # Boundary and Initial Conditions - 'KG': 68, # Kyrgyzstan - 'AZ': 69, # Azerbaijan - 'ATX': 70, # EMEP-external Remaining North-East Atlantic Ocean - 'RUX': 71, # EMEP-external part of Russian Federation - 'RS': 72, # Serbia - 'SRB': 72, # Serbia (HJM) - 'KOS': 72, # Kosovo (HJM) - 'ME': 73, # Montenegro - 'MNE':73, # Montenegro (HJM) - 'RFE': 74, # Rest of Russian Federation in the extended EMEP domain - 'KZE': 75, # Rest of Kazakhstan in the extended EMEP domain - 'UZO': 76, # Uzbekistan in the former official EMEP domain - 'TMO': 77, # Turkmenistan in the former official EMEP domain - 'UZE': 78, # Rest of Uzbekistan in the extended EMEP domain - 'TME': 79, # Rest of Turkmenistan in the extended EMEP domain - 'CAS': 80, # Caspian Sea - 'TJ': 81, # Tajikistan - 'ARO': 82, # Aral Lake in the former official EMEP domain - 'ARE': 83, # Rest of Aral Lake in the extended EMEP domain - 'ASM': 84, # Modified Remaining Asian Areas in the former official EMEP domain - 'ASE': 85, # Remaining Asian Areas in the extended EMEP domain - 'AOE': 86, # Arctic Ocean in the extended EMEP domain - 'RFX': 87, # Extended EMEP External Part of Russian Federation - 'ASX': 88, # Extended EMEP External Part of Asia - 'PAX': 89, # Extended EMEP External Part of Pacific Ocean - 'AOX': 90, # Extended EMEP External Part of Arctic Ocean - 'NAX': 91, # Extended EMEP External Part of North Africa - 'KZT': 92, # Kazakhstan - 'RUE': 93, # Russian Federation in the extended EMEP domain (RU+RFE+RUX) - 'UZ': 94, # Uzbekistan - 'TM': 95, # Turkmenistan - 'AST': 96, # Asian areas in the extended EMEP domain (ASM+ASE+ARO+ARE+CAS) - 'FYU': 99, # Former Yugoslavia - 'BEF': 301, # Belgium (Flanders) - 'BA2': 302, # Baltic Sea EU Cargo o12m - 'BA3': 303, # Baltic Sea ROW Cargo o12m - 'BA4': 304, # Baltic Sea EU Cargo i12m - 'BA5': 305, # Baltic Sea ROW Cargo i12m - 'BA6': 306, # Baltic Sea EU Ferry o12m - 'BA7': 307, # Baltic Sea ROW Ferry o12m - 'BA8': 308, # Baltic Sea EU Ferry i12m - 'BA9': 309, # Baltic Sea ROW Ferry i12m - 'NO2': 312, # North Sea EU Cargo o12m - 'NO3': 313, # North Sea ROW Cargo o12m - 'NO4': 314, # North Sea EU Cargo i12m - 'NO5': 315, # North Sea ROW Cargo i12m - 'NO6': 316, # North Sea EU Ferry o12m - 'NO7': 317, # North Sea ROW Ferry o12m - 'NO8': 318, # North Sea EU Ferry i12m - 'NO9': 319, # North Sea ROW Ferry i12m - 'AT2': 322, # Remaining North-East Atlantic Ocean EU Cargo 'o1': 2 # m - 'AT3': 323, # Remaining North-East Atlantic Ocean ROW Cargo 'o1': 2 # m - 'AT4': 324, # Remaining North-East Atlantic Ocean EU Cargo 'i1': 2 # m - 'AT5': 325, # Remaining North-East Atlantic Ocean ROW Cargo 'i1': 2 # m - 'AT6': 326, # Remaining North-East Atlantic Ocean EU Ferry 'o1': 2 # m - 'AT7': 327, # Remaining North-East Atlantic Ocean ROW Ferry 'o1': 2 # m - 'AT8': 328, # Remaining North-East Atlantic Ocean EU Ferry 'i1': 2 # m - 'AT9': 329, # Remaining North-East Atlantic Ocean ROW Ferry 'i1': 2 # m - 'ME2': 332, # Mediterranean Sea EU Cargo o12m - 'ME3': 333, # Mediterranean Sea ROW Cargo o12m - 'ME4': 334, # Mediterranean Sea EU Cargo i12m - 'ME5': 335, # Mediterranean Sea ROW Cargo i12m - 'ME6': 336, # Mediterranean Sea EU Ferry o12m - 'ME7': 337, # Mediterranean Sea ROW Ferry o12m - 'ME8': 338, # Mediterranean Sea EU Ferry i12m - 'ME9': 339, # Mediterranean Sea ROW Ferry i12m - 'BL2': 342, # Black Sea EU Cargo o12m - 'BL3': 343, # Black Sea ROW Cargo o12m - 'BL4': 344, # Black Sea EU Cargo i12m - 'BL5': 345, # Black Sea ROW Cargo i12m - 'BL6': 346, # Black Sea EU Ferry o12m - 'BL7': 347, # Black Sea ROW Ferry o12m - 'BL8': 348, # Black Sea EU Ferry i12m - 'BL9': 349, # Black Sea ROW Ferry i12m - 'GL': 601, # Greenland + "AL": 1, # Albania + "AT": 2, # Austria + "BE": 3, # Belgium + "BG": 4, # Bulgaria + "FCS": 5, # Former Czechoslovakia + "DK": 6, # Denmark + "FI": 7, # Finland + "FR": 8, # France + "FGD": 9, # Former German Democratic Republic + "FFR": 10, # Former Federal Republic of Germany + "GR": 11, # Greece + "HU": 12, # Hungary + "IS": 13, # Iceland + "IE": 14, # Ireland + "IT": 15, # Italy + "LU": 16, # Luxembourg + "NL": 17, # Netherlands + "NO": 18, # Norway + "PL": 19, # Poland + "PT": 20, # Portugal + "RO": 21, # Romania + "ES": 22, # Spain + "AD": 22, # Andorra (assigned to Spain) + "SE": 23, # Sweden + "CH": 24, # Switzerland + "TR": 25, # Turkey + "FSU": 26, # Former USSR + "GB": 27, # United Kingdom + "VOL": 28, # Volcanic emissions + "REM": 29, # Remaining land Areas + "BAS": 30, # Baltic Sea + "NOS": 31, # North Sea + "ATL": 32, # Remaining North-East Atlantic Ocean + "MED": 33, # Mediterranean Sea + "BLS": 34, # Black Sea + "NAT": 35, # Natural marine emissions + "RUO": 36, # Kola & Karelia + "RUP": 37, # St.Petersburg & Novgorod-Pskov + "RUA": 38, # Kaliningrad + "BY": 39, # Belarus + "UA": 40, # Ukraine + "MD": 41, # Republic of Moldova + "RUR": 42, # Rest of the Russian Federation + "EE": 43, # Estonia + "LV": 44, # Latvia + "LT": 45, # Lithuania + "CZ": 46, # Czech Republic + "SK": 47, # Slovakia + "SI": 48, # Slovenia + "HR": 49, # Croatia + "BA": 50, # Bosnia and Herzegovina + "CS": 51, # Serbia and Montenegro + "MK": 52, # The former Yugoslav Republic of Macedonia + "KZ": 53, # Kazakhstan in the former official EMEP domain + "GE": 54, # Georgia + "CY": 55, # Cyprus + "CYN": 55, # Cyprus (HJM) + "AM": 56, # Armenia + "MT": 57, # Malta + "ASI": 58, # Remaining Asian areas + "LI": 59, # Liechtenstein + "DE": 60, # Germany + "RU": 61, # Russian Federation in the former official EMEP domain + "MC": 62, # Monaco + "NOA": 63, # North Africa + "MAR": 63, # Maroko (HJM) + "TUN": 63, # Tunisia + "DZA": 63, # Algeria (HJM) + "EU": 64, # European Community + "US": 65, # United States + "CA": 66, # Canada + "BIC": 67, # Boundary and Initial Conditions + "KG": 68, # Kyrgyzstan + "AZ": 69, # Azerbaijan + "ATX": 70, # EMEP-external Remaining North-East Atlantic Ocean + "RUX": 71, # EMEP-external part of Russian Federation + "RS": 72, # Serbia + "SRB": 72, # Serbia (HJM) + "KOS": 72, # Kosovo (HJM) + "ME": 73, # Montenegro + "MNE": 73, # Montenegro (HJM) + "RFE": 74, # Rest of Russian Federation in the extended EMEP domain + "KZE": 75, # Rest of Kazakhstan in the extended EMEP domain + "UZO": 76, # Uzbekistan in the former official EMEP domain + "TMO": 77, # Turkmenistan in the former official EMEP domain + "UZE": 78, # Rest of Uzbekistan in the extended EMEP domain + "TME": 79, # Rest of Turkmenistan in the extended EMEP domain + "CAS": 80, # Caspian Sea + "TJ": 81, # Tajikistan + "ARO": 82, # Aral Lake in the former official EMEP domain + "ARE": 83, # Rest of Aral Lake in the extended EMEP domain + "ASM": 84, # Modified Remaining Asian Areas in the former official EMEP domain + "ASE": 85, # Remaining Asian Areas in the extended EMEP domain + "AOE": 86, # Arctic Ocean in the extended EMEP domain + "RFX": 87, # Extended EMEP External Part of Russian Federation + "ASX": 88, # Extended EMEP External Part of Asia + "PAX": 89, # Extended EMEP External Part of Pacific Ocean + "AOX": 90, # Extended EMEP External Part of Arctic Ocean + "NAX": 91, # Extended EMEP External Part of North Africa + "KZT": 92, # Kazakhstan + "RUE": 93, # Russian Federation in the extended EMEP domain (RU+RFE+RUX) + "UZ": 94, # Uzbekistan + "TM": 95, # Turkmenistan + "AST": 96, # Asian areas in the extended EMEP domain (ASM+ASE+ARO+ARE+CAS) + "FYU": 99, # Former Yugoslavia + "BEF": 301, # Belgium (Flanders) + "BA2": 302, # Baltic Sea EU Cargo o12m + "BA3": 303, # Baltic Sea ROW Cargo o12m + "BA4": 304, # Baltic Sea EU Cargo i12m + "BA5": 305, # Baltic Sea ROW Cargo i12m + "BA6": 306, # Baltic Sea EU Ferry o12m + "BA7": 307, # Baltic Sea ROW Ferry o12m + "BA8": 308, # Baltic Sea EU Ferry i12m + "BA9": 309, # Baltic Sea ROW Ferry i12m + "NO2": 312, # North Sea EU Cargo o12m + "NO3": 313, # North Sea ROW Cargo o12m + "NO4": 314, # North Sea EU Cargo i12m + "NO5": 315, # North Sea ROW Cargo i12m + "NO6": 316, # North Sea EU Ferry o12m + "NO7": 317, # North Sea ROW Ferry o12m + "NO8": 318, # North Sea EU Ferry i12m + "NO9": 319, # North Sea ROW Ferry i12m + "AT2": 322, # Remaining North-East Atlantic Ocean EU Cargo 'o1': 2 # m + "AT3": 323, # Remaining North-East Atlantic Ocean ROW Cargo 'o1': 2 # m + "AT4": 324, # Remaining North-East Atlantic Ocean EU Cargo 'i1': 2 # m + "AT5": 325, # Remaining North-East Atlantic Ocean ROW Cargo 'i1': 2 # m + "AT6": 326, # Remaining North-East Atlantic Ocean EU Ferry 'o1': 2 # m + "AT7": 327, # Remaining North-East Atlantic Ocean ROW Ferry 'o1': 2 # m + "AT8": 328, # Remaining North-East Atlantic Ocean EU Ferry 'i1': 2 # m + "AT9": 329, # Remaining North-East Atlantic Ocean ROW Ferry 'i1': 2 # m + "ME2": 332, # Mediterranean Sea EU Cargo o12m + "ME3": 333, # Mediterranean Sea ROW Cargo o12m + "ME4": 334, # Mediterranean Sea EU Cargo i12m + "ME5": 335, # Mediterranean Sea ROW Cargo i12m + "ME6": 336, # Mediterranean Sea EU Ferry o12m + "ME7": 337, # Mediterranean Sea ROW Ferry o12m + "ME8": 338, # Mediterranean Sea EU Ferry i12m + "ME9": 339, # Mediterranean Sea ROW Ferry i12m + "BL2": 342, # Black Sea EU Cargo o12m + "BL3": 343, # Black Sea ROW Cargo o12m + "BL4": 344, # Black Sea EU Cargo i12m + "BL5": 345, # Black Sea ROW Cargo i12m + "BL6": 346, # Black Sea EU Ferry o12m + "BL7": 347, # Black Sea ROW Ferry o12m + "BL8": 348, # Black Sea EU Ferry i12m + "BL9": 349, # Black Sea ROW Ferry i12m + "GL": 601, # Greenland } diff --git a/grid.py b/grid.py index ba5e8aa..a0747b0 100644 --- a/grid.py +++ b/grid.py @@ -10,15 +10,17 @@ import math from math import pi, sin, cos, asin, atan2 -def read_rotpole(filename): - """ Read rotated pole coordinates from COSMO filename. """ - - with netCDF4.Dataset(filename) as nc: - rotpole = nc.variables['rotated_pole'] - pollon = rotpole.getncattr('grid_north_pole_longitude') - pollat = rotpole.getncattr('grid_north_pole_latitude') - return pollon, pollat +def read_rotpole(filename): + """ Read rotated pole coordinates from COSMO filename. """ + + with netCDF4.Dataset(filename) as nc: + rotpole = nc.variables["rotated_pole"] + pollon = rotpole.getncattr("grid_north_pole_longitude") + pollat = rotpole.getncattr("grid_north_pole_latitude") + + return pollon, pollat + def rotated_grid_transform(position, direction, north_pole): """ @@ -34,12 +36,11 @@ def rotated_grid_transform(position, direction, north_pole): lon = (lon * pi) / 180.0 lat = (lat * pi) / 180.0 - NP_lon = north_pole[0] + NP_lon = north_pole[0] NP_lat = north_pole[1] - theta = 90 - NP_lat # Rotation around y-axis - phi = NP_lon - 180 # Rotation around z-axis - + theta = 90 - NP_lat # Rotation around y-axis + phi = NP_lon - 180 # Rotation around z-axis # Convert degrees to radians phi = (phi * pi) / 180.0 @@ -50,23 +51,37 @@ def rotated_grid_transform(position, direction, north_pole): y = sin(lon) * cos(lat) z = sin(lat) - if direction == 1: # Regular -> Rotated - - x_new = cos(theta) * cos(phi) * x + cos(theta) * sin(phi) * y + sin(theta) * z + if direction == 1: # Regular -> Rotated + + x_new = ( + cos(theta) * cos(phi) * x + + cos(theta) * sin(phi) * y + + sin(theta) * z + ) y_new = -sin(phi) * x + cos(phi) * y - z_new = -sin(theta) * cos(phi) * x - sin(theta) * sin(phi) * y + cos(theta) * z - - elif direction == 2: # Rotated -> Regular - + z_new = ( + -sin(theta) * cos(phi) * x + - sin(theta) * sin(phi) * y + + cos(theta) * z + ) + + elif direction == 2: # Rotated -> Regular + phi = -phi theta = -theta - x_new = cos(theta) * cos(phi) * x + sin(phi) * y + sin(theta) * cos(phi) * z - y_new = -cos(theta) * sin(phi) * x + cos(phi) * y - sin(theta) * sin(phi) * z + x_new = ( + cos(theta) * cos(phi) * x + sin(phi) * y + sin(theta) * cos(phi) * z + ) + y_new = ( + -cos(theta) * sin(phi) * x + + cos(phi) * y + - sin(theta) * sin(phi) * z + ) z_new = -sin(theta) * x + cos(theta) * z - + else: - raise Exception('Invalid direction, value must be either 1 or 2.') + raise Exception("Invalid direction, value must be either 1 or 2.") # Convert cartesian back to spherical coordinates lon_new = atan2(y_new, x_new) diff --git a/hourly_emissions/produce_hourly_emi.py b/hourly_emissions/produce_hourly_emi.py index 58804be..7727191 100644 --- a/hourly_emissions/produce_hourly_emi.py +++ b/hourly_emissions/produce_hourly_emi.py @@ -28,7 +28,7 @@ def daterange(start_date, end_date): Consecutive dates, starting from start_date and ending the day before end_date. """ - for n in range(int((end_date - start_date).days)+1): + for n in range(int((end_date - start_date).days) + 1): yield start_date + datetime.timedelta(n) @@ -105,7 +105,7 @@ def extract_to_grid(grid_to_index, country_vals): return res -class CountryToGridTranslator(): +class CountryToGridTranslator: """A wrapper for extract_to_grid that stores the grid_to_index matrix. The job of this class could be done more elegantly by a lambda, but that doesn't work with Pool (lambdas are not pickleable). @@ -116,6 +116,7 @@ class CountryToGridTranslator(): >>> val_on_emigrid = CountryToGridTranslator(emigrid_to_index) """ + def __init__(self, grid_to_index): self.grid_to_index = grid_to_index @@ -154,42 +155,50 @@ def write_metadata(outfile, org_file, emi_file, ver_file, variables): variables : list(str) List of variable-names to be created """ - rlat = emi_file.dimensions['rlat'].size - rlon = emi_file.dimensions['rlon'].size - level = ver_file.dimensions['level'].size - - outfile.createDimension('time') - outfile.createDimension('rlon', rlon) - outfile.createDimension('rlat', rlat) - outfile.createDimension('bnds', 2) - outfile.createDimension('level', level) - - var_srcfile = [('time', org_file), - ('rotated_pole', org_file), - ('level', org_file), - ('level_bnds', org_file), - ('rlon', emi_file), - ('rlat', emi_file)] + rlat = emi_file.dimensions["rlat"].size + rlon = emi_file.dimensions["rlon"].size + level = ver_file.dimensions["level"].size + + outfile.createDimension("time") + outfile.createDimension("rlon", rlon) + outfile.createDimension("rlat", rlat) + outfile.createDimension("bnds", 2) + outfile.createDimension("level", level) + + var_srcfile = [ + ("time", org_file), + ("rotated_pole", org_file), + ("level", org_file), + ("level_bnds", org_file), + ("rlon", emi_file), + ("rlat", emi_file), + ] for varname, src_file in var_srcfile: - outfile.createVariable(varname=varname, - datatype=src_file[varname].datatype, - dimensions=src_file[varname].dimensions) + outfile.createVariable( + varname=varname, + datatype=src_file[varname].datatype, + dimensions=src_file[varname].dimensions, + ) outfile[varname].setncatts(src_file[varname].__dict__) - outfile['time'][:] = org_file['time'][:] - outfile['level'][:] = ver_file['layer_mid'][:] - outfile['level_bnds'][:] = ( - np.array([ver_file['layer_bot'][:], ver_file['layer_top'][:]])) - outfile['rlat'][:] = emi_file['rlat'][:] - outfile['rlon'][:] = emi_file['rlon'][:] + outfile["time"][:] = org_file["time"][:] + outfile["level"][:] = ver_file["layer_mid"][:] + outfile["level_bnds"][:] = np.array( + [ver_file["layer_bot"][:], ver_file["layer_top"][:]] + ) + outfile["rlat"][:] = emi_file["rlat"][:] + outfile["rlon"][:] = emi_file["rlon"][:] for varname in variables: - outfile.createVariable(varname=varname, - datatype='float32', - dimensions=('time', 'level', 'rlat', 'rlon')) + outfile.createVariable( + varname=varname, + datatype="float32", + dimensions=("time", "level", "rlat", "rlon"), + ) outfile[varname].units = "kg m-2 s-1" + def extract_matrices(infile, var_list, indices, transform=None): """Extract the array specified by indices for each variable in var_list from infile. @@ -306,42 +315,42 @@ def process_day(date, path_template, lists, matrices, datasets): """ print(date.strftime("Processing %x...")) - with Dataset(datasets['emi_path']) as emi,\ - Dataset(datasets['org_path']) as org,\ - Dataset(datasets['ver_path']) as ver: - rlat = emi.dimensions['rlat'].size - rlon = emi.dimensions['rlon'].size - levels = ver.dimensions['level'].size + with Dataset(datasets["emi_path"]) as emi, Dataset( + datasets["org_path"] + ) as org, Dataset(datasets["ver_path"]) as ver: + rlat = emi.dimensions["rlat"].size + rlon = emi.dimensions["rlon"].size + levels = ver.dimensions["level"].size for hour in range(24): day_hour = datetime.datetime.combine(date, datetime.time(hour)) of_path = day_hour.strftime(path_template) with Dataset(of_path, "w") as of: - write_metadata(outfile=of, - org_file=org, - emi_file=emi, - ver_file=ver, - variables=lists['variables']) - - for v, var in enumerate(lists['variables']): + write_metadata( + outfile=of, + org_file=org, + emi_file=emi, + ver_file=ver, + variables=lists["variables"], + ) + + for v, var in enumerate(lists["variables"]): oae_vals = np.zeros((levels, rlat, rlon)) - for cat, tp, vp in zip(lists['cats'][v], - lists['tps'][v], - lists['vps'][v]): - emi_mat = matrices['emi_mats'][cat] - hod_mat = matrices['hod_mats'][hour][tp] - dow_mat = matrices['dow_mats'][tp] - moy_mat = matrices['moy_mats'][tp] - ver_mat = matrices['ver_mats'][vp] + for cat, tp, vp in zip( + lists["cats"][v], lists["tps"][v], lists["vps"][v] + ): + emi_mat = matrices["emi_mats"][cat] + hod_mat = matrices["hod_mats"][hour][tp] + dow_mat = matrices["dow_mats"][tp] + moy_mat = matrices["moy_mats"][tp] + ver_mat = matrices["ver_mats"][vp] # Add dimensions so numpy can broadcast ver_mat = np.reshape(ver_mat, (levels, 1, 1)) - oae_vals += (emi_mat * - hod_mat * - dow_mat * - moy_mat * - ver_mat) + oae_vals += ( + emi_mat * hod_mat * dow_mat * moy_mat * ver_mat + ) # Careful, automatic reshaping! of[var][0, :] = oae_vals @@ -354,11 +363,22 @@ def process_day_one_arg(arg): return process_day(*arg) -def generate_arguments(start_date, end_date, - emi_path, org_path, ver_path, - hod_path, dow_path, moy_path, - var_list, catlist, tplist, vplist, - output_path, output_prefix): +def generate_arguments( + start_date, + end_date, + emi_path, + org_path, + ver_path, + hod_path, + dow_path, + moy_path, + var_list, + catlist, + tplist, + vplist, + output_path, + output_prefix, +): """Prepare the arguments for process_day() (mainly extract the relevant data from netcdf-files) and yield them as tuples. @@ -417,17 +437,16 @@ def generate_arguments(start_date, end_date, # Create grid-country-mapping print("Creating gridpoint -> index mapping...") with Dataset(emi_path) as emi, Dataset(moy_path) as moy: - mapping_result = pool.apply_async(country_id_mapping, - args=(moy['country'][:], - emi['country_ids'][:])) + mapping_result = pool.apply_async( + country_id_mapping, + args=(moy["country"][:], emi["country_ids"][:]), + ) # Time- and (grid-country-mapping)-independent data - args_indep = [(emi_path, - catlist, - np.s_[:]), - (ver_path, - vplist, - np.s_[:])] + args_indep = [ + (emi_path, catlist, np.s_[:]), + (ver_path, vplist, np.s_[:]), + ] print("Extracting average emissions and vertical profiles...") res_indep = pool.starmap(extract_matrices, args_indep) @@ -441,26 +460,21 @@ def generate_arguments(start_date, end_date, print("... finished gridpoint -> index mapping") # Time- and (grid-country-mapping)-dependent data - assert start_date.month == end_date.month, ("Adjust the script to " - "prepare data for multiple" - " months.") + assert start_date.month == end_date.month, ( + "Adjust the script to " "prepare data for multiple" " months." + ) month_id = start_date.month - 1 - args_dep = [(moy_path, - tplist, - np.s_[month_id, :], - val_on_emigrid)] + args_dep = [(moy_path, tplist, np.s_[month_id, :], val_on_emigrid)] for day in range(7): # always extract the whole week - args_dep.append(tuple([dow_path, - tplist, - np.s_[day, :], - val_on_emigrid])) + args_dep.append( + tuple([dow_path, tplist, np.s_[day, :], val_on_emigrid]) + ) for hour in range(24): - args_dep.append(tuple([hod_path, - tplist, - np.s_[hour, :], - val_on_emigrid])) + args_dep.append( + tuple([hod_path, tplist, np.s_[hour, :], val_on_emigrid]) + ) print("Extracting month, day and hour profiles...") # List of dicts (containing 2D arrays) is not ideal for cache @@ -471,22 +485,28 @@ def generate_arguments(start_date, end_date, hod_mats = res_dep[8:] print("... finished temporal profiles") - lists = {'variables': var_list, - 'cats': catlist, - 'tps': tplist, - 'vps': vplist} - datasets = {'org_path': org_path, - 'emi_path': emi_path, - 'ver_path': ver_path} + lists = { + "variables": var_list, + "cats": catlist, + "tps": tplist, + "vps": vplist, + } + datasets = { + "org_path": org_path, + "emi_path": emi_path, + "ver_path": ver_path, + } stop = time.time() print("Finished extracting data in " + str(stop - start)) for day_date in daterange(start_date, end_date): - matrices = {'emi_mats': emi_mats, - 'hod_mats': hod_mats, - 'dow_mats': dow_mats[day_date.weekday()], - 'moy_mats': moy_mats, - 'ver_mats': ver_mats} + matrices = { + "emi_mats": emi_mats, + "hod_mats": hod_mats, + "dow_mats": dow_mats[day_date.weekday()], + "moy_mats": moy_mats, + "ver_mats": ver_mats, + } yield (day_date, path_template, lists, matrices, datasets) @@ -499,144 +519,201 @@ def create_lists(): for (s, nfr) in [(s, nfr) for s in tracers for nfr in categories]: var_list.append(s + "_" + nfr + "_E") - catlist_prelim = ( + catlist_prelim = [ + ["CO2_A_AREA", "CO2_A_POINT"], # for CO2_A_E + ["CO2_B_AREA", "CO2_B_POINT"], # for CO2_B_E + ["CO2_C_AREA"], # for CO2_C_E + ["CO2_F_AREA"], # for CO2_F_E + [ + "CO2_D_AREA", + "CO2_D_POINT", + "CO2_E_AREA", + "CO2_G_AREA", + "CO2_H_AREA", + "CO2_H_POINT", + "CO2_I_AREA", + "CO2_J_AREA", + "CO2_J_POINT", + "CO2_K_AREA", + "CO2_L_AREA", + ], # for CO2_O_E + [ + "CO2_A_AREA", + "CO2_A_POINT", + "CO2_B_AREA", + "CO2_B_POINT", + "CO2_C_AREA", + "CO2_F_AREA", + "CO2_D_AREA", + "CO2_D_POINT", + "CO2_E_AREA", + "CO2_G_AREA", + "CO2_H_AREA", + "CO2_H_POINT", + "CO2_I_AREA", + "CO2_J_AREA", + "CO2_J_POINT", + "CO2_K_AREA", + "CO2_L_AREA", + ], # for CO2_ALL_E + ["CO_A_AREA", "CO_A_POINT"], # for CO_A_E + ["CO_B_AREA", "CO_B_POINT"], # for CO_B_E + ["CO_C_AREA"], # for CO_C_E + ["CO_F_AREA"], # for CO_F_E + [ + "CO_D_AREA", + "CO_D_POINT", + "CO_E_AREA", + "CO_G_AREA", + "CO_H_AREA", + "CO_H_POINT", + "CO_I_AREA", + "CO_J_AREA", + "CO_J_POINT", + "CO_K_AREA", + "CO_L_AREA", + ], # for CO_O_E [ - ["CO2_A_AREA", "CO2_A_POINT"], # for CO2_A_E - ["CO2_B_AREA", "CO2_B_POINT"], # for CO2_B_E - ["CO2_C_AREA"], # for CO2_C_E - ["CO2_F_AREA"], # for CO2_F_E - ["CO2_D_AREA", "CO2_D_POINT", - "CO2_E_AREA", - "CO2_G_AREA", - "CO2_H_AREA", "CO2_H_POINT", - "CO2_I_AREA", - "CO2_J_AREA", "CO2_J_POINT", - "CO2_K_AREA", - "CO2_L_AREA"], # for CO2_O_E - ["CO2_A_AREA", "CO2_A_POINT", - "CO2_B_AREA", "CO2_B_POINT", - "CO2_C_AREA", - "CO2_F_AREA", - "CO2_D_AREA", "CO2_D_POINT", - "CO2_E_AREA", - "CO2_G_AREA", - "CO2_H_AREA", "CO2_H_POINT", - "CO2_I_AREA", - "CO2_J_AREA", "CO2_J_POINT", - "CO2_K_AREA", - "CO2_L_AREA"], # for CO2_ALL_E - ["CO_A_AREA", "CO_A_POINT"], # for CO_A_E - ["CO_B_AREA", "CO_B_POINT"], # for CO_B_E - ["CO_C_AREA"], # for CO_C_E - ["CO_F_AREA"], # for CO_F_E - ["CO_D_AREA", "CO_D_POINT", - "CO_E_AREA", - "CO_G_AREA", - "CO_H_AREA", "CO_H_POINT", - "CO_I_AREA", - "CO_J_AREA", "CO_J_POINT", - "CO_K_AREA", - "CO_L_AREA"], # for CO_O_E - ["CO_A_AREA", "CO_A_POINT", - "CO_B_AREA", "CO_B_POINT", - "CO_C_AREA", - "CO_F_AREA", - "CO_D_AREA", "CO_D_POINT", - "CO_E_AREA", - "CO_G_AREA", - "CO_H_AREA", "CO_H_POINT", - "CO_I_AREA", - "CO_J_AREA", "CO_J_POINT", - "CO_K_AREA", - "CO_L_AREA"], # for CO_ALL_E - ["CH4_A_AREA", "CH4_A_POINT"], # for CH4_A_E - ["CH4_B_AREA", "CH4_B_POINT"], # for CH4_B_E - ["CH4_C_AREA"], # for CH4_C_E - ["CH4_F_AREA"], # for CH4_F_E - ["CH4_D_AREA", "CH4_D_POINT", - "CH4_E_AREA", - "CH4_G_AREA", - "CH4_H_AREA", "CH4_H_POINT", - "CH4_I_AREA", - "CH4_J_AREA", "CH4_J_POINT", - "CO2_K_AREA", - "CO2_L_AREA"], # for CH4_O_E - ["CH4_A_AREA", "CH4_A_POINT", - "CH4_B_AREA", "CH4_B_POINT", - "CH4_C_AREA", - "CH4_F_AREA", - "CH4_D_AREA", "CH4_D_POINT", - "CH4_E_AREA", - "CH4_G_AREA", - "CH4_H_AREA", "CH4_H_POINT", - "CH4_I_AREA", - "CH4_J_AREA", "CH4_J_POINT", - "CO2_K_AREA", - "CO2_L_AREA"] # for CH4_ALL_E - ]) - - tplist_prelim = ( + "CO_A_AREA", + "CO_A_POINT", + "CO_B_AREA", + "CO_B_POINT", + "CO_C_AREA", + "CO_F_AREA", + "CO_D_AREA", + "CO_D_POINT", + "CO_E_AREA", + "CO_G_AREA", + "CO_H_AREA", + "CO_H_POINT", + "CO_I_AREA", + "CO_J_AREA", + "CO_J_POINT", + "CO_K_AREA", + "CO_L_AREA", + ], # for CO_ALL_E + ["CH4_A_AREA", "CH4_A_POINT"], # for CH4_A_E + ["CH4_B_AREA", "CH4_B_POINT"], # for CH4_B_E + ["CH4_C_AREA"], # for CH4_C_E + ["CH4_F_AREA"], # for CH4_F_E [ - ['GNFR_A', 'GNFR_A'], # for s_A_E - ['GNFR_B', 'GNFR_B'], # for s_B_E - ['GNFR_C'], # for s_C_E - ['GNFR_F'], # for s_F_E - ['GNFR_D', 'GNFR_D', - 'GNFR_E', - 'GNFR_G', - 'GNFR_H', 'GNFR_H', - 'GNFR_I', - 'GNFR_J', 'GNFR_J', - 'GNFR_K', - 'GNFR_L'], # for s_O_E - ['GNFR_A', 'GNFR_A', - 'GNFR_B', 'GNFR_B', - 'GNFR_C', - 'GNFR_F', - 'GNFR_D', 'GNFR_D', - 'GNFR_E', - 'GNFR_G', - 'GNFR_H', 'GNFR_H', - 'GNFR_I', - 'GNFR_J', 'GNFR_J', - 'GNFR_K', - 'GNFR_L'] # for s_ALL_E - ]) + "CH4_D_AREA", + "CH4_D_POINT", + "CH4_E_AREA", + "CH4_G_AREA", + "CH4_H_AREA", + "CH4_H_POINT", + "CH4_I_AREA", + "CH4_J_AREA", + "CH4_J_POINT", + "CO2_K_AREA", + "CO2_L_AREA", + ], # for CH4_O_E + [ + "CH4_A_AREA", + "CH4_A_POINT", + "CH4_B_AREA", + "CH4_B_POINT", + "CH4_C_AREA", + "CH4_F_AREA", + "CH4_D_AREA", + "CH4_D_POINT", + "CH4_E_AREA", + "CH4_G_AREA", + "CH4_H_AREA", + "CH4_H_POINT", + "CH4_I_AREA", + "CH4_J_AREA", + "CH4_J_POINT", + "CO2_K_AREA", + "CO2_L_AREA", + ], # for CH4_ALL_E + ] + + tplist_prelim = [ + ["GNFR_A", "GNFR_A"], # for s_A_E + ["GNFR_B", "GNFR_B"], # for s_B_E + ["GNFR_C"], # for s_C_E + ["GNFR_F"], # for s_F_E + [ + "GNFR_D", + "GNFR_D", + "GNFR_E", + "GNFR_G", + "GNFR_H", + "GNFR_H", + "GNFR_I", + "GNFR_J", + "GNFR_J", + "GNFR_K", + "GNFR_L", + ], # for s_O_E + [ + "GNFR_A", + "GNFR_A", + "GNFR_B", + "GNFR_B", + "GNFR_C", + "GNFR_F", + "GNFR_D", + "GNFR_D", + "GNFR_E", + "GNFR_G", + "GNFR_H", + "GNFR_H", + "GNFR_I", + "GNFR_J", + "GNFR_J", + "GNFR_K", + "GNFR_L", + ], # for s_ALL_E + ] """The vertical profile is only applied to area sources. All area sources have emissions at the floor level. As such, their profiles are of the shape [1,0,0,...], like GNFR_L""" - - vplist_prelim = ( + + vplist_prelim = [ + ["GNFR_L", "GNFR_A"], # for s_A_E + ["GNFR_L", "GNFR_B"], # for s_B_E + ["GNFR_L"], # for s_C_E + ["GNFR_L"], # for s_F_E + [ + "GNFR_L", + "GNFR_D", + "GNFR_L", + "GNFR_L", + "GNFR_L", + "GNFR_H", + "GNFR_L", + "GNFR_L", + "GNFR_J", + "GNFR_L", + "GNFR_L", + ], # for s_O_E [ - ['GNFR_L', 'GNFR_A'], # for s_A_E - ['GNFR_L', 'GNFR_B'], # for s_B_E - ['GNFR_L'], # for s_C_E - ['GNFR_L'], # for s_F_E - ['GNFR_L', 'GNFR_D', - 'GNFR_L', - 'GNFR_L', - 'GNFR_L', 'GNFR_H', - 'GNFR_L', - 'GNFR_L', 'GNFR_J', - 'GNFR_L', - 'GNFR_L'], # for s_O_E - ['GNFR_L', 'GNFR_A', - 'GNFR_L', 'GNFR_B', - 'GNFR_L', - 'GNFR_L', - 'GNFR_L', 'GNFR_D', - 'GNFR_L', - 'GNFR_L', - 'GNFR_L', 'GNFR_H', - 'GNFR_L', - 'GNFR_L', 'GNFR_J', - 'GNFR_L', - 'GNFR_L'] # for s_ALL_E - ]) + "GNFR_L", + "GNFR_A", + "GNFR_L", + "GNFR_B", + "GNFR_L", + "GNFR_L", + "GNFR_L", + "GNFR_D", + "GNFR_L", + "GNFR_L", + "GNFR_L", + "GNFR_H", + "GNFR_L", + "GNFR_L", + "GNFR_J", + "GNFR_L", + "GNFR_L", + ], # for s_ALL_E + ] tplist_prelim *= 3 - vplist_prelim *= 3 #tplist_prelim + vplist_prelim *= 3 # tplist_prelim # Make sure catlist, tplist, vplist have the same shape catlist, tplist, vplist = [], [], [] @@ -644,9 +721,9 @@ def create_lists(): subcat = [] subtp = [] subvp = [] - for cat, tp, vp in zip(catlist_prelim[v], - tplist_prelim[v], - vplist_prelim[v]): + for cat, tp, vp in zip( + catlist_prelim[v], tplist_prelim[v], vplist_prelim[v] + ): subcat.append(cat) subtp.append(tp) subvp.append(vp) @@ -657,34 +734,41 @@ def create_lists(): return var_list, catlist, tplist, vplist -def main(path_emi, path_org, output_path, output_name, prof_path, - start_date, end_date): +def main( + path_emi, + path_org, + output_path, + output_name, + prof_path, + start_date, + end_date, +): point1 = time.time() var_list, catlist, tplist, vplist = create_lists() # Prepare arguments - args = generate_arguments(start_date=start_date, - end_date=end_date, - emi_path=path_emi, - org_path=path_org, - ver_path=prof_path + "vertical_profiles.nc", - hod_path=prof_path + "hourofday.nc", - dow_path=prof_path + "dayofweek.nc", - moy_path=prof_path + "monthofyear.nc", - var_list=var_list, - catlist=catlist, - tplist=tplist, - vplist=vplist, - output_path=output_path, - output_prefix=output_name) + args = generate_arguments( + start_date=start_date, + end_date=end_date, + emi_path=path_emi, + org_path=path_org, + ver_path=prof_path + "vertical_profiles.nc", + hod_path=prof_path + "hourofday.nc", + dow_path=prof_path + "dayofweek.nc", + moy_path=prof_path + "monthofyear.nc", + var_list=var_list, + catlist=catlist, + tplist=tplist, + vplist=vplist, + output_path=output_path, + output_prefix=output_name, + ) with Pool(14) as pool: # 2 weeks in parallel # Use imap for smaller memory requirements (arguments are only built # when they are needed) - res = pool.imap(func=process_day_one_arg, - iterable=args, - chunksize=1) + res = pool.imap(func=process_day_one_arg, iterable=args, chunksize=1) # Iterate through all results to prevent the script from returning # before imap is done @@ -695,21 +779,27 @@ def main(path_emi, path_org, output_path, output_name, prof_path, print("Overall runtime: {} s".format(point2 - point1)) -if __name__ == '__main__': +if __name__ == "__main__": # CHE - path_emi = "../testdata/CHE_TNO_offline/emis_2015_Europe.nc" #"./../emis_2015_Europe.nc" - path_org = ("../testdata/hourly_emi_brd/CO2_CO_NOX_Berlin-coarse_2015010110.nc") + path_emi = ( + "../testdata/CHE_TNO_offline/emis_2015_Europe.nc" + ) # "./../emis_2015_Europe.nc" + path_org = ( + "../testdata/hourly_emi_brd/CO2_CO_NOX_Berlin-coarse_2015010110.nc" + ) output_path = "./output_CHE/" output_name = "Europe_CHE_" prof_path = "./input_profiles_CHE/" start_date = datetime.date(2015, 1, 1) - end_date = datetime.date(2015, 1, 9) #included - - main(path_emi=path_emi, - path_org=path_org, - output_path=output_path, - output_name=output_name, - prof_path=prof_path, - start_date=start_date, - end_date=end_date) + end_date = datetime.date(2015, 1, 9) # included + + main( + path_emi=path_emi, + path_org=path_org, + output_path=output_path, + output_name=output_name, + prof_path=prof_path, + start_date=start_date, + end_date=end_date, + ) diff --git a/main_ch.py b/main_ch.py index b12a9b1..a0ae287 100644 --- a/main_ch.py +++ b/main_ch.py @@ -2,30 +2,33 @@ from glob import glob import numpy as np -#species = ['CO2', 'CO', 'CH4'] -species = ['CO2'] -gnfr_cat = [ "A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M" ] +# species = ['CO2', 'CO', 'CH4'] +species = ["CO2"] +gnfr_cat = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M"] + def get_ch_emi(filename): """read in meteotest swiss inventory data output: emi_trans[lon][lat] (emission data per species per category) """ - no_data= -9999 + no_data = -9999 emi = np.loadtxt(filename, skiprows=6) - - emi_new = np.empty((np.shape(emi)[0],np.shape(emi)[1])) # create the same dimention as emi + + emi_new = np.empty( + (np.shape(emi)[0], np.shape(emi)[1]) + ) # create the same dimention as emi for i in range(np.shape(emi)[0]): - for j in range(np.shape(emi)[1]): - emi_new[i,j]=emi[np.shape(emi)[0]-1-i,j] - if emi_new[i,j]==no_data: - emi_new[i,j]=0 - + for j in range(np.shape(emi)[1]): + emi_new[i, j] = emi[np.shape(emi)[0] - 1 - i, j] + if emi_new[i, j] == no_data: + emi_new[i, j] = 0 + emi_trans = np.transpose(emi_new) - - + return emi_trans + def main(cfg_path): """ The main script for processing Swiss inventories. Takes a configuration file as input""" @@ -36,80 +39,87 @@ def main(cfg_path): """Load or compute the country mask""" country_mask = get_country_mask(cfg) country_mask = np.transpose(country_mask) - mask = country_mask != country_codes['CH'] - print('Only use data inside "CH" according to country_mask ' - '(country code %d)' % country_codes['CH']) + mask = country_mask != country_codes["CH"] + print( + 'Only use data inside "CH" according to country_mask ' + "(country code %d)" % country_codes["CH"] + ) """Set names for longitude and latitude""" - lonname = "rlon"; latname="rlat" - if (cfg.pollon==180 or cfg.pollon==0) and cfg.pollat==90: - lonname = "lon"; latname="lat" + lonname = "rlon" + latname = "rlat" + if (cfg.pollon == 180 or cfg.pollon == 0) and cfg.pollat == 90: + lonname = "lon" + latname = "lat" """Load or compute the interpolation map""" - interpolation = get_interpolation(cfg, None, inv_name=cfg.origin, - filename='mapping_' + cfg.origin + '.npy') - + interpolation = get_interpolation( + cfg, + None, + inv_name=cfg.origin, + filename="mapping_" + cfg.origin + ".npy", + ) """Starts writing out the output file""" - output_path = os.path.join(cfg.output_path, - "emis_" + str(cfg.year) + "_" + cfg.gridname + ".nc") - with nc.Dataset(output_path,"w") as out: - prepare_output_file(cfg,out,country_mask) + output_path = os.path.join( + cfg.output_path, "emis_" + str(cfg.year) + "_" + cfg.gridname + ".nc" + ) + with nc.Dataset(output_path, "w") as out: + prepare_output_file(cfg, out, country_mask) """ Swiss inventory specific (works with MeteoTest, MaiolicaCH4, and CarboCountCO2) """ total_flux = {} for var in cfg.species: - total_flux[var] = np.zeros((cfg.ny,cfg.nx)) - - for cat in cfg.ch_cat: - for var in cfg.species: - if cfg.origin == 'meteotest': - constfile = os.path.join(cfg.input_path, - ''.join(['e',cat.lower(),'15_', - var.lower(),'*']) - ) + total_flux[var] = np.zeros((cfg.ny, cfg.nx)) + + for cat in cfg.ch_cat: + for var in cfg.species: + if cfg.origin == "meteotest": + constfile = os.path.join( + cfg.input_path, + "".join(["e", cat.lower(), "15_", var.lower(), "*"]), + ) out_var_name = var + "_" + cat - elif cfg.origin == 'carbocount': - constfile = os.path.join(cfg.input_path, - ''.join([cat.lower(),'10_', - '*_kg.txt']) - ) + elif cfg.origin == "carbocount": + constfile = os.path.join( + cfg.input_path, + "".join([cat.lower(), "10_", "*_kg.txt"]), + ) out_var_name = var + "_" + cfg.mapping[cat] - elif cfg.origin == 'maiolica': - constfile = os.path.join(cfg.input_path, - ''.join(['ch4_',cat.lower(), - '.txt']) - ) + elif cfg.origin == "maiolica": + constfile = os.path.join( + cfg.input_path, "".join(["ch4_", cat.lower(), ".txt"]) + ) out_var_name = var + "_" + cfg.mapping[cat] else: print("Wrong origin in the config file.") raise ValueError - emi= np.zeros((cfg.ch_xn,cfg.ch_yn)) - for filename in sorted(glob(constfile)): + emi = np.zeros((cfg.ch_xn, cfg.ch_yn)) + for filename in sorted(glob(constfile)): print(filename) - emi += get_ch_emi(filename) #(lon,lat) + emi += get_ch_emi(filename) # (lon,lat) start = time.time() - out_var = np.zeros((cfg.ny,cfg.nx)) + out_var = np.zeros((cfg.ny, cfg.nx)) for lon in range(np.shape(emi)[0]): for lat in range(np.shape(emi)[1]): - for (x,y,r) in interpolation[lon,lat]: - out_var[y,x] += emi[lon,lat]*r + for (x, y, r) in interpolation[lon, lat]: + out_var[y, x] += emi[lon, lat] * r end = time.time() - print("it takes ",end-start,"sec") + print("it takes ", end - start, "sec") """calculate the areas (m^^2) of the COSMO grid""" - cosmo_area = 1./gridbox_area(cfg) + cosmo_area = 1.0 / gridbox_area(cfg) """unit conversion""" - if cfg.origin == 'carbocount': + if cfg.origin == "carbocount": """convert unit from kg.year-1.cell-1 to kg.m-2.s-1""" - out_var *= cosmo_area.T*convfac + out_var *= cosmo_area.T * convfac else: """convert unit from g.year-1.cell-1 to kg.m-2.s-1""" - out_var *= cosmo_area.T*convfac/1000 + out_var *= cosmo_area.T * convfac / 1000 """mask values outside CH for consistency""" out_var[mask] = 0 @@ -118,7 +128,7 @@ def main(cfg_path): out_var[out_var < 0] = 0 if not out_var_name in out.variables.keys(): - out.createVariable(out_var_name,float,(latname,lonname)) + out.createVariable(out_var_name, float, (latname, lonname)) if lonname == "rlon" and latname == "rlat": out[out_var_name].grid_mapping = "rotated_pole" out[out_var_name].units = "kg m-2 s-1" @@ -128,27 +138,26 @@ def main(cfg_path): total_flux[var] += out_var - """Calcluate total emission/flux per species""" for s in cfg.species: - out.createVariable(s,float,(latname,lonname)) + out.createVariable(s, float, (latname, lonname)) out[s].units = "kg m-2 s-1" if lonname == "rlon" and latname == "rlat": out[s].grid_mapping = "rotated_pole" out[s][:] = total_flux[s] - + """Create dummy variables for merging inventories""" for s in species: if not s in out.variables.keys(): - out.createVariable(s,float,(latname,lonname)) + out.createVariable(s, float, (latname, lonname)) if lonname == "rlon" and latname == "rlat": out[s].grid_mapping = "rotated_pole" out[s].units = "kg m-2 s-1" out[s][:] = 0 for gnfr in gnfr_cat: - varname = s + '_' + gnfr + varname = s + "_" + gnfr if not varname in out.variables.keys(): - out.createVariable(varname,float,(latname,lonname)) + out.createVariable(varname, float, (latname, lonname)) if lonname == "rlon" and latname == "rlat": out[varname].grid_mapping = "rotated_pole" out[varname].units = "kg m-2 s-1" diff --git a/main_edgar.py b/main_edgar.py index e3ddc19..171607a 100644 --- a/main_edgar.py +++ b/main_edgar.py @@ -1,6 +1,7 @@ from make_online_emissions import * from glob import glob + def get_lat_lon(filename): with open(filename) as f: lat = [] @@ -12,27 +13,33 @@ def get_lat_lon(filename): lon.append(float(data[1])) except ValueError: continue - return lat,lon + return lat, lon + -def get_emi(filename,cfg): - lon_dim,lat_dim,lon_var,lat_var = get_dim_var(None,"edgar",cfg) - emi = np.zeros((lon_dim,lat_dim)) +def get_emi(filename, cfg): + lon_dim, lat_dim, lon_var, lat_var = get_dim_var(None, "edgar", cfg) + emi = np.zeros((lon_dim, lat_dim)) with open(filename) as f: for l in f: data = l.split(";") try: lat = float(data[0]) lon = float(data[1]) - if (lat <=lat_var[-1]) and (lon <= lon_var[-1]) and \ - (lat >=lat_var[0]) and (lon >=lon_var[0]): - emi_lon = round((lon - cfg.edgar_xmin)/cfg.edgar_dx) - emi_lat = round((lat - cfg.edgar_ymin)/cfg.edgar_dy) - - emi[emi_lon,emi_lat] = float(data[2]) + if ( + (lat <= lat_var[-1]) + and (lon <= lon_var[-1]) + and (lat >= lat_var[0]) + and (lon >= lon_var[0]) + ): + emi_lon = round((lon - cfg.edgar_xmin) / cfg.edgar_dx) + emi_lat = round((lat - cfg.edgar_ymin) / cfg.edgar_dy) + + emi[emi_lon, emi_lat] = float(data[2]) except ValueError: continue return emi + def main(cfg_path): """ The main script for processing TNO inventory. Takes a configuration file as input""" @@ -44,52 +51,50 @@ def main(cfg_path): country_mask = get_country_mask(cfg) """Load or compute the interpolation map""" - interpolation = get_interpolation(cfg,None,inv_name="edgar") + interpolation = get_interpolation(cfg, None, inv_name="edgar") """Starts writing out the output file""" - output_path = cfg.output_path+"emis_"+str(cfg.year)+"_"+cfg.gridname+".nc" - with nc.Dataset(output_path,"w") as out: - prepare_output_file(cfg,out,country_mask) - + output_path = ( + cfg.output_path + "emis_" + str(cfg.year) + "_" + cfg.gridname + ".nc" + ) + with nc.Dataset(output_path, "w") as out: + prepare_output_file(cfg, out, country_mask) """ EDGAR specific""" - out_var = np.zeros((cfg.ny,cfg.nx)) #sum of all sources + out_var = np.zeros((cfg.ny, cfg.nx)) # sum of all sources for cat in cfg.edgar_cat: - path = os.path.join(cfg.input_path,cat) - files = glob(path+"/*_2015_*") - if len(files)>1 or len(files)<0: + path = os.path.join(cfg.input_path, cat) + files = glob(path + "/*_2015_*") + if len(files) > 1 or len(files) < 0: print("There are too many or too few files") print(files) else: filename = files[0] print(filename) - + start = time.time() - #lat,lon = get_lat_lon(filename) - emi = get_emi(filename,cfg) #(lon,lat) + # lat,lon = get_lat_lon(filename) + emi = get_emi(filename, cfg) # (lon,lat) for lon in range(emi.shape[0]): for lat in range(emi.shape[1]): - for (x,y,r) in interpolation[lon,lat]: + for (x, y, r) in interpolation[lon, lat]: # EDGAR inventory is in tons per grid cell - out_var[y,x] += emi[lon,lat]*r + out_var[y, x] += emi[lon, lat] * r end = time.time() - print("it takes ",end-start,"sec") + print("it takes ", end - start, "sec") """convert unit from ton.year-1.cell-1 to kg.m-2.s-1""" - + """calculate the areas (m^^2) of the COSMO grid""" - cosmo_area = 1./gridbox_area(cfg) - out_var *= cosmo_area.T*convfac*1000 + cosmo_area = 1.0 / gridbox_area(cfg) + out_var *= cosmo_area.T * convfac * 1000 - out_var_name = cfg.species+"_EDGAR" - out.createVariable(out_var_name,float,("rlat","rlon")) + out_var_name = cfg.species + "_EDGAR" + out.createVariable(out_var_name, float, ("rlat", "rlon")) out[out_var_name].units = "kg m-2 s-1" out[out_var_name][:] = out_var - - if __name__ == "__main__": cfg_name = sys.argv[1] main("./config_" + cfg_name) - diff --git a/main_tno.py b/main_tno.py index 1312f22..56bd569 100644 --- a/main_tno.py +++ b/main_tno.py @@ -2,18 +2,19 @@ import country_code -#species = ['CO2', 'CO', 'CH4'] -species = ['CO2'] -gnfr_cat = [ "A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M" ] +# species = ['CO2', 'CO', 'CH4'] +species = ["CO2"] +gnfr_cat = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M"] exclude_CH = False only_CH = False -def interpolate_tno_to_cosmo_grid(tno,cfg): - return interpolate_to_cosmo_grid(tno,"tno",cfg) +def interpolate_tno_to_cosmo_grid(tno, cfg): + return interpolate_to_cosmo_grid(tno, "tno", cfg) -def interpolate_tno_to_cosmo_point(source,tno,cfg): + +def interpolate_tno_to_cosmo_point(source, tno, cfg): """This function returns the indices of the cosmo grid cell that contains the point source input : - source : The index of the point source from the TNO inventory @@ -21,12 +22,12 @@ def interpolate_tno_to_cosmo_point(source,tno,cfg): - cfg : the configuration file output : - (cosmo_indx,cosmo_indy) : the indices of the cosmo grid cell containing the source -""" +""" lon_source = tno["longitude_source"][source] lat_source = tno["latitude_source"][source] - - return interpolate_to_cosmo_point(lat_source,lon_source,cfg) - + + return interpolate_to_cosmo_point(lat_source, lon_source, cfg) + # transform = ccrs.RotatedPole(pole_longitude=cfg.pollon, pole_latitude=cfg.pollat) # point = transform.transform_point(lon_source,lat_source,ccrs.PlateCarree()) @@ -36,8 +37,7 @@ def interpolate_tno_to_cosmo_point(source,tno,cfg): # return (cosmo_indx,cosmo_indy) - -def var_name(s,snap,cat_kind): +def var_name(s, snap, cat_kind): """Returns the name of a variable for a given species and snap input : - s : species name ("CO2", "CH4" ...) @@ -46,27 +46,28 @@ def var_name(s,snap,cat_kind): output : - returns a string which concatenate the species with the category number """ - out_var_name = s+"_" - if cat_kind=="SNAP": - if snap==70: + out_var_name = s + "_" + if cat_kind == "SNAP": + if snap == 70: out_var_name += "07_" else: - if snap>9: - out_var_name += str(snap)+"_" + if snap > 9: + out_var_name += str(snap) + "_" else: - out_var_name += "0"+str(snap)+"_" - elif cat_kind=="NFR": + out_var_name += "0" + str(snap) + "_" + elif cat_kind == "NFR": out_var_name += snap - elif cat_kind=="NFR_BAU": - out_var_name += snap+"_BAU_" - elif cat_kind=="NFR_CC": - out_var_name += snap+"_CC_" + elif cat_kind == "NFR_BAU": + out_var_name += snap + "_BAU_" + elif cat_kind == "NFR_CC": + out_var_name += snap + "_CC_" else: print("Wrong cat_kind in the config file. Must be SNAP or NFR") raise ValueError return out_var_name -def tno_file(species,cfg): + +def tno_file(species, cfg): """Returns the path to the TNO file for a given species. For some inputs, the CO2 is in a different file than the other species for instance. There are two parameters in the config file to distinguish them. @@ -74,7 +75,7 @@ def tno_file(species,cfg): """ # get the TNO inventory - if species=="CO2": + if species == "CO2": tno_path = cfg.tnoCamsPath else: tno_path = cfg.tnoMACCIIIPath @@ -85,16 +86,15 @@ def tno_file(species,cfg): first_year = 2000 last_year = 2011 - if cfg.year>last_year: + if cfg.year > last_year: base_year = str(last_year) - elif cfg.year=0 and indx=0 and indy= 0 + and indx < cfg.nx + and indy >= 0 + and indy < cfg.ny + ): + out_var_point[indy, indx] += var[i] end = time.time() - print("it takes ",end-start,"sec") - ## TO DO : + print("it takes ", end - start, "sec") + ## TO DO : ## - Add the factor from 2011 to 2015 - """convert unit from kg.year-1.cell-1 to kg.m-2.s-1""" """calculate the areas (m^^2) of the COSMO grid""" - cosmo_area = 1./gridbox_area(cfg) - out_var_point*= cosmo_area.T*convfac - out_var_area *= cosmo_area.T*convfac - out_var_name = var_name(s,snap,cfg.cat_kind) - - for (t,sel,out_var) in zip(["AREA","POINT"], - [selection_snap_area,selection_snap_point], - [out_var_area,out_var_point]): + cosmo_area = 1.0 / gridbox_area(cfg) + out_var_point *= cosmo_area.T * convfac + out_var_area *= cosmo_area.T * convfac + out_var_name = var_name(s, snap, cfg.cat_kind) + + for (t, sel, out_var) in zip( + ["AREA", "POINT"], + [selection_snap_area, selection_snap_point], + [out_var_area, out_var_point], + ): if sel.any(): if exclude_CH or only_CH: out_var[mask] = 0 if not out_var_name in out.variables.keys(): - out.createVariable(out_var_name,float, - (latname,lonname)) + out.createVariable( + out_var_name, float, (latname, lonname) + ) out[out_var_name].units = "kg m-2 s-1" if lonname == "rlon" and latname == "rlat": - out[out_var_name].grid_mapping = "rotated_pole" + out[ + out_var_name + ].grid_mapping = "rotated_pole" out[out_var_name][:] = out_var else: out[out_var_name][:] += out_var total_flux[s] += out_var - """Calcluate total emission/flux per species""" for s in species_list: - out.createVariable(s,float,(latname,lonname)) + out.createVariable(s, float, (latname, lonname)) out[s].units = "kg m-2 s-1" if lonname == "rlon" and latname == "rlat": out[s].grid_mapping = "rotated_pole" @@ -249,22 +280,21 @@ def main(cfg_path): """Create dummy variables for merging inventories""" for s in species: if not s in out.variables.keys(): - out.createVariable(s,float,(latname,lonname)) + out.createVariable(s, float, (latname, lonname)) if lonname == "rlon" and latname == "rlat": out[s].grid_mapping = "rotated_pole" out[s].units = "kg m-2 s-1" out[s][:] = 0 for gnfr in gnfr_cat: - varname = s + '_' + gnfr + varname = s + "_" + gnfr if not varname in out.variables.keys(): - out.createVariable(varname,float,(latname,lonname)) + out.createVariable(varname, float, (latname, lonname)) if lonname == "rlon" and latname == "rlat": out[varname].grid_mapping = "rotated_pole" out[varname].units = "kg m-2 s-1" out[varname][:] = 0 - + if __name__ == "__main__": cfg_name = sys.argv[1] main("./config_" + cfg_name) - diff --git a/main_tno_example.py b/main_tno_example.py index 3ce4965..ac55df8 100644 --- a/main_tno_example.py +++ b/main_tno_example.py @@ -1,7 +1,7 @@ from make_online_emissions import * -def var_name(s,cat,cat_kind): +def var_name(s, cat, cat_kind): """Returns the name of a variable for a given species and cat input : - s : species name ("CO2", "CH4" ...) @@ -10,14 +10,14 @@ def var_name(s,cat,cat_kind): output : - returns a string which concatenates the species with the category number """ - out_var_name = s+"_" - if cat_kind=="SNAP": - if cat>9: - out_var_name += str(cat)+"_" + out_var_name = s + "_" + if cat_kind == "SNAP": + if cat > 9: + out_var_name += str(cat) + "_" else: - out_var_name += "0"+str(cat)+"_" - elif cat_kind=="NFR": - out_var_name += cat +"_" + out_var_name += "0" + str(cat) + "_" + elif cat_kind == "NFR": + out_var_name += cat + "_" else: print("Wrong cat_kind in the config file. Must be SNAP or NFR") raise ValueError @@ -25,104 +25,128 @@ def var_name(s,cat,cat_kind): return out_var_name - - def main(cfg_path): """ The main script for processing TNO inventory. Takes a configuration file as input""" """Load the configuration file""" cfg = load_cfg(cfg_path) - + """Load or compute the country mask""" country_mask = get_country_mask(cfg) """Set names for longitude and latitude""" - if (cfg.pollon==180 or cfg.pollon==0) and cfg.pollat==90: - lonname = "lon"; latname="lat" - print('Non-rotated grid: pollon = %f, pollat = %f' % (cfg.pollon, cfg.pollat)) + if (cfg.pollon == 180 or cfg.pollon == 0) and cfg.pollat == 90: + lonname = "lon" + latname = "lat" + print( + "Non-rotated grid: pollon = %f, pollat = %f" + % (cfg.pollon, cfg.pollat) + ) else: - lonname = "rlon"; latname="rlat" - print('Rotated grid: pollon = %f, pollat = %f' % (cfg.pollon, cfg.pollat)) - + lonname = "rlon" + latname = "rlat" + print( + "Rotated grid: pollon = %f, pollat = %f" % (cfg.pollon, cfg.pollat) + ) + """Starts writing out the output file""" output_path = cfg.output_path - output_file = os.path.join(cfg.output_path,cfg.output_name) - with nc.Dataset(output_file,"w") as out: - prepare_output_file(cfg,out,country_mask) + output_file = os.path.join(cfg.output_path, cfg.output_name) + with nc.Dataset(output_file, "w") as out: + prepare_output_file(cfg, out, country_mask) - """Load or compute the interpolation maps""" + """Load or compute the interpolation maps""" with nc.Dataset(cfg.tnofile) as tno: - interpolation = get_interpolation(cfg,tno) + interpolation = get_interpolation(cfg, tno) - """From here onward, quite specific for TNO""" + """From here onward, quite specific for TNO""" """Mask corresponding to the area/point sources""" - selection_area = tno["source_type_index"][:]==1 - selection_point = tno["source_type_index"][:]==2 + selection_area = tno["source_type_index"][:] == 1 + selection_point = tno["source_type_index"][:] == 2 """Area of the COSMO grid cells""" - cosmo_area = 1./gridbox_area(cfg) + cosmo_area = 1.0 / gridbox_area(cfg) for cat in cfg.output_cat: """In emission_category_index, we have the index of the category, starting with 1.""" """mask corresponding to the given category""" selection_cat = np.array( - [tno["emission_category_index"][:] == cfg.tno_cat.index(cat)+1]) + [ + tno["emission_category_index"][:] + == cfg.tno_cat.index(cat) + 1 + ] + ) """mask corresponding to the given category for area/point""" - selection_cat_area = np.array([selection_cat.any(0),selection_area]).all(0) - selection_cat_point = np.array([selection_cat.any(0),selection_point]).all(0) + selection_cat_area = np.array( + [selection_cat.any(0), selection_area] + ).all(0) + selection_cat_point = np.array( + [selection_cat.any(0), selection_point] + ).all(0) species_list = cfg.species for s in species_list: - print("Species",s,"Category",cat) - out_var_area = np.zeros((cfg.ny,cfg.nx)) - out_var_point = np.zeros((cfg.ny,cfg.nx)) + print("Species", s, "Category", cat) + out_var_area = np.zeros((cfg.ny, cfg.nx)) + out_var_point = np.zeros((cfg.ny, cfg.nx)) var = tno[s.lower()][:] start = time.time() - for (i,source) in enumerate(var): + for (i, source) in enumerate(var): if selection_cat_area[i]: - lon_ind = tno["longitude_index"][i]-1 - lat_ind = tno["latitude_index"][i]-1 - for (x,y,r) in interpolation[lon_ind,lat_ind]: - out_var_area[y,x]+=var[i]*r + lon_ind = tno["longitude_index"][i] - 1 + lat_ind = tno["latitude_index"][i] - 1 + for (x, y, r) in interpolation[lon_ind, lat_ind]: + out_var_area[y, x] += var[i] * r if selection_cat_point[i]: - (indx,indy) = interpolate_to_cosmo_point( + (indx, indy) = interpolate_to_cosmo_point( tno["latitude_source"][i], - tno["longitude_source"][i], - cfg) - if indx>=0 and indx=0 and indy= 0 + and indx < cfg.nx + and indy >= 0 + and indy < cfg.ny + ): + out_var_point[indy, indx] += var[i] end = time.time() - print("it takes ",end-start,"sec") + print("it takes ", end - start, "sec") """convert unit from kg.year-1.cell-1 to kg.m-2.s-1""" """calculate the areas (m^^2) of the COSMO grid""" - out_var_point*= cosmo_area.T*convfac - out_var_area *= cosmo_area.T*convfac - - out_var_name = var_name(s,cat,cfg.cat_kind) - for (t,sel,out_var) in zip(["AREA","POINT"], - [selection_cat_area,selection_cat_point], - [out_var_area,out_var_point]): + out_var_point *= cosmo_area.T * convfac + out_var_area *= cosmo_area.T * convfac + + out_var_name = var_name(s, cat, cfg.cat_kind) + for (t, sel, out_var) in zip( + ["AREA", "POINT"], + [selection_cat_area, selection_cat_point], + [out_var_area, out_var_point], + ): if sel.any(): - out.createVariable(out_var_name+t,float,(latname,lonname)) - out[out_var_name+t].units = "kg m-2 s-1" + out.createVariable( + out_var_name + t, float, (latname, lonname) + ) + out[out_var_name + t].units = "kg m-2 s-1" if lonname == "rlon" and latname == "rlat": - out[out_var_name+t].grid_mapping = "rotated_pole" - out[out_var_name+t][:] = out_var + out[ + out_var_name + t + ].grid_mapping = "rotated_pole" + out[out_var_name + t][:] = out_var + - if __name__ == "__main__": try: config = sys.argv[1] except IndexError: config = "./config_tno" main(config) - diff --git a/main_vprm.py b/main_vprm.py index 9684423..0f00fe9 100644 --- a/main_vprm.py +++ b/main_vprm.py @@ -8,87 +8,105 @@ import os import importlib import config_VPRM_5km as cfg -from datetime import datetime,timedelta +from datetime import datetime, timedelta from multiprocessing import Pool - -def doit(i,interp_tot,cosmo_area): - inidate = datetime(2015,1,1)+timedelta(days=i) - print("taking care of ",inidate.strftime("%Y%m%d")) - for s,sname in zip(["gpp_","ra_"],["CO2_GPP_F","CO2_RA_F"]): +def doit(i, interp_tot, cosmo_area): + inidate = datetime(2015, 1, 1) + timedelta(days=i) + print("taking care of ", inidate.strftime("%Y%m%d")) + for s, sname in zip(["gpp_", "ra_"], ["CO2_GPP_F", "CO2_RA_F"]): for h in range(24): time_ = inidate + timedelta(hours=h) time_str = time_.strftime("%Y%m%d%H") - print(" ",s,time_str) + print(" ", s, time_str) + + # output_path = "/scratch/snx3000/haussaij/VPRM/int2lm/"+s+"_"+time_str+".nc" + filename_out = s + "_" + time_str + ".nc" - #output_path = "/scratch/snx3000/haussaij/VPRM/int2lm/"+s+"_"+time_str+".nc" - filename_out = s+"_"+time_str+".nc" - - output_path = "/scratch/snx3000/haussaij/VPRM/output/"+filename_out + output_path = ( + "/scratch/snx3000/haussaij/VPRM/output/" + filename_out + ) if os.path.exists(output_path): - print('exists:',output_path) + print("exists:", output_path) continue - data_path = "/scratch/snx3000/haussaij/VPRM/input/"+ inidate.strftime("%Y%m%d")+"/"+s+"_"+time_str+".nc" + data_path = ( + "/scratch/snx3000/haussaij/VPRM/input/" + + inidate.strftime("%Y%m%d") + + "/" + + s + + "_" + + time_str + + ".nc" + ) data = Dataset(data_path)[sname][:] - out_var_area = np.zeros((cfg.ny,cfg.nx)) + out_var_area = np.zeros((cfg.ny, cfg.nx)) for lon in range(data.shape[1]): for lat in range(data.shape[0]): - for (x,y,r) in interp_tot[lat,lon]: - #vprm_dat is in umol/m^2/s. Each cell is exactly 1km^2 - out_var_area[y,x] += data[lat,lon]*r*cfg.tno_dx*cfg.tno_dy #now out_var_area is in umol.cell-1.s-1 - - with Dataset(output_path,"w") as outf: - prepare_output_file(cfg,outf) + for (x, y, r) in interp_tot[lat, lon]: + # vprm_dat is in umol/m^2/s. Each cell is exactly 1km^2 + out_var_area[y, x] += ( + data[lat, lon] * r * cfg.tno_dx * cfg.tno_dy + ) # now out_var_area is in umol.cell-1.s-1 + + with Dataset(output_path, "w") as outf: + prepare_output_file(cfg, outf) # Add the time variable - #example_file = "/store/empa/em05/haussaij/CHE/input_che/vprm_old/vprm_smartcarb/processed/"+filename_out - #with Dataset(example_file) as ex : - - outf.createDimension("time") - outf.createVariable(varname="time", - datatype='float64',#ex["time"].datatype, - dimensions=())#ex["time"].dimensions) - outf['time'].units = 'seconds since 2015-01-01 00:00' - outf['time'].calendar = 'proleptic_gregorian' - #outf["time"].setncatts(ex["time"].__dict__) + # example_file = "/store/empa/em05/haussaij/CHE/input_che/vprm_old/vprm_smartcarb/processed/"+filename_out + # with Dataset(example_file) as ex : + outf.createDimension("time") + outf.createVariable( + varname="time", + datatype="float64", # ex["time"].datatype, + dimensions=(), + ) # ex["time"].dimensions) + outf["time"].units = "seconds since 2015-01-01 00:00" + outf["time"].calendar = "proleptic_gregorian" + # outf["time"].setncatts(ex["time"].__dict__) """convert unit from umol.cell-1.s-1 to kg.m-2.s-1""" """calculate the areas (m^^2) of the COSMO grid""" m_co2 = 44.01 - out_var_area *= cosmo_area.T*1e-9*m_co2 + out_var_area *= cosmo_area.T * 1e-9 * m_co2 out_var_name = sname - outf.createVariable(out_var_name,float,("rlat","rlon")) + outf.createVariable(out_var_name, float, ("rlat", "rlon")) outf[out_var_name].units = "kg m-2 s-1" outf[out_var_name][:] = out_var_area - def get_interp_tot(cfg): - for n,i in enumerate([6,2,4,1,3,5]):#[1,3,5,6,2,4]): # Weird order, following the README of Julia - input_path = cfg.vprm_path %str(i) + for n, i in enumerate( + [6, 2, 4, 1, 3, 5] + ): # [1,3,5,6,2,4]): # Weird order, following the README of Julia + input_path = cfg.vprm_path % str(i) with Dataset(input_path) as inf: - interpolation = get_interpolation(cfg,inf,filename="EU"+str(i)+"_mapping.npy",inv_name="vprm") - - if n==0: - x,y = inf["lat"][:].shape - interp_tot = np.empty((x*2,y*3),dtype=object) - - if n<3: - interp_tot[:x,n*y:(n+1)*y] = interpolation.T + interpolation = get_interpolation( + cfg, + inf, + filename="EU" + str(i) + "_mapping.npy", + inv_name="vprm", + ) + + if n == 0: + x, y = inf["lat"][:].shape + interp_tot = np.empty((x * 2, y * 3), dtype=object) + + if n < 3: + interp_tot[:x, n * y : (n + 1) * y] = interpolation.T else: - interp_tot[x:,(n-3)*y:(n-2)*y] = interpolation.T + interp_tot[x:, (n - 3) * y : (n - 2) * y] = interpolation.T return interp_tot + def get_interp_5km(cfg): input_path = cfg.vprm_path with Dataset(input_path) as inf: - interpolation = get_interpolation(cfg,inf,inv_name="vprm") - - return interpolation.T + interpolation = get_interpolation(cfg, inf, inv_name="vprm") + return interpolation.T def main(cfg_path): @@ -101,13 +119,14 @@ def main(cfg_path): else: interp_tot = get_interp_tot(cfg) - cosmo_area = 1./gridbox_area(cfg) + cosmo_area = 1.0 / gridbox_area(cfg) with Pool(36) as pool: - pool.starmap(doit,[(i,interp_tot,cosmo_area) for i in range(10)])#range(1,10)]) - + pool.starmap( + doit, [(i, interp_tot, cosmo_area) for i in range(10)] + ) # range(1,10)]) + -if __name__=="__main__": +if __name__ == "__main__": cfg_name = sys.argv[1] main("./config_" + cfg_name) - diff --git a/make_online_emissions.py b/make_online_emissions.py index e388b9a..555f9c2 100644 --- a/make_online_emissions.py +++ b/make_online_emissions.py @@ -17,38 +17,49 @@ # constants to convert from kg/yr to kg/s/m2 day_per_yr = 365.25 -hr_per_yr = day_per_yr * 24. +hr_per_yr = day_per_yr * 24.0 sec_per_day = 86400 -convfac = 1./(day_per_yr*sec_per_day) +convfac = 1.0 / (day_per_yr * sec_per_day) npool = 18 + def load_cfg(cfg_path): """Load config file""" try: sys.path.append(os.path.dirname(os.path.realpath(cfg_path))) cfg = importlib.import_module(os.path.basename(cfg_path)) except IndexError: - print('ERROR: no config file provided!') + print("ERROR: no config file provided!") sys.exit(1) except ImportError: - print('ERROR: failed to import config module "%s"!' % os.path.basename(cfg_path)) + print( + 'ERROR: failed to import config module "%s"!' + % os.path.basename(cfg_path) + ) sys.exit(1) return cfg + def gridbox_area(cfg): """Calculate 2D array of the areas (m^^2) of the output COSMO grid""" - radius=6375000. #the earth radius in meters - deg2rad=np.pi/180. - dlat = cfg.dy*deg2rad - dlon = cfg.dx*deg2rad + radius = 6375000.0 # the earth radius in meters + deg2rad = np.pi / 180.0 + dlat = cfg.dy * deg2rad + dlon = cfg.dx * deg2rad """Box area at equator""" - dd=2.*pow(radius,2)*dlon*np.sin(0.5*dlat) - areas = np.array([[dd*np.cos(deg2rad*cfg.ymin+j*dlat) for j in range(cfg.ny)] for foo in range(cfg.nx)]) - return areas + dd = 2.0 * pow(radius, 2) * dlon * np.sin(0.5 * dlat) + areas = np.array( + [ + [dd * np.cos(deg2rad * cfg.ymin + j * dlat) for j in range(cfg.ny)] + for foo in range(cfg.nx) + ] + ) + return areas + -def prepare_output_file(cfg,out,country_mask=[]): +def prepare_output_file(cfg, out, country_mask=[]): """Starts writing out the output file : - Dimensions and variables for longitude and latitude - Rotated pole @@ -59,58 +70,59 @@ def prepare_output_file(cfg,out,country_mask=[]): - country_mask : the country_mask that has already been calculated prior outputs: None - """ + """ """Create the dimensions and the rotated pole""" - if (cfg.pollon==180 or cfg.pollon==0) and cfg.pollat==90: - lonname = "lon"; latname="lat" + if (cfg.pollon == 180 or cfg.pollon == 0) and cfg.pollat == 90: + lonname = "lon" + latname = "lat" rotated = False else: - lonname = "rlon"; latname="rlat" + lonname = "rlon" + latname = "rlat" rotated = True - out.createVariable("rotated_pole",str) + out.createVariable("rotated_pole", str) out["rotated_pole"].grid_mapping_name = "rotated_latitude_longitude" out["rotated_pole"].grid_north_pole_latitude = cfg.pollat out["rotated_pole"].grid_north_pole_longitude = cfg.pollon - out["rotated_pole"].north_pole_grid_longitude = 0. + out["rotated_pole"].north_pole_grid_longitude = 0.0 + + out.createDimension(lonname, cfg.nx) + out.createDimension(latname, cfg.ny) - out.createDimension(lonname,cfg.nx) - out.createDimension(latname,cfg.ny) - """Create the variable associated to the dimensions""" - out.createVariable(lonname,"float32",lonname) + out.createVariable(lonname, "float32", lonname) out[lonname].axis = "X" out[lonname].units = "degrees_east" out[lonname].standard_name = "longitude" - lon_range = np.arange(cfg.xmin,cfg.xmin+cfg.dx*cfg.nx,cfg.dx) + lon_range = np.arange(cfg.xmin, cfg.xmin + cfg.dx * cfg.nx, cfg.dx) if len(lon_range) == cfg.nx + 1: lon_range = lon_range[:-1] out[lonname][:] = lon_range - out.createVariable(latname,"float32",latname) + out.createVariable(latname, "float32", latname) out[latname].axis = "Y" out[latname].units = "degrees_north" out[latname].standard_name = "latitude" - lat_range = np.arange(cfg.ymin,cfg.ymin+cfg.dy*cfg.ny,cfg.dy) + lat_range = np.arange(cfg.ymin, cfg.ymin + cfg.dy * cfg.ny, cfg.dy) if len(lat_range) == cfg.ny + 1: lat_range = lat_range[:-1] - out[latname][:] = lat_range + out[latname][:] = lat_range """Create the variable associated with the country_mask""" if len(country_mask): mask_name = "country_ids" - out.createVariable(mask_name,"short",(latname,lonname)) + out.createVariable(mask_name, "short", (latname, lonname)) if rotated: out[mask_name].grid_mapping = "rotated_pole" out[mask_name].long_name = "EMEP_country_code" out[mask_name][:] = country_mask - ################################## ## Regarding the country mask ## ################################## -def check_country(country,points): +def check_country(country, points): """For a given country, return if the grid cell defined by points is within the country input : - country @@ -118,15 +130,18 @@ def check_country(country,points): output : - True if the grid cell is within the country """ - bounds = country.bounds #(minx, miny, maxx, maxy) - if ((bounds[0]>max([k[0] for k in points])) or - (bounds[2]max([k[1] for k in points])) or - (bounds[3] max([k[0] for k in points])) + or (bounds[2] < min([k[0] for k in points])) + or (bounds[1] > max([k[1] for k in points])) + or (bounds[3] < min([k[1] for k in points])) + ): return False else: return True + def compute_country_mask(cfg): """Returns the name of the country for each cosmo grid cell. If for a given grid cell, no country is found (Ocean for example), the country code 0 is assigned. @@ -140,8 +155,9 @@ def compute_country_mask(cfg): print("Creating the country mask") natural_earth = True if natural_earth: - shpfilename = shpreader.natural_earth(resolution='110m', - category='cultural', name='admin_0_countries') + shpfilename = shpreader.natural_earth( + resolution="110m", category="cultural", name="admin_0_countries" + ) iso3 = "ADM0_A3" else: shpfilename = "/usr/local/exelis/idl82/resource/maps/shape/cntry08.shp" @@ -149,11 +165,10 @@ def compute_country_mask(cfg): reader = shpreader.Reader(shpfilename) - country_mask = np.empty((cfg.nx,cfg.ny)) - + country_mask = np.empty((cfg.nx, cfg.ny)) - cosmo_xlocs = np.arange(cfg.xmin, cfg.xmin+cfg.dx*cfg.nx, cfg.dx) - cosmo_ylocs = np.arange(cfg.ymin, cfg.ymin+cfg.dy*cfg.ny, cfg.dy) + cosmo_xlocs = np.arange(cfg.xmin, cfg.xmin + cfg.dx * cfg.nx, cfg.dx) + cosmo_ylocs = np.arange(cfg.ymin, cfg.ymin + cfg.dy * cfg.ny, cfg.dy) """ Be careful with numpy.arange(). Floating point numbers are not exactly represented. Thus, the length of the generated list could have one entry @@ -165,24 +180,30 @@ def compute_country_mask(cfg): if len(cosmo_ylocs) == cfg.ny + 1: cosmo_ylocs = cosmo_ylocs[:-1] - transform = ccrs.RotatedPole(pole_longitude=cfg.pollon, pole_latitude=cfg.pollat) - incr=0 - no_country_code=[] + transform = ccrs.RotatedPole( + pole_longitude=cfg.pollon, pole_latitude=cfg.pollat + ) + incr = 0 + no_country_code = [] european = [] non_euro = [] for country in reader.records(): - if country.attributes["CONTINENT"]=="Europe": + if country.attributes["CONTINENT"] == "Europe": european.append(country) else: non_euro.append(country) - for (a,x) in enumerate(cosmo_xlocs): - for (b,y) in enumerate(cosmo_ylocs): + for (a, x) in enumerate(cosmo_xlocs): + for (b, y) in enumerate(cosmo_ylocs): """Progress bar""" - incr+=1 - sys.stdout.write('\r') - sys.stdout.write(" {:.1f}%".format((100/((cfg.nx*cfg.ny)-1)*((a*cfg.ny)+b)))) + incr += 1 + sys.stdout.write("\r") + sys.stdout.write( + " {:.1f}%".format( + (100 / ((cfg.nx * cfg.ny) - 1) * ((a * cfg.ny) + b)) + ) + ) sys.stdout.flush() mask = [] @@ -192,79 +213,94 @@ def compute_country_mask(cfg): # cosmo_cell_x = [x+cfg.dx,x+cfg.dx,x,x] # cosmo_cell_y = [y+cfg.dy,y,y,y+cfg.dy] # Or the center of the cell - cosmo_cell_x = np.array([x+cfg.dx/2,x+cfg.dx/2,x-cfg.dx/2,x-cfg.dx/2]) - cosmo_cell_y = np.array([y+cfg.dy/2,y-cfg.dy/2,y-cfg.dy/2,y+cfg.dy/2]) - - points = ccrs.PlateCarree().transform_points(transform,cosmo_cell_x,cosmo_cell_y) + cosmo_cell_x = np.array( + [x + cfg.dx / 2, x + cfg.dx / 2, x - cfg.dx / 2, x - cfg.dx / 2] + ) + cosmo_cell_y = np.array( + [y + cfg.dy / 2, y - cfg.dy / 2, y - cfg.dy / 2, y + cfg.dy / 2] + ) + + points = ccrs.PlateCarree().transform_points( + transform, cosmo_cell_x, cosmo_cell_y + ) polygon_cosmo = Polygon(points) """To be faster, only check european countries at first""" - for country in european:#reader.records(): - if check_country(country,points): + for country in european: # reader.records(): + if check_country(country, points): # if x+cfg.dxbounds[2] or y>bounds[3]: # continue if polygon_cosmo.intersects(country.geometry): mask.append(country.attributes[iso3]) """If not found among the european countries, check elsewhere""" - if len(mask)==0: - for country in non_euro:#reader.records(): - if check_country(country,points): + if len(mask) == 0: + for country in non_euro: # reader.records(): + if check_country(country, points): # if x+cfg.dxbounds[2] or y>bounds[3]: # continue if polygon_cosmo.intersects(country.geometry): mask.append(country.attributes[iso3]) - """If more than one country, assign the one which has the greatest area""" - if len(mask)>1: + if len(mask) > 1: area = 0 - for country in [rec for rec in reader.records() if rec.attributes[iso3] in mask]: + for country in [ + rec + for rec in reader.records() + if rec.attributes[iso3] in mask + ]: new_area = polygon_cosmo.intersection(country.geometry).area - if area 0: - print("The following countries were found, but didn't have a corresponding code") + end = time.time() + print("it took", end - start, "seconds") + if len(no_country_code) > 0: + print( + "The following countries were found, but didn't have a corresponding code" + ) print(set(no_country_code)) - np.save(os.path.join(cfg.output_path,"country_mask.npy"),country_mask) + np.save(os.path.join(cfg.output_path, "country_mask.npy"), country_mask) return country_mask def get_country_mask(cfg): """Calculate the country mask""" add_country_mask = True - cmask_path = os.path.join(cfg.output_path,"country_mask.npy") - cmask_path_nc = os.path.join(cfg.output_path,"country_mask.nc") + cmask_path = os.path.join(cfg.output_path, "country_mask.npy") + cmask_path_nc = os.path.join(cfg.output_path, "country_mask.nc") if os.path.isfile(cmask_path_nc): - print("Do you want to use the country mask found in %s ?" % cmask_path_nc) + print( + "Do you want to use the country mask found in %s ?" % cmask_path_nc + ) s = input("[y]/n \n") - if s=='y' or s=='': + if s == "y" or s == "": with nc.Dataset(cmask_path_nc, "r") as inf: - return inf.variables['country_mask'][:].T + return inf.variables["country_mask"][:].T elif os.path.isfile(cmask_path): - print("Do you wanna overwite the country mask found in %s ?" % cmask_path) + print( + "Do you wanna overwite the country mask found in %s ?" % cmask_path + ) s = input("y/[n] \n") - add_country_mask = (s=="y") + add_country_mask = s == "y" if add_country_mask: country_mask = compute_country_mask(cfg) @@ -273,10 +309,13 @@ def get_country_mask(cfg): return country_mask + ################################### ## Regarding the point sources ## ################################### -def interpolate_to_cosmo_point(lat_source,lon_source,cfg,proj=ccrs.PlateCarree()): +def interpolate_to_cosmo_point( + lat_source, lon_source, cfg, proj=ccrs.PlateCarree() +): """This function returns the indices of the cosmo grid cell that contains the point source input : - lat_source : The latitude of the point source @@ -285,22 +324,23 @@ def interpolate_to_cosmo_point(lat_source,lon_source,cfg,proj=ccrs.PlateCarree() - proj : The cartopy projection of the lat/lon of the point source output : - (cosmo_indx,cosmo_indy) : the indices of the cosmo grid cell containing the source -""" - - transform = ccrs.RotatedPole(pole_longitude=cfg.pollon, pole_latitude=cfg.pollat) - point = transform.transform_point(lon_source,lat_source,proj) +""" - cosmo_indx = int(np.floor((point[0]-cfg.xmin)/cfg.dx)) - cosmo_indy = int(np.floor((point[1]-cfg.ymin)/cfg.dy)) + transform = ccrs.RotatedPole( + pole_longitude=cfg.pollon, pole_latitude=cfg.pollat + ) + point = transform.transform_point(lon_source, lat_source, proj) - return (cosmo_indx,cosmo_indy) + cosmo_indx = int(np.floor((point[0] - cfg.xmin) / cfg.dx)) + cosmo_indy = int(np.floor((point[1] - cfg.ymin) / cfg.dy)) + return (cosmo_indx, cosmo_indy) ################################## ## Regarding the area sources ## ################################## -def interpolate_single_cell(cfg,points): +def interpolate_single_cell(cfg, points): """This function determines which COSMO cell coincides with a TNO cell input : - cfg : the configuration file @@ -317,8 +357,12 @@ def interpolate_single_cell(cfg,points): """ """Information about the cosmo grid""" - cosmo_xlocs = np.arange(cfg["xmin"], cfg["xmin"]+cfg["dx"]*cfg["nx"], cfg["dx"]) - cosmo_ylocs = np.arange(cfg["ymin"], cfg["ymin"]+cfg["dy"]*cfg["ny"], cfg["dy"]) + cosmo_xlocs = np.arange( + cfg["xmin"], cfg["xmin"] + cfg["dx"] * cfg["nx"], cfg["dx"] + ) + cosmo_ylocs = np.arange( + cfg["ymin"], cfg["ymin"] + cfg["dy"] * cfg["ny"], cfg["dy"] + ) """ Be careful with numpy.arange(). Floating point numbers are not exactly @@ -333,7 +377,7 @@ def interpolate_single_cell(cfg,points): """This is the interpolation that will be returned""" """Initialization""" - mapping=[] + mapping = [] """Make a polygon out of the points, and get its area.""" """Note : the unit/meaning of the area is in degree^2, but it doesn't matter since we only want the ratio. we assume the area is on a flat earth.""" @@ -342,98 +386,123 @@ def interpolate_single_cell(cfg,points): area_tno = polygon_tno.area """Find the cosmo cells that intersect it""" - for (a,x) in enumerate(cosmo_xlocs): + for (a, x) in enumerate(cosmo_xlocs): """Get the corners of the cosmo cell""" - cosmo_cell_x = [x+cfg["dx"]/2,x+cfg["dx"]/2,x-cfg["dx"]/2,x-cfg["dx"]/2] - if (min(cosmo_cell_x)>max([k[0] for k in points])) or (max(cosmo_cell_x) max([k[0] for k in points])) or ( + max(cosmo_cell_x) < min([k[0] for k in points]) + ): continue - for (b,y) in enumerate(cosmo_ylocs): - cosmo_cell_y = [y+cfg["dy"]/2,y-cfg["dy"]/2,y-cfg["dy"]/2,y+cfg["dy"]/2] - - if (min(cosmo_cell_y)>max([k[1] for k in points])) or (max(cosmo_cell_y) max([k[1] for k in points])) or ( + max(cosmo_cell_y) < min([k[1] for k in points]) + ): continue - points_cosmo = [k for k in zip(cosmo_cell_x,cosmo_cell_y)] + points_cosmo = [k for k in zip(cosmo_cell_x, cosmo_cell_y)] polygon_cosmo = Polygon(points_cosmo) if polygon_cosmo.intersects(polygon_tno): inter = polygon_cosmo.intersection(polygon_tno) - mapping.append((a,b,inter.area/area_tno)) + mapping.append((a, b, inter.area / area_tno)) return mapping -def cell_corners(lon_var,lat_var,inv_name,i,j,cfg): +def cell_corners(lon_var, lat_var, inv_name, i, j, cfg): if inv_name == "tno": x_tno = lon_var[i] y_tno = lat_var[j] - cell_x= np.array([ - x_tno+cfg.tno_dx/2, - x_tno+cfg.tno_dx/2, - x_tno-cfg.tno_dx/2, - x_tno-cfg.tno_dx/2]) - cell_y= np.array([ - y_tno+cfg.tno_dy/2, - y_tno-cfg.tno_dy/2, - y_tno-cfg.tno_dy/2, - y_tno+cfg.tno_dy/2]) + cell_x = np.array( + [ + x_tno + cfg.tno_dx / 2, + x_tno + cfg.tno_dx / 2, + x_tno - cfg.tno_dx / 2, + x_tno - cfg.tno_dx / 2, + ] + ) + cell_y = np.array( + [ + y_tno + cfg.tno_dy / 2, + y_tno - cfg.tno_dy / 2, + y_tno - cfg.tno_dy / 2, + y_tno + cfg.tno_dy / 2, + ] + ) proj = ccrs.PlateCarree() elif inv_name == "vprm": - globe = ccrs.Globe(ellipse=None,semimajor_axis = 6370000,semiminor_axis = 6370000) + globe = ccrs.Globe( + ellipse=None, semimajor_axis=6370000, semiminor_axis=6370000 + ) lambert = ccrs.LambertConformal( - central_longitude = 12.5, central_latitude = 51.604, - standard_parallels=[51.604],globe=globe) - - center_lambert = lambert.transform_point(lon_var[j,i],lat_var[j,i],ccrs.PlateCarree()) - cell_x = np.array([ - center_lambert[0]+cfg.tno_dx/2, - center_lambert[0]+cfg.tno_dx/2, - center_lambert[0]-cfg.tno_dx/2, - center_lambert[0]-cfg.tno_dx/2]) - cell_y = np.array([ - center_lambert[1]+cfg.tno_dy/2, - center_lambert[1]-cfg.tno_dy/2, - center_lambert[1]-cfg.tno_dy/2, - center_lambert[1]+cfg.tno_dy/2]) + central_longitude=12.5, + central_latitude=51.604, + standard_parallels=[51.604], + globe=globe, + ) + + center_lambert = lambert.transform_point( + lon_var[j, i], lat_var[j, i], ccrs.PlateCarree() + ) + cell_x = np.array( + [ + center_lambert[0] + cfg.tno_dx / 2, + center_lambert[0] + cfg.tno_dx / 2, + center_lambert[0] - cfg.tno_dx / 2, + center_lambert[0] - cfg.tno_dx / 2, + ] + ) + cell_y = np.array( + [ + center_lambert[1] + cfg.tno_dy / 2, + center_lambert[1] - cfg.tno_dy / 2, + center_lambert[1] - cfg.tno_dy / 2, + center_lambert[1] + cfg.tno_dy / 2, + ] + ) proj = lambert elif inv_name == "edgar": x_tno = lon_var[i] y_tno = lat_var[j] - cell_x= np.array([ - x_tno+cfg.edgar_dx, - x_tno+cfg.edgar_dx, - x_tno, - x_tno]) - cell_y= np.array([ - y_tno+cfg.edgar_dy, - y_tno-cfg.edgar_dy, - y_tno, - y_tno]) + cell_x = np.array( + [x_tno + cfg.edgar_dx, x_tno + cfg.edgar_dx, x_tno, x_tno] + ) + cell_y = np.array( + [y_tno + cfg.edgar_dy, y_tno - cfg.edgar_dy, y_tno, y_tno] + ) + proj = ccrs.PlateCarree() + elif ( + inv_name == "meteotest" + or inv_name == "maiolica" + or inv_name == "carbocount" + ): + x1_ch, y1_ch = swiss2wgs84(lat_var[j], lon_var[i]) # i-lon, j-lat + x2_ch, y2_ch = swiss2wgs84(lat_var[j] + 200, lon_var[i] + 200) + cell_x = np.array([x2_ch, x2_ch, x1_ch, x1_ch]) + cell_y = np.array([y2_ch, y1_ch, y1_ch, y2_ch]) proj = ccrs.PlateCarree() - elif inv_name == 'meteotest' or inv_name == 'maiolica' or \ - inv_name == 'carbocount': - x1_ch, y1_ch = swiss2wgs84(lat_var[j],lon_var[i]) # i-lon, j-lat - x2_ch, y2_ch = swiss2wgs84(lat_var[j]+200,lon_var[i]+200) - cell_x= np.array([ - x2_ch, - x2_ch, - x1_ch, - x1_ch]) - cell_y= np.array([ - y2_ch, - y1_ch, - y1_ch, - y2_ch]) - proj = ccrs.PlateCarree() else: - print("Inventory %s is not supported yet. Consider defining your own or using tno or vprm." % inv_name) - - + print( + "Inventory %s is not supported yet. Consider defining your own or using tno or vprm." + % inv_name + ) - return cell_x,cell_y,proj + return cell_x, cell_y, proj -def get_dim_var(inv,inv_name,cfg): +def get_dim_var(inv, inv_name, cfg): if inv_name == "tno": lon_dim = inv.dimensions["longitude"].size lat_dim = inv.dimensions["latitude"].size @@ -445,80 +514,103 @@ def get_dim_var(inv,inv_name,cfg): lon_var = inv["lon"][:] lat_var = inv["lat"][:] elif inv_name == "edgar": - lon_var = np.arange(cfg.edgar_xmin,cfg.edgar_xmax,cfg.edgar_dx) - lat_var = np.arange(cfg.edgar_ymin,cfg.edgar_ymax,cfg.edgar_dy) + lon_var = np.arange(cfg.edgar_xmin, cfg.edgar_xmax, cfg.edgar_dx) + lat_var = np.arange(cfg.edgar_ymin, cfg.edgar_ymax, cfg.edgar_dy) lon_dim = len(lon_var) lat_dim = len(lat_var) - elif inv_name == 'meteotest' or inv_name == 'maiolica' or \ - inv_name == 'carbocount': - lon_var = np.array( [ cfg.ch_xll+i*cfg.ch_cell for i in range(0,cfg.ch_xn) ] ) - lat_var = np.array( [ cfg.ch_yll+i*cfg.ch_cell for i in range(0,cfg.ch_yn) ] ) + elif ( + inv_name == "meteotest" + or inv_name == "maiolica" + or inv_name == "carbocount" + ): + lon_var = np.array( + [cfg.ch_xll + i * cfg.ch_cell for i in range(0, cfg.ch_xn)] + ) + lat_var = np.array( + [cfg.ch_yll + i * cfg.ch_cell for i in range(0, cfg.ch_yn)] + ) lon_dim = np.shape(lon_var)[0] lat_dim = np.shape(lat_var)[0] else: - print("Inventory %s is not supported yet. Consider defining your own or using tno or vprm." % inv_name) + print( + "Inventory %s is not supported yet. Consider defining your own or using tno or vprm." + % inv_name + ) - return lon_dim,lat_dim,lon_var,lat_var + return lon_dim, lat_dim, lon_var, lat_var -def interpolate_to_cosmo_grid(tno,inv_name,cfg,filename): - print("Retrieving the interpolation between the cosmo and the inventory grids") +def interpolate_to_cosmo_grid(tno, inv_name, cfg, filename): + print( + "Retrieving the interpolation between the cosmo and the inventory grids" + ) start = time.time() - - transform = ccrs.RotatedPole(pole_longitude=cfg.pollon, pole_latitude=cfg.pollat) - lon_dim,lat_dim,lon_var,lat_var = get_dim_var(tno,inv_name,cfg) + transform = ccrs.RotatedPole( + pole_longitude=cfg.pollon, pole_latitude=cfg.pollat + ) + + lon_dim, lat_dim, lon_var, lat_var = get_dim_var(tno, inv_name, cfg) """This is the interpolation that will be returned""" - mapping = np.empty((lon_dim,lat_dim),dtype=object) + mapping = np.empty((lon_dim, lat_dim), dtype=object) # Annoying ... """Transform the cfg to a dictionary""" - config = dict({ - "dx" : cfg.dx, - "dy" : cfg.dy, - "nx" : cfg.nx, - "ny" : cfg.ny, - "xmin" : cfg.xmin, - "ymin" : cfg.ymin}) - + config = dict( + { + "dx": cfg.dx, + "dy": cfg.dy, + "nx": cfg.nx, + "ny": cfg.ny, + "xmin": cfg.xmin, + "ymin": cfg.ymin, + } + ) with Pool(npool) as pool: for i in range(lon_dim): - print("ongoing :",i) + print("ongoing :", i) points = [] for j in range(lat_dim): - tno_cell_x,tno_cell_y,proj = cell_corners(lon_var,lat_var,inv_name,i,j,cfg) - points.append(transform.transform_points(proj,tno_cell_x,tno_cell_y)) - - mapping[i,:] = pool.starmap(interpolate_single_cell,[(config,points[j]) for j in range(lat_dim)]) + tno_cell_x, tno_cell_y, proj = cell_corners( + lon_var, lat_var, inv_name, i, j, cfg + ) + points.append( + transform.transform_points(proj, tno_cell_x, tno_cell_y) + ) + + mapping[i, :] = pool.starmap( + interpolate_single_cell, + [(config, points[j]) for j in range(lat_dim)], + ) end = time.time() print("\nInterpolation is over") - print("it took ",end-start,"seconds") + print("it took ", end - start, "seconds") - np.save(os.path.join(cfg.output_path,filename),mapping) + np.save(os.path.join(cfg.output_path, filename), mapping) return mapping -def get_interpolation(cfg,tno,inv_name = "tno",filename="mapping.npy"): +def get_interpolation(cfg, tno, inv_name="tno", filename="mapping.npy"): """retrieve the interpolation between the tno and cosmo grids.""" - make_interpolate = True - mapping_path = os.path.join(cfg.output_path,filename) + make_interpolate = True + mapping_path = os.path.join(cfg.output_path, filename) if os.path.isfile(mapping_path): print("Do you wanna overwite the mapping found in %s ?" % mapping_path) s = input("y/[n] \n") - make_interpolate = (s=="y") - + make_interpolate = s == "y" + if make_interpolate: - interpolation = interpolate_to_cosmo_grid(tno,inv_name,cfg,filename) + interpolation = interpolate_to_cosmo_grid(tno, inv_name, cfg, filename) else: interpolation = np.load(mapping_path) return interpolation -def swiss2wgs84(x,y): +def swiss2wgs84(x, y): """ Convert Swiss LV03 coordinates (x easting and y northing) to WGS 84 based on swisstopo approximated soluation (0.1" accuracy). @@ -529,21 +621,20 @@ def swiss2wgs84(x,y): y = (y - 600000.0) / 1000000.0 lon = ( - 2.6779094 - + 4.728982 * y - + 0.791484 * y * x - + 0.1306 * y * x**2 - - 0.0436 * y**3 + 2.6779094 + + 4.728982 * y + + 0.791484 * y * x + + 0.1306 * y * x ** 2 + - 0.0436 * y ** 3 ) / 0.36 lat = ( - 16.9023892 - + 3.238272 * x - - 0.270978 * y**2 - - 0.002528 * x**2 - - 0.0447 * y**2 * x - - 0.0140 * x**3 + 16.9023892 + + 3.238272 * x + - 0.270978 * y ** 2 + - 0.002528 * x ** 2 + - 0.0447 * y ** 2 * x + - 0.0140 * x ** 3 ) / 0.36 return lon, lat - diff --git a/merge_inventories.py b/merge_inventories.py index 2df78fa..7a7b7f0 100644 --- a/merge_inventories.py +++ b/merge_inventories.py @@ -19,13 +19,15 @@ import grid inpath = "/newhome/jae/gitlab/online-emission-processing/testdata/" -base_inventory = os.path.join(inpath, 'emis_2018_tno_CO2_CO_CH4_1km.nc') -outfile = os.path.join(inpath, 'merged_inventories.nc') +base_inventory = os.path.join(inpath, "emis_2018_tno_CO2_CO_CH4_1km.nc") +outfile = os.path.join(inpath, "merged_inventories.nc") + +list_inventories = [ + "emis_2018_maiolica_CH4_1km.nc", + "emis_2018_carbocount_CO2_1km.nc", + "emis_2018_meteotest_CO_1km.nc", +] -list_inventories = ['emis_2018_maiolica_CH4_1km.nc', - 'emis_2018_carbocount_CO2_1km.nc', - 'emis_2018_meteotest_CO_1km.nc' - ] def main(cfg_path): """ The main script for processing TNO inventory. @@ -37,7 +39,7 @@ def main(cfg_path): """Load or compute the country mask""" country_mask = get_country_mask(cfg) country_mask = np.transpose(country_mask) - mask = country_mask == country_codes['CH'] + mask = country_mask == country_codes["CH"] sum_co2_tno = 0 sum_co2_merged = 0 @@ -45,29 +47,35 @@ def main(cfg_path): sum_co_merged = 0 sum_ch4_tno = 0 sum_ch4_merged = 0 - with nc.Dataset(base_inventory, "r") as src, nc.Dataset(outfile, "w") as dst: + with nc.Dataset(base_inventory, "r") as src, nc.Dataset( + outfile, "w" + ) as dst: print(base_inventory) # Copy global attributes all at once via dictionary dst.setncatts(src.__dict__) - # Copy all dimensions + # Copy all dimensions for name, dimension in src.dimensions.items(): dst.createDimension(name, len(dimension)) # Copy all variables for name, variable in src.variables.items(): print(name) - x = dst.createVariable(name, variable.datatype, - variable.dimensions, - zlib=True, complevel=9) + x = dst.createVariable( + name, + variable.datatype, + variable.dimensions, + zlib=True, + complevel=9, + ) # Copy variable attributes all at once via dictionary dst[name].setncatts(src[name].__dict__) - if 'rotated_pole' not in name: + if "rotated_pole" not in name: dst[name][:] = src[name][:] - sum_co2_tno = np.sum(src['CO2'][:]) - sum_co_tno = np.sum(src['CO'][:]) - sum_ch4_tno = np.sum(src['CH4'][:]) + sum_co2_tno = np.sum(src["CO2"][:]) + sum_co_tno = np.sum(src["CO"][:]) + sum_ch4_tno = np.sum(src["CH4"][:]) # Merge inventories for fname in list_inventories: @@ -75,27 +83,33 @@ def main(cfg_path): print(infile) with nc.Dataset(infile, "r") as inv: for name, variable in src.variables.items(): - if ('CO2' in name and 'carbocount' in infile) or \ - ('CH4' in name and 'maiolica' in infile) or \ - (('CO_' in name or name == 'CO') and 'meteotest' in infile): - print('Overwriting variable %s' % name) + if ( + ("CO2" in name and "carbocount" in infile) + or ("CH4" in name and "maiolica" in infile) + or ( + ("CO_" in name or name == "CO") + and "meteotest" in infile + ) + ): + print("Overwriting variable %s" % name) var_dst = dst.variables[name][:] - np.copyto(var_dst, inv.variables[name][:], - where=mask[:]) + np.copyto( + var_dst, inv.variables[name][:], where=mask[:] + ) dst.variables[name][:] = var_dst - sum_co2_merged = np.sum(dst['CO2'][:]) - sum_co_merged = np.sum(dst['CO'][:]) - sum_ch4_merged = np.sum(dst['CH4'][:]) + sum_co2_merged = np.sum(dst["CO2"][:]) + sum_co_merged = np.sum(dst["CO"][:]) + sum_ch4_merged = np.sum(dst["CH4"][:]) + + print("sum_co2_tno", sum_co2_tno) + print("sum_co2_merged", sum_co2_merged) + print("sum_co_tno", sum_co_tno) + print("sum_co_merged", sum_co_merged) + print("sum_ch4_tno", sum_ch4_tno) + print("sum_ch4_merged", sum_ch4_merged) - print('sum_co2_tno',sum_co2_tno) - print('sum_co2_merged',sum_co2_merged) - print('sum_co_tno',sum_co_tno) - print('sum_co_merged',sum_co_merged) - print('sum_ch4_tno',sum_ch4_tno) - print('sum_ch4_merged',sum_ch4_merged) -if __name__ == '__main__': +if __name__ == "__main__": cfg_name = sys.argv[1] main("./config_" + cfg_name) - diff --git a/profiles/temporal_profiles.py b/profiles/temporal_profiles.py index 7f35340..6c80c82 100644 --- a/profiles/temporal_profiles.py +++ b/profiles/temporal_profiles.py @@ -40,8 +40,8 @@ # variables: # float varname1(monthofyear, country) ; # varname1:units = "1" ; -# varname1:long_name = "monthly scaling factor for 12 months"; -# varname1:comment = "first month is Jan, last month is Dec"; +# varname1:long_name = "monthly scaling factor for 12 months"; +# varname1:comment = "first month is Jan, last month is Dec"; # [...] # short countryID(country) ; # countryID:long_name = "EMEP country code" ; @@ -66,33 +66,34 @@ import numpy as np import netCDF4 import pytz -from datetime import datetime +from datetime import datetime from country_code import country_codes as cc """Weekly and annual profiles are availble for these tracers""" -#TRACERS = ['CO', 'CO2', 'NH3', 'NMVOC', 'NO', 'NO2', 'NOx', 'P10', 'P25', 'SOx'] -TRACERS = ['CO2', 'CO', 'CH4'] +# TRACERS = ['CO', 'CO2', 'NH3', 'NMVOC', 'NO', 'NO2', 'NOx', 'P10', 'P25', 'SOx'] +TRACERS = ["CO2", "CO", "CH4"] only_ones = False """Set output type (normal, CH0, outCH0)""" -output_type = 'normal' -#output_type = 'CH0' -#output_type = 'outCH0' +output_type = "normal" +# output_type = 'CH0' +# output_type = 'outCH0' -def permute_cycle_tz(tz,cycle): +def permute_cycle_tz(tz, cycle): # tz in the shape "+0400" - shift = int(int(tz[1:3])*((tz[0]=="+")-0.5)*2) + shift = int(int(tz[1:3]) * ((tz[0] == "+") - 0.5) * 2) # for now I don't consider half hours for the time zone try: - answer = [cycle[shift-i] for i in range(len(cycle),0,-1)] + answer = [cycle[shift - i] for i in range(len(cycle), 0, -1)] except IndexError: - answer = [cycle[shift-i+24] for i in range(len(cycle),0,-1)] + answer = [cycle[shift - i + 24] for i in range(len(cycle), 0, -1)] return answer + def get_country_tz(countries): country_exception = dict( FGD="+0100", @@ -100,58 +101,60 @@ def get_country_tz(countries): RUO="+0400", RUP="+0300", RUA="+0200", - RUR="+0300", # Moscau - RU ="+0300", #Moscau + RUR="+0300", # Moscau + RU="+0300", # Moscau YUG="+0100", - CS ="+0100", + CS="+0100", NOA="+0100", - PT ="+0000", #Portugal mainland, not azores - KZ ="+0500", #West of Kazakhstan + PT="+0000", # Portugal mainland, not azores + KZ="+0500", # West of Kazakhstan ) all_tz = dict() for country in countries: try: - country_names = [name for name,code in cc.items() if (code==country)] - country_name="" + country_names = [ + name for name, code in cc.items() if (code == country) + ] + country_name = "" # Try find the name of the country, with 2 character for name in country_names: - if len(name)==2: + if len(name) == 2: country_name = name # If it's not there, maybe it's one of the exceptions - if country_name=="": + if country_name == "": for name in country_names: try: - print("exception : ",name,country_exception[name]) + print("exception : ", name, country_exception[name]) country_name = name continue except KeyError: pass # If I still didn't find it, then nevermind - if country_name=="": + if country_name == "": raise KeyError except KeyError: - print(country_names,"not a proper country name in pytz") + print(country_names, "not a proper country name in pytz") continue # For countries with several time zones, they are listed in the exception and one tz is assigned. if country_name in country_exception: - all_tz[country]= country_exception[country_name] - else: + all_tz[country] = country_exception[country_name] + else: zones = pytz.country_timezones[country_name] - + fmt = "%z" - zero = pytz.utc.localize(datetime(2015, 1, 1, 0, 0, 0)) + zero = pytz.utc.localize(datetime(2015, 1, 1, 0, 0, 0)) for zone in zones: loc_tz = pytz.timezone(zone) loc_zero = zero.astimezone(loc_tz) hour = loc_zero.strftime("%z") - all_tz[country]= hour + all_tz[country] = hour return all_tz def read_daily_profiles(path): - filename = os.path.join(path, 'HourlyFac.dat') + filename = os.path.join(path, "HourlyFac.dat") snaps = [] data = {} @@ -160,23 +163,24 @@ def read_daily_profiles(path): for line in profile_file: values = line.split() snap = values[0].strip() - data[snap] = np.array(values[1:], 'f4') + data[snap] = np.array(values[1:], "f4") return snaps, data + def read_temporal_profile_simple(path): started = False data = [] snaps = [] - with open(path,encoding="latin") as prof: + with open(path, encoding="latin") as prof: for l in prof: splitted = l.split(";") - if splitted[0]=="1": - started=True + if splitted[0] == "1": + started = True # if splitted[0]=="11": # started=False if started: - if splitted[1]=="F1": + if splitted[1] == "F1": snaps.append("F") elif "F" in splitted[1]: continue @@ -184,9 +188,10 @@ def read_temporal_profile_simple(path): snaps.append(splitted[1]) data.append([float(i) for i in splitted[3:]]) - return snaps,np.array(data) + return snaps, np.array(data) -def read_temporal_profiles(tracer, kind, path='timeprofiles'): + +def read_temporal_profiles(tracer, kind, path="timeprofiles"): """\ Read temporal profiles for given `tracer` for 'weekly' or 'annual' profiles. @@ -195,14 +200,14 @@ def read_temporal_profiles(tracer, kind, path='timeprofiles'): countries = [] snaps = [] - if kind == 'weekly': - filename = os.path.join(path, 'DailyFac_%s.dat' % tracer) - elif kind == 'annual': - filename = os.path.join(path, 'MonthlyFac_%s.dat' % tracer) + if kind == "weekly": + filename = os.path.join(path, "DailyFac_%s.dat" % tracer) + elif kind == "annual": + filename = os.path.join(path, "MonthlyFac_%s.dat" % tracer) else: raise ValueError("kind has to be 'weekly' or 'annual' not '%s'" % kind) - with open(filename, 'r') as profile_file: + with open(filename, "r") as profile_file: for line in profile_file: values = line.split() country, snap = int(values[0]), str(values[1]) @@ -210,13 +215,12 @@ def read_temporal_profiles(tracer, kind, path='timeprofiles'): countries.append(country) snaps.append(snap) - data[country, snap] = np.array(values[2:], 'f4') + data[country, snap] = np.array(values[2:], "f4") return list(set(countries)), list(set(snaps)), data - -def read_all_data(tracers, path='timeprofiles'): +def read_all_data(tracers, path="timeprofiles"): daily_profiles = {} weekly_profiles = {} @@ -230,88 +234,102 @@ def read_all_data(tracers, path='timeprofiles'): for tracer in tracers: # weekly - c,s,d = read_temporal_profiles(tracer, 'weekly', path) + c, s, d = read_temporal_profiles(tracer, "weekly", path) weekly_profiles[tracer] = d countries += c snaps += s # weekly - c,s,d = read_temporal_profiles(tracer, 'annual', path) + c, s, d = read_temporal_profiles(tracer, "annual", path) annual_profiles[tracer] = d countries += c snaps += s - return sorted(set(countries)), sorted(set(snaps)), daily_profiles, weekly_profiles, annual_profiles + return ( + sorted(set(countries)), + sorted(set(snaps)), + daily_profiles, + weekly_profiles, + annual_profiles, + ) -def create_netcdf(path,countries): +def create_netcdf(path, countries): profile_names = ["hourofday", "dayofweek", "monthofyear", "hourofyear"] - var_names = ["hourofday", "dayofweek", "monthofyear", "hourofyear"] - if output_type == 'normal': + var_names = ["hourofday", "dayofweek", "monthofyear", "hourofyear"] + if output_type == "normal": pass - elif output_type == 'CH0': - profile_names[0] = 'hourofday_CH0' - elif output_type == 'outCH0': - profile_names[0] = 'hourofday_outCH0' + elif output_type == "CH0": + profile_names[0] = "hourofday_CH0" + elif output_type == "outCH0": + profile_names[0] = "hourofday_outCH0" else: - raise ValueError("Variable 'output_type' has to be 'normal' or 'CH0' " - "or 'outCH0', not '%s'" % output_type) + raise ValueError( + "Variable 'output_type' has to be 'normal' or 'CH0' " + "or 'outCH0', not '%s'" % output_type + ) - for (profile, var, size) in zip(profile_names, var_names, [24,7,12,8784]): - filename = os.path.join(path,profile + ".nc") - - with netCDF4.Dataset(filename, 'w') as nc: + for (profile, var, size) in zip( + profile_names, var_names, [24, 7, 12, 8784] + ): + filename = os.path.join(path, profile + ".nc") + + with netCDF4.Dataset(filename, "w") as nc: # global attributes (add input data) - nc.setncatts({ - 'DESCRIPTION': 'Temporal profiles for emissions', - 'DATAORIGIN': 'TNO time profiles', - 'CREATOR': 'Jean-Matthieu Haussaire', - 'EMAIL': 'jean-matthieu.haussaire@empa.ch', - 'AFFILIATION': 'Empa Duebendorf, Switzerland', - 'DATE CREATED': time.ctime(time.time()), - }) - - # create dimensions + nc.setncatts( + { + "DESCRIPTION": "Temporal profiles for emissions", + "DATAORIGIN": "TNO time profiles", + "CREATOR": "Jean-Matthieu Haussaire", + "EMAIL": "jean-matthieu.haussaire@empa.ch", + "AFFILIATION": "Empa Duebendorf, Switzerland", + "DATE CREATED": time.ctime(time.time()), + } + ) + + # create dimensions nc.createDimension(var, size=size) - nc.createDimension('country', size=len(countries)) + nc.createDimension("country", size=len(countries)) + + nc_cid = nc.createVariable("country", "i2", ("country")) + nc_cid[:] = np.array(countries, "i2") + nc_cid.long_name = "EMEP country code" - nc_cid = nc.createVariable('country', 'i2', ('country')) - nc_cid[:] = np.array(countries, 'i2') - nc_cid.long_name = 'EMEP country code' - def write_single_variable(path, profile, values, tracer, snap): filename = os.path.join(path, profile + ".nc") - if 'hourofday' in profile: - descr = 'diurnal scaling factor' - varname = 'hourofday' - if profile == "dayofweek": - descr ='day-of-week scaling factor' - varname = 'dayofweek' + if "hourofday" in profile: + descr = "diurnal scaling factor" + varname = "hourofday" + if profile == "dayofweek": + descr = "day-of-week scaling factor" + varname = "dayofweek" if profile == "monthofyear": - descr = 'month-of-year scaling factor' - varname = 'monthofyear' + descr = "month-of-year scaling factor" + varname = "monthofyear" if profile == "hourofyear": - descr = 'hour-of-year scaling factor' - varname = 'hourofyear' - - - with netCDF4.Dataset(filename, 'a') as nc: + descr = "hour-of-year scaling factor" + varname = "hourofyear" + + with netCDF4.Dataset(filename, "a") as nc: # create variables and attributes # if profile == "hourofday" or profile == "hourofyear": # nc_var = nc.createVariable(tracer+"_"+snap, 'f4', (profile)) # else: - nc_var = nc.createVariable(tracer+"_"+snap, 'f4', (varname, 'country')) + nc_var = nc.createVariable( + tracer + "_" + snap, "f4", (varname, "country") + ) if complex_profile: - nc_var.long_name = "%s for species %s and SNAP %s" % (descr,tracer,snap) + nc_var.long_name = "%s for species %s and SNAP %s" % ( + descr, + tracer, + snap, + ) else: - nc_var.long_name = "%s for GNFR %s" % (descr,snap) - nc_var.units = '1' - nc_var[:]= values - - - + nc_var.long_name = "%s for GNFR %s" % (descr, snap) + nc_var.units = "1" + nc_var[:] = values def main_complex(path): @@ -320,63 +338,74 @@ def main_complex(path): # read all data countries, snaps, daily, weekly, annual = read_all_data(TRACERS) - countries = [0]+countries + countries = [0] + countries n_countries = len(countries) - create_netcdf(path,countries) + create_netcdf(path, countries) # day of week and month of year dow = np.ones((7, n_countries)) moy = np.ones((12, n_countries)) hod = np.ones((24, n_countries)) - + country_tz = get_country_tz(countries) - - - for (tracer,snap) in itertools.product(TRACERS, snaps): + + for (tracer, snap) in itertools.product(TRACERS, snaps): if not only_ones: - for i, country in enumerate(countries): + for i, country in enumerate(countries): try: - hod[:,i] = permute_cycle_tz(country_tz[country],daily[snap]) + hod[:, i] = permute_cycle_tz( + country_tz[country], daily[snap] + ) except KeyError: pass try: - dow[:,i] = weekly[tracer][country, snap] + dow[:, i] = weekly[tracer][country, snap] if mean: - dow[:5,i] = np.ones(5)*weekly[tracer][country, snap][:5].mean() + dow[:5, i] = ( + np.ones(5) + * weekly[tracer][country, snap][:5].mean() + ) except KeyError: pass try: - moy[:,i] = annual[tracer][country, snap] + moy[:, i] = annual[tracer][country, snap] except KeyError: pass - write_single_variable(path,"hourofday",hod,tracer,snap) - write_single_variable(path,"dayofweek",dow,tracer,snap) - write_single_variable(path,"monthofyear",moy,tracer,snap) + write_single_variable(path, "hourofday", hod, tracer, snap) + write_single_variable(path, "dayofweek", dow, tracer, snap) + write_single_variable(path, "monthofyear", moy, tracer, snap) # TODO: hour of year hours_in_year = np.arange(0, 8784) + def main_simple(path): mean = False - #snaps = ["A","B","C","D","E","F","G","H","I","J"] + # snaps = ["A","B","C","D","E","F","G","H","I","J"] countries = np.arange(74) # remove all countries not known, and not worth an exception - countries = np.delete(countries, - [5,26,28,29,30,31,32,33,34,35,58,64,67,70,71]) + countries = np.delete( + countries, [5, 26, 28, 29, 30, 31, 32, 33, 34, 35, 58, 64, 67, 70, 71] + ) n_countries = len(countries) country_tz = get_country_tz(countries) + create_netcdf(path, countries) - create_netcdf(path,countries) - - snaps,daily = read_temporal_profile_simple("CHE_input/timeprofiles-hour-in-day_GNFR.csv") - snaps,weekly = read_temporal_profile_simple("CHE_input/timeprofiles-day-in-week_GNFR.csv") - snaps,monthly = read_temporal_profile_simple("CHE_input/timeprofiles-month-in-year_GNFR.csv") + snaps, daily = read_temporal_profile_simple( + "CHE_input/timeprofiles-hour-in-day_GNFR.csv" + ) + snaps, weekly = read_temporal_profile_simple( + "CHE_input/timeprofiles-day-in-week_GNFR.csv" + ) + snaps, monthly = read_temporal_profile_simple( + "CHE_input/timeprofiles-month-in-year_GNFR.csv" + ) print(snaps) @@ -385,51 +414,52 @@ def main_simple(path): moy = np.ones((12, n_countries)) hod = np.ones((24, n_countries)) - for snap_ind,snap in enumerate(snaps): - for i, country in enumerate(countries): + for snap_ind, snap in enumerate(snaps): + for i, country in enumerate(countries): try: - hod[:,i] = permute_cycle_tz(country_tz[country],daily[snap_ind,:]) + hod[:, i] = permute_cycle_tz( + country_tz[country], daily[snap_ind, :] + ) except KeyError: pass try: - dow[:,i] = weekly[snap_ind,:] + dow[:, i] = weekly[snap_ind, :] if mean: - dow[:5,i] = np.ones(5)*weekly[snap_ind,:5].mean() + dow[:5, i] = np.ones(5) * weekly[snap_ind, :5].mean() except KeyError: pass try: - moy[:,i] = monthly[snap_ind] + moy[:, i] = monthly[snap_ind] except KeyError: pass - # Use hourofday profile to mask CH if needed - if output_type == 'normal': - write_single_variable(path, "hourofday" ,hod, "GNFR", snap) - elif output_type == 'CH0': - hod[:,23] = 0. + if output_type == "normal": + write_single_variable(path, "hourofday", hod, "GNFR", snap) + elif output_type == "CH0": + hod[:, 23] = 0.0 write_single_variable(path, "hourofday_CH0", hod, "GNFR", snap) - elif output_type == 'outCH0': + elif output_type == "outCH0": for i, country in enumerate(countries): - if i != 23: hod[:,i] = 0. + if i != 23: + hod[:, i] = 0.0 write_single_variable(path, "hourofday_outCH0", hod, "GNFR", snap) else: - raise ValueError("Variable 'output_type' has to be 'normal' or 'CH0' " - "or 'outCH0', not '%s'" % output_type) + raise ValueError( + "Variable 'output_type' has to be 'normal' or 'CH0' " + "or 'outCH0', not '%s'" % output_type + ) # Write other profiles - write_single_variable(path,"dayofweek",dow,"GNFR",snap) - write_single_variable(path,"monthofyear",moy,"GNFR",snap) - - + write_single_variable(path, "dayofweek", dow, "GNFR", snap) + write_single_variable(path, "monthofyear", moy, "GNFR", snap) + if __name__ == "__main__": - complex_profile=False + complex_profile = False if complex_profile: main_complex("./output") else: main_simple("./output") - - diff --git a/profiles/vertical_profiles.py b/profiles/vertical_profiles.py index b11f132..ed801e6 100644 --- a/profiles/vertical_profiles.py +++ b/profiles/vertical_profiles.py @@ -22,54 +22,57 @@ # SNAP-*:units = "1" ; # SNAP-*:long_name = "vertical scale factor for sources of SNAP-* category"; -def get_all_levels(levels): #levels are the top of the layer + +def get_all_levels(levels): # levels are the top of the layer layer_top = levels - layer_bot = [0]+levels[:-1] - layer_mid = [(i+j)/2 for (i,j) in zip(layer_top,layer_bot)] - - return layer_bot,layer_mid,layer_top - + layer_bot = [0] + levels[:-1] + layer_mid = [(i + j) / 2 for (i, j) in zip(layer_top, layer_bot)] + + return layer_bot, layer_mid, layer_top def write_netcdf(filename, categories, cat_name, levels, scale_factors): - layer_bot,layer_mid,layer_top = get_all_levels(levels) + layer_bot, layer_mid, layer_top = get_all_levels(levels) - with netCDF4.Dataset(filename, 'w') as nc: + with netCDF4.Dataset(filename, "w") as nc: # global attributes (add input data) - nc.setncatts({ - 'DESCRIPTION': 'Vertical profiles for emissions', - 'DATAORIGIN': 'based on profiles developed for COST-728 action', - 'CREATOR': 'Jean-Matthieu Haussaire', - 'EMAIL': 'jean-matthieu.haussaire@empa.ch', - 'AFFILIATION': 'Empa Duebendorf, Switzerland', - 'DATE CREATED': time.ctime(time.time()), - }) - - # create dimensions - nc.createDimension('level', size=len(levels)) - - + nc.setncatts( + { + "DESCRIPTION": "Vertical profiles for emissions", + "DATAORIGIN": "based on profiles developed for COST-728 action", + "CREATOR": "Jean-Matthieu Haussaire", + "EMAIL": "jean-matthieu.haussaire@empa.ch", + "AFFILIATION": "Empa Duebendorf, Switzerland", + "DATE CREATED": time.ctime(time.time()), + } + ) + + # create dimensions + nc.createDimension("level", size=len(levels)) + # create variables and attributes - nc_bot = nc.createVariable('layer_bot', 'f4', ('level')) - nc_bot.long_name = 'bottom of layer above ground' - nc_bot.units = 'm' - nc_bot[:] = layer_bot - - nc_mid = nc.createVariable('layer_mid', 'f4', ('level')) - nc_mid.long_name = 'middle of layer above ground' - nc_mid.units = 'm' - nc_mid[:] = layer_mid - - nc_top = nc.createVariable('layer_top', 'f4', ('level')) - nc_top.long_name = 'top of layer above ground' - nc_top.units = 'm' - nc_top[:] = layer_top - - for (i,cat) in enumerate(categories): - nc_sca = nc.createVariable(cat_name+cat, 'f4', ('level')) - nc_sca.long_name = 'vertical scale factor for sources of %s category' % cat - nc_sca.units = '1' + nc_bot = nc.createVariable("layer_bot", "f4", ("level")) + nc_bot.long_name = "bottom of layer above ground" + nc_bot.units = "m" + nc_bot[:] = layer_bot + + nc_mid = nc.createVariable("layer_mid", "f4", ("level")) + nc_mid.long_name = "middle of layer above ground" + nc_mid.units = "m" + nc_mid[:] = layer_mid + + nc_top = nc.createVariable("layer_top", "f4", ("level")) + nc_top.long_name = "top of layer above ground" + nc_top.units = "m" + nc_top[:] = layer_top + + for (i, cat) in enumerate(categories): + nc_sca = nc.createVariable(cat_name + cat, "f4", ("level")) + nc_sca.long_name = ( + "vertical scale factor for sources of %s category" % cat + ) + nc_sca.units = "1" nc_sca[:] = scale_factors[i] @@ -77,48 +80,46 @@ def read_profiles(filename, nlevel=16): levels = [] categories = [] profiles = [] - + with open(filename) as profile_file: levels = [int(i) for i in profile_file.readline().split("\t")[1:]] - - all_sevens= [] - for line in profile_file: + + all_sevens = [] + for line in profile_file: values = line.split() cat = values[0] - profile = values[1:] - if cat=="34": + profile = values[1:] + if cat == "34": categories.append("3") categories.append("4") profiles.append(profile) profiles.append(profile) - elif cat=="7.01": + elif cat == "7.01": categories.append("7") all_sevens.append([float(i) for i in profile]) elif "7." in cat: all_sevens.append([float(i) for i in profile]) - elif cat=="F1": + elif cat == "F1": categories.append("F") all_sevens.append([float(i) for i in profile]) elif "F" in cat: all_sevens.append([float(i) for i in profile]) else: categories.append(cat) - if len(all_sevens)>0: + if len(all_sevens) > 0: profiles.append(np.array(all_sevens).mean(0)) - all_sevens=[] + all_sevens = [] profiles.append(profile) - - - return categories, np.array(profiles, 'f4'),levels - + return categories, np.array(profiles, "f4"), levels def main(filename): - categories, profiles,levels = read_profiles('vert_profiles/vert_prof_che_gnfr.dat') + categories, profiles, levels = read_profiles( + "vert_profiles/vert_prof_che_gnfr.dat" + ) write_netcdf(filename, categories, "GNFR_", levels, profiles) - if __name__ == "__main__": - main('output/vertical_profiles.nc') + main("output/vertical_profiles.nc")