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/conf/hermes_prace.conf b/conf/hermes_prace.conf new file mode 100644 index 0000000000000000000000000000000000000000..79ddc8068fda2cad13d1e5218bfa3c130bb6c100 --- /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 3101a2fccdef07220c15f2b448d56b2971c302d3..29e45519417dd45d8653e067c070a6ff48d15506 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..4f9395652cb1910bf4a145db98c394ef71161f32 --- /dev/null +++ b/hermesv3_gr/modules/emision_inventories/point_gfas_emission_inventory.py @@ -0,0 +1,435 @@ +#!/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'}) + + # 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'][:] + 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..162b8ba16ce1374ae759b47a365d016e46f063a8 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,49 @@ 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