From cc374bc61473b3be4a9aa82e8267fd12d31ed621 Mon Sep 17 00:00:00 2001 From: mguevara Date: Fri, 9 Nov 2018 11:34:36 +0100 Subject: [PATCH 01/10] Edited README.md file --- README.md | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index b393f1f..03a4695 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,25 @@ -# HERMESv3 Global/Regional +# HERMESv3_GR - High-Elective Resolution Modelling Emission System version 3 – Global_Regional [![Codacy Badge](https://api.codacy.com/project/badge/Grade/34fc5d6c803444178034b99dd28c7e3c)](https://www.codacy.com/app/carlestena/hermesv3_gr?utm_source=earth.bsc.es&utm_medium=referral&utm_content=gitlab/es/hermesv3_gr&utm_campaign=Badge_Grade) + +HERMESv3_GR is a highly customizable emission processing system that calculates emissions from different sources, regions and pollutants over a user-specified global or regional model grid and fulfills the requirements of several CTMs. + +The main characteristics of HERMESv3_GR are as follows: +- User-defined grid and choice between different map projections: Emissions can be computed on any global or regional domain with a regular lat-lon, rotated lat-lon, mercator or lambert conformal conic projection. +- Choice between different emission inventories: the emission data library of HERMESv3_GR includes current state-of-the-art global and regional inventories that cover different sources (anthropogenic, biomass burning, volcanoes), pollutants (ozone precursor gases, acidifying gases and primary particulates) and base years (past, present and future). Moreover, country-specific scaling and masking factors defined by the user can be applied to the base inventories in order to combine and adjust them. +- Choice between different vertical, temporal and speciation profiles: HERMESv3_GR includes a dataset of profiles reported by the literature, but it also allows the user to add its own weight factors for any pollutant sector and specie. Additionally, the model is able to combine base inventories with gridded temporal profiles, which can be of importance for those pollutant sectors whose temporal variation is not uniform across the space due to local influences (e.g. residential combustion emissions and temperature). +- Choice between different CTMs: The generated emission files can be used as input for the CMAQ, WRF-CHEM and NMMB-MONARCH chemical transport models. +- Choice between different chemical mechanisms: base pollutants can be mapped to several gas-phase and aerosol chemical mechanism, including CB05, CB05e51, RADM2, AERO5, AERO6 and MADE/SORGAM. All these mechanisms are widely used in the air quality modelling community. +- Parallel implementation: The emission core module of HERMESv3_GR is parallelized using a map partitioning strategy, which allows decreasing the execution time and memory consumption of the model. This feature can be of importance when using the model in operational air quality forecasting systems, for which the simulations need to be completed within the required time constraints. + +## Availability + +HERMESv3_GR is distributed free of charge under the licence [GNU GPL v3.0](https://www.gnu.org/licenses/quick-guide-gplv3.html). + +## Citing + +Guevara, M., Tena, C., Porquet, M., Jorba, O., Pérez García-Pando, C., 2018. HERMESv3_GR: A stand-alone multiscale emission processing system. Extended abstract of the 17th Annual CMAS Conference, Chapell Hill, NC, October 22-24. + +## Support + +Due to our limited time and resources, the developing team can not guarantee a regular support. Nevertheless, we will be happy to provide advice for new users. +Questions should be send to [Marc Guevara](marc.guevara@bsc.es) and [Carles Tena](carles.tena@bsc.es) \ No newline at end of file -- GitLab From 0ba5377a29bac6ead748884e6ad1f06d42bbda54 Mon Sep 17 00:00:00 2001 From: mguevara Date: Fri, 9 Nov 2018 11:38:22 +0100 Subject: [PATCH 02/10] Edited README.md file --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index 03a4695..4aa2182 100644 --- a/README.md +++ b/README.md @@ -11,6 +11,8 @@ The main characteristics of HERMESv3_GR are as follows: - Choice between different chemical mechanisms: base pollutants can be mapped to several gas-phase and aerosol chemical mechanism, including CB05, CB05e51, RADM2, AERO5, AERO6 and MADE/SORGAM. All these mechanisms are widely used in the air quality modelling community. - Parallel implementation: The emission core module of HERMESv3_GR is parallelized using a map partitioning strategy, which allows decreasing the execution time and memory consumption of the model. This feature can be of importance when using the model in operational air quality forecasting systems, for which the simulations need to be completed within the required time constraints. +HERMESv3_GR is an in-house model fully developed by the [Barcelona Supercomputing Center](https://www.bsc.es/). + ## Availability HERMESv3_GR is distributed free of charge under the licence [GNU GPL v3.0](https://www.gnu.org/licenses/quick-guide-gplv3.html). -- GitLab From 2c9de82f04d5649e5f4a1399e90f5f6615076fe7 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Fri, 30 Nov 2018 11:29:23 +0100 Subject: [PATCH 03/10] BUg corrected --- hermesv3_gr/modules/writing/writer_monarch.py | 4 ++-- run_test.py | 18 ++++++++++++++++++ 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/hermesv3_gr/modules/writing/writer_monarch.py b/hermesv3_gr/modules/writing/writer_monarch.py index 3321b06..e501926 100644 --- a/hermesv3_gr/modules/writing/writer_monarch.py +++ b/hermesv3_gr/modules/writing/writer_monarch.py @@ -369,8 +369,8 @@ class WriterMonarch(Writer): # Rotated pole mapping = netcdf.createVariable('rotated_pole', 'c') mapping.grid_mapping_name = 'rotated_latitude_longitude' - mapping.grid_north_pole_latitude = self.grid.new_pole_latitude_degrees - mapping.grid_north_pole_longitude = 90 - self.grid.new_pole_longitude_degrees + mapping.grid_north_pole_latitude = 90 - self.grid.new_pole_latitude_degrees + mapping.grid_north_pole_longitude = self.grid.new_pole_longitude_degrees elif LambertConformalConic: # CRS mapping = netcdf.createVariable('Lambert_conformal', 'i') diff --git a/run_test.py b/run_test.py index 6eb2926..711c425 100644 --- a/run_test.py +++ b/run_test.py @@ -1,4 +1,22 @@ +#!/usr/bin/env python # coding=utf-8 + +# Copyright 2018 Earth Sciences Department, BSC-CNS +# +# This file is part of HERMESv3_GR. +# +# HERMESv3_GR is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# HERMESv3_GR is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with HERMESv3_GR. If not, see . """Script to run the tests for EarthDiagnostics and generate the code coverage report""" import os -- GitLab From a32554f5c6a96997146f7ff453e7886d2415f013 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Wed, 2 Jan 2019 15:19:36 +0100 Subject: [PATCH 04/10] to check --- preproc/GFAS_hourly_Parameters.csv | 27 +++ preproc/gfas12_h_preproc.py | 282 +++++++++++++++++++++++++++++ 2 files changed, 309 insertions(+) create mode 100644 preproc/GFAS_hourly_Parameters.csv create mode 100755 preproc/gfas12_h_preproc.py diff --git a/preproc/GFAS_hourly_Parameters.csv b/preproc/GFAS_hourly_Parameters.csv new file mode 100644 index 0000000..f2497eb --- /dev/null +++ b/preproc/GFAS_hourly_Parameters.csv @@ -0,0 +1,27 @@ +Name;Short_Name;Units;id +Carbon Monoxide;co;kg m-2 s-1;param81.210.192 +83;83;kg m-2 s-1;param83.210.192 +Nitrogen Oxides expressed as monoxide nitrogen;nox_no;kg m-2 s-1;param85.210.192 +Particulate Matter PM2.5;pm25;kg m-2 s-1;param87.210.192 +Organic Carbon;oc;kg m-2 s-1;param90.210.192 +Black Carbon;bc;kg m-2 s-1;param91.210.192 +Sulfur Dioxide;so2;kg m-2 s-1;param102.210.192 +Methanol;ch3oh;kg m-2 s-1;param103.210.192 +Higher Alkanes;hialkanes;kg m-2 s-1;param112.210.192 +Formaldehyde;ch2o;kg m-2 s-1;param113.210.192 +Acetaldehyde;c2h4o;kg m-2 s-1;param114.210.192 +Acetone;c3h6o;kg m-2 s-1;param115.210.192 +Ammonia;nh3;kg m-2 s-1;param116.210.192 +Dimethyl Sulfide;c2h6s;kg m-2 s-1;param117.210.192 +Ethane;c2h6;kg m-2 s-1;param118.210.192 +Toluene;c7h8;kg m-2 s-1;param231.210.192 +Benzene;c6h6;kg m-2 s-1;param232.210.192 +Xylene;c8h10;kg m-2 s-1;param233.210.192 +Butenes;c4h8;kg m-2 s-1;param234.210.192 +Pentenes;c5h10;kg m-2 s-1;param235.210.192 +Hexene;c6h12;kg m-2 s-1;param236.210.192 +Octene;c8h16;kg m-2 s-1;param237.210.192 +Butanes;c4h10;kg m-2 s-1;param238.210.192 +Pentanes;c5h12;kg m-2 s-1;param239.210.192 +Hexanes;c6h14;kg m-2 s-1;param240.210.192 +Heptane;c7h16;kg m-2 s-1;param241.210.192 diff --git a/preproc/gfas12_h_preproc.py b/preproc/gfas12_h_preproc.py new file mode 100755 index 0000000..9cc93fe --- /dev/null +++ b/preproc/gfas12_h_preproc.py @@ -0,0 +1,282 @@ +#!/usr/bin/env python + +# Copyright 2018 Earth Sciences Department, BSC-CNS +# +# This file is part of HERMESv3_GR. +# +# HERMESv3_GR is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# HERMESv3_GR is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with HERMESv3_GR. If not, see . + + +import os +from netCDF4 import Dataset +import cf_units +import pandas as pd +import datetime +from datetime import datetime, timedelta + +# ============== README ====================== +""" +downloading website: http://apps.ecmwf.int/datasets/data/cams-gfas/ +reference: https://www.biogeosciences.net/9/527/2012/ +Besides citing HERMESv3_GR, users must also acknowledge the use of the corresponding emission inventories in their works +""" + +# ============== CONFIGURATION PARAMETERS ====================== +INPUT_PATH = '/esarchive/recon/ecmwf/gfas/original_files/ga_mc_sfc_gfas_ecmf/' +INPUT_NAME = 'gfas_hourly_.grb' +OUTPUT_PATH = '/esarchive/recon/ecmwf/gfas' + +STARTING_DATE = datetime(year=2018, month=11, day=01) +ENDIND_DATE = datetime(year=2018, month=11, day=01) + +PARAMETERS_FILE = '/esarchive/recon/ecmwf/gfas/original_files/ga_mc_sfc_gfas_ecmf/GFAS_hourly_Parameters.csv' +# ============================================================== + + +def create_bounds(coords, number_vertices=2): + """ + Calculate the vertices coordinates. + + :param coords: Coordinates in degrees (latitude or longitude) + :type coords: numpy.array + + :param number_vertices: Non mandatory parameter that informs the number of vertices that must have the boundaries. + (by default 2) + :type number_vertices: int + + :return: Array with as many elements as vertices for each value of coords. + :rtype: numpy.array + """ + import numpy as np + + interval = coords[1] - coords[0] + + coords_left = coords - interval / 2 + coords_right = coords + interval / 2 + if number_vertices == 2: + bound_coords = np.dstack((coords_left, coords_right)) + elif number_vertices == 4: + bound_coords = np.dstack((coords_left, coords_right, coords_right, coords_left)) + else: + raise ValueError('The number of vertices of the boudaries must be 2 or 4') + + return bound_coords + + +def get_grid_area(filename): + """ + Calculate the area for each cell of the grid using CDO + + :param filename: Path to the file to calculate the cell area + :type filename: str + + :return: Area of each cell of the grid. + :rtype: numpy.array + """ + from cdo import Cdo + from netCDF4 import Dataset + + cdo = Cdo() + s = cdo.gridarea(input=filename) + nc_aux = Dataset(s, mode='r') + grid_area = nc_aux.variables['cell_area'][:] + nc_aux.close() + + return grid_area + + +def write_netcdf(output_name_path, data_list, center_lats, center_lons, grid_cell_area, date): + """ + Write a NetCDF with the given information. + + :param output_name_path: Complete path to the output NetCDF to be stored. + :type output_name_path: str + + :param data_list + + :param center_lats: Latitudes of the center of each cell. + :type center_lats: numpy.array + + :param center_lons: Longitudes of the center of each cell. + :type center_lons: numpy.array + + :param grid_cell_area: Area of each cell of the grid. + :type grid_cell_area: numpy.array + + :param date: Date of the current netCDF. + :type date: datetime.datetime + + """ + print output_name_path + # Creating NetCDF & Dimensions + nc_output = Dataset(output_name_path, mode='w', format="NETCDF4") + nc_output.createDimension('nv', 2) + nc_output.createDimension('lon', center_lons.shape[0]) + nc_output.createDimension('lat', center_lats.shape[0]) + nc_output.createDimension('time', None) + + # TIME + time = nc_output.createVariable('time', 'd', ('time',), zlib=True) + # time.units = "{0} since {1}".format(tstep_units, global_atts['Start_DateTime'].strftime('%Y-%m-%d %H:%M:%S')) + time.units = "hours since {0}".format(date.strftime('%Y-%m-%d %H:%M:%S')) + time.standard_name = "time" + time.calendar = "gregorian" + time.long_name = "time" + time[:] = [0] + + # LATITUDE + lat = nc_output.createVariable('lat', 'f', ('lat',), zlib=True) + lat.bounds = "lat_bnds" + lat.units = "degrees_north" + lat.axis = "Y" + lat.long_name = "latitude" + lat.standard_name = "latitude" + lat[:] = center_lats + + lat_bnds = nc_output.createVariable('lat_bnds', 'f', ('lat', 'nv',), zlib=True) + lat_bnds[:] = create_bounds(center_lats) + + # LONGITUDE + lon = nc_output.createVariable('lon', 'f', ('lon',), zlib=True) + lon.bounds = "lon_bnds" + lon.units = "degrees_east" + lon.axis = "X" + lon.long_name = "longitude" + lon.standard_name = "longitude" + lon[:] = center_lons + + lon_bnds = nc_output.createVariable('lon_bnds', 'f', ('lon', 'nv',), zlib=True) + lon_bnds[:] = create_bounds(center_lons) + + for var in data_list: + # VARIABLE + nc_var = nc_output.createVariable(var['name'], 'f', ('time', 'lat', 'lon',), zlib=True) + nc_var.units = var['units'].symbol + nc_var.long_name = var['long_name'] + nc_var.coordinates = 'lat lon' + nc_var.grid_mapping = 'crs' + nc_var.cell_measures = 'area: cell_area' + nc_var[:] = var['data'] + + # CELL AREA + cell_area = nc_output.createVariable('cell_area', 'f', ('lat', 'lon',)) + cell_area.long_name = "area of the grid cell" + cell_area.standard_name = "area" + cell_area.units = "m2" + cell_area[:] = grid_cell_area + + # CRS + crs = nc_output.createVariable('crs', 'i') + crs.grid_mapping_name = "latitude_longitude" + crs.semi_major_axis = 6371000.0 + crs.inverse_flattening = 0 + + nc_output.setncattr('title', 'GFASv1.2 inventory') + nc_output.setncattr('Conventions', 'CF-1.6', ) + nc_output.setncattr('institution', 'ECMWF', ) + nc_output.setncattr('source', 'GFAS', ) + nc_output.setncattr('history', 'Re-writing of the GFAS input to follow the CF-1.6 conventions;\n' + + '2019-01-02: Added boundaries;\n' + + '2019-01-02: Added global attributes;\n' + + '2019-01-02: Re-naming pollutant;\n' + + '2019-01-02: Added cell_area variable;\n' + + '2019-01-02: Added new varaibles;\n') + nc_output.setncattr('references', 'downloading website: http://apps.ecmwf.int/datasets/data/cams-gfas/\n' + + 'reference: https://www.biogeosciences.net/9/527/2012/', ) + nc_output.setncattr('comment', 'Re-writing done by Carles Tena (carles.tena@bsc.es) from the BSC-CNS ' + + '(Barcelona Supercomputing Center)', ) + + nc_output.close() + return True + + +def do_transformation(input_file, date, output_dir, variables_list): + """ + Transform the original file into a NEtCDF file that follows the conventions. + + :param input_file: + :type input_file: str + + :param date: Date of the file to do the transformation. + :type date: datetime.datetime + + :param output_dir: Path where have to be stored the output file. + :type output_dir: str + + :param variables_list: LIst of dictionaries with the information of each variable of the output files. + :type variables_list: list + """ + from cdo import Cdo + cdo = Cdo() + + nc_temp = cdo.copy(input=input_file, options='-R -r -f nc4c -z zip_4') + + nc_in = Dataset(nc_temp, mode='r') + + cell_area = get_grid_area(nc_temp) + + lats = nc_in.variables['lat'][:] + lons = nc_in.variables['lon'][:] + + for variable in variables_list: + variable['data'] = nc_in.variables[variable['original_name']][:] + + nc_in.close() + + out_path_aux = os.path.join(output_dir, '1hourly', 'multivar') + if not os.path.exists(out_path_aux): + os.makedirs(out_path_aux) + out_path_aux = os.path.join(out_path_aux, 'ga_{0}.nc'.format(date.strftime('%Y%m%d'))) + write_netcdf(out_path_aux, variables_list, lats, lons, cell_area, date) + + return True + + +def do_var_list(variables_file): + """ + Create the List of dictionaries + + :param variables_file: CSV file with the information of each variable + :type variables_file: str + + :return: Dictionaries list with the information of each variable. + :rtype: list + """ + df = pd.read_csv(variables_file, sep=';') + list_aux = [] + for i, element in df.iterrows(): + # print element + dict_aux = { + 'original_name': str(element.id), + 'name': element['Short_Name'], + 'long_name': element['Name'], + 'units': cf_units.Unit(element['Units']), + } + list_aux.append(dict_aux) + return list_aux + + +if __name__ == '__main__': + + var_list = do_var_list(PARAMETERS_FILE) + + date_aux = STARTING_DATE + while date_aux <= ENDIND_DATE: + f = os.path.join(INPUT_PATH, INPUT_NAME.replace('', date_aux.strftime('%Y%m%d'))) + if os.path.isfile(f): + do_transformation(f, date_aux, OUTPUT_PATH, var_list) + else: + print 'ERROR: file {0} not found'.format(f) + + date_aux = date_aux + timedelta(days=1) -- GitLab From 1cfcef1349dcfa5d06e2a0d8310d716b9f3ff04e Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Fri, 4 Jan 2019 09:57:34 +0100 Subject: [PATCH 05/10] Added param83 to csv --- preproc/GFAS_hourly_Parameters.csv | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/preproc/GFAS_hourly_Parameters.csv b/preproc/GFAS_hourly_Parameters.csv index f2497eb..7b8f414 100644 --- a/preproc/GFAS_hourly_Parameters.csv +++ b/preproc/GFAS_hourly_Parameters.csv @@ -1,6 +1,6 @@ Name;Short_Name;Units;id Carbon Monoxide;co;kg m-2 s-1;param81.210.192 -83;83;kg m-2 s-1;param83.210.192 +Non-Methane Hydro-Carbons;nmhc;kg m-2 s-1;param83.210.192 Nitrogen Oxides expressed as monoxide nitrogen;nox_no;kg m-2 s-1;param85.210.192 Particulate Matter PM2.5;pm25;kg m-2 s-1;param87.210.192 Organic Carbon;oc;kg m-2 s-1;param90.210.192 -- GitLab From d55cfa5db7d54827c1d927ce87282f45b52a1e19 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Wed, 9 Jan 2019 19:03:11 +0100 Subject: [PATCH 06/10] Doing Point distribution for GFAS --- conf/EI_configuration.csv | 25 +- conf/hermes_prace.conf | 79 +++++ .../emision_inventories/emission_inventory.py | 47 ++- .../point_gfas_emission_inventory.py | 275 ++++++++++++++++++ hermesv3_gr/modules/grids/grid.py | 3 +- hermesv3_gr/modules/grids/grid_global.py | 6 +- hermesv3_gr/modules/vertical/vertical_gfas.py | 8 + 7 files changed, 412 insertions(+), 31 deletions(-) create mode 100644 conf/hermes_prace.conf create mode 100755 hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py diff --git a/conf/EI_configuration.csv b/conf/EI_configuration.csv index 2dee5eb..a25b3a5 100644 --- a/conf/EI_configuration.csv +++ b/conf/EI_configuration.csv @@ -1,15 +1,16 @@ ei;sector;ref_year;active;factor_mask;regrid_mask;pollutants;path;frequency;source_type;p_vertical;p_month;p_day;p_hour;p_speciation;comment -GFASv12;;2015;1;;;co,nox_no,pm25,oc,bc,so2,ch3oh,c2h5oh,c3h8,c2h4,c3h6,c5h8,terpenes,hialkenes,hialkanes,ch2o,c2h4o,c3h6o,nh3,c2h6s,c2h6,c7h8,c6h6,c8h10,c4h8,c5h10,c6h12,c8h16,c4h10,c5h12,c6h14,c7h16;/ecmwf/gfas/daily_mean;daily;area;method=sovief,approach=uniform;;;H001;E001; -HTAPv2;energy;2010;1;;;co,nox_no2,pm10,pm25,oc,bc,so2,nh3,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc21,voc22,voc23,voc24;/jrc/htapv2/monthly_mean;monthly;area;V001;;D002;H002;E002; -HTAPv2;industry;2010;1;;;co,nox_no2,pm10,pm25,oc,bc,so2,nh3,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc18,voc19,voc20,voc21,voc22,voc23,voc24;/jrc/htapv2/monthly_mean;monthly;area;V002;;D003;H004;E003; -HTAPv2;residential;2010;1;;;co,nox_no2,pm10,pm25,oc,bc,so2,nh3,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc19,voc21,voc22,voc23,voc24;/jrc/htapv2/monthly_mean;monthly;area;;;D003;H003;E004; -HTAPv2;transport;2010;1;;;co,nox_no2,pm10,pm25,oc,bc,so2,nh3,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc21,voc22,voc23;/jrc/htapv2/monthly_mean;monthly;area;;;D005;weekday=H006, saturday=H009, sunday=H010;E005; -HTAPv2;agriculture;2010;1;;;nh3;/jrc/htapv2/monthly_mean;monthly;area;;;D001;H007;E006; -HTAPv2;air_lto;2010;1;;;co,nox_no2,pm10,pm25,oc,bc,so2,voc02,voc03,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc17,voc21,voc22,voc23;/jrc/htapv2/yearly_mean;yearly;area;V003;M001;D001;H001;E007; -HTAPv2;air_cds;2010;1;;;co,nox_no2,pm10,pm25,oc,bc,so2,voc02,voc03,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc17,voc21,voc22,voc23;/jrc/htapv2/yearly_mean;yearly;area;V004;M001;D001;H001;E007; -HTAPv2;air_crs;2010;1;;;co,nox_no2,pm10,pm25,oc,bc,so2,voc02,voc03,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc17,voc21,voc22,voc23;/jrc/htapv2/yearly_mean;yearly;area;V005;M001;D001;H001;E007; -HTAPv2;ships;2010;1;;;co,nox_no2,pm10,pm25,oc,bc,so2,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc18,voc19,voc21,voc22,voc23,voc24;/jrc/htapv2/yearly_mean;yearly;area;;M001;D001;H001;E008; -wiedinmyer;;2010;1;;;bc,c2h2,c2h4,c3h6,c6h6,ch2o,ch3cooh,ch3oh,co,hcl,nh3,nox_no,oc,pm10,pm25,so2;/ucar/wiedinmyer/yearly_mean;yearly;area;;M001;D001;H008;E009; +GFASv12;;2015;1;;;co;/ecmwf/gfas/daily_mean;daily;point;method=sovief,approach=uniform;;;;E001; +GFASv12;;2015;0;;;co,nox_no,pm25,oc,bc,so2,ch3oh,c2h5oh,c3h8,c2h4,c3h6,c5h8,terpenes,hialkenes,hialkanes,ch2o,c2h4o,c3h6o,nh3,c2h6s,c2h6,c7h8,c6h6,c8h10,c4h8,c5h10,c6h12,c8h16,c4h10,c5h12,c6h14,c7h16;/ecmwf/gfas/daily_mean;daily;area;method=sovief,approach=uniform;;;H001;E001; +HTAPv2;energy;2010;0;;;co,nox_no2,pm10,pm25,oc,bc,so2,nh3,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc21,voc22,voc23,voc24;/jrc/htapv2/monthly_mean;monthly;area;V001;;D002;H002;E002; +HTAPv2;industry;2010;0;;;co,nox_no2,pm10,pm25,oc,bc,so2,nh3,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc18,voc19,voc20,voc21,voc22,voc23,voc24;/jrc/htapv2/monthly_mean;monthly;area;V002;;D003;H004;E003; +HTAPv2;residential;2010;0;;;co,nox_no2,pm10,pm25,oc,bc,so2,nh3,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc19,voc21,voc22,voc23,voc24;/jrc/htapv2/monthly_mean;monthly;area;;;D003;H003;E004; +HTAPv2;transport;2010;0;;;co,nox_no2,pm10,pm25,oc,bc,so2,nh3,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc21,voc22,voc23;/jrc/htapv2/monthly_mean;monthly;area;;;D005;weekday=H006, saturday=H009, sunday=H010;E005; +HTAPv2;agriculture;2010;0;;;nh3;/jrc/htapv2/monthly_mean;monthly;area;;;D001;H007;E006; +HTAPv2;air_lto;2010;0;;;co,nox_no2,pm10,pm25,oc,bc,so2,voc02,voc03,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc17,voc21,voc22,voc23;/jrc/htapv2/yearly_mean;yearly;area;V003;M001;D001;H001;E007; +HTAPv2;air_cds;2010;0;;;co,nox_no2,pm10,pm25,oc,bc,so2,voc02,voc03,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc17,voc21,voc22,voc23;/jrc/htapv2/yearly_mean;yearly;area;V004;M001;D001;H001;E007; +HTAPv2;air_crs;2010;0;;;co,nox_no2,pm10,pm25,oc,bc,so2,voc02,voc03,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc17,voc21,voc22,voc23;/jrc/htapv2/yearly_mean;yearly;area;V005;M001;D001;H001;E007; +HTAPv2;ships;2010;0;;;co,nox_no2,pm10,pm25,oc,bc,so2,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc18,voc19,voc21,voc22,voc23,voc24;/jrc/htapv2/yearly_mean;yearly;area;;M001;D001;H001;E008; +wiedinmyer;;2010;0;;;bc,c2h2,c2h4,c3h6,c6h6,ch2o,ch3cooh,ch3oh,co,hcl,nh3,nox_no,oc,pm10,pm25,so2;/ucar/wiedinmyer/yearly_mean;yearly;area;;M001;D001;H008;E009; TNO_MACC-III;snap1;2011;0;;;co,nox_no2,so2,nh3,pm10,pm25,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc19,voc21,voc22,voc23,voc24;/tno/tno_macc_iii/yearly_mean/;yearly;area;V001;M002;D002;H002;E010; TNO_MACC-III;snap2;2011;0;;;co,nox_no2,so2,nh3,pm10,pm25,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc19,voc21,voc22,voc23,voc24;/tno/tno_macc_iii/yearly_mean/;yearly;area;;M003;D003;H003;E011; TNO_MACC-III;snap34;2011;0;;;co,nox_no2,so2,nh3,pm10,pm25,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc19,voc21,voc22,voc23,voc24;/tno/tno_macc_iii/yearly_mean/;yearly;area;V002;M004;D003;H004;E012; @@ -86,7 +87,7 @@ EMEP;i_offroad;2015;0;;;co,nox_no2,pm10,pm25,so2,nmvoc,nh3;/ceip/emep EMEP;j_waste;2015;0;;;co,nox_no2,pm10,pm25,so2,nmvoc,nh3;/ceip/emepv18/yearly_mean;yearly;area;;M001;D001;H001;E083; EMEP;k_agrilivestock;2015;0;;;co,nox_no2,pm10,pm25,so2,nmvoc,nh3;/ceip/emepv18/yearly_mean;yearly;area;;M009;D001;H007;E084; EMEP;l_agriother;2015;0;;;co,nox_no2,pm10,pm25,so2,nmvoc,nh3;/ceip/emepv18/yearly_mean;yearly;area;;M009;D001;H007;E085; -carn;;2015;1;;;so2;/mtu/carnetal/yearly_mean;yearly;point;;M001;D001;H001;E086; +carn;;2015;0;;;so2;/mtu/carnetal/yearly_mean;yearly;point;;M001;D001;H001;E086; CEDS;agriculture;2014;0;;;nox_no2,nh3;/jgcri/ceds/monthly_mean;monthly;area;;;D001;H007;E087; CEDS;energy;2014;0;;;co,nox_no2,oc,bc,so2,nh3,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc18,voc19,voc20,voc21,voc22,voc23,voc24;/jgcri/ceds/monthly_mean;monthly;area;V001;;D002;H002;E088; CEDS;industrial;2014;0;;;co,nox_no2,oc,bc,so2,nh3,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc18,voc19,voc20,voc21,voc22,voc23,voc24;/jgcri/ceds/monthly_mean;monthly;area;V002;;D003;H004;E089; diff --git a/conf/hermes_prace.conf b/conf/hermes_prace.conf new file mode 100644 index 0000000..79ddc80 --- /dev/null +++ b/conf/hermes_prace.conf @@ -0,0 +1,79 @@ +[GENERAL] +log_level = 3 +input_dir = /gpfs/projects/pr1ehu00/pr1ehu03/hermesv3_gr +data_path = /gpfs/scratch/bsc32/bsc32538/HERMES_data +output_dir = /gpfs/projects/pr1ehu00/pr1ehu03/hermesv3_gr_OUT +output_name = HERMESv3_GFAS_.nc +start_date = 2016/05/01 00:00:00 +# ***** end_date = start_date [DEFAULT] ***** +# end_date = 2017/01/02 00:00:00 +# ***** output_timestep_type = [hourly, daily, monthly, yearly] ***** +output_timestep_type = hourly +output_timestep_num = 24 +output_timestep_freq = 1 + +[DOMAIN] + +# ***** output_model = [MONARCH, CMAQ, WRF_CHEM] ***** +output_model = MONARCH +# output_model = CMAQ +# output_model = WRF_CHEM +output_attributes = /data/global_attributes.csv + +# ***** domain_type=[global, lcc, rotated, mercator] ***** +domain_type = global +# domain_type = lcc +# domain_type = rotated +# domain_type = mercator +vertical_description = /data/profiles/vertical/Benchmark_15layers_vertical_description.csv +auxiliar_files_path = /data/auxiliar_files/_ + +# if domain_type == global: + inc_lat = 1. + inc_lon = 1.40625 + +# if domain_type == rotated: + #centre_lat = 35 + #centre_lon = 20 + #west_boundary = -51 + #south_boundary = -35 + #inc_rlat = 0.1 + #inc_rlon = 0.1 + +# if domain_type == lcc: + #lat_1 = 37 + #lat_2 = 43 + #lon_0 = -3 + #lat_0 = 40 + #nx = 478 + #ny = 398 + #inc_x = 12000 + #inc_y = 12000 + #x_0 = -2131849.000 + #y_0 = -2073137.875 + +# if domain_type == mercator: + #lat_ts = -1.5 + #lon_0 = -18 + #nx = 210 + #ny = 236 + #inc_x = 50000 + #inc_y = 50000 + #x_0 = -126017.5 + #y_0 = -5407460 + + +[EMISSION_INVENTORY_CONFIGURATION] + +cross_table = /conf/EI_configuration.csv + +[EMISSION_INVENTORY_PROFILES] + +p_vertical = /data/profiles/vertical/Vertical_profile.csv +p_month = /data/profiles/temporal/TemporalProfile_Monthly.csv +p_day = /data/profiles/temporal/TemporalProfile_Daily.csv +p_hour = /data/profiles/temporal/TemporalProfile_Hourly.csv +p_speciation = /data/profiles/speciation/Speciation_profile_cb05_aero5_MONARCH.csv + +molecular_weights = /data/profiles/speciation/MolecularWeights.csv +world_info = /data/profiles/temporal/tz_world_country_iso3166.csv diff --git a/hermesv3_gr/modules/emision_inventories/emission_inventory.py b/hermesv3_gr/modules/emision_inventories/emission_inventory.py index 3101a2f..89eee3a 100644 --- a/hermesv3_gr/modules/emision_inventories/emission_inventory.py +++ b/hermesv3_gr/modules/emision_inventories/emission_inventory.py @@ -207,8 +207,12 @@ class EmissionInventory(object): if self.source_type == 'area': extension = 'nc' + elif self.source_type == 'point': - extension = 'csv' + if self.inventory_name == 'GFASv12': + extension = 'nc' + else: + extension = 'csv' else: settings.write_log('ERROR: Check the .err file to get more info.') if settings.rank == 0: @@ -227,8 +231,9 @@ class EmissionInventory(object): elif self.input_frequency == 'monthly': file_name = "{0}_{1}{2}.{3}".format(pollutant, self.reference_year, self.date.strftime("%m"), extension) elif self.input_frequency == 'daily': + print pollutant, self.reference_year, self.date.strftime("%m"), self.date.strftime("%d"), extension file_name = "{0}_{1}{2}{3}.{4}".format(pollutant, self.reference_year, self.date.strftime("%m"), - self.date.strftime("%D"), extension) + self.date.strftime("%d"), extension) else: settings.write_log('ERROR: Check the .err file to get more info.') if settings.rank == 0: @@ -281,6 +286,7 @@ class EmissionInventory(object): """ import pandas as pd import re + from point_gfas_emission_inventory import PointGfasEmissionInventory from gfas_emission_inventory import GfasEmissionInventory from point_source_emission_inventory import PointSourceEmissionInventory @@ -345,19 +351,30 @@ class EmissionInventory(object): p_hour=emission_inventory.p_hour, p_speciation=emission_inventory.p_speciation)) elif emission_inventory.source_type == 'point': - emission_inventory_list.append( - PointSourceEmissionInventory(options, grid, date, emission_inventory.ei, - emission_inventory.source_type, emission_inventory.sector, pollutants, - emission_inventory_path, - emission_inventory.frequency, vertical_output_profile, - reference_year=emission_inventory.ref_year, - factors=emission_inventory.factor_mask, - regrid_mask=emission_inventory.regrid_mask, - p_vertical=emission_inventory.p_vertical, - p_month=p_month, - p_day=emission_inventory.p_day, - p_hour=emission_inventory.p_hour, - p_speciation=emission_inventory.p_speciation)) + if emission_inventory.ei == 'GFASv12': + emission_inventory_list.append( + PointGfasEmissionInventory( + options, grid, date, emission_inventory.ei, emission_inventory.source_type, + emission_inventory.sector, pollutants, emission_inventory_path, + emission_inventory.frequency, vertical_output_profile, + reference_year=emission_inventory.ref_year, factors=emission_inventory.factor_mask, + regrid_mask=emission_inventory.regrid_mask, p_vertical=emission_inventory.p_vertical, + p_month=p_month, p_day=emission_inventory.p_day, p_hour=emission_inventory.p_hour, + p_speciation=emission_inventory.p_speciation)) + else: + emission_inventory_list.append( + PointSourceEmissionInventory(options, grid, date, emission_inventory.ei, + emission_inventory.source_type, emission_inventory.sector, pollutants, + emission_inventory_path, + emission_inventory.frequency, vertical_output_profile, + reference_year=emission_inventory.ref_year, + factors=emission_inventory.factor_mask, + regrid_mask=emission_inventory.regrid_mask, + p_vertical=emission_inventory.p_vertical, + p_month=p_month, + p_day=emission_inventory.p_day, + p_hour=emission_inventory.p_hour, + p_speciation=emission_inventory.p_speciation)) else: settings.write_log('ERROR: Check the .err file to get more info.') if settings.rank == 0: diff --git a/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py b/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py new file mode 100755 index 0000000..248b5af --- /dev/null +++ b/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py @@ -0,0 +1,275 @@ +#!/usr/bin/env python + +# Copyright 2018 Earth Sciences Department, BSC-CNS +# +# This file is part of HERMESv3_GR. +# +# HERMESv3_GR is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# HERMESv3_GR is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with HERMESv3_GR. If not, see . + +import sys +import os +import timeit + +import hermesv3_gr.config.settings as settings +from emission_inventory import EmissionInventory + + +class PointGfasEmissionInventory(EmissionInventory): + """ + Class that defines the content and the methodology for the GFAS emission inventories + + :param current_date: Date of the day to simulate. + :type current_date: datetime.datetime + + :param inventory_name: Name of the inventory to use. + :type inventory_name: str + + :param sector: Name of the sector of the inventory to use. + :type sector: str + + :param pollutants: List of the pollutant name to take into account. + :type pollutants: list of str + + :param frequency: Frequency of the inputs. [yearly, monthly, daily] + :type frequency: str + + :param reference_year: year of reference of the information of the dataset. + :type reference_year: int + + :param factors: NOT IMPLEMENTED YET + :type factors: NOT IMPLEMENTED YET + + :param regrid_mask: NOT IMPLEMENTED YET + :type regrid_mask: NOT IMPLEMENTED YET + + :param p_vertical: ID of the vertical profile to use. + :type p_vertical: str + + :param p_month: ID of the temporal monthly profile to use. + :type p_month: str + + :param p_day: ID of the temporal daily profile to use. + :type p_day: str + + :param p_hour: ID of the temporal hourly profile to use. + :type p_hour: str + + :param p_speciation: ID of the speciation profile to use. + :type p_speciation: str + """ + + def __init__(self, options, grid, current_date, inventory_name, source_type, sector, pollutants, inputs_path, + frequency, vertical_output_profile, + reference_year=2010, factors=None, regrid_mask=None, p_vertical=None, p_month=None, p_day=None, + p_hour=None, p_speciation=None): + from hermesv3_gr.modules.vertical.vertical_gfas import GfasVerticalDistribution + + st_time = timeit.default_timer() + settings.write_log('\t\tCreating GFAS emission inventory as point source.', level=3) + super(PointGfasEmissionInventory, self).__init__( + options, grid, current_date, inventory_name, source_type, sector, pollutants, inputs_path, frequency, + vertical_output_profile, + reference_year=reference_year, factors=factors, regrid_mask=regrid_mask, p_vertical=None, + p_month=p_month, p_day=p_day, p_hour=p_hour, p_speciation=p_speciation) + + # self.approach = self.get_approach(p_vertical) + self.method = self.get_method(p_vertical) + + # self.altitude = self.get_altitude() + + self.vertical = GfasVerticalDistribution(vertical_output_profile, self.get_approach(p_vertical), self.get_altitude()) + + settings.write_time('PointGfasEmissionInventory', 'Init', timeit.default_timer() - st_time, level=3) + + def __str__(self): + string = "PointGfasEmissionInventory:\n" + string += "\t self.method = {0}\n".format(self.method) + string += "\t self.vertical = {0}\n".format(self.vertical) + + return string + + def get_input_path(self, pollutant=None, extension='nc'): + """ + Completes the path of the NetCDF that contains the needed information of the given pollutant. + + :param pollutant: Name of the pollutant of the NetCDF. + :type pollutant: str + + :param extension: Extension of the input file. + :type: str + + :return: Full path of the needed NetCDF. + :rtype: str + """ + st_time = timeit.default_timer() + + file_path = os.path.join(self.inputs_path, 'multivar', 'ga_{0}.{1}'.format( + self.date.strftime('%Y%m%d'), extension)) + + # Checking input file + if not os.path.exists(file_path): + settings.write_log('ERROR: Check the .err file to get more info.') + if settings.rank == 0: + raise IOError('ERROR: File {0} not found.'.format(file_path)) + sys.exit(1) + + settings.write_time('PointGfasEmissionInventory', 'get_input_path', timeit.default_timer() - st_time, level=3) + + return file_path + + def get_altitude(self): + """ + Extract the altitude values depending on the choosen method. + + :return: Array with the alittude of each fire. + :rtype: numpy.array + """ + from hermesv3_gr.tools.netcdf_tools import extract_vars + + st_time = timeit.default_timer() + + if self.method == 'sovief': + alt_var = 'apt' + elif self.method == 'prm': + alt_var = 'mami' + else: + alt_var = None + + print "ERROR: Only 'sovief' and 'prm' methods are accepted." + + [alt] = extract_vars(self.get_input_path(), [alt_var]) + + alt = alt['data'] + + settings.write_time('PointGfasEmissionInventory', 'get_altitude', timeit.default_timer() - st_time, level=3) + return alt + + @ staticmethod + def get_approach(p_vertical): + """ + Extract the given approach value. + + :return: Approach value + :rtype: str + """ + import re + + st_time = timeit.default_timer() + + return_value = None + aux_list = re.split(', |,| , | ,', p_vertical) + for element in aux_list: + aux_value = re.split('=| =|= | = ', element) + if aux_value[0] == 'approach': + return_value = aux_value[1] + + settings.write_time('PointGfasEmissionInventory', 'get_approach', timeit.default_timer() - st_time, level=3) + + return return_value + + @ staticmethod + def get_method(p_vertical): + """ + Extract the given method value. + + :return: Method value + :rtype: str + """ + import re + + st_time = timeit.default_timer() + + return_value = None + aux_list = re.split(', |,| , | ,', p_vertical) + for element in aux_list: + aux_value = re.split('=| =|= | = ', element) + if aux_value[0] == 'method': + return_value = aux_value[1] + + settings.write_time('PointGfasEmissionInventory', 'get_method', timeit.default_timer() - st_time, level=3) + + return return_value + + def do_vertical_allocation(self, values): + """ + Allocates the fire emissions on their top level. + + :param values: 2D array with the fire emissions + :type values: numpy.array + + :return: Emissions already allocated on the top altitude of each fire. + :rtype: numpy.array + """ + print 'do_vertical_allocation' + sys.exit() + st_time = timeit.default_timer() + + return_value = self.vertical.do_vertical_interpolation_allocation(values, self.altitude) + + settings.write_time('PointGfasEmissionInventory', 'do_vertical_allocation', timeit.default_timer() - st_time, + level=3) + + return return_value + + def do_regrid(self): + import pandas as pd + import geopandas as gpd + from shapely.geometry import Point + import numpy as np + from netCDF4 import Dataset + + st_time = timeit.default_timer() + settings.write_log("\tAllocating GFAS as point sources on grid:", level=2) + + + netcdf = Dataset(self.pollutant_dicts[0]['path'], mode='r') + if settings.rank == 0: + gdf = gpd.GeoDataFrame(crs={'init': 'epsg:4326'}) + + gdf['src_index'] = np.where(self.vertical.altitude.flatten() > 0)[0] + gdf['altitude'] = self.vertical.altitude.flatten()[gdf['src_index']] + + lat_1d = netcdf.variables['lat'][:] + lon_1d = netcdf.variables['lon'][:] + # 1D to 2D + lats = np.array([lat_1d] * len(lon_1d)).T + lons = np.array([lon_1d] * len(lat_1d)) + + gdf['geometry'] = [Point(xy) for xy in zip(lons.flatten()[gdf['src_index']], + lats.flatten()[gdf['src_index']])] + grid_shp = self.grid.to_shapefile() + + gdf = gpd.sjoin(gdf, grid_shp.to_crs(gdf.crs), how='left') + + gdf = np.array_split(gdf, settings.size) + else: + gdf = None + + gdf = settings.comm.scatter(gdf, root=0) + + for num, pollutant in enumerate(self.pollutant_dicts): + settings.write_log('\t\tPollutant {0} ({1}/{2})'.format( + pollutant['name'], num + 1, len(self.pollutant_dicts)), level=3) + print ('\t\tPollutant {0} ({1}/{2})'.format(pollutant['name'], num + 1, len(self.pollutant_dicts))) + gdf[pollutant['name']] = netcdf.variables[pollutant['name']][:].flatten()[gdf['src_index']] + print gdf['FID'] + + netcdf.close() + self.vertical = None + + settings.write_time('PointGfasEmissionInventory', 'do_regrid', timeit.default_timer() - st_time, level=2) + print 'regrid done' + sys.exit() +if __name__ == "__main__": + pass diff --git a/hermesv3_gr/modules/grids/grid.py b/hermesv3_gr/modules/grids/grid.py index 0c424d1..27a173a 100644 --- a/hermesv3_gr/modules/grids/grid.py +++ b/hermesv3_gr/modules/grids/grid.py @@ -250,7 +250,8 @@ class Grid(object): regular_latlon=True) # Calculate the cell area of the auxiliary NetCDF file - self.cell_area = self.get_cell_area() + # TODO restore + # self.cell_area = self.get_cell_area() # Re-writes the NetCDF adding the cell area write_netcdf(self.coords_netcdf_file, self.center_latitudes, self.center_longitudes, diff --git a/hermesv3_gr/modules/grids/grid_global.py b/hermesv3_gr/modules/grids/grid_global.py index ebbd97e..b00f654 100644 --- a/hermesv3_gr/modules/grids/grid_global.py +++ b/hermesv3_gr/modules/grids/grid_global.py @@ -85,9 +85,9 @@ class GlobalGrid(Grid): self.shape = (timestep_num, len(self.vertical_description), self.x_upper_bound-self.x_lower_bound, self.y_upper_bound-self.y_lower_bound) - - self.cell_area = self.get_cell_area()[self.x_lower_bound:self.x_upper_bound, - self.y_lower_bound:self.y_upper_bound] + # TODO restore + # self.cell_area = self.get_cell_area()[self.x_lower_bound:self.x_upper_bound, + # self.y_lower_bound:self.y_upper_bound] settings.write_time('GlobalGrid', 'Init', timeit.default_timer() - st_time, level=1) diff --git a/hermesv3_gr/modules/vertical/vertical_gfas.py b/hermesv3_gr/modules/vertical/vertical_gfas.py index 600ab7f..9c17a36 100644 --- a/hermesv3_gr/modules/vertical/vertical_gfas.py +++ b/hermesv3_gr/modules/vertical/vertical_gfas.py @@ -42,6 +42,14 @@ class GfasVerticalDistribution(VerticalDistribution): settings.write_time('GfasVerticalDistribution', 'Init', timeit.default_timer() - st_time, level=3) + def __str__(self): + string = "GFAS Vertical distribution:\n" + string += "\t self.approach = {0}\n".format(self.approach) + string += "\t self.output_heights = {0}\n".format(self.output_heights) + string += "\t self.altitude = {0}\n".format(self.altitude) + + return string + @staticmethod def calculate_widths(heights_list): """ -- GitLab From fe648a9f0f7e8efe9837125f761ef0fdf86650e7 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Fri, 11 Jan 2019 15:32:47 +0100 Subject: [PATCH 07/10] Updated GFAS csv preproc for the new variables. to check --- .../emision_inventories/emission_inventory.py | 1 + .../point_gfas_emission_inventory.py | 56 +++++++++++++++++-- hermesv3_gr/modules/grids/grid.py | 4 +- hermesv3_gr/modules/grids/grid_global.py | 6 +- hermesv3_gr/modules/vertical/vertical_gfas.py | 10 +++- 5 files changed, 67 insertions(+), 10 deletions(-) diff --git a/hermesv3_gr/modules/emision_inventories/emission_inventory.py b/hermesv3_gr/modules/emision_inventories/emission_inventory.py index 89eee3a..29e4551 100644 --- a/hermesv3_gr/modules/emision_inventories/emission_inventory.py +++ b/hermesv3_gr/modules/emision_inventories/emission_inventory.py @@ -106,6 +106,7 @@ class EmissionInventory(object): world_mask_file=os.path.join(os.path.dirname(options.auxiliar_files_path), '{0}_WorldMask.nc'.format(inventory_name))) + self.input_pollutants = pollutants self.pollutant_dicts = self.create_pollutants_dicts(pollutants) self.masking.check_regrid_mask(self.pollutant_dicts[0]['path']) diff --git a/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py b/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py index 248b5af..019a135 100755 --- a/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py +++ b/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py @@ -232,7 +232,6 @@ class PointGfasEmissionInventory(EmissionInventory): st_time = timeit.default_timer() settings.write_log("\tAllocating GFAS as point sources on grid:", level=2) - netcdf = Dataset(self.pollutant_dicts[0]['path'], mode='r') if settings.rank == 0: gdf = gpd.GeoDataFrame(crs={'init': 'epsg:4326'}) @@ -242,6 +241,7 @@ class PointGfasEmissionInventory(EmissionInventory): lat_1d = netcdf.variables['lat'][:] lon_1d = netcdf.variables['lon'][:] + lon_1d[lon_1d > 180] -= 360 # 1D to 2D lats = np.array([lat_1d] * len(lon_1d)).T lons = np.array([lon_1d] * len(lat_1d)) @@ -250,7 +250,7 @@ class PointGfasEmissionInventory(EmissionInventory): lats.flatten()[gdf['src_index']])] grid_shp = self.grid.to_shapefile() - gdf = gpd.sjoin(gdf, grid_shp.to_crs(gdf.crs), how='left') + gdf = gpd.sjoin(gdf.to_crs(grid_shp.crs), grid_shp, how='left') gdf = np.array_split(gdf, settings.size) else: @@ -263,13 +263,61 @@ class PointGfasEmissionInventory(EmissionInventory): pollutant['name'], num + 1, len(self.pollutant_dicts)), level=3) print ('\t\tPollutant {0} ({1}/{2})'.format(pollutant['name'], num + 1, len(self.pollutant_dicts))) gdf[pollutant['name']] = netcdf.variables[pollutant['name']][:].flatten()[gdf['src_index']] - print gdf['FID'] netcdf.close() - self.vertical = None settings.write_time('PointGfasEmissionInventory', 'do_regrid', timeit.default_timer() - st_time, level=2) print 'regrid done' + del gdf['src_index'], gdf['index_right'] + self.emissions = gdf + + return True + + def calculate_altitudes(self, vertical_description_path): + """ + Calculate the number layer to allocate the point source. + + :param vertical_description_path: Path to the file that contains the vertical description + :type vertical_description_path: str + + :return: True + :rtype: bool + """ + import pandas as pd + + st_time = timeit.default_timer() + settings.write_log("\t\tCalculating vertical allocation.", level=3) + df = pd.read_csv(vertical_description_path, sep=';') + + self.emissions['layer'] = None + for i, line in df.iterrows(): + self.emissions.loc[self.emissions['altitude'] <= line['height_magl'], 'layer'] = line['Ilayer'] - 1 + self.emissions.loc[self.emissions['altitude'] <= line['height_magl'], 'altitude'] = None + del self.emissions['altitude'] + + self.emissions = self.emissions.groupby(['FID', 'layer']).sum() + self.emissions.reset_index(inplace=True) + + self.emissions = self.vertical.distribute_vertically(self.emissions) + + settings.write_time('PointGfasEmissionInventory', 'calculate_altitudes', timeit.default_timer() - st_time, + level=2) + print 'calculate_altitudes' + sys.exit() + return True + + def point_source_by_cell(self): + """ + Sums the different emissions that are allocated in the same cell and layer. + + :return: None + """ + st_time = timeit.default_timer() sys.exit() + settings.write_time('PointGfasEmissionInventory', 'Init', timeit.default_timer() - st_time, level=3) + + return None + + if __name__ == "__main__": pass diff --git a/hermesv3_gr/modules/grids/grid.py b/hermesv3_gr/modules/grids/grid.py index 27a173a..fe15d96 100644 --- a/hermesv3_gr/modules/grids/grid.py +++ b/hermesv3_gr/modules/grids/grid.py @@ -250,8 +250,8 @@ class Grid(object): regular_latlon=True) # Calculate the cell area of the auxiliary NetCDF file - # TODO restore - # self.cell_area = self.get_cell_area() + + self.cell_area = self.get_cell_area() # Re-writes the NetCDF adding the cell area write_netcdf(self.coords_netcdf_file, self.center_latitudes, self.center_longitudes, diff --git a/hermesv3_gr/modules/grids/grid_global.py b/hermesv3_gr/modules/grids/grid_global.py index b00f654..ebbd97e 100644 --- a/hermesv3_gr/modules/grids/grid_global.py +++ b/hermesv3_gr/modules/grids/grid_global.py @@ -85,9 +85,9 @@ class GlobalGrid(Grid): self.shape = (timestep_num, len(self.vertical_description), self.x_upper_bound-self.x_lower_bound, self.y_upper_bound-self.y_lower_bound) - # TODO restore - # self.cell_area = self.get_cell_area()[self.x_lower_bound:self.x_upper_bound, - # self.y_lower_bound:self.y_upper_bound] + + self.cell_area = self.get_cell_area()[self.x_lower_bound:self.x_upper_bound, + self.y_lower_bound:self.y_upper_bound] settings.write_time('GlobalGrid', 'Init', timeit.default_timer() - st_time, level=1) diff --git a/hermesv3_gr/modules/vertical/vertical_gfas.py b/hermesv3_gr/modules/vertical/vertical_gfas.py index 9c17a36..dd78175 100644 --- a/hermesv3_gr/modules/vertical/vertical_gfas.py +++ b/hermesv3_gr/modules/vertical/vertical_gfas.py @@ -17,7 +17,7 @@ # You should have received a copy of the GNU General Public License # along with HERMESv3_GR. If not, see . - +import sys import timeit import hermesv3_gr.config.settings as settings from vertical import VerticalDistribution @@ -180,6 +180,14 @@ class GfasVerticalDistribution(VerticalDistribution): level=3) return fire_list + def distribute_vertically(self, emissions_df): + print emissions_df + + vert_emissions = [] + + for layer, emis in emissions_df.groupby('layer'): + print layer + if __name__ == '__main__': pass -- GitLab From 4e949b0c1b844365d4fb30c978c40fdf8d2839e3 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Mon, 14 Jan 2019 18:17:26 +0100 Subject: [PATCH 08/10] GFAS as point source emission done. also implemented the Regrid Mask and the factor mask options (working in serial and in parallel) --- conf/EI_configuration.csv | 25 ++-- .../point_gfas_emission_inventory.py | 138 ++++++++++++++++-- .../point_source_emission_inventory.py | 2 +- hermesv3_gr/modules/vertical/vertical_gfas.py | 43 +++++- 4 files changed, 177 insertions(+), 31 deletions(-) diff --git a/conf/EI_configuration.csv b/conf/EI_configuration.csv index a25b3a5..2d33a26 100644 --- a/conf/EI_configuration.csv +++ b/conf/EI_configuration.csv @@ -1,16 +1,15 @@ ei;sector;ref_year;active;factor_mask;regrid_mask;pollutants;path;frequency;source_type;p_vertical;p_month;p_day;p_hour;p_speciation;comment -GFASv12;;2015;1;;;co;/ecmwf/gfas/daily_mean;daily;point;method=sovief,approach=uniform;;;;E001; -GFASv12;;2015;0;;;co,nox_no,pm25,oc,bc,so2,ch3oh,c2h5oh,c3h8,c2h4,c3h6,c5h8,terpenes,hialkenes,hialkanes,ch2o,c2h4o,c3h6o,nh3,c2h6s,c2h6,c7h8,c6h6,c8h10,c4h8,c5h10,c6h12,c8h16,c4h10,c5h12,c6h14,c7h16;/ecmwf/gfas/daily_mean;daily;area;method=sovief,approach=uniform;;;H001;E001; -HTAPv2;energy;2010;0;;;co,nox_no2,pm10,pm25,oc,bc,so2,nh3,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc21,voc22,voc23,voc24;/jrc/htapv2/monthly_mean;monthly;area;V001;;D002;H002;E002; -HTAPv2;industry;2010;0;;;co,nox_no2,pm10,pm25,oc,bc,so2,nh3,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc18,voc19,voc20,voc21,voc22,voc23,voc24;/jrc/htapv2/monthly_mean;monthly;area;V002;;D003;H004;E003; -HTAPv2;residential;2010;0;;;co,nox_no2,pm10,pm25,oc,bc,so2,nh3,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc19,voc21,voc22,voc23,voc24;/jrc/htapv2/monthly_mean;monthly;area;;;D003;H003;E004; -HTAPv2;transport;2010;0;;;co,nox_no2,pm10,pm25,oc,bc,so2,nh3,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc21,voc22,voc23;/jrc/htapv2/monthly_mean;monthly;area;;;D005;weekday=H006, saturday=H009, sunday=H010;E005; -HTAPv2;agriculture;2010;0;;;nh3;/jrc/htapv2/monthly_mean;monthly;area;;;D001;H007;E006; -HTAPv2;air_lto;2010;0;;;co,nox_no2,pm10,pm25,oc,bc,so2,voc02,voc03,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc17,voc21,voc22,voc23;/jrc/htapv2/yearly_mean;yearly;area;V003;M001;D001;H001;E007; -HTAPv2;air_cds;2010;0;;;co,nox_no2,pm10,pm25,oc,bc,so2,voc02,voc03,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc17,voc21,voc22,voc23;/jrc/htapv2/yearly_mean;yearly;area;V004;M001;D001;H001;E007; -HTAPv2;air_crs;2010;0;;;co,nox_no2,pm10,pm25,oc,bc,so2,voc02,voc03,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc17,voc21,voc22,voc23;/jrc/htapv2/yearly_mean;yearly;area;V005;M001;D001;H001;E007; -HTAPv2;ships;2010;0;;;co,nox_no2,pm10,pm25,oc,bc,so2,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc18,voc19,voc21,voc22,voc23,voc24;/jrc/htapv2/yearly_mean;yearly;area;;M001;D001;H001;E008; -wiedinmyer;;2010;0;;;bc,c2h2,c2h4,c3h6,c6h6,ch2o,ch3cooh,ch3oh,co,hcl,nh3,nox_no,oc,pm10,pm25,so2;/ucar/wiedinmyer/yearly_mean;yearly;area;;M001;D001;H008;E009; +GFASv12;;2015;1;;;co,nox_no,pm25,oc,bc,so2,ch3oh,c2h5oh,c3h8,c2h4,c3h6,c5h8,terpenes,hialkenes,hialkanes,ch2o,c2h4o,c3h6o,nh3,c2h6s,c2h6,c7h8,c6h6,c8h10,c4h8,c5h10,c6h12,c8h16,c4h10,c5h12,c6h14,c7h16;/ecmwf/gfas/daily_mean;daily;point;method=sovief,approach=uniform;;;H001;E001; +HTAPv2;energy;2010;1;;;co,nox_no2,pm10,pm25,oc,bc,so2,nh3,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc21,voc22,voc23,voc24;/jrc/htapv2/monthly_mean;monthly;area;V001;;D002;H002;E002; +HTAPv2;industry;2010;1;;;co,nox_no2,pm10,pm25,oc,bc,so2,nh3,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc18,voc19,voc20,voc21,voc22,voc23,voc24;/jrc/htapv2/monthly_mean;monthly;area;V002;;D003;H004;E003; +HTAPv2;residential;2010;1;;;co,nox_no2,pm10,pm25,oc,bc,so2,nh3,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc19,voc21,voc22,voc23,voc24;/jrc/htapv2/monthly_mean;monthly;area;;;D003;H003;E004; +HTAPv2;transport;2010;1;;;co,nox_no2,pm10,pm25,oc,bc,so2,nh3,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc21,voc22,voc23;/jrc/htapv2/monthly_mean;monthly;area;;;D005;weekday=H006, saturday=H009, sunday=H010;E005; +HTAPv2;agriculture;2010;1;;;nh3;/jrc/htapv2/monthly_mean;monthly;area;;;D001;H007;E006; +HTAPv2;air_lto;2010;1;;;co,nox_no2,pm10,pm25,oc,bc,so2,voc02,voc03,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc17,voc21,voc22,voc23;/jrc/htapv2/yearly_mean;yearly;area;V003;M001;D001;H001;E007; +HTAPv2;air_cds;2010;1;;;co,nox_no2,pm10,pm25,oc,bc,so2,voc02,voc03,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc17,voc21,voc22,voc23;/jrc/htapv2/yearly_mean;yearly;area;V004;M001;D001;H001;E007; +HTAPv2;air_crs;2010;1;;;co,nox_no2,pm10,pm25,oc,bc,so2,voc02,voc03,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc17,voc21,voc22,voc23;/jrc/htapv2/yearly_mean;yearly;area;V005;M001;D001;H001;E007; +HTAPv2;ships;2010;1;;;co,nox_no2,pm10,pm25,oc,bc,so2,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc18,voc19,voc21,voc22,voc23,voc24;/jrc/htapv2/yearly_mean;yearly;area;;M001;D001;H001;E008; +wiedinmyer;;2010;1;;;bc,c2h2,c2h4,c3h6,c6h6,ch2o,ch3cooh,ch3oh,co,hcl,nh3,nox_no,oc,pm10,pm25,so2;/ucar/wiedinmyer/yearly_mean;yearly;area;;M001;D001;H008;E009; TNO_MACC-III;snap1;2011;0;;;co,nox_no2,so2,nh3,pm10,pm25,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc19,voc21,voc22,voc23,voc24;/tno/tno_macc_iii/yearly_mean/;yearly;area;V001;M002;D002;H002;E010; TNO_MACC-III;snap2;2011;0;;;co,nox_no2,so2,nh3,pm10,pm25,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc19,voc21,voc22,voc23,voc24;/tno/tno_macc_iii/yearly_mean/;yearly;area;;M003;D003;H003;E011; TNO_MACC-III;snap34;2011;0;;;co,nox_no2,so2,nh3,pm10,pm25,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc19,voc21,voc22,voc23,voc24;/tno/tno_macc_iii/yearly_mean/;yearly;area;V002;M004;D003;H004;E012; @@ -87,7 +86,7 @@ EMEP;i_offroad;2015;0;;;co,nox_no2,pm10,pm25,so2,nmvoc,nh3;/ceip/emep EMEP;j_waste;2015;0;;;co,nox_no2,pm10,pm25,so2,nmvoc,nh3;/ceip/emepv18/yearly_mean;yearly;area;;M001;D001;H001;E083; EMEP;k_agrilivestock;2015;0;;;co,nox_no2,pm10,pm25,so2,nmvoc,nh3;/ceip/emepv18/yearly_mean;yearly;area;;M009;D001;H007;E084; EMEP;l_agriother;2015;0;;;co,nox_no2,pm10,pm25,so2,nmvoc,nh3;/ceip/emepv18/yearly_mean;yearly;area;;M009;D001;H007;E085; -carn;;2015;0;;;so2;/mtu/carnetal/yearly_mean;yearly;point;;M001;D001;H001;E086; +carn;;2015;1;;;so2;/mtu/carnetal/yearly_mean;yearly;point;;M001;D001;H001;E086; CEDS;agriculture;2014;0;;;nox_no2,nh3;/jgcri/ceds/monthly_mean;monthly;area;;;D001;H007;E087; CEDS;energy;2014;0;;;co,nox_no2,oc,bc,so2,nh3,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc18,voc19,voc20,voc21,voc22,voc23,voc24;/jgcri/ceds/monthly_mean;monthly;area;V001;;D002;H002;E088; CEDS;industrial;2014;0;;;co,nox_no2,oc,bc,so2,nh3,voc01,voc02,voc03,voc04,voc05,voc06,voc07,voc08,voc09,voc12,voc13,voc14,voc15,voc16,voc17,voc18,voc19,voc20,voc21,voc22,voc23,voc24;/jgcri/ceds/monthly_mean;monthly;area;V002;;D003;H004;E089; diff --git a/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py b/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py index 019a135..4f93956 100755 --- a/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py +++ b/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py @@ -20,6 +20,7 @@ import sys import os import timeit +import warnings import hermesv3_gr.config.settings as settings from emission_inventory import EmissionInventory @@ -232,11 +233,94 @@ class PointGfasEmissionInventory(EmissionInventory): st_time = timeit.default_timer() settings.write_log("\tAllocating GFAS as point sources on grid:", level=2) + netcdf = Dataset(self.pollutant_dicts[0]['path'], mode='r') + gdf = gpd.GeoDataFrame(crs={'init': 'epsg:4326'}) + + # gdf['src_index'] = np.where(self.vertical.altitude.flatten() > 0)[0] + # print self.input_pollutants[0] + first_var = netcdf.variables[self.input_pollutants[0]][:] + # first_var = netcdf.variables['bc'][:] + + # print first_var + gdf['src_index'] = np.where(first_var.flatten() > 0)[0] + + # Masking + if self.masking.regrid_mask is not None: + gdf = gdf.loc[gdf['src_index'].isin(np.where(self.masking.regrid_mask.flatten() > 0)[0]), ] + + gdf['altitude'] = self.vertical.altitude.flatten()[gdf['src_index']] + + lat_1d = netcdf.variables['lat'][:] + lon_1d = netcdf.variables['lon'][:] + lon_1d[lon_1d > 180] -= 360 + # 1D to 2D + lats = np.array([lat_1d] * len(lon_1d)).T + lons = np.array([lon_1d] * len(lat_1d)) + + gdf['geometry'] = [Point(xy) for xy in zip(lons.flatten()[gdf['src_index']], + lats.flatten()[gdf['src_index']])] + gdf['src_area'] = netcdf.variables['cell_area'][:].flatten()[gdf['src_index']] + grid_shp = self.grid.to_shapefile(full_grid=False) + + temp_coords = Dataset(os.path.join(self.grid.temporal_path, 'temporal_coords.nc'), mode='r') + cell_area = temp_coords.variables['cell_area'][:].flatten() + grid_shp['dst_area'] = cell_area[grid_shp['FID']] + # grid_shp['dst_area'] = grid_shp.to_crs({'init': 'epsg:3857'}).area + # print grid_shp.crs['units'] + # sys.exit() + + gdf = gpd.sjoin(gdf.to_crs(grid_shp.crs), grid_shp, how='inner') + + for num, pollutant in enumerate(self.pollutant_dicts): + settings.write_log('\t\tPollutant {0} ({1}/{2})'.format( + pollutant['name'], num + 1, len(self.pollutant_dicts)), level=3) + + aux = netcdf.variables[pollutant['name']][:].flatten()[gdf['src_index']] + + if self.masking.scale_mask is not None: + aux = aux * self.masking.scale_mask.flatten()[gdf['src_index']] + + gdf[pollutant['name']] = aux * gdf['src_area'] + # print 'masa {0}: {1} '.format(pollutant['name'], gdf[pollutant['name']].sum()) + gdf[pollutant['name']] = (aux / gdf['dst_area']) * netcdf.variables['cell_area'][:].flatten()[ + gdf['src_index']] + # print netcdf.variables['bc'][:].sum() + + netcdf.close() + + settings.write_time('PointGfasEmissionInventory', 'do_regrid', timeit.default_timer() - st_time, level=2) + del gdf['src_index'], gdf['index_right'] + self.emissions = gdf + + return True + + def do_regrid_balanced_NOT_USED(self): + import pandas as pd + import geopandas as gpd + from shapely.geometry import Point + import numpy as np + from netCDF4 import Dataset + + st_time = timeit.default_timer() + settings.write_log("\tAllocating GFAS as point sources on grid:", level=2) + netcdf = Dataset(self.pollutant_dicts[0]['path'], mode='r') if settings.rank == 0: gdf = gpd.GeoDataFrame(crs={'init': 'epsg:4326'}) - gdf['src_index'] = np.where(self.vertical.altitude.flatten() > 0)[0] + # gdf['src_index'] = np.where(self.vertical.altitude.flatten() > 0)[0] + # print self.input_pollutants[0] + first_var = netcdf.variables[self.input_pollutants[0]][:] + # first_var = netcdf.variables['bc'][:] + + # print first_var + gdf['src_index'] = np.where(first_var.flatten() > 0)[0] + + # if self.masking.factors_mask_values is not None: + # print self.masking.factors_mask_values + # sys.exit() + + gdf['altitude'] = self.vertical.altitude.flatten()[gdf['src_index']] lat_1d = netcdf.variables['lat'][:] @@ -248,10 +332,18 @@ class PointGfasEmissionInventory(EmissionInventory): gdf['geometry'] = [Point(xy) for xy in zip(lons.flatten()[gdf['src_index']], lats.flatten()[gdf['src_index']])] + gdf['src_area'] = netcdf.variables['cell_area'][:].flatten()[gdf['src_index']] grid_shp = self.grid.to_shapefile() - gdf = gpd.sjoin(gdf.to_crs(grid_shp.crs), grid_shp, how='left') + temp_coords = Dataset(os.path.join(self.grid.temporal_path, 'temporal_coords.nc'), mode = 'r') + cell_area = temp_coords.variables['cell_area'][:].flatten() + grid_shp['dst_area'] = cell_area[grid_shp['FID']] + # grid_shp['dst_area'] = grid_shp.to_crs({'init': 'epsg:3857'}).area + # print grid_shp.crs['units'] + # sys.exit() + gdf = gpd.sjoin(gdf.to_crs(grid_shp.crs), grid_shp, how='inner') + print gdf gdf = np.array_split(gdf, settings.size) else: gdf = None @@ -262,7 +354,11 @@ class PointGfasEmissionInventory(EmissionInventory): settings.write_log('\t\tPollutant {0} ({1}/{2})'.format( pollutant['name'], num + 1, len(self.pollutant_dicts)), level=3) print ('\t\tPollutant {0} ({1}/{2})'.format(pollutant['name'], num + 1, len(self.pollutant_dicts))) - gdf[pollutant['name']] = netcdf.variables[pollutant['name']][:].flatten()[gdf['src_index']] + aux = netcdf.variables[pollutant['name']][:].flatten()[gdf['src_index']] + gdf[pollutant['name']] = aux * gdf['src_area'] + # print 'masa {0}: {1} '.format(pollutant['name'], gdf[pollutant['name']].sum()) + gdf[pollutant['name']] = (aux / gdf['dst_area']) * netcdf.variables['cell_area'][:].flatten()[gdf['src_index']] + # print netcdf.variables['bc'][:].sum() netcdf.close() @@ -288,22 +384,32 @@ class PointGfasEmissionInventory(EmissionInventory): st_time = timeit.default_timer() settings.write_log("\t\tCalculating vertical allocation.", level=3) df = pd.read_csv(vertical_description_path, sep=';') + if self.vertical.approach == 'surface': + self.emissions['altitude'] = 0 self.emissions['layer'] = None for i, line in df.iterrows(): - self.emissions.loc[self.emissions['altitude'] <= line['height_magl'], 'layer'] = line['Ilayer'] - 1 + layer = line['Ilayer'] - 1 + self.emissions.loc[self.emissions['altitude'] <= line['height_magl'], 'layer'] = layer self.emissions.loc[self.emissions['altitude'] <= line['height_magl'], 'altitude'] = None + + # Fires with higher altitudes than the max layer limit goes to the last layer + if len(self.emissions[~self.emissions['altitude'].isna()]) > 0: + self.emissions.loc[~self.emissions['altitude'].isna(), 'layer'] = layer + warnings.warn('WARNING: One or more fires have an altitude of fire emission injection higher than the top' + + ' layer of the model defined in the {0} file'.format(vertical_description_path)) + + print self.emissions.loc[~self.emissions['altitude'].isna()] del self.emissions['altitude'] self.emissions = self.emissions.groupby(['FID', 'layer']).sum() self.emissions.reset_index(inplace=True) - - self.emissions = self.vertical.distribute_vertically(self.emissions) + # print '---> Rank {0}, {1}, {2}'.format(settings.rank, self.emissions, len(self.emissions)) + # sys.exit() + self.emissions = self.vertical.distribute_vertically(self.emissions, self.input_pollutants) settings.write_time('PointGfasEmissionInventory', 'calculate_altitudes', timeit.default_timer() - st_time, level=2) - print 'calculate_altitudes' - sys.exit() return True def point_source_by_cell(self): @@ -312,11 +418,17 @@ class PointGfasEmissionInventory(EmissionInventory): :return: None """ - st_time = timeit.default_timer() - sys.exit() - settings.write_time('PointGfasEmissionInventory', 'Init', timeit.default_timer() - st_time, level=3) - - return None + emissions = self.emissions + self.location = emissions.loc[:, ['FID', 'layer']] + + self.emissions = [] + for input_pollutant in self.input_pollutants: + dict_aux = { + 'name': input_pollutant, + 'units': '...', + 'data': emissions.loc[:, input_pollutant].values + } + self.emissions.append(dict_aux) if __name__ == "__main__": diff --git a/hermesv3_gr/modules/emision_inventories/point_source_emission_inventory.py b/hermesv3_gr/modules/emision_inventories/point_source_emission_inventory.py index 8c22e11..a66f65b 100755 --- a/hermesv3_gr/modules/emision_inventories/point_source_emission_inventory.py +++ b/hermesv3_gr/modules/emision_inventories/point_source_emission_inventory.py @@ -121,9 +121,9 @@ class PointSourceEmissionInventory(EmissionInventory): df = gpd.GeoDataFrame(df.loc[:, ['Emis', 'Alt_Injection']], crs=self.crs, geometry=geometry) df = df.to_crs(grid_shape.crs) + df = gpd.sjoin(df, grid_shape, how="inner", op='intersects') - # Drops duplicates when the point source is on the boundary of the cell df = df[~df.index.duplicated(keep='first')] if self.location is None: diff --git a/hermesv3_gr/modules/vertical/vertical_gfas.py b/hermesv3_gr/modules/vertical/vertical_gfas.py index dd78175..162b8ba 100644 --- a/hermesv3_gr/modules/vertical/vertical_gfas.py +++ b/hermesv3_gr/modules/vertical/vertical_gfas.py @@ -180,14 +180,49 @@ class GfasVerticalDistribution(VerticalDistribution): level=3) return fire_list - def distribute_vertically(self, emissions_df): - print emissions_df + def calculate_weight_layer_dict(self, layer): + weight_layer_dict = {x: None for x in xrange(layer + 1)} + if self.approach == '50_top': + weight_layer_dict[layer] = 0.5 + to_distribute = 0.5 + layer = layer - 1 + elif self.approach == 'uniform': + to_distribute = 1. - vert_emissions = [] + total_width = self.output_heights[layer] + + previous_height = 0 + for i, height in enumerate(self.output_heights[0:layer + 1]): + partial_width = height - previous_height + weight_layer_dict[i] = to_distribute * partial_width / total_width + + previous_height = height + + return weight_layer_dict + def distribute_vertically(self, emissions_df, input_pollutant_list): + import pandas as pd + + vert_emissions = [] for layer, emis in emissions_df.groupby('layer'): - print layer + if layer == 0: + vert_emissions.append(emis) + else: + weight_layer_dict = self.calculate_weight_layer_dict(layer) + for layer, weight in weight_layer_dict.iteritems(): + aux_emis = emis.copy() + aux_emis.loc[:, 'layer'] = layer + aux_emis.loc[:, input_pollutant_list] = aux_emis[input_pollutant_list].multiply(weight) + + vert_emissions.append(aux_emis) + if len(vert_emissions) == 0: + vert_emissions = emissions_df + else: + vert_emissions = pd.concat(vert_emissions) + vert_emissions = vert_emissions.groupby(['FID', 'layer']).sum() + vert_emissions.reset_index(inplace=True) + return vert_emissions if __name__ == '__main__': pass -- GitLab From ad2dafbfc1bb92fa406340f30f133f670342b7a2 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Mon, 18 Feb 2019 17:12:05 +0100 Subject: [PATCH 09/10] Correcting PEP8 errors --- conf/hermes_prace.conf | 79 ------------------- .../emision_inventories/emission_inventory.py | 4 +- .../point_gfas_emission_inventory.py | 15 +--- hermesv3_gr/modules/vertical/vertical_gfas.py | 1 + 4 files changed, 7 insertions(+), 92 deletions(-) delete mode 100644 conf/hermes_prace.conf diff --git a/conf/hermes_prace.conf b/conf/hermes_prace.conf deleted file mode 100644 index 79ddc80..0000000 --- a/conf/hermes_prace.conf +++ /dev/null @@ -1,79 +0,0 @@ -[GENERAL] -log_level = 3 -input_dir = /gpfs/projects/pr1ehu00/pr1ehu03/hermesv3_gr -data_path = /gpfs/scratch/bsc32/bsc32538/HERMES_data -output_dir = /gpfs/projects/pr1ehu00/pr1ehu03/hermesv3_gr_OUT -output_name = HERMESv3_GFAS_.nc -start_date = 2016/05/01 00:00:00 -# ***** end_date = start_date [DEFAULT] ***** -# end_date = 2017/01/02 00:00:00 -# ***** output_timestep_type = [hourly, daily, monthly, yearly] ***** -output_timestep_type = hourly -output_timestep_num = 24 -output_timestep_freq = 1 - -[DOMAIN] - -# ***** output_model = [MONARCH, CMAQ, WRF_CHEM] ***** -output_model = MONARCH -# output_model = CMAQ -# output_model = WRF_CHEM -output_attributes = /data/global_attributes.csv - -# ***** domain_type=[global, lcc, rotated, mercator] ***** -domain_type = global -# domain_type = lcc -# domain_type = rotated -# domain_type = mercator -vertical_description = /data/profiles/vertical/Benchmark_15layers_vertical_description.csv -auxiliar_files_path = /data/auxiliar_files/_ - -# if domain_type == global: - inc_lat = 1. - inc_lon = 1.40625 - -# if domain_type == rotated: - #centre_lat = 35 - #centre_lon = 20 - #west_boundary = -51 - #south_boundary = -35 - #inc_rlat = 0.1 - #inc_rlon = 0.1 - -# if domain_type == lcc: - #lat_1 = 37 - #lat_2 = 43 - #lon_0 = -3 - #lat_0 = 40 - #nx = 478 - #ny = 398 - #inc_x = 12000 - #inc_y = 12000 - #x_0 = -2131849.000 - #y_0 = -2073137.875 - -# if domain_type == mercator: - #lat_ts = -1.5 - #lon_0 = -18 - #nx = 210 - #ny = 236 - #inc_x = 50000 - #inc_y = 50000 - #x_0 = -126017.5 - #y_0 = -5407460 - - -[EMISSION_INVENTORY_CONFIGURATION] - -cross_table = /conf/EI_configuration.csv - -[EMISSION_INVENTORY_PROFILES] - -p_vertical = /data/profiles/vertical/Vertical_profile.csv -p_month = /data/profiles/temporal/TemporalProfile_Monthly.csv -p_day = /data/profiles/temporal/TemporalProfile_Daily.csv -p_hour = /data/profiles/temporal/TemporalProfile_Hourly.csv -p_speciation = /data/profiles/speciation/Speciation_profile_cb05_aero5_MONARCH.csv - -molecular_weights = /data/profiles/speciation/MolecularWeights.csv -world_info = /data/profiles/temporal/tz_world_country_iso3166.csv diff --git a/hermesv3_gr/modules/emision_inventories/emission_inventory.py b/hermesv3_gr/modules/emision_inventories/emission_inventory.py index 29e4551..0c3df52 100644 --- a/hermesv3_gr/modules/emision_inventories/emission_inventory.py +++ b/hermesv3_gr/modules/emision_inventories/emission_inventory.py @@ -365,8 +365,8 @@ class EmissionInventory(object): else: emission_inventory_list.append( PointSourceEmissionInventory(options, grid, date, emission_inventory.ei, - emission_inventory.source_type, emission_inventory.sector, pollutants, - emission_inventory_path, + emission_inventory.source_type, emission_inventory.sector, + pollutants, emission_inventory_path, emission_inventory.frequency, vertical_output_profile, reference_year=emission_inventory.ref_year, factors=emission_inventory.factor_mask, diff --git a/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py b/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py index 4f93956..3c675e9 100755 --- a/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py +++ b/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py @@ -89,7 +89,8 @@ class PointGfasEmissionInventory(EmissionInventory): # self.altitude = self.get_altitude() - self.vertical = GfasVerticalDistribution(vertical_output_profile, self.get_approach(p_vertical), self.get_altitude()) + self.vertical = GfasVerticalDistribution(vertical_output_profile, self.get_approach(p_vertical), + self.get_altitude()) settings.write_time('PointGfasEmissionInventory', 'Init', timeit.default_timer() - st_time, level=3) @@ -308,19 +309,10 @@ class PointGfasEmissionInventory(EmissionInventory): if settings.rank == 0: gdf = gpd.GeoDataFrame(crs={'init': 'epsg:4326'}) - # gdf['src_index'] = np.where(self.vertical.altitude.flatten() > 0)[0] - # print self.input_pollutants[0] first_var = netcdf.variables[self.input_pollutants[0]][:] - # first_var = netcdf.variables['bc'][:] - # print first_var gdf['src_index'] = np.where(first_var.flatten() > 0)[0] - # if self.masking.factors_mask_values is not None: - # print self.masking.factors_mask_values - # sys.exit() - - gdf['altitude'] = self.vertical.altitude.flatten()[gdf['src_index']] lat_1d = netcdf.variables['lat'][:] @@ -357,7 +349,8 @@ class PointGfasEmissionInventory(EmissionInventory): aux = netcdf.variables[pollutant['name']][:].flatten()[gdf['src_index']] gdf[pollutant['name']] = aux * gdf['src_area'] # print 'masa {0}: {1} '.format(pollutant['name'], gdf[pollutant['name']].sum()) - gdf[pollutant['name']] = (aux / gdf['dst_area']) * netcdf.variables['cell_area'][:].flatten()[gdf['src_index']] + gdf[pollutant['name']] = (aux / gdf['dst_area']) * netcdf.variables['cell_area'][:].flatten()[ + gdf['src_index']] # print netcdf.variables['bc'][:].sum() netcdf.close() diff --git a/hermesv3_gr/modules/vertical/vertical_gfas.py b/hermesv3_gr/modules/vertical/vertical_gfas.py index 162b8ba..b845208 100644 --- a/hermesv3_gr/modules/vertical/vertical_gfas.py +++ b/hermesv3_gr/modules/vertical/vertical_gfas.py @@ -224,5 +224,6 @@ class GfasVerticalDistribution(VerticalDistribution): vert_emissions.reset_index(inplace=True) return vert_emissions + if __name__ == '__main__': pass -- GitLab From 43c3f4a502875da0dc384e473c0cdea350ef52c5 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Mon, 18 Feb 2019 17:14:57 +0100 Subject: [PATCH 10/10] Correcting PEP8 errors --- .../emision_inventories/point_gfas_emission_inventory.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py b/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py index 3c675e9..35663df 100755 --- a/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py +++ b/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py @@ -327,7 +327,7 @@ class PointGfasEmissionInventory(EmissionInventory): gdf['src_area'] = netcdf.variables['cell_area'][:].flatten()[gdf['src_index']] grid_shp = self.grid.to_shapefile() - temp_coords = Dataset(os.path.join(self.grid.temporal_path, 'temporal_coords.nc'), mode = 'r') + temp_coords = Dataset(os.path.join(self.grid.temporal_path, 'temporal_coords.nc'), mode='r') cell_area = temp_coords.variables['cell_area'][:].flatten() grid_shp['dst_area'] = cell_area[grid_shp['FID']] # grid_shp['dst_area'] = grid_shp.to_crs({'init': 'epsg:3857'}).area -- GitLab