From c84eccb3cc1faae4a313ba04438c4207cbbd040c Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Mon, 26 Aug 2019 14:27:14 +0200 Subject: [PATCH 01/10] Updating changelog --- CHANGELOG | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/CHANGELOG b/CHANGELOG index 8238cd1..222e0e1 100755 --- a/CHANGELOG +++ b/CHANGELOG @@ -16,6 +16,19 @@ - Sector Manager: 1. Aviation sector + 2. Shipping port sector + 3. Livestock sector + 4. Crop operations sector + 5. Crop fertilizers sector + 6. Agricultural machinery sector + 7. Residential combustion sector + 8. Recreational boats sector + 9. Point sources sector + 10. Road traffic sector + 11. Traffic area (evaporative & small cities) sector - Writing options: - 1. Default writer \ No newline at end of file + 1. Default writer + 2. CMAQ writer + 3. MONARCH writer + 4. WRF-Chem writer \ No newline at end of file -- GitLab From c7c8d595353e6b6a426bd1e83ffcaefa86de95f4 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Mon, 26 Aug 2019 18:56:31 +0200 Subject: [PATCH 02/10] -Agricultural machinery: corrected parallel bug --- hermesv3_bu/sectors/agricultural_machinery_sector.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hermesv3_bu/sectors/agricultural_machinery_sector.py b/hermesv3_bu/sectors/agricultural_machinery_sector.py index 12a7be3..78c9a6e 100755 --- a/hermesv3_bu/sectors/agricultural_machinery_sector.py +++ b/hermesv3_bu/sectors/agricultural_machinery_sector.py @@ -35,7 +35,7 @@ class AgriculturalMachinerySector(AgriculturalSector): self.crop_machinery_nuts3 = self.read_profiles(crop_machinery_nuts3) self.crop_distribution = self.get_crop_distribution_by_nut( - self.crop_distribution, machinery_distibution_nut_shapefile_path, nut_code='nuts3_id') + self.crop_distribution, machinery_distibution_nut_shapefile_path, nut_code='nuts3_id') self.months = self.get_date_array_by_month() @@ -86,7 +86,7 @@ class AgriculturalMachinerySector(AgriculturalSector): crop_distribution.drop(columns=self.crop_list, inplace=True) crop_distribution.rename(columns={nut_code: 'NUT_code'}, inplace=True) - IoShapefile(self.comm).write_shapefile_serial(crop_distribution, crop_distribution_nut_path) + IoShapefile(self.comm).write_shapefile_parallel(crop_distribution, crop_distribution_nut_path) else: crop_distribution = IoShapefile(self.comm).read_shapefile(crop_distribution_nut_path) -- GitLab From 438250fd1c05cfa598df011fcbff16585f75e0af Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Mon, 26 Aug 2019 18:57:32 +0200 Subject: [PATCH 03/10] -Traffic added check that checks that all the profiles that are used in the road links are in the profiles files --- conf/hermes.conf | 22 +++++------ hermesv3_bu/sectors/traffic_sector.py | 53 ++++++++++++++++++++++++++- 2 files changed, 63 insertions(+), 12 deletions(-) diff --git a/conf/hermes.conf b/conf/hermes.conf index ca1c47c..d4ae43b 100755 --- a/conf/hermes.conf +++ b/conf/hermes.conf @@ -3,20 +3,20 @@ log_level = 3 input_dir = /home/Earth/ctena/Models/hermesv3_bu_data data_path = /esarchive/recon output_dir = /scratch/Earth/HERMESv3_BU_OUT -output_name = HERMESv3__tr_area.nc +output_name = HERMESv3__traffic.nc emission_summary = 0 start_date = 2016/11/29 00:00:00 # ----- end_date = start_date [DEFAULT] ----- # end_date = 2010/01/01 00:00:00 output_timestep_num = 24 auxiliary_files_path = /scratch/Earth/HERMESv3_BU_aux/_ -erase_auxiliary_files = 1 +erase_auxiliary_files = 0 [DOMAIN] # domain_type=[lcc, rotated, mercator, regular] -domain_type = lcc +domain_type = regular # output_type=[MONARCH, CMAQ, WRF_CHEM, DEFAULT] output_model = DEFAULT output_attributes = /writing/global_attributes_WRF-Chem.csv @@ -107,12 +107,12 @@ vertical_description = /profiles/vertical/MONARCH_Global_48layers_ver #y_0 = -5407460 # if domain_type == regular: - lat_orig = 41.1 - lon_orig = 1.8 - inc_lat = 0.1 - inc_lon = 0.1 - n_lat = 10 - n_lon = 10 + lat_orig = 40.5 + lon_orig = 0.0 + inc_lat = 0.05 + inc_lon = 0.05 + n_lat = 50 + n_lon = 70 [CLIPPING] @@ -147,8 +147,8 @@ agricultural_machinery_processors = 0 residential_processors = 0 recreational_boats_processors = 0 point_sources_processors = 0 -traffic_processors = 0 -traffic_area_processors = 1 +traffic_processors = 1 +traffic_area_processors = 0 [SHAPEFILES] diff --git a/hermesv3_bu/sectors/traffic_sector.py b/hermesv3_bu/sectors/traffic_sector.py index db584b1..a8b2f96 100755 --- a/hermesv3_bu/sectors/traffic_sector.py +++ b/hermesv3_bu/sectors/traffic_sector.py @@ -89,7 +89,7 @@ class TrafficSector(Sector): self.hourly_profiles = self.read_all_hourly_profiles(hourly_mean_profiles_path, hourly_weekday_profiles_path, hourly_saturday_profiles_path, hourly_sunday_profiles_path) - + self.check_profiles() self.expanded = self.expand_road_links() del self.fleet_compo, self.speed_hourly, self.monthly_profiles, self.weekly_profiles, self.hourly_profiles @@ -103,6 +103,57 @@ class TrafficSector(Sector): self.logger.write_time_log('TrafficSector', '__init__', timeit.default_timer() - spent_time) + def check_profiles(self): + spent_time = timeit.default_timer() + # Checking speed profiles IDs + links_speed = set(np.unique(np.concatenate([ + np.unique(self.road_links['sp_hour_su'].dropna().values), + np.unique(self.road_links['sp_hour_mo'].dropna().values), + np.unique(self.road_links['sp_hour_tu'].dropna().values), + np.unique(self.road_links['sp_hour_we'].dropna().values), + np.unique(self.road_links['sp_hour_th'].dropna().values), + np.unique(self.road_links['sp_hour_fr'].dropna().values), + np.unique(self.road_links['sp_hour_sa'].dropna().values), + ]))) + # The '0' speed profile means that we don't know the speed profile and it will be replaced by a flat profile + speed = set(np.unique(np.concatenate([self.speed_hourly.index.values, [0]]))) + + speed_res = links_speed - speed + if len(speed_res) > 0: + raise ValueError("The following speed profile IDs reported in the road links shapefile do not appear " + + "in the hourly speed profiles file. {0}".format(speed_res)) + + # Checking monthly profiles IDs + links_month = set(np.unique(self.road_links['aadt_m_mn'].dropna().values)) + month = set(self.monthly_profiles.index.values) + month_res = links_month - month + if len(month_res) > 0: + raise ValueError("The following monthly profile IDs reported in the road links shapefile do not appear " + + "in the monthly profiles file. {0}".format(month_res)) + + # Checking weekly profiles IDs + links_week = set(np.unique(self.road_links['aadt_week'].dropna().values)) + week = set(self.weekly_profiles.index.values) + week_res = links_week - week + if len(week_res) > 0: + raise ValueError("The following weekly profile IDs reported in the road links shapefile do not appear " + + "in the weekly profiles file. {0}".format(week_res)) + + # Checking hourly profiles IDs + links_hour = set(np.unique(np.concatenate([ + np.unique(self.road_links['aadt_h_mn'].dropna().values), + np.unique(self.road_links['aadt_h_wd'].dropna().values), + np.unique(self.road_links['aadt_h_sat'].dropna().values), + np.unique(self.road_links['aadt_h_sun'].dropna().values), + ]))) + hour = set(self.hourly_profiles.index.values) + hour_res = links_hour - hour + if len(hour_res) > 0: + raise ValueError("The following hourly profile IDs reported in the road links shapefile do not appear " + + "in the hourly profiles file. {0}".format(hour_res)) + + self.logger.write_time_log('TrafficSector', 'check_profiles', timeit.default_timer() - spent_time) + def read_all_hourly_profiles(self, hourly_mean_profiles_path, hourly_weekday_profiles_path, hourly_saturday_profiles_path, hourly_sunday_profiles_path): hourly_profiles = pd.concat([self.read_hourly_profiles(hourly_mean_profiles_path), -- GitLab From 24a18b7f69d81eca27c1a9a6319176f433f8ebb3 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Tue, 27 Aug 2019 11:28:02 +0200 Subject: [PATCH 04/10] -Traffic: Forced to espg: 3035 the rotated and regular grids to calculate the link length in meters - All sectors: Pass Grid object instead grid shapefile --- conf/hermes.conf | 58 ++++++++----------- .../agricultural_crop_fertilizers_sector.py | 12 ++-- .../agricultural_crop_operations_sector.py | 4 +- .../sectors/agricultural_machinery_sector.py | 4 +- hermesv3_bu/sectors/agricultural_sector.py | 14 ++--- hermesv3_bu/sectors/aviation_sector.py | 22 +++---- hermesv3_bu/sectors/livestock_sector.py | 16 ++--- hermesv3_bu/sectors/point_source_sector.py | 8 +-- .../sectors/recreational_boats_sector.py | 12 ++-- hermesv3_bu/sectors/residential_sector.py | 12 ++-- hermesv3_bu/sectors/sector.py | 9 +-- hermesv3_bu/sectors/sector_manager.py | 22 +++---- hermesv3_bu/sectors/shipping_port_sector.py | 10 ++-- hermesv3_bu/sectors/traffic_area_sector.py | 14 ++--- hermesv3_bu/sectors/traffic_sector.py | 19 +++--- 15 files changed, 117 insertions(+), 119 deletions(-) diff --git a/conf/hermes.conf b/conf/hermes.conf index d4ae43b..e11afec 100755 --- a/conf/hermes.conf +++ b/conf/hermes.conf @@ -2,8 +2,8 @@ log_level = 3 input_dir = /home/Earth/ctena/Models/hermesv3_bu_data data_path = /esarchive/recon -output_dir = /scratch/Earth/HERMESv3_BU_OUT -output_name = HERMESv3__traffic.nc +output_dir = /scratch/Earth/HERMESv3_BU_OUT/serial +output_name = HERMESv3__residential.nc emission_summary = 0 start_date = 2016/11/29 00:00:00 # ----- end_date = start_date [DEFAULT] ----- @@ -16,7 +16,7 @@ erase_auxiliary_files = 0 [DOMAIN] # domain_type=[lcc, rotated, mercator, regular] -domain_type = regular +domain_type = lcc # output_type=[MONARCH, CMAQ, WRF_CHEM, DEFAULT] output_model = DEFAULT output_attributes = /writing/global_attributes_WRF-Chem.csv @@ -31,12 +31,12 @@ vertical_description = /profiles/vertical/MONARCH_Global_48layers_ver # inc_rlat = 0.2 # inc_rlon = 0.2 - centre_lat = 40.5 - centre_lon = -3.5 - west_boundary = -7.0 - south_boundary = -7.0 - inc_rlat = 0.4 - inc_rlon = 0.4 + # centre_lat = 40.5 + # centre_lon = -3.5 + # west_boundary = -7.0 + # south_boundary = -7.0 + # inc_rlat = 0.4 + # inc_rlon = 0.4 # if domain_type == lcc: @@ -46,14 +46,6 @@ vertical_description = /profiles/vertical/MONARCH_Global_48layers_ver lon_0 = -3 lat_0 = 40 - # TEST - #nx = 30 - #ny = 30 - #inc_x = 10000 - #inc_y = 10000 - #x_0 = 253151.59375 - #y_0 = 43862.90625 - # CATALUNYA #nx = 278 #ny = 298 @@ -97,22 +89,22 @@ vertical_description = /profiles/vertical/MONARCH_Global_48layers_ver #y_0 = -20137.891 # 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 + # lat_ts = -1.5 + # lon_0 = -18 + # nx = 10 + # ny = 10 + # inc_x = 50000 + # inc_y = 50000 + # x_0 = -126017.5 + # y_0 = -5407460 # if domain_type == regular: - lat_orig = 40.5 - lon_orig = 0.0 - inc_lat = 0.05 - inc_lon = 0.05 - n_lat = 50 - n_lon = 70 + # lat_orig = 40.5 + # lon_orig = 0.0 + # inc_lat = 0.05 + # inc_lon = 0.05 + # n_lat = 50 + # n_lon = 70 [CLIPPING] @@ -144,10 +136,10 @@ livestock_processors = 0 crop_operations_processors = 0 crop_fertilizers_processors = 0 agricultural_machinery_processors = 0 -residential_processors = 0 +residential_processors = 1 recreational_boats_processors = 0 point_sources_processors = 0 -traffic_processors = 1 +traffic_processors = 0 traffic_area_processors = 0 diff --git a/hermesv3_bu/sectors/agricultural_crop_fertilizers_sector.py b/hermesv3_bu/sectors/agricultural_crop_fertilizers_sector.py index 924d70d..b291236 100755 --- a/hermesv3_bu/sectors/agricultural_crop_fertilizers_sector.py +++ b/hermesv3_bu/sectors/agricultural_crop_fertilizers_sector.py @@ -15,7 +15,7 @@ formula = True class AgriculturalCropFertilizersSector(AgriculturalSector): - def __init__(self, comm_agr, comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, + def __init__(self, comm_agr, comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, crop_list, nut_shapefile, land_uses_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path, landuse_by_nut, crop_by_nut, crop_from_landuse_path, cultivated_ratio, fertilizer_rate, crop_f_parameter, crop_f_fertilizers, gridded_ph, gridded_cec, @@ -24,7 +24,7 @@ class AgriculturalCropFertilizersSector(AgriculturalSector): spent_time = timeit.default_timer() logger.write_log('===== AGRICULTURAL CROP FERTILIZERS SECTOR =====') super(AgriculturalCropFertilizersSector, self).__init__( - comm_agr, comm, logger, auxiliary_dir, grid_shp, clip, date_array, nut_shapefile, source_pollutants, + comm_agr, comm, logger, auxiliary_dir, grid, clip, date_array, nut_shapefile, source_pollutants, vertical_levels, crop_list, land_uses_path, landuse_by_nut, crop_by_nut, crop_from_landuse_path, None, None, None, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path) @@ -132,10 +132,10 @@ class AgriculturalCropFertilizersSector(AgriculturalSector): def to_dst_resolution(self, src_shapefile, value): spent_time = timeit.default_timer() - intersection = self.spatial_overlays(src_shapefile.to_crs(self.grid_shp.crs).reset_index(), - self.grid_shp.reset_index()) + intersection = self.spatial_overlays(src_shapefile.to_crs(self.grid.shapefile.crs).reset_index(), + self.grid.shapefile.reset_index()) intersection['area'] = intersection.geometry.area - dst_shapefile = self.grid_shp.reset_index().copy() + dst_shapefile = self.grid.shapefile.reset_index().copy() dst_shapefile['involved_area'] = intersection.groupby('FID')['area'].sum() intersection_with_dst_areas = pd.merge(intersection, dst_shapefile.loc[:, ['FID', 'involved_area']], how='left', on='FID') @@ -155,7 +155,7 @@ class AgriculturalCropFertilizersSector(AgriculturalSector): def to_dst_resolution_parallel(self, src_shapefile, index, value): spent_time = timeit.default_timer() - grid_shp = self.grid_shp.loc[index, :].copy() + grid_shp = self.grid.shapefile.loc[index, :].copy() src_shapefile = self.comm.bcast(src_shapefile, root=0) src_shapefile = src_shapefile.to_crs(grid_shp.crs) src_shapefile = src_shapefile[src_shapefile.within(grid_shp.unary_union)] diff --git a/hermesv3_bu/sectors/agricultural_crop_operations_sector.py b/hermesv3_bu/sectors/agricultural_crop_operations_sector.py index 9930f47..f456687 100755 --- a/hermesv3_bu/sectors/agricultural_crop_operations_sector.py +++ b/hermesv3_bu/sectors/agricultural_crop_operations_sector.py @@ -11,7 +11,7 @@ from hermesv3_bu.logger.log import Log class AgriculturalCropOperationsSector(AgriculturalSector): - def __init__(self, comm_agr, comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, + def __init__(self, comm_agr, comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, crop_list, nut_shapefile_path, land_uses_path, ef_dir, monthly_profiles_path, weekly_profiles_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path, landuse_by_nut, crop_by_nut, crop_from_landuse_path): @@ -90,7 +90,7 @@ class AgriculturalCropOperationsSector(AgriculturalSector): spent_time = timeit.default_timer() logger.write_log('===== AGRICULTURAL CROP OPERATIONS SECTOR =====') super(AgriculturalCropOperationsSector, self).__init__( - comm_agr, comm, logger, auxiliary_dir, grid_shp, clip, date_array, nut_shapefile_path, source_pollutants, + comm_agr, comm, logger, auxiliary_dir, grid, clip, date_array, nut_shapefile_path, source_pollutants, vertical_levels, crop_list, land_uses_path, landuse_by_nut, crop_by_nut, crop_from_landuse_path, ef_dir, monthly_profiles_path, weekly_profiles_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path) diff --git a/hermesv3_bu/sectors/agricultural_machinery_sector.py b/hermesv3_bu/sectors/agricultural_machinery_sector.py index 78c9a6e..7fb03bf 100755 --- a/hermesv3_bu/sectors/agricultural_machinery_sector.py +++ b/hermesv3_bu/sectors/agricultural_machinery_sector.py @@ -15,7 +15,7 @@ from hermesv3_bu.logger.log import Log class AgriculturalMachinerySector(AgriculturalSector): - def __init__(self, comm_agr, comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, + def __init__(self, comm_agr, comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, crop_list, nut_shapefile, machinery_list, land_uses_path, ef_files_dir, monthly_profiles_path, weekly_profiles_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path, landuse_by_nut, crop_by_nut, crop_from_landuse_path, @@ -26,7 +26,7 @@ class AgriculturalMachinerySector(AgriculturalSector): logger.write_log('===== AGRICULTURAL MACHINERY SECTOR =====') super(AgriculturalMachinerySector, self).__init__( - comm_agr, comm, logger, auxiliary_dir, grid_shp, clip, date_array, nut_shapefile, source_pollutants, + comm_agr, comm, logger, auxiliary_dir, grid, clip, date_array, nut_shapefile, source_pollutants, vertical_levels, crop_list, land_uses_path, landuse_by_nut, crop_by_nut, crop_from_landuse_path, ef_files_dir, monthly_profiles_path, weekly_profiles_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path) diff --git a/hermesv3_bu/sectors/agricultural_sector.py b/hermesv3_bu/sectors/agricultural_sector.py index 00f1973..31118ee 100755 --- a/hermesv3_bu/sectors/agricultural_sector.py +++ b/hermesv3_bu/sectors/agricultural_sector.py @@ -18,7 +18,7 @@ from pandas import DataFrame class AgriculturalSector(Sector): - def __init__(self, comm_agr, comm, logger, auxiliary_dir, grid_shp, clip, date_array, nut_shapefile, + def __init__(self, comm_agr, comm, logger, auxiliary_dir, grid, clip, date_array, nut_shapefile, source_pollutants, vertical_levels, crop_list, land_uses_path, land_use_by_nut, crop_by_nut, crop_from_landuse_path, ef_files_dir, monthly_profiles_path, weekly_profiles_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path): @@ -110,7 +110,7 @@ class AgriculturalSector(Sector): spent_time = timeit.default_timer() super(AgriculturalSector, self).__init__( - comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, monthly_profiles_path, weekly_profiles_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path) @@ -128,7 +128,7 @@ class AgriculturalSector(Sector): def involved_grid_cells(self, src_shp): spent_time = timeit.default_timer() - grid_shp = IoShapefile(self.comm).split_shapefile(self.grid_shp) + grid_shp = IoShapefile(self.comm).split_shapefile(self.grid.shapefile) src_union = src_shp.to_crs(grid_shp.crs).geometry.unary_union grid_shp = grid_shp.loc[grid_shp.intersects(src_union), :] @@ -418,10 +418,10 @@ class AgriculturalSector(Sector): spent_time = timeit.default_timer() crop_list = list(np.setdiff1d(crop_distribution.columns.values, ['NUT', 'geometry'])) - crop_distribution = crop_distribution.to_crs(self.grid_shp.crs) + crop_distribution = crop_distribution.to_crs(self.grid.shapefile.crs) crop_distribution['src_inter_fraction'] = crop_distribution.geometry.area - crop_distribution = self.spatial_overlays(crop_distribution.reset_index(), self.grid_shp.reset_index(), + crop_distribution = self.spatial_overlays(crop_distribution.reset_index(), self.grid.shapefile.reset_index(), how='intersection') crop_distribution['src_inter_fraction'] = \ crop_distribution.geometry.area / crop_distribution['src_inter_fraction'] @@ -431,8 +431,8 @@ class AgriculturalSector(Sector): crop_distribution = crop_distribution.loc[:, crop_list + ['FID']].groupby('FID').sum() - crop_distribution = gpd.GeoDataFrame(crop_distribution, crs=self.grid_shp.crs, - geometry=self.grid_shp.loc[crop_distribution.index, 'geometry']) + crop_distribution = gpd.GeoDataFrame(crop_distribution, crs=self.grid.shapefile.crs, + geometry=self.grid.shapefile.loc[crop_distribution.index, 'geometry']) crop_distribution.reset_index(inplace=True) crop_distribution.set_index('FID', inplace=True) diff --git a/hermesv3_bu/sectors/aviation_sector.py b/hermesv3_bu/sectors/aviation_sector.py index ab7c605..064a1f9 100755 --- a/hermesv3_bu/sectors/aviation_sector.py +++ b/hermesv3_bu/sectors/aviation_sector.py @@ -34,7 +34,7 @@ class AviationSector(Sector): - Taxi in - Post-taxi in """ - def __init__(self, comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + def __init__(self, comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, airport_list, plane_list, airport_shapefile_path, airport_runways_shapefile_path, airport_runways_corners_shapefile_path, airport_trajectories_shapefile_path, operations_path, planes_path, times_path, ef_dir, weekly_profiles_path, hourly_profiles_path, speciation_map_path, @@ -126,7 +126,7 @@ class AviationSector(Sector): spent_time = timeit.default_timer() logger.write_log('===== AVIATION SECTOR =====') super(AviationSector, self).__init__( - comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, None, + comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, None, weekly_profiles_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path) @@ -363,8 +363,8 @@ class AviationSector(Sector): spent_time = timeit.default_timer() if self.comm.Get_rank() == 0: airport_shapefile = airport_shapefile.reset_index() - airport_shapefile = gpd.sjoin(airport_shapefile.to_crs(self.grid_shp.crs), - self.clip.shapefile.to_crs(self.grid_shp.crs), how='inner', op='intersects') + airport_shapefile = gpd.sjoin(airport_shapefile.to_crs(self.grid.shapefile.crs), + self.clip.shapefile.to_crs(self.grid.shapefile.crs), how='inner', op='intersects') shp_airport_list = list(np.unique(airport_shapefile['airport_id'].values)) @@ -423,9 +423,9 @@ class AviationSector(Sector): airport_shapefile = airport_shapefile.loc[self.airport_list_full, :].copy() if not os.path.exists(os.path.dirname(airport_distribution_path)): os.makedirs(os.path.dirname(airport_distribution_path)) - airport_shapefile.to_crs(self.grid_shp.crs, inplace=True) + airport_shapefile.to_crs(self.grid.shapefile.crs, inplace=True) airport_shapefile['area'] = airport_shapefile.area - airport_distribution = self.spatial_overlays(airport_shapefile, self.grid_shp.reset_index(), + airport_distribution = self.spatial_overlays(airport_shapefile, self.grid.shapefile.reset_index(), how='intersection') airport_distribution['fraction'] = airport_distribution.area / airport_distribution['area'] airport_distribution.drop(columns=['idx2', 'area', 'geometry', 'cons'], inplace=True) @@ -482,13 +482,13 @@ class AviationSector(Sector): if self.comm.rank == 0: runway_shapefile['{0}_f'.format(phase_type)] = runway_shapefile.groupby('airport_id').apply(normalize) - runway_shapefile.to_crs(self.grid_shp.crs, inplace=True) + runway_shapefile.to_crs(self.grid.shapefile.crs, inplace=True) runway_shapefile['length'] = runway_shapefile.length # duplicating each runway by involved cell - runway_shapefile = gpd.sjoin(runway_shapefile.reset_index(), self.grid_shp.reset_index(), how="inner", + runway_shapefile = gpd.sjoin(runway_shapefile.reset_index(), self.grid.shapefile.reset_index(), how="inner", op='intersects') # Adding cell geometry - runway_shapefile = runway_shapefile.merge(self.grid_shp.reset_index(), on='FID', how='left') + runway_shapefile = runway_shapefile.merge(self.grid.shapefile.reset_index(), on='FID', how='left') # Intersection between line (roadway) and polygon (cell) # runway_shapefile['geometry'] = runway_shapefile.apply(do_intersection, axis=1) runway_shapefile['mini_length'] = runway_shapefile.apply(get_intersection_length, axis=1) @@ -594,8 +594,8 @@ class AviationSector(Sector): trajectories_distr.reset_index(inplace=True) # HORIZONTAL DISTRIBUTION - aux_grid = self.grid_shp.to_crs(trajectories_distr.crs).reset_index() - # trajectories_distr.to_crs(self.grid_shp.crs, inplace=True) + aux_grid = self.grid.shapefile.to_crs(trajectories_distr.crs).reset_index() + # trajectories_distr.to_crs(self.grid.shapefile.crs, inplace=True) # duplicating each runway by involved cell trajectories_distr = gpd.sjoin(trajectories_distr, aux_grid, how="inner", op='intersects') # Adding cell geometry diff --git a/hermesv3_bu/sectors/livestock_sector.py b/hermesv3_bu/sectors/livestock_sector.py index 8bf264d..bede311 100755 --- a/hermesv3_bu/sectors/livestock_sector.py +++ b/hermesv3_bu/sectors/livestock_sector.py @@ -22,7 +22,7 @@ class LivestockSector(Sector): """ Class that contains all the information and methods to calculate the livestock emissions. """ - def __init__(self, comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + def __init__(self, comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, animal_list, gridded_livestock_path, correction_split_factors_path, temperature_dir, wind_speed_dir, denominator_yearly_factor_dir, ef_dir, monthly_profiles_path, weekly_profiles_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path, @@ -154,7 +154,7 @@ class LivestockSector(Sector): spent_time = timeit.default_timer() logger.write_log('===== LIVESTOCK SECTOR =====') super(LivestockSector, self).__init__( - comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, monthly_profiles_path, weekly_profiles_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path) @@ -315,14 +315,14 @@ class LivestockSector(Sector): """ spent_time = timeit.default_timer() self.logger.write_log('\t\tCreating animal shapefile into destiny resolution', message_level=3) - self.grid_shp.reset_index(inplace=True) + self.grid.shapefile.reset_index(inplace=True) # Changing coordinates system to the grid one - animal_distribution.to_crs(self.grid_shp.crs, inplace=True) + animal_distribution.to_crs(self.grid.shapefile.crs, inplace=True) # Getting src area animal_distribution['src_inter_fraction'] = animal_distribution.geometry.area # Making the intersection between the src distribution and the destiny grid - animal_distribution = self.spatial_overlays(animal_distribution.reset_index(), self.grid_shp, + animal_distribution = self.spatial_overlays(animal_distribution.reset_index(), self.grid.shapefile, how='intersection') # Getting proportion of intersection in the src cell (src_area/portion_area) animal_distribution['src_inter_fraction'] = \ @@ -333,10 +333,10 @@ class LivestockSector(Sector): # Sum by destiny cell animal_distribution = animal_distribution.loc[:, self.animal_list + ['FID']].groupby('FID').sum() - self.grid_shp.set_index('FID', drop=False, inplace=True) + self.grid.shapefile.set_index('FID', drop=False, inplace=True) # Adding geometry and coordinates system from the destiny grid shapefile - animal_distribution = gpd.GeoDataFrame(animal_distribution, crs=self.grid_shp.crs, - geometry=self.grid_shp.loc[animal_distribution.index, 'geometry']) + animal_distribution = gpd.GeoDataFrame(animal_distribution, crs=self.grid.shapefile.crs, + geometry=self.grid.shapefile.loc[animal_distribution.index, 'geometry']) animal_distribution.reset_index(inplace=True) self.logger.write_time_log('LivestockSector', 'animals_shapefile_to_dst_resolution', timeit.default_timer() - spent_time) diff --git a/hermesv3_bu/sectors/point_source_sector.py b/hermesv3_bu/sectors/point_source_sector.py index a13d765..b9555a9 100755 --- a/hermesv3_bu/sectors/point_source_sector.py +++ b/hermesv3_bu/sectors/point_source_sector.py @@ -46,14 +46,14 @@ class PointSourceSector(Sector): :param sector_list: List os sectors (SNAPS) to take into account. 01, 03, 04, 09 :type sector_list: list """ - def __init__(self, comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + def __init__(self, comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, catalog_path, monthly_profiles_path, weekly_profiles_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, sector_list, measured_emission_path, molecular_weights_path, plume_rise=False, plume_rise_pahts=None): spent_time = timeit.default_timer() super(PointSourceSector, self).__init__( - comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, monthly_profiles_path, weekly_profiles_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path) @@ -856,9 +856,9 @@ class PointSourceSector(Sector): def point_source_to_fid(self, catalog): catalog.reset_index(inplace=True) - catalog = catalog.to_crs(self.grid_shp.crs) + catalog = catalog.to_crs(self.grid.shapefile.crs) - catalog = gpd.sjoin(catalog, self.grid_shp.reset_index(), how="inner", op='intersects') + catalog = gpd.sjoin(catalog, self.grid.shapefile.reset_index(), how="inner", op='intersects') # Drops duplicates when the point source is on the boundary of the cell catalog = catalog[~catalog.index.duplicated(keep='first')] diff --git a/hermesv3_bu/sectors/recreational_boats_sector.py b/hermesv3_bu/sectors/recreational_boats_sector.py index bbad747..462c7d6 100755 --- a/hermesv3_bu/sectors/recreational_boats_sector.py +++ b/hermesv3_bu/sectors/recreational_boats_sector.py @@ -14,14 +14,14 @@ from hermesv3_bu.logger.log import Log class RecreationalBoatsSector(Sector): - def __init__(self, comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + def __init__(self, comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, boat_list, density_map_path, boats_data_path, ef_file_path, monthly_profiles_path, weekly_profiles_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path): spent_time = timeit.default_timer() super(RecreationalBoatsSector, self).__init__( - comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, monthly_profiles_path, weekly_profiles_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path) @@ -42,9 +42,9 @@ class RecreationalBoatsSector(Sector): src_density_map = IoRaster(self.comm).to_shapefile_serie(density_map_path, nodata=0) src_density_map = src_density_map.loc[src_density_map['data'] > 0] src_density_map['data'] = src_density_map['data'] / src_density_map['data'].sum() - src_density_map.to_crs(self.grid_shp.crs, inplace=True) + src_density_map.to_crs(self.grid.shapefile.crs, inplace=True) src_density_map['src_inter_fraction'] = src_density_map.area - src_density_map = self.spatial_overlays(src_density_map, self.grid_shp.reset_index(), + src_density_map = self.spatial_overlays(src_density_map, self.grid.shapefile.reset_index(), how='intersection') src_density_map['src_inter_fraction'] = src_density_map.area / src_density_map['src_inter_fraction'] @@ -52,8 +52,8 @@ class RecreationalBoatsSector(Sector): axis="index") src_density_map = src_density_map.loc[:, ['FID', 'data']].groupby('FID').sum() - src_density_map = gpd.GeoDataFrame(src_density_map, crs=self.grid_shp.crs, - geometry=self.grid_shp.loc[src_density_map.index, 'geometry']) + src_density_map = gpd.GeoDataFrame(src_density_map, crs=self.grid.shapefile.crs, + geometry=self.grid.shapefile.loc[src_density_map.index, 'geometry']) src_density_map.reset_index(inplace=True) IoShapefile(self.comm).write_shapefile_serial(src_density_map, density_map_auxpath) diff --git a/hermesv3_bu/sectors/residential_sector.py b/hermesv3_bu/sectors/residential_sector.py index 495b5e7..4cf8316 100755 --- a/hermesv3_bu/sectors/residential_sector.py +++ b/hermesv3_bu/sectors/residential_sector.py @@ -16,7 +16,7 @@ from hermesv3_bu.logger.log import Log class ResidentialSector(Sector): - def __init__(self, comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + def __init__(self, comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, fuel_list, prov_shapefile, ccaa_shapefile, population_density_map, population_type_map, population_type_nuts2, population_type_nuts3, energy_consumption_nuts3, energy_consumption_nuts2, residential_spatial_proxies, residential_ef_files_path, @@ -25,7 +25,7 @@ class ResidentialSector(Sector): spent_time = timeit.default_timer() super(ResidentialSector, self).__init__( - comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, None, None, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path) @@ -124,11 +124,11 @@ class ResidentialSector(Sector): def to_dst_resolution(self, src_distribution): spent_time = timeit.default_timer() - src_distribution.to_crs(self.grid_shp.crs, inplace=True) + src_distribution.to_crs(self.grid.shapefile.crs, inplace=True) # src_distribution.reset_index().to_file( # os.path.join(self.auxiliary_dir, 'residential', 'fuel_distribution_src.shp')) src_distribution['src_inter_fraction'] = src_distribution.geometry.area - src_distribution = self.spatial_overlays(src_distribution, self.grid_shp.reset_index(), how='intersection') + src_distribution = self.spatial_overlays(src_distribution, self.grid.shapefile.reset_index(), how='intersection') # src_distribution.reset_index().to_file( # os.path.join(self.auxiliary_dir, 'residential', 'fuel_distribution_raw.shp')) src_distribution['src_inter_fraction'] = src_distribution.geometry.area / src_distribution[ @@ -138,8 +138,8 @@ class ResidentialSector(Sector): src_distribution["src_inter_fraction"], axis="index") src_distribution = src_distribution.loc[:, self.fuel_list + ['FID']].groupby('FID').sum() - src_distribution = gpd.GeoDataFrame(src_distribution, crs=self.grid_shp.crs, - geometry=self.grid_shp.loc[src_distribution.index, 'geometry']) + src_distribution = gpd.GeoDataFrame(src_distribution, crs=self.grid.shapefile.crs, + geometry=self.grid.shapefile.loc[src_distribution.index, 'geometry']) src_distribution.reset_index(inplace=True) self.logger.write_time_log('ResidentialSector', 'to_dst_resolution', timeit.default_timer() - spent_time) diff --git a/hermesv3_bu/sectors/sector.py b/hermesv3_bu/sectors/sector.py index 8aae0b5..f49e379 100755 --- a/hermesv3_bu/sectors/sector.py +++ b/hermesv3_bu/sectors/sector.py @@ -9,11 +9,12 @@ import pandas as pd import geopandas as gpd from geopandas import GeoDataFrame from mpi4py import MPI +from hermesv3_bu.grids.grid import Grid class Sector(object): - def __init__(self, comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + def __init__(self, comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, monthly_profiles_path, weekly_profiles_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path): """ @@ -29,8 +30,8 @@ class Sector(object): created yet. :type auxiliary_dir: str - :param grid_shp: Shapefile with the grid horizontal distribution. - :type grid_shp: GeoDataFrame + :param grid: Grid object + :type grid: Grid :param date_array: List of datetimes. :type date_array: list(datetime.datetime, ...) @@ -74,7 +75,7 @@ class Sector(object): self.comm = comm self.logger = logger self.auxiliary_dir = auxiliary_dir - self.grid_shp = grid_shp + self.grid = grid self.clip = clip self.date_array = date_array self.source_pollutants = source_pollutants diff --git a/hermesv3_bu/sectors/sector_manager.py b/hermesv3_bu/sectors/sector_manager.py index 048eeb6..be6fe1b 100755 --- a/hermesv3_bu/sectors/sector_manager.py +++ b/hermesv3_bu/sectors/sector_manager.py @@ -38,7 +38,7 @@ class SectorManager(object): from hermesv3_bu.sectors.aviation_sector import AviationSector self.sector = AviationSector( comm_world.Split(color, sector_procs.index(comm_world.Get_rank())), self.logger, - arguments.auxiliary_files_path, grid.shapefile, clip, date_array, + arguments.auxiliary_files_path, grid, clip, date_array, arguments.aviation_source_pollutants, grid.vertical_desctiption, arguments.airport_list, arguments.plane_list, arguments.airport_shapefile_path, arguments.airport_runways_shapefile_path, arguments.airport_runways_corners_shapefile_path, arguments.airport_trajectories_shapefile_path, @@ -50,7 +50,7 @@ class SectorManager(object): from hermesv3_bu.sectors.shipping_port_sector import ShippingPortSector self.sector = ShippingPortSector( comm_world.Split(color, sector_procs.index(comm_world.Get_rank())), self.logger, - arguments.auxiliary_files_path, grid.shapefile, clip, date_array, + arguments.auxiliary_files_path, grid, clip, date_array, arguments.shipping_port_source_pollutants, grid.vertical_desctiption, arguments.vessel_list, arguments.port_list, arguments.hoteling_shapefile_path, arguments.maneuvering_shapefile_path, arguments.shipping_port_ef_path, arguments.shipping_port_engine_percent_path, @@ -63,7 +63,7 @@ class SectorManager(object): from hermesv3_bu.sectors.livestock_sector import LivestockSector self.sector = LivestockSector( comm_world.Split(color, sector_procs.index(comm_world.Get_rank())), self.logger, - arguments.auxiliary_files_path, grid.shapefile, clip, date_array, + arguments.auxiliary_files_path, grid, clip, date_array, arguments.livestock_source_pollutants, grid.vertical_desctiption, arguments.animal_list, arguments.gridded_livestock, arguments.correction_split_factors, arguments.temperature_daily_files_path, arguments.wind_speed_daily_files_path, @@ -78,7 +78,7 @@ class SectorManager(object): comm_agr = comm_world.Split(agg_color, agg_procs.index(comm_world.Get_rank())) comm = comm_agr.Split(color, sector_procs.index(comm_world.Get_rank())) self.sector = AgriculturalCropOperationsSector( - comm_agr, comm, logger, arguments.auxiliary_files_path, grid.shapefile, clip, date_array, + comm_agr, comm, logger, arguments.auxiliary_files_path, grid, clip, date_array, arguments.crop_operations_source_pollutants, grid.vertical_desctiption, arguments.crop_operations_list, arguments.nuts2_shapefile, arguments.land_uses_path, arguments.crop_operations_ef_files_dir, @@ -93,7 +93,7 @@ class SectorManager(object): comm_agr = comm_world.Split(agg_color, agg_procs.index(comm_world.Get_rank())) comm = comm_agr.Split(color, sector_procs.index(comm_world.Get_rank())) self.sector = AgriculturalCropFertilizersSector( - comm_agr, comm, logger, arguments.auxiliary_files_path, grid.shapefile, clip, date_array, + comm_agr, comm, logger, arguments.auxiliary_files_path, grid, clip, date_array, arguments.crop_fertilizers_source_pollutants, grid.vertical_desctiption, arguments.crop_fertilizers_list, arguments.nuts2_shapefile, arguments.land_uses_path, arguments.crop_fertilizers_hourly_profiles, arguments.speciation_map, @@ -111,7 +111,7 @@ class SectorManager(object): comm_agr = comm_world.Split(agg_color, agg_procs.index(comm_world.Get_rank())) comm = comm_agr.Split(color, sector_procs.index(comm_world.Get_rank())) self.sector = AgriculturalMachinerySector( - comm_agr, comm, logger, arguments.auxiliary_files_path, grid.shapefile, clip, date_array, + comm_agr, comm, logger, arguments.auxiliary_files_path, grid, clip, date_array, arguments.crop_machinery_source_pollutants, grid.vertical_desctiption, arguments.crop_machinery_list, arguments.nuts2_shapefile, arguments.machinery_list, arguments.land_uses_path, arguments.crop_machinery_ef_path, @@ -128,7 +128,7 @@ class SectorManager(object): from hermesv3_bu.sectors.residential_sector import ResidentialSector self.sector = ResidentialSector( comm_world.Split(color, sector_procs.index(comm_world.Get_rank())), self.logger, - arguments.auxiliary_files_path, grid.shapefile, clip, date_array, + arguments.auxiliary_files_path, grid, clip, date_array, arguments.residential_source_pollutants, grid.vertical_desctiption, arguments.fuel_list, arguments.nuts3_shapefile, arguments.nuts2_shapefile, arguments.population_density_map, arguments.population_type_map, arguments.population_type_nuts2, arguments.population_type_nuts3, @@ -142,7 +142,7 @@ class SectorManager(object): from hermesv3_bu.sectors.recreational_boats_sector import RecreationalBoatsSector self.sector = RecreationalBoatsSector( comm_world.Split(color, sector_procs.index(comm_world.Get_rank())), self.logger, - arguments.auxiliary_files_path, grid.shapefile, clip, date_array, + arguments.auxiliary_files_path, grid, clip, date_array, arguments.recreational_boats_source_pollutants, grid.vertical_desctiption, arguments.recreational_boats_list, arguments.recreational_boats_density_map, arguments.recreational_boats_by_type, arguments.recreational_boats_ef_path, @@ -154,7 +154,7 @@ class SectorManager(object): from hermesv3_bu.sectors.point_source_sector import PointSourceSector self.sector = PointSourceSector( comm_world.Split(color, sector_procs.index(comm_world.Get_rank())), self.logger, - arguments.auxiliary_files_path, grid.shapefile, clip, date_array, + arguments.auxiliary_files_path, grid, clip, date_array, arguments.point_source_pollutants, grid.vertical_desctiption, arguments.point_source_catalog, arguments.point_source_monthly_profiles, arguments.point_source_weekly_profiles, arguments.point_source_hourly_profiles, @@ -175,7 +175,7 @@ class SectorManager(object): from hermesv3_bu.sectors.traffic_sector import TrafficSector self.sector = TrafficSector( comm_world.Split(color, sector_procs.index(comm_world.Get_rank())), self.logger, - arguments.auxiliary_files_path, grid.shapefile, clip, date_array, arguments.traffic_pollutants, + arguments.auxiliary_files_path, grid, clip, date_array, arguments.traffic_pollutants, grid.vertical_desctiption, arguments.road_link_path, arguments.fleet_compo_path, arguments.traffic_speed_hourly_path, arguments.traffic_monthly_profiles, arguments.traffic_weekly_profiles, arguments.traffic_hourly_profiles_mean, @@ -193,7 +193,7 @@ class SectorManager(object): from hermesv3_bu.sectors.traffic_area_sector import TrafficAreaSector self.sector = TrafficAreaSector( comm_world.Split(color, sector_procs.index(comm_world.Get_rank())), self.logger, - arguments.auxiliary_files_path, grid.shapefile, clip, date_array, arguments.traffic_area_pollutants, + arguments.auxiliary_files_path, grid, clip, date_array, arguments.traffic_area_pollutants, grid.vertical_desctiption, arguments.population_density_map, arguments.speciation_map, arguments.molecular_weights, arguments.do_evaporative, arguments.traffic_area_gas_path, arguments.population_nuts3, arguments.nuts3_shapefile, diff --git a/hermesv3_bu/sectors/shipping_port_sector.py b/hermesv3_bu/sectors/shipping_port_sector.py index 316959c..0e316c2 100755 --- a/hermesv3_bu/sectors/shipping_port_sector.py +++ b/hermesv3_bu/sectors/shipping_port_sector.py @@ -10,13 +10,13 @@ from hermesv3_bu.io_server.io_shapefile import IoShapefile class ShippingPortSector(Sector): - def __init__(self, comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + def __init__(self, comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, vessel_list, port_list, hoteling_shapefile_path, maneuvering_shapefile_path, ef_dir, engine_percent_path, tonnage_path, load_factor_path, power_path, monthly_profiles_path, weekly_profiles_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path): """ - Initialise the Shipping port sectopr class + Initialise the Shipping port sector class :param comm: Communicator for the sector calculation. :type comm: MPI.COMM @@ -83,7 +83,7 @@ class ShippingPortSector(Sector): logger.write_log('===== SHIPPING PORT SECTOR =====') super(ShippingPortSector, self).__init__( - comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, monthly_profiles_path, weekly_profiles_path, hourly_profiles_path, speciation_map_path, speciation_profiles_path, molecular_weights_path) @@ -579,9 +579,9 @@ class ShippingPortSector(Sector): dataframe.reset_index(inplace=True) dataframe.drop(columns=['code'], inplace=True) - dataframe.to_crs(self.grid_shp.crs, inplace=True) + dataframe.to_crs(self.grid.shapefile.crs, inplace=True) dataframe['src_inter_fraction'] = dataframe.geometry.area - dataframe = self.spatial_overlays(dataframe, self.grid_shp, how='intersection') + dataframe = self.spatial_overlays(dataframe, self.grid.shapefile, how='intersection') dataframe['src_inter_fraction'] = dataframe.geometry.area / dataframe['src_inter_fraction'] dataframe[self.source_pollutants] = dataframe[self.source_pollutants].multiply(dataframe["src_inter_fraction"], diff --git a/hermesv3_bu/sectors/traffic_area_sector.py b/hermesv3_bu/sectors/traffic_area_sector.py index 64db3da..3047954 100755 --- a/hermesv3_bu/sectors/traffic_area_sector.py +++ b/hermesv3_bu/sectors/traffic_area_sector.py @@ -14,7 +14,7 @@ pmc_list = ['pmc', 'PMC'] class TrafficAreaSector(Sector): - def __init__(self, comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + def __init__(self, comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, population_tif_path, speciation_map_path, molecular_weights_path, do_evaporative, gasoline_path, total_pop_by_prov, nuts_shapefile, speciation_profiles_evaporative, evaporative_ef_file, temperature_dir, @@ -24,7 +24,7 @@ class TrafficAreaSector(Sector): logger.write_log('===== TRAFFIC AREA SECTOR =====') super(TrafficAreaSector, self).__init__( - comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, None, None, None, speciation_map_path, None, molecular_weights_path) self.do_evaporative = do_evaporative @@ -152,14 +152,14 @@ class TrafficAreaSector(Sector): if not os.path.exists(pop_nut_cell_path): - pop_by_nut = pop_by_nut.to_crs(self.grid_shp.crs) + pop_by_nut = pop_by_nut.to_crs(self.grid.shapefile.crs) # del pop_by_nut[''] pop_by_nut['area_in'] = pop_by_nut.geometry.area # df = gpd.overlay(pop_by_nut, grid_shp, how='intersection') - df = self.spatial_overlays(pop_by_nut, self.grid_shp.reset_index(), how='intersection') + df = self.spatial_overlays(pop_by_nut, self.grid.shapefile.reset_index(), how='intersection') - df.crs = self.grid_shp.crs + df.crs = self.grid.shapefile.crs df.loc[:, 'data'] = df['data'] * (df.geometry.area / df['area_in']) del pop_by_nut['area_in'] if write_file: @@ -198,7 +198,7 @@ class TrafficAreaSector(Sector): aux_df = df.loc[:, ['FID'] + vehicle_type_list].groupby('FID').sum() aux_df.loc[:, 'FID'] = aux_df.index - geom = self.grid_shp.loc[aux_df.index, 'geometry'] + geom = self.grid.shapefile.loc[aux_df.index, 'geometry'] df = gpd.GeoDataFrame(aux_df, geometry=geom, crs=pop_nut_cell.crs) IoShapefile(self.comm).write_shapefile_serial(df, veh_by_cell_path) @@ -384,7 +384,7 @@ class TrafficAreaSector(Sector): p_names = small_cities.columns.values - aux_grid = self.grid_shp.loc[small_cities.index.values, :].reset_index().copy() + aux_grid = self.grid.shapefile.loc[small_cities.index.values, :].reset_index().copy() aux_grid = self.add_timezone(aux_grid) aux_grid.set_index('FID', inplace=True) diff --git a/hermesv3_bu/sectors/traffic_sector.py b/hermesv3_bu/sectors/traffic_sector.py index a8b2f96..2fb9e4d 100755 --- a/hermesv3_bu/sectors/traffic_sector.py +++ b/hermesv3_bu/sectors/traffic_sector.py @@ -21,6 +21,7 @@ libc.malloc_trim(0) MIN_RAIN = 0.254 # After USEPA (2011) RECOVERY_RATIO = 0.0872 # After Amato et al. (2012) +FINAL_PROJ = {'init': 'epsg:3035'} aerosols = ['oc', 'ec', 'pno3', 'pso4', 'pmfine', 'pmc', 'poa', 'poc', 'pec', 'pcl', 'pnh4', 'pna', 'pmg', 'pk', 'pca', @@ -43,7 +44,7 @@ class TrafficSector(Sector): relative to the timesteps. """ - def __init__(self, comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + def __init__(self, comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, road_link_path, fleet_compo_path, speed_hourly_path, monthly_profiles_path, weekly_profiles_path, hourly_mean_profiles_path, hourly_weekday_profiles_path, hourly_saturday_profiles_path, hourly_sunday_profiles_path, ef_common_path, vehicle_list=None, load=0.5, speciation_map_path=None, @@ -55,7 +56,7 @@ class TrafficSector(Sector): spent_time = timeit.default_timer() logger.write_log('===== TRAFFIC SECTOR =====') super(TrafficSector, self).__init__( - comm, logger, auxiliary_dir, grid_shp, clip, date_array, source_pollutants, vertical_levels, + comm, logger, auxiliary_dir, grid, clip, date_array, source_pollutants, vertical_levels, monthly_profiles_path, weekly_profiles_path, None, speciation_map_path, None, molecular_weights_path) self.resuspension_correction = resuspension_correction @@ -363,6 +364,7 @@ class TrafficSector(Sector): self.logger.write_time_log('TrafficSector', 'read_road_links', timeit.default_timer() - spent_time) libc.malloc_trim(0) + return df def read_ef(self, emission_type, pollutant_name): @@ -1253,14 +1255,18 @@ class TrafficSector(Sector): if not os.path.exists(self.link_to_grid_csv): link_emissions_aux = link_emissions.loc[link_emissions['tstep'] == 0, :] - link_emissions_aux = link_emissions_aux.to_crs(self.grid_shp.crs) + if self.grid.grid_type in ['Lambert Conformal Conic', 'Mercator']: + grid_aux = self.grid.shapefile + else: + grid_aux = self.grid.shapefile.to_crs(FINAL_PROJ) + + link_emissions_aux = link_emissions_aux.to_crs(grid_aux.shapefile.crs) - link_emissions_aux = gpd.sjoin(link_emissions_aux, self.grid_shp.reset_index(), - how="inner", op='intersects') + link_emissions_aux = gpd.sjoin(link_emissions_aux, grid_aux.reset_index(), how="inner", op='intersects') link_emissions_aux = link_emissions_aux.loc[:, ['Link_ID', 'geometry', 'FID']] - link_emissions_aux = link_emissions_aux.merge(self.grid_shp.reset_index().loc[:, ['FID', 'geometry']], + link_emissions_aux = link_emissions_aux.merge(grid_aux.reset_index().loc[:, ['FID', 'geometry']], on='FID', how='left') length_list = [] @@ -1366,7 +1372,6 @@ 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: - gpd.GeoDataFrame().to_file df_in.to_file(os.path.join(self.output_dir, 'roads.shp')) count = 0 -- GitLab From 1addb23332ce624e59100025dab158d9094026a7 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Tue, 27 Aug 2019 11:28:35 +0200 Subject: [PATCH 05/10] -Traffic: Forced to espg: 3035 the rotated and regular grids to calculate the link length in meters - All sectors: Pass Grid object instead grid shapefile --- hermesv3_bu/sectors/traffic_sector.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hermesv3_bu/sectors/traffic_sector.py b/hermesv3_bu/sectors/traffic_sector.py index 2fb9e4d..5008f7f 100755 --- a/hermesv3_bu/sectors/traffic_sector.py +++ b/hermesv3_bu/sectors/traffic_sector.py @@ -1259,7 +1259,7 @@ class TrafficSector(Sector): grid_aux = self.grid.shapefile else: grid_aux = self.grid.shapefile.to_crs(FINAL_PROJ) - + link_emissions_aux = link_emissions_aux.to_crs(grid_aux.shapefile.crs) link_emissions_aux = gpd.sjoin(link_emissions_aux, grid_aux.reset_index(), how="inner", op='intersects') -- GitLab From 10adbc8d9ec3a34917e5d2370d87b4ec02389392 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Tue, 27 Aug 2019 11:36:54 +0200 Subject: [PATCH 06/10] PEP8 Corrections --- conf/hermes.conf | 6 +++--- hermesv3_bu/sectors/aviation_sector.py | 7 ++++--- hermesv3_bu/sectors/residential_sector.py | 3 ++- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/conf/hermes.conf b/conf/hermes.conf index e11afec..bd0ccbe 100755 --- a/conf/hermes.conf +++ b/conf/hermes.conf @@ -3,7 +3,7 @@ log_level = 3 input_dir = /home/Earth/ctena/Models/hermesv3_bu_data data_path = /esarchive/recon output_dir = /scratch/Earth/HERMESv3_BU_OUT/serial -output_name = HERMESv3__residential.nc +output_name = HERMESv3__point_sources.nc emission_summary = 0 start_date = 2016/11/29 00:00:00 # ----- end_date = start_date [DEFAULT] ----- @@ -136,9 +136,9 @@ livestock_processors = 0 crop_operations_processors = 0 crop_fertilizers_processors = 0 agricultural_machinery_processors = 0 -residential_processors = 1 +residential_processors = 0 recreational_boats_processors = 0 -point_sources_processors = 0 +point_sources_processors = 1 traffic_processors = 0 traffic_area_processors = 0 diff --git a/hermesv3_bu/sectors/aviation_sector.py b/hermesv3_bu/sectors/aviation_sector.py index 064a1f9..399efd3 100755 --- a/hermesv3_bu/sectors/aviation_sector.py +++ b/hermesv3_bu/sectors/aviation_sector.py @@ -364,7 +364,8 @@ class AviationSector(Sector): if self.comm.Get_rank() == 0: airport_shapefile = airport_shapefile.reset_index() airport_shapefile = gpd.sjoin(airport_shapefile.to_crs(self.grid.shapefile.crs), - self.clip.shapefile.to_crs(self.grid.shapefile.crs), how='inner', op='intersects') + self.clip.shapefile.to_crs(self.grid.shapefile.crs), how='inner', + op='intersects') shp_airport_list = list(np.unique(airport_shapefile['airport_id'].values)) @@ -485,8 +486,8 @@ class AviationSector(Sector): runway_shapefile.to_crs(self.grid.shapefile.crs, inplace=True) runway_shapefile['length'] = runway_shapefile.length # duplicating each runway by involved cell - runway_shapefile = gpd.sjoin(runway_shapefile.reset_index(), self.grid.shapefile.reset_index(), how="inner", - op='intersects') + runway_shapefile = gpd.sjoin(runway_shapefile.reset_index(), self.grid.shapefile.reset_index(), + how="inner", op='intersects') # Adding cell geometry runway_shapefile = runway_shapefile.merge(self.grid.shapefile.reset_index(), on='FID', how='left') # Intersection between line (roadway) and polygon (cell) diff --git a/hermesv3_bu/sectors/residential_sector.py b/hermesv3_bu/sectors/residential_sector.py index 4cf8316..17f6904 100755 --- a/hermesv3_bu/sectors/residential_sector.py +++ b/hermesv3_bu/sectors/residential_sector.py @@ -128,7 +128,8 @@ class ResidentialSector(Sector): # src_distribution.reset_index().to_file( # os.path.join(self.auxiliary_dir, 'residential', 'fuel_distribution_src.shp')) src_distribution['src_inter_fraction'] = src_distribution.geometry.area - src_distribution = self.spatial_overlays(src_distribution, self.grid.shapefile.reset_index(), how='intersection') + src_distribution = self.spatial_overlays(src_distribution, self.grid.shapefile.reset_index(), + how='intersection') # src_distribution.reset_index().to_file( # os.path.join(self.auxiliary_dir, 'residential', 'fuel_distribution_raw.shp')) src_distribution['src_inter_fraction'] = src_distribution.geometry.area / src_distribution[ -- GitLab From 1f9dbba8a61a2d0df8ebfd9a9189780347071a46 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Tue, 27 Aug 2019 12:36:34 +0200 Subject: [PATCH 07/10] PEP8 Corrections --- conf/hermes.conf | 8 ++++---- hermesv3_bu/sectors/traffic_sector.py | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/conf/hermes.conf b/conf/hermes.conf index bd0ccbe..d84e469 100755 --- a/conf/hermes.conf +++ b/conf/hermes.conf @@ -3,7 +3,7 @@ log_level = 3 input_dir = /home/Earth/ctena/Models/hermesv3_bu_data data_path = /esarchive/recon output_dir = /scratch/Earth/HERMESv3_BU_OUT/serial -output_name = HERMESv3__point_sources.nc +output_name = HERMESv3__traffic.nc emission_summary = 0 start_date = 2016/11/29 00:00:00 # ----- end_date = start_date [DEFAULT] ----- @@ -138,8 +138,8 @@ crop_fertilizers_processors = 0 agricultural_machinery_processors = 0 residential_processors = 0 recreational_boats_processors = 0 -point_sources_processors = 1 -traffic_processors = 0 +point_sources_processors = 0 +traffic_processors = 1 traffic_area_processors = 0 @@ -312,7 +312,7 @@ resuspension_correction = 1 write_rline = 0 traffic_pollutants = nox_no2, nh3, co, so2, pm, voc, ch4 -# vehicle_types = PCD_13 PCD_14 PCD_15 PCG_25 PCG_26 PCG_27 +vehicle_types = PCD_13 PCD_14 PCD_15 PCG_25 PCG_26 PCG_27 # load = [0, 0.5, 1] load = 0.5 road_link_path = /traffic/road_links/2015/road_links_2015.shp diff --git a/hermesv3_bu/sectors/traffic_sector.py b/hermesv3_bu/sectors/traffic_sector.py index 5008f7f..5bc53aa 100755 --- a/hermesv3_bu/sectors/traffic_sector.py +++ b/hermesv3_bu/sectors/traffic_sector.py @@ -1260,7 +1260,7 @@ class TrafficSector(Sector): else: grid_aux = self.grid.shapefile.to_crs(FINAL_PROJ) - link_emissions_aux = link_emissions_aux.to_crs(grid_aux.shapefile.crs) + link_emissions_aux = link_emissions_aux.to_crs(grid_aux.crs) link_emissions_aux = gpd.sjoin(link_emissions_aux, grid_aux.reset_index(), how="inner", op='intersects') -- GitLab From a8a2238345a852953021e3d7aa3b98ba8c798126 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Tue, 27 Aug 2019 14:34:00 +0200 Subject: [PATCH 08/10] - Traffic: Use metric projection for rotated amd regular projections --- conf/hermes.conf | 6 +++--- hermesv3_bu/sectors/traffic_sector.py | 4 +++- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/conf/hermes.conf b/conf/hermes.conf index d84e469..683290f 100755 --- a/conf/hermes.conf +++ b/conf/hermes.conf @@ -3,7 +3,7 @@ log_level = 3 input_dir = /home/Earth/ctena/Models/hermesv3_bu_data data_path = /esarchive/recon output_dir = /scratch/Earth/HERMESv3_BU_OUT/serial -output_name = HERMESv3__traffic.nc +output_name = HERMESv3__traffic_area.nc emission_summary = 0 start_date = 2016/11/29 00:00:00 # ----- end_date = start_date [DEFAULT] ----- @@ -139,8 +139,8 @@ agricultural_machinery_processors = 0 residential_processors = 0 recreational_boats_processors = 0 point_sources_processors = 0 -traffic_processors = 1 -traffic_area_processors = 0 +traffic_processors = 0 +traffic_area_processors = 1 [SHAPEFILES] diff --git a/hermesv3_bu/sectors/traffic_sector.py b/hermesv3_bu/sectors/traffic_sector.py index 5bc53aa..73be4f9 100755 --- a/hermesv3_bu/sectors/traffic_sector.py +++ b/hermesv3_bu/sectors/traffic_sector.py @@ -21,7 +21,7 @@ libc.malloc_trim(0) MIN_RAIN = 0.254 # After USEPA (2011) RECOVERY_RATIO = 0.0872 # After Amato et al. (2012) -FINAL_PROJ = {'init': 'epsg:3035'} +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', @@ -1258,6 +1258,8 @@ class TrafficSector(Sector): if self.grid.grid_type in ['Lambert Conformal Conic', 'Mercator']: grid_aux = self.grid.shapefile else: + # For REGULAR and ROTATED grids, shapefile projection is trasnfromed to a metric projected coordinate + # system to derive the length in km. grid_aux = self.grid.shapefile.to_crs(FINAL_PROJ) link_emissions_aux = link_emissions_aux.to_crs(grid_aux.crs) -- GitLab From d00e0604fc287b33bc9ff21e57a10b20aec7a9f8 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Tue, 27 Aug 2019 14:34:34 +0200 Subject: [PATCH 09/10] - Traffic: Use metric projection for rotated amd regular projections comment --- hermesv3_bu/sectors/traffic_sector.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hermesv3_bu/sectors/traffic_sector.py b/hermesv3_bu/sectors/traffic_sector.py index 73be4f9..e98843e 100755 --- a/hermesv3_bu/sectors/traffic_sector.py +++ b/hermesv3_bu/sectors/traffic_sector.py @@ -1258,7 +1258,7 @@ class TrafficSector(Sector): if self.grid.grid_type in ['Lambert Conformal Conic', 'Mercator']: grid_aux = self.grid.shapefile else: - # For REGULAR and ROTATED grids, shapefile projection is trasnfromed to a metric projected coordinate + # For REGULAR and ROTATED grids, shapefile projection is transfromed to a metric projected coordinate # system to derive the length in km. grid_aux = self.grid.shapefile.to_crs(FINAL_PROJ) -- GitLab From c29eba665cc624608fc7cfa83870b73d76c625cf Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Tue, 27 Aug 2019 14:34:57 +0200 Subject: [PATCH 10/10] - Traffic: Use metric projection for rotated amd regular projections comment --- hermesv3_bu/sectors/traffic_sector.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hermesv3_bu/sectors/traffic_sector.py b/hermesv3_bu/sectors/traffic_sector.py index e98843e..41fcb2f 100755 --- a/hermesv3_bu/sectors/traffic_sector.py +++ b/hermesv3_bu/sectors/traffic_sector.py @@ -1258,7 +1258,7 @@ class TrafficSector(Sector): if self.grid.grid_type in ['Lambert Conformal Conic', 'Mercator']: grid_aux = self.grid.shapefile else: - # For REGULAR and ROTATED grids, shapefile projection is transfromed to a metric projected coordinate + # For REGULAR and ROTATED grids, shapefile projection is transformed to a metric projected coordinate # system to derive the length in km. grid_aux = self.grid.shapefile.to_crs(FINAL_PROJ) -- GitLab