diff --git a/CHANGELOG b/CHANGELOG index 81ff9542506469619bc3a6375abbc03a25003d1e..9a76d15812d2cd918101180c2eb5f194f640d161 100755 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,3 +1,8 @@ +1.0.6 + 2022/02/18 + + - Add traffic scenario shapefile + 1.0.5 2020/12/09 diff --git a/conf/hermes.conf b/conf/hermes.conf index 189868ec06cb2e9817c0ca24b5a7f8898a447a33..37fd0691eb2301f81b19bb86554c602ee3dfce1e 100755 --- a/conf/hermes.conf +++ b/conf/hermes.conf @@ -330,6 +330,7 @@ traffic_speciation_profile_tyre = /profiles/speciation/traffic/tyre_b traffic_speciation_profile_road = /profiles/speciation/traffic/road_base.csv traffic_speciation_profile_brake = /profiles/speciation/traffic/brake_base.csv traffic_speciation_profile_resuspension = /profiles/speciation/traffic/resuspension_base.csv +traffic_scenario = /traffic/shp_ZPE [TRAFFIC AREA SECTOR] traffic_area_pollutants = nox_no2,nmvoc,so2,co,nh3,pm10,pm25 diff --git a/environment.yml b/environment.yml index 584e815d9fbd3e5843f104f8f1f0bbab6c8866e2..a32aac5fd301c024bf78134856ba1a520a84c12d 100755 --- a/environment.yml +++ b/environment.yml @@ -7,7 +7,7 @@ channels: - anaconda dependencies: - - python = 2 + - python = 3 - numpy - netcdf4 >= 1.3.1 - python-cdo >= 1.3.6 diff --git a/hermesv3_bu/__init__.py b/hermesv3_bu/__init__.py index 92192eed4fc88eaee9fdbafbc8e2139c43ba00e2..382021f30ff6ce35872aa0de4ebe2a43f50c21d7 100755 --- a/hermesv3_bu/__init__.py +++ b/hermesv3_bu/__init__.py @@ -1 +1 @@ -__version__ = "1.0.4" +__version__ = "1.0.6" diff --git a/hermesv3_bu/config/config.py b/hermesv3_bu/config/config.py index f9141b84c8773cedba4d9373abdf59d830009cd3..9e6e07a08831c81fdc8f80d91316a912a4ab8350 100755 --- a/hermesv3_bu/config/config.py +++ b/hermesv3_bu/config/config.py @@ -5,6 +5,7 @@ from configargparse import ArgParser import os from mpi4py import MPI from hermesv3_bu.tools.checker import error_exit +from hermesv3_bu import __version__ class Config(ArgParser): @@ -40,6 +41,7 @@ class Config(ArgParser): p = ArgParser() p.add_argument('-c', '--my-config', required=False, is_config_file=True, help='Path to the configuration file.') + p.add_argument('--version', '-V', action='version', version="%(prog)s " + __version__) # ===== GENERAL ===== p.add_argument('--log_level', required=True, help='Level of detail of the running process information.', type=int, choices=[1, 2, 3]) @@ -625,6 +627,8 @@ class Config(ArgParser): p.add_argument('--traffic_speciation_profile_resuspension', required=False, help="Defines the path to the CSV file that contains the speciation profiles for the " + "resuspension emissions.") + p.add_argument('--traffic_scenario', required=False, default=None, + help="Defines the path to the shapefile that contains the traffic scenario by input pollutant.") # ***** TRAFFIC AREA SECTOR ***** p.add_argument('--traffic_area_pollutants', required=False, @@ -725,6 +729,9 @@ class Config(ArgParser): self.comm.Barrier() self.create_dir(arguments.auxiliary_files_path) + if arguments.clipping == "None": + arguments.clipping = None + # Booleans arguments.do_traffic = arguments.traffic_processors > 0 arguments.do_traffic_area = arguments.traffic_area_processors > 0 @@ -794,6 +801,8 @@ class Config(ArgParser): # Traffic lists arguments.traffic_pollutants = self._parse_list(arguments.traffic_pollutants) + if arguments.vehicle_types == "All": + arguments.vehicle_types = None arguments.vehicle_types = self._parse_list(arguments.vehicle_types) # Traffic area bools diff --git a/hermesv3_bu/io_server/io_shapefile.py b/hermesv3_bu/io_server/io_shapefile.py index 59ece649cf221a86df7b49557988fe4114f3affc..bdc397f7040375e05c8810240569c0c85b3c530a 100755 --- a/hermesv3_bu/io_server/io_shapefile.py +++ b/hermesv3_bu/io_server/io_shapefile.py @@ -88,6 +88,33 @@ class IoShapefile(IoServer): return data + def read_shapefile_broadcast(self, path, rank=0, crs=None): + """ + The master process reads the shapefile and broadcast it. + + :param path: Path to the shapefile + :type path: str + + :param rank: Master rank + :type rank: int + + :param crs: Projection desired. + :type crs: None, str + + :return: Shapefile + :rtype: GeoDataFrame + """ + if self.comm.Get_rank() == rank: + data = self.read_shapefile_serial(path) + if crs is not None: + data = data.to_crs(crs) + else: + data = None + + data = self.comm.bcast(data, root=0) + + return data + def split_shapefile(self, data, rank=0): """ diff --git a/hermesv3_bu/sectors/sector_manager.py b/hermesv3_bu/sectors/sector_manager.py index b726d3229b1de9baf3e42333a15446b68c3298d2..575730a6a150eb3c9f30855a7ad815235f3d8d6d 100755 --- a/hermesv3_bu/sectors/sector_manager.py +++ b/hermesv3_bu/sectors/sector_manager.py @@ -188,7 +188,7 @@ class SectorManager(object): arguments.temperature_hourly_files_path, arguments.output_dir, arguments.molecular_weights, arguments.resuspension_correction, arguments.precipitation_files_path, arguments.do_hot, arguments.do_cold, arguments.do_tyre_wear, arguments.do_brake_wear, arguments.do_road_wear, - arguments.do_resuspension, arguments.write_rline) + arguments.do_resuspension, arguments.write_rline, arguments.traffic_scenario) elif sector == 'traffic_area' and comm_world.Get_rank() in sector_procs: from hermesv3_bu.sectors.traffic_area_sector import TrafficAreaSector diff --git a/hermesv3_bu/sectors/traffic_sector.py b/hermesv3_bu/sectors/traffic_sector.py index 35e1e3015a240673397b673901880a2721add382..f2736069556c5750f3930707639798e3cdc9e08b 100755 --- a/hermesv3_bu/sectors/traffic_sector.py +++ b/hermesv3_bu/sectors/traffic_sector.py @@ -17,7 +17,6 @@ from hermesv3_bu.io_server.io_shapefile import IoShapefile from hermesv3_bu.tools.checker import check_files, error_exit import gc -from memory_profiler import profile from ctypes import cdll, CDLL cdll.LoadLibrary("libc.so.6") libc = CDLL("libc.so.6") @@ -33,7 +32,7 @@ FINAL_PROJ = {'init': 'epsg:3035'} # https://epsg.io/3035 ETRS89 / LAEA Europe aerosols = ['oc', 'ec', 'pno3', 'pso4', 'pmfine', 'pmc', 'poa', 'poc', 'pec', 'pcl', 'pnh4', 'pna', 'pmg', 'pk', 'pca', 'pncom', 'pfe', 'pal', 'psi', 'pti', 'pmn', 'ph2o', 'pmothr'] pmc_list = ['pmc', 'PMC'] -rline_shp = False +rline_shp = True class TrafficSector(Sector): @@ -56,7 +55,7 @@ class TrafficSector(Sector): hot_cold_speciation=None, tyre_speciation=None, road_speciation=None, brake_speciation=None, resuspension_speciation=None, temp_common_path=None, output_dir=None, molecular_weights_path=None, resuspension_correction=True, precipitation_path=None, do_hot=True, do_cold=True, do_tyre_wear=True, - do_brake_wear=True, do_road_wear=True, do_resuspension=True, write_rline=False): + do_brake_wear=True, do_road_wear=True, do_resuspension=True, write_rline=False, traffic_scenario=None): spent_time = timeit.default_timer() logger.write_log('===== TRAFFIC SECTOR =====') @@ -136,6 +135,9 @@ class TrafficSector(Sector): hourly_saturday_profiles_path, hourly_sunday_profiles_path) self.check_profiles() self.expanded = self.expand_road_links() + if traffic_scenario in ['=', '', ' ', 'None']: + traffic_scenario = None + self.scenario = self.get_scenario(traffic_scenario) del self.fleet_compo, self.speed_hourly, self.monthly_profiles, self.weekly_profiles, self.hourly_profiles libc.malloc_trim(0) @@ -167,6 +169,30 @@ class TrafficSector(Sector): libc.malloc_trim(0) gc.collect() + def get_scenario(self, scenario_path): + if scenario_path in ['', 'None']: + scenario_path = None + if scenario_path is not None and not os.path.exists(scenario_path): + msg = "ERROR!!! " + msg += "Traffic scenario file '{0}' not found!".format(scenario_path) + error_exit(msg) + + if scenario_path is not None: + self.logger.write_log('\t\tGetting emission scenario', message_level=2) + scenario_shp = IoShapefile(self.comm).read_shapefile_broadcast(scenario_path, crs=self.road_links.crs) + + scenario = gpd.sjoin(GeoDataFrame(index=self.road_links.index, geometry=self.road_links.geometry.centroid, + crs=self.road_links.crs), scenario_shp, how='left') + scenario = scenario.loc[:, scenario_shp.columns] + scenario = pd.DataFrame(scenario.drop(columns='geometry'), index=pd.Index(scenario.index, name='Link_ID')) + + scenario.fillna(1, inplace=True) + + else: + scenario = None + + return scenario + def check_profiles(self): spent_time = timeit.default_timer() # Checking speed profiles IDs @@ -512,8 +538,8 @@ class TrafficSector(Sector): # Checks that the splited DataFrames contain the full DataFrame if (len(df_code_slope_road) + len(df_code_slope) + len(df_code_road) + len(df_code)) != len(df): - # TODO check that error - error_exit('ERROR in blablavbla') + error_exit('ERROR in the Emission Factor file {0}. Check the Road.Slope and Mode info'.format( + ef_path)) return df_code_slope_road, df_code_slope, df_code_road, df_code elif emission_type == 'cold' or emission_type == 'tyre' or emission_type == 'road' or \ @@ -1356,6 +1382,18 @@ class TrafficSector(Sector): self.logger.write_time_log('TrafficSector', 'speciate_traffic', timeit.default_timer() - spent_time) return df_out + def apply_scenario(self, emissions): + emissions.index = emissions.index.set_levels(emissions.index.levels[0].astype(int), level=0) + + self.logger.write_log('\t\tApplying emission scenario', message_level=2) + emis_aux = emissions.join(self.scenario, on='Link_ID', rsuffix='_f') + + for pollutant in self.scenario.columns: + if pollutant in emissions.columns: + self.logger.write_log('\t\t\tApplying emission scenario for {0}'.format(pollutant), message_level=3) + emissions.loc[emis_aux.index, pollutant] *= emis_aux.loc[:, '{0}_f'.format(pollutant)] + return emissions + def calculate_emissions(self): spent_time = timeit.default_timer() self.logger.write_log('\tCalculating Road traffic emissions', message_level=1) @@ -1419,6 +1457,9 @@ class TrafficSector(Sector): libc.malloc_trim(0) df_accum.set_index(['Link_ID', 'tstep'], inplace=True) + if self.scenario is not None: + df_accum = self.apply_scenario(df_accum) + if self.write_rline: self.write_rline_output(df_accum.copy()) self.logger.write_log('\t\tRoad link emissions to grid.', message_level=2) @@ -1557,7 +1598,8 @@ class TrafficSector(Sector): df_in = df_in.to_crs({u'units': u'm', u'no_defs': True, u'ellps': u'intl', u'proj': u'utm', u'zone': 31}) if rline_shp: - df_in.to_file(os.path.join(self.output_dir, 'roads.shp')) + if not os.path.exists(os.path.join(self.output_dir, 'roads.shp')): + df_in.to_file(os.path.join(self.output_dir, 'roads.shp')) count = 0 for i, line in df_in.iterrows(): @@ -1604,7 +1646,8 @@ class TrafficSector(Sector): df_out.set_index('Link_ID', inplace=True) df_out.sort_index(inplace=True) - df_out.to_csv(os.path.join(self.output_dir, 'roads.txt'), index=False, sep=' ') + if not os.path.exists(os.path.join(self.output_dir, 'roads.txt')): + df_out.to_csv(os.path.join(self.output_dir, 'roads.txt'), index=False, sep=' ') self.comm.Barrier() self.logger.write_log('\t\tTraffic emissions calculated', message_level=2) self.logger.write_time_log('TrafficSector', 'write_rline_roadlinks', timeit.default_timer() - spent_time) diff --git a/hermesv3_bu/writer/default_writer.py b/hermesv3_bu/writer/default_writer.py index c791c32f808ebac7891b83e47c9744ced78eda98..f0799e9a5dd0f72db302a755dd8b7b6239143116 100755 --- a/hermesv3_bu/writer/default_writer.py +++ b/hermesv3_bu/writer/default_writer.py @@ -222,14 +222,14 @@ class DefaultWriter(Writer): var.units = self.pollutant_info.loc[var_name, 'units'] var.missing_value = -999.0 var.coordinates = 'lat lon' - # if self.grid.grid_type == 'Regular Lat-Lon': - # var.grid_mapping = 'Latitude_Longitude' - # elif self.grid.grid_type == 'Lambert Conformal Conic': - # var.grid_mapping = 'Lambert_Conformal' - # elif self.grid.grid_type in ['Rotated', 'Rotated_nested']: - # var.grid_mapping = 'rotated_pole' - # elif self.grid.grid_type == 'Mercator': - # var.grid_mapping = 'mercator' + if self.grid.grid_type == 'Regular Lat-Lon': + var.grid_mapping = 'Latitude_Longitude' + elif self.grid.grid_type == 'Lambert Conformal Conic': + var.grid_mapping = 'Lambert_Conformal' + elif self.grid.grid_type in ['Rotated', 'Rotated_nested']: + var.grid_mapping = 'rotated_pole' + elif self.grid.grid_type == 'Mercator': + var.grid_mapping = 'mercator' # ========== METADATA ========== self.logger.write_log('\tCreating NetCDF metadata', message_level=2) @@ -242,14 +242,14 @@ class DefaultWriter(Writer): mapping.semi_major_axis = 6371000.0 mapping.inverse_flattening = 0 - # elif self.grid.grid_type == 'Lambert Conformal Conic': - # mapping = netcdf.createVariable('Lambert_Conformal', 'i') - # mapping.grid_mapping_name = "lambert_conformal_conic" - # mapping.standard_parallel = [self.grid.attributes['lat_2'], self.grid.attributes['lat_1']] - # mapping.longitude_of_central_meridian = self.grid.attributes['lon_0'] - # mapping.latitude_of_projection_origin = self.grid.attributes['lat_0'] - # mapping.false_easting = self.grid.attributes['x_0'] - # mapping.false_northing = self.grid.attributes['y_0'] + elif self.grid.grid_type == 'Lambert Conformal Conic': + mapping = netcdf.createVariable('Lambert_Conformal', 'i') + mapping.grid_mapping_name = "lambert_conformal_conic" + mapping.standard_parallel = [self.grid.attributes['lat_2'], self.grid.attributes['lat_1']] + mapping.longitude_of_central_meridian = self.grid.attributes['lon_0'] + mapping.latitude_of_projection_origin = self.grid.attributes['lat_0'] + mapping.false_easting = self.grid.attributes['x_0'] + mapping.false_northing = self.grid.attributes['y_0'] elif self.grid.grid_type in ['Rotated', 'Rotated_nested']: mapping = netcdf.createVariable('rotated_pole', 'c')