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 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)