diff --git a/README.md b/README.md index b393f1f503d737afb63e16c85f268c9cbc7ee09c..4aa21822ba27d1726eee9823c03800eafd878a21 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,27 @@ -# 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. + +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). + +## 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 diff --git a/conf/EI_configuration.csv b/conf/EI_configuration.csv index 2dee5ebe57b90359c0438a57d6271f2afae63bf5..2d33a264e3341b27ac46a24a79970d8664eaac22 100644 --- a/conf/EI_configuration.csv +++ b/conf/EI_configuration.csv @@ -1,5 +1,5 @@ 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; +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; diff --git a/hermesv3_gr/modules/emision_inventories/emission_inventory.py b/hermesv3_gr/modules/emision_inventories/emission_inventory.py index 3101a2fccdef07220c15f2b448d56b2971c302d3..0c3df52470998c04e3c611ed282cff29992e0426 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']) @@ -207,8 +208,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 +232,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 +287,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 +352,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 0000000000000000000000000000000000000000..35663df73328033b9d423ba3195146e731750c49 --- /dev/null +++ b/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py @@ -0,0 +1,428 @@ +#!/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 warnings + +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') + 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'}) + + first_var = netcdf.variables[self.input_pollutants[0]][:] + + gdf['src_index'] = np.where(first_var.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() + + 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 + + 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))) + 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() + + 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=';') + if self.vertical.approach == 'surface': + self.emissions['altitude'] = 0 + + self.emissions['layer'] = None + for i, line in df.iterrows(): + 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) + # 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) + return True + + def point_source_by_cell(self): + """ + Sums the different emissions that are allocated in the same cell and layer. + + :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__": + pass 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 8c22e114aa539794cd8a57109f9cb83820ab6e0d..a66f65b638c0665d6d174262d298f27a46eb6cf4 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/grids/grid.py b/hermesv3_gr/modules/grids/grid.py index 0c424d1ecd90aa27f068164d514c93216535a649..fe15d965030915faaf3fe781bb0df1b9964a2d3f 100644 --- a/hermesv3_gr/modules/grids/grid.py +++ b/hermesv3_gr/modules/grids/grid.py @@ -250,6 +250,7 @@ class Grid(object): regular_latlon=True) # Calculate the cell area of the auxiliary NetCDF file + self.cell_area = self.get_cell_area() # Re-writes the NetCDF adding the cell area diff --git a/hermesv3_gr/modules/vertical/vertical_gfas.py b/hermesv3_gr/modules/vertical/vertical_gfas.py index 600ab7f6fef3559a064f0f37ea209847ada65260..b845208ad25749d2392741840b10322b7c9f2796 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 @@ -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): """ @@ -172,6 +180,50 @@ class GfasVerticalDistribution(VerticalDistribution): level=3) return fire_list + 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. + + 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'): + 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 diff --git a/hermesv3_gr/modules/writing/writer_monarch.py b/hermesv3_gr/modules/writing/writer_monarch.py index 3321b06c0f0945da74668f8e8d50d6c6bb503663..e5019262657d9e68302d2e049212a64ed015dafc 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/preproc/GFAS_hourly_Parameters.csv b/preproc/GFAS_hourly_Parameters.csv new file mode 100644 index 0000000000000000000000000000000000000000..7b8f414d2ed167ab10262cdfea3a90e8826b63b4 --- /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 +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 +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 0000000000000000000000000000000000000000..9cc93fef153e21e8666355ab350c815145f1cf56 --- /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) diff --git a/run_test.py b/run_test.py index 6eb2926793f63a9c62eee79c06273d12f95fab82..711c425bbcb719bc46c277b18d1f4f754326274b 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