From 2fd02bfd82d5fcddfce029f238fb9813f502c35d Mon Sep 17 00:00:00 2001 From: ctena Date: Wed, 9 Mar 2022 17:10:59 +0100 Subject: [PATCH 1/4] Added plume_rise_output option --- CHANGELOG | 5 +++++ conf/hermes.conf | 1 + hermesv3_bu/__init__.py | 2 +- hermesv3_bu/config/config.py | 6 ++++++ hermesv3_bu/sectors/point_source_sector.py | 17 ++++++++++++++++- hermesv3_bu/sectors/sector_manager.py | 6 +++++- 6 files changed, 34 insertions(+), 3 deletions(-) diff --git a/CHANGELOG b/CHANGELOG index 9a76d15..ddbbbb6 100755 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,3 +1,8 @@ +1.0.7 + 2022/03/09 + + - Point Source plume rise output as CSV added + 1.0.6 2022/02/18 diff --git a/conf/hermes.conf b/conf/hermes.conf index 37fd069..e0be66e 100755 --- a/conf/hermes.conf +++ b/conf/hermes.conf @@ -293,6 +293,7 @@ recreational_boats_speciation_profiles = /profiles/speciation/recreat [POINT SOURCES] point_source_pollutants = nox_no2,nmvoc,so2,co,nh3,pm10,pm25,ch4,n2o,co2 plume_rise = False +plume_rise_output = True # point_source_snaps = 09 point_source_catalog = /point_sources/Maestra_focos_2015_plume_rise.shp point_source_monthly_profiles = /profiles/temporal/point_sources/monthly_profiles.csv diff --git a/hermesv3_bu/__init__.py b/hermesv3_bu/__init__.py index 382021f..9e604c0 100755 --- a/hermesv3_bu/__init__.py +++ b/hermesv3_bu/__init__.py @@ -1 +1 @@ -__version__ = "1.0.6" +__version__ = "1.0.7" diff --git a/hermesv3_bu/config/config.py b/hermesv3_bu/config/config.py index 9e6e07a..2176f95 100755 --- a/hermesv3_bu/config/config.py +++ b/hermesv3_bu/config/config.py @@ -534,6 +534,8 @@ class Config(ArgParser): help="List of pollutants considered for the calculation of the point sources sector.") p.add_argument('--plume_rise', required=False, help="Boolean that defines if the plume rise algorithm is activated or not.") + p.add_argument('--plume_rise_output', required=False, type=str, default='False', + help="Boolean that defines if you want to produce the plume-rise output as CSV.") p.add_argument('--point_source_snaps', required=False, help="Defines the SNAP source categories considered during the emission calculation " + "[01, 03, 04, 09].") @@ -784,6 +786,10 @@ class Config(ArgParser): # Point Source bools arguments.plume_rise = self._parse_bool(arguments.plume_rise) + if arguments.plume_rise: + arguments.plume_rise_output = self._parse_bool(arguments.plume_rise_output) + else: + arguments.plume_rise_output = False # Point Source lists arguments.point_source_pollutants = self._parse_list(arguments.point_source_pollutants) diff --git a/hermesv3_bu/sectors/point_source_sector.py b/hermesv3_bu/sectors/point_source_sector.py index 6a404ba..3089cfe 100755 --- a/hermesv3_bu/sectors/point_source_sector.py +++ b/hermesv3_bu/sectors/point_source_sector.py @@ -52,7 +52,7 @@ class PointSourceSector(Sector): 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): + molecular_weights_path, plume_rise=False, plume_rise_filename=False, plume_rise_pahts=None): spent_time = timeit.default_timer() logger.write_log('===== POINT SOURCES SECTOR =====') check_files( @@ -64,6 +64,7 @@ class PointSourceSector(Sector): speciation_profiles_path, molecular_weights_path) self.plume_rise = plume_rise + self.plume_rise_filename = plume_rise_filename self.catalog = self.read_catalog_shapefile(catalog_path, sector_list) self.check_catalog() self.catalog_measured = self.read_catalog_for_measured_emissions(catalog_path, sector_list) @@ -755,12 +756,26 @@ class PointSourceSector(Sector): # Step 4: Plume rise catalog['h_top'] = (1.5 * catalog['Ah']) + catalog['Height'] catalog['h_bot'] = (0.5 * catalog['Ah']) + catalog['Height'] + if self.plume_rise_filename is not None: + self.write_plume_rise_output(catalog) catalog.drop(columns=['Height', 'Diameter', 'Speed', 'Temp', 'date_utc', 'temp_sfc', 'friction_v', 'pbl', 'obukhov_len', 'temp_top', 'wSpeed_top', 'Fb', 'S', 'Ah'], inplace=True) self.logger.write_time_log('PointSourceSector', 'get_plume_rise_top_bot', timeit.default_timer() - spent_time) return catalog + def write_plume_rise_output(self, catalog): + catalog_aux = catalog.drop( + columns=['index', 'Height', 'Diameter', 'Speed', 'Temp', 'P_month', 'P_week', 'P_hour', 'P_spec', + 'geometry', 'FID', 'month', 'weekday', 'hour', 'date_as_date', 'temp_sfc', 'friction_v', 'pbl', + 'obukhov_len', 'temp_top', 'wSpeed_top', 'Fb', 'S', ]).copy() + catalog_aux = self.comm.gather(catalog_aux, root=0) + + if self.comm.Get_rank() == 0: + catalog_aux = pd.concat(catalog_aux) + catalog_aux.to_csv(self.plume_rise_filename) + return True + def set_layer(self, catalog): spent_time = timeit.default_timer() diff --git a/hermesv3_bu/sectors/sector_manager.py b/hermesv3_bu/sectors/sector_manager.py index 575730a..5158230 100755 --- a/hermesv3_bu/sectors/sector_manager.py +++ b/hermesv3_bu/sectors/sector_manager.py @@ -153,6 +153,10 @@ class SectorManager(object): elif sector == 'point_sources' and comm_world.Get_rank() in sector_procs: from hermesv3_bu.sectors.point_source_sector import PointSourceSector + if arguments.plume_rise_output and arguments.plume_rise: + plume_rise_out_file = arguments.output_name.replace('.nc', '_ps_plumerise.csv') + else: + plume_rise_out_file = None self.sector = PointSourceSector( comm_world.Split(color, sector_procs.index(comm_world.Get_rank())), self.__logger, arguments.auxiliary_files_path, grid, clip, date_array, @@ -161,7 +165,7 @@ class SectorManager(object): arguments.point_source_weekly_profiles, arguments.point_source_hourly_profiles, arguments.speciation_map, arguments.point_source_speciation_profiles, arguments.point_source_snaps, arguments.point_source_measured_emissions, arguments.molecular_weights, - plume_rise=arguments.plume_rise, plume_rise_pahts={ + plume_rise=arguments.plume_rise, plume_rise_filename=plume_rise_out_file, plume_rise_pahts={ 'friction_velocity_dir': arguments.friction_velocity_dir, 'pblh_dir': arguments.pblh_dir, 'obukhov_length_dir': arguments.obukhov_length_dir, -- GitLab From 55659bedb6de3d03731899cfb9d465e44a7abb88 Mon Sep 17 00:00:00 2001 From: ctena Date: Thu, 17 Mar 2022 12:56:57 +0100 Subject: [PATCH 2/4] Re-formating PlumeRise vestribution output --- hermesv3_bu/sectors/point_source_sector.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/hermesv3_bu/sectors/point_source_sector.py b/hermesv3_bu/sectors/point_source_sector.py index 3089cfe..9a253af 100755 --- a/hermesv3_bu/sectors/point_source_sector.py +++ b/hermesv3_bu/sectors/point_source_sector.py @@ -766,14 +766,15 @@ class PointSourceSector(Sector): def write_plume_rise_output(self, catalog): catalog_aux = catalog.drop( - columns=['index', 'Height', 'Diameter', 'Speed', 'Temp', 'P_month', 'P_week', 'P_hour', 'P_spec', + columns=['index', 'step', 'Diameter', 'Speed', 'Temp', 'P_month', 'P_week', 'P_hour', 'P_spec', 'geometry', 'FID', 'month', 'weekday', 'hour', 'date_as_date', 'temp_sfc', 'friction_v', 'pbl', 'obukhov_len', 'temp_top', 'wSpeed_top', 'Fb', 'S', ]).copy() + catalog_aux.rename(columns={'Height': 'h_stck'}, inplace=True) catalog_aux = self.comm.gather(catalog_aux, root=0) if self.comm.Get_rank() == 0: catalog_aux = pd.concat(catalog_aux) - catalog_aux.to_csv(self.plume_rise_filename) + catalog_aux.to_csv(self.plume_rise_filename, index=False) return True def set_layer(self, catalog): -- GitLab From 341b6b285275a663c7f06a910f7002b2158cb68c Mon Sep 17 00:00:00 2001 From: ctena Date: Fri, 18 Mar 2022 16:01:40 +0100 Subject: [PATCH 3/4] Re-formating PlumeRise meteo paths with patterns --- conf/hermes.conf | 20 +++---- hermesv3_bu/config/config.py | 28 ++++----- hermesv3_bu/sectors/point_source_sector.py | 67 ++++++++++------------ hermesv3_bu/sectors/sector.py | 7 +++ hermesv3_bu/sectors/sector_manager.py | 20 +++---- hermesv3_bu/tools/checker.py | 1 + 6 files changed, 72 insertions(+), 71 deletions(-) diff --git a/conf/hermes.conf b/conf/hermes.conf index 37fd069..c40608e 100755 --- a/conf/hermes.conf +++ b/conf/hermes.conf @@ -163,16 +163,16 @@ temperature_daily_files_path = /ecmwf/era5/original_files/reorder/dai wind_speed_daily_files_path = /ecmwf/era5/original_files/reorder/daily_mean/sfcWind/ precipitation_files_path = /ecmwf/era5/original_files/reorder/1hourly/prlr/ -temperature_4d_dir = /esarchive/exp/monarch/a1wd/regional/hourly/t -temperature_sfc_dir = /esarchive/exp/monarch/a1wd/regional/hourly/t2 -u_wind_speed_4d_dir = /esarchive/exp/monarch/a1wd/regional/hourly/u -v_wind_speed_4d_dir = /esarchive/exp/monarch/a1wd/regional/hourly/v -u10_wind_speed_dir = /esarchive/exp/monarch/a1wd/regional/hourly/u10 -v10_wind_speed_dir = /esarchive/exp/monarch/a1wd/regional/hourly/v10 -friction_velocity_dir = /esarchive/exp/monarch/a1wd/regional/hourly/ustar -pblh_dir = /esarchive/exp/monarch/a1wd/regional/hourly/mixed_layer_height -obukhov_length_dir = /esarchive/exp/monarch/a1wd/regional/hourly/rmol -layer_thickness_dir = /esarchive/exp/monarch/a1wd/regional/hourly/layer_thickness +temperature_4d_path = /esarchive/exp/monarch/a1wd/regional/hourly/t/t_00.nc +temperature_sfc_path = /esarchive/exp/monarch/a1wd/regional/hourly/t2/t2_00.nc +u_wind_speed_4d_path = /esarchive/exp/monarch/a1wd/regional/hourly/u/u_00.nc +v_wind_speed_4d_path = /esarchive/exp/monarch/a1wd/regional/hourly/v/v_00.nc +u10_wind_speed_path = /esarchive/exp/monarch/a1wd/regional/hourly/u10/u10_00.nc +v10_wind_speed_path = /esarchive/exp/monarch/a1wd/regional/hourly/v10/v10_00.nc +friction_velocity_path = /esarchive/exp/monarch/a1wd/regional/hourly/ustar/ustar_00.nc +pblh_path = /esarchive/exp/monarch/a1wd/regional/hourly/mixed_layer_height/mixed_layer_height_00.nc +obukhov_length_path = /esarchive/exp/monarch/a1wd/regional/hourly/rmol/rmol_00.nc +layer_thickness_path = /esarchive/exp/monarch/a1wd/regional/hourly/layer_thickness/layer_thickness_00.nc [AVIATION SECTOR] diff --git a/hermesv3_bu/config/config.py b/hermesv3_bu/config/config.py index 9e6e07a..913b6ff 100755 --- a/hermesv3_bu/config/config.py +++ b/hermesv3_bu/config/config.py @@ -234,33 +234,33 @@ class Config(ArgParser): "to make a polygon or nothing to use the default clip: domain extension") # ===== METEO PATHS ===== - p.add_argument('--temperature_hourly_files_path', required=False, type=str, default='True', + p.add_argument('--temperature_hourly_files_path', required=False, type=str, default=None, help="Defines the path to the NetCDF files containing hourly mean 2m temperature data.") - p.add_argument('--temperature_daily_files_path', required=False, type=str, default='True', + p.add_argument('--temperature_daily_files_path', required=False, type=str, default=None, help="Defines the path to the NetCDF files containing daily mean 2m temperature data.") - p.add_argument('--wind_speed_daily_files_path', required=False, type=str, default='True', + p.add_argument('--wind_speed_daily_files_path', required=False, type=str, default=None, help="Defines the path to the NetCDF files containing daily mean 10m wind speed data.") - p.add_argument('--precipitation_files_path', required=False, type=str, default='True', + p.add_argument('--precipitation_files_path', required=False, type=str, default=None, help="Defines the path to the NetCDF files containing hourly mean precipitation data.") - p.add_argument('--temperature_4d_dir', required=False, type=str, default='True', + p.add_argument('--temperature_4d_path', required=False, type=str, default=None, help="Defines the path to the NetCDF files containing hourly mean 4D temperature data.") - p.add_argument('--temperature_sfc_dir', required=False, type=str, default='True', + p.add_argument('--temperature_sfc_path', required=False, type=str, default=None, help="Defines the path to the NetCDF files containing hourly mean surface temperature data.") - p.add_argument('--u_wind_speed_4d_dir', required=False, type=str, default='True', + p.add_argument('--u_wind_speed_4d_path', required=False, type=str, default=None, help="Defines the path to the NetCDF files containing hourly mean 4D U wind component data.") - p.add_argument('--v_wind_speed_4d_dir', required=False, type=str, default='True', + p.add_argument('--v_wind_speed_4d_path', required=False, type=str, default=None, help="Defines the path to the NetCDF files containing hourly mean 4D V wind component data.") - p.add_argument('--u10_wind_speed_dir', required=False, type=str, default='True', + p.add_argument('--u10_wind_speed_path', required=False, type=str, default=None, help="Defines the path to the NetCDF files containing daily mean 10m U wind component data.") - p.add_argument('--v10_wind_speed_dir', required=False, type=str, default='True', + p.add_argument('--v10_wind_speed_path', required=False, type=str, default=None, help="Defines the path to the NetCDF files containing daily mean 10m V wind component data.") - p.add_argument('--friction_velocity_dir', required=False, type=str, default='True', + p.add_argument('--friction_velocity_path', required=False, type=str, default=None, help="Defines the path to the NetCDF files containing hourly mean 4D friction velocity data.") - p.add_argument('--pblh_dir', required=False, type=str, default='True', + p.add_argument('--pblh_path', required=False, type=str, default=None, help="Defines the path to the NetCDF files containing hourly mean PBL height data.") - p.add_argument('--obukhov_length_dir', required=False, type=str, default='True', + p.add_argument('--obukhov_length_path', required=False, type=str, default=None, help="Defines the path to the NetCDF files containing hourly mean Obukhov length data.") - p.add_argument('--layer_thickness_dir', required=False, type=str, default='True', + p.add_argument('--layer_thickness_path', required=False, type=str, default=None, help="Defines the path to the NetCDF files containing hourly mean 4D layer thickness data.") # ***** AVIATION SECTOR ***** diff --git a/hermesv3_bu/sectors/point_source_sector.py b/hermesv3_bu/sectors/point_source_sector.py index 6a404ba..0bbae01 100755 --- a/hermesv3_bu/sectors/point_source_sector.py +++ b/hermesv3_bu/sectors/point_source_sector.py @@ -468,10 +468,9 @@ class PointSourceSector(Sector): return dataframe[['X', 'Y']] def get_plumerise_meteo(self, catalog): - def get_sfc_value(dataframe, dir_path, var_name): + def get_sfc_value(dataframe, filepath, var_name): from netCDF4 import Dataset, num2date - nc_path = os.path.join(dir_path, - '{0}_{1}.nc'.format(var_name, dataframe.name.replace(hour=0).strftime("%Y%m%d%H"))) + nc_path = self.parse_netcdf_path(filepath, date_aux=dataframe.name) check_files(nc_path) netcdf = Dataset(nc_path, mode='r') # time_index @@ -492,10 +491,9 @@ class PointSourceSector(Sector): return dataframe[[var_name]] - def get_layers(dataframe, dir_path, var_name): + def get_layers(dataframe, filepath, var_name): from netCDF4 import Dataset, num2date - nc_path = os.path.join(dir_path, - '{0}_{1}.nc'.format(var_name, dataframe.name.replace(hour=0).strftime("%Y%m%d%H"))) + nc_path = self.parse_netcdf_path(filepath, date_aux=dataframe.name) check_files(nc_path) netcdf = Dataset(nc_path, mode='r') # time_index @@ -528,12 +526,11 @@ class PointSourceSector(Sector): return dataframe[['layers']] - def get_temp_top(dataframe, dir_path, var_name): + def get_temp_top(dataframe, filepath, var_name): from netCDF4 import Dataset, num2date from scipy.interpolate import interp1d as interpolate - nc_path = os.path.join(dir_path, - '{0}_{1}.nc'.format(var_name, dataframe.name.replace(hour=0).strftime("%Y%m%d%H"))) + nc_path = self.parse_netcdf_path(filepath, date_aux=dataframe.name) check_files(nc_path) netcdf = Dataset(nc_path, mode='r') # time_index @@ -572,8 +569,7 @@ class PointSourceSector(Sector): def get_wind_speed_10m(dataframe, u_dir_path, v_dir_path, u_var_name, v_var_name): from netCDF4 import Dataset, num2date # === u10 === - u10_nc_path = os.path.join( - u_dir_path, '{0}_{1}.nc'.format(u_var_name, dataframe.name.replace(hour=0).strftime("%Y%m%d%H"))) + u10_nc_path = self.parse_netcdf_path(u_dir_path, date_aux=dataframe.name) check_files(u10_nc_path) u10_netcdf = Dataset(u10_nc_path, mode='r') # time_index @@ -593,8 +589,7 @@ class PointSourceSector(Sector): dataframe['u10'] = var[dataframe['Y'], dataframe['X']] # === v10 === - v10_nc_path = os.path.join( - v_dir_path, '{0}_{1}.nc'.format(v_var_name, dataframe.name.replace(hour=0).strftime("%Y%m%d%H"))) + v10_nc_path = self.parse_netcdf_path(v_dir_path, date_aux=dataframe.name) check_files(v10_nc_path) v10_netcdf = Dataset(v10_nc_path, mode='r') @@ -614,8 +609,7 @@ class PointSourceSector(Sector): from netCDF4 import Dataset, num2date from scipy.interpolate import interp1d as interpolate # === u10 === - u10_nc_path = os.path.join( - u_dir_path, '{0}_{1}.nc'.format(u_var_name, dataframe.name.replace(hour=0).strftime("%Y%m%d%H"))) + u10_nc_path = self.parse_netcdf_path(u_dir_path, date_aux=dataframe.name) check_files(u10_nc_path) u10_netcdf = Dataset(u10_nc_path, mode='r') # time_index @@ -639,8 +633,7 @@ class PointSourceSector(Sector): dataframe['u_{0}'.format(i)] = t_lay # === v10 === - v10_nc_path = os.path.join( - v_dir_path, '{0}_{1}.nc'.format(v_var_name, dataframe.name.replace(hour=0).strftime("%Y%m%d%H"))) + v10_nc_path = self.parse_netcdf_path(v_dir_path, date_aux=dataframe.name) check_files(v10_nc_path) v10_netcdf = Dataset(v10_nc_path, mode='r') @@ -669,49 +662,49 @@ class PointSourceSector(Sector): # TODO Use IoNetCDF spent_time = timeit.default_timer() - - meteo_xy = self.get_meteo_xy(catalog.groupby('Code').first(), os.path.join( - self.plume_rise_pahts['temperature_sfc_dir'], - 't2_{0}.nc'.format(self.date_array[0].replace(hour=0).strftime("%Y%m%d%H")))) + print(self.parse_netcdf_path(self.plume_rise_pahts['temperature_sfc_path'], date_aux=self.date_array[0])) + meteo_xy = self.get_meteo_xy(catalog.groupby('Code').first(), + self.parse_netcdf_path(self.plume_rise_pahts['temperature_sfc_path'], + date_aux=self.date_array[0])) catalog = catalog.merge(meteo_xy, left_index=True, right_index=True) # ===== 3D Meteo variables ===== # Adding stc_temp - self.logger.write_log('\t\tGetting temperature from {0}'.format(self.plume_rise_pahts['temperature_sfc_dir']), + self.logger.write_log('\t\tGetting temperature from {0}'.format(self.plume_rise_pahts['temperature_sfc_path']), message_level=3) catalog['temp_sfc'] = catalog.groupby('date_utc')['X', 'Y'].apply( - lambda x: get_sfc_value(x, self.plume_rise_pahts['temperature_sfc_dir'], 't2')) + lambda x: get_sfc_value(x, self.plume_rise_pahts['temperature_sfc_path'], 't2')) self.logger.write_log('\t\tGetting friction velocity from {0}'.format( - self.plume_rise_pahts['friction_velocity_dir']), message_level=3) + self.plume_rise_pahts['friction_velocity_path']), message_level=3) catalog['friction_v'] = catalog.groupby('date_utc')['X', 'Y'].apply( - lambda x: get_sfc_value(x, self.plume_rise_pahts['friction_velocity_dir'], 'ustar')) + lambda x: get_sfc_value(x, self.plume_rise_pahts['friction_velocity_path'], 'ustar')) self.logger.write_log('\t\tGetting PBL height from {0}'.format( - self.plume_rise_pahts['pblh_dir']), message_level=3) + self.plume_rise_pahts['pblh_path']), message_level=3) catalog['pbl'] = catalog.groupby('date_utc')['X', 'Y'].apply( - lambda x: get_sfc_value(x, self.plume_rise_pahts['pblh_dir'], 'mixed_layer_height')) + lambda x: get_sfc_value(x, self.plume_rise_pahts['pblh_path'], 'mixed_layer_height')) self.logger.write_log('\t\tGetting obukhov length from {0}'.format( - self.plume_rise_pahts['obukhov_length_dir']), message_level=3) + self.plume_rise_pahts['obukhov_length_path']), message_level=3) catalog['obukhov_len'] = catalog.groupby('date_utc')['X', 'Y'].apply( - lambda x: get_sfc_value(x, self.plume_rise_pahts['obukhov_length_dir'], 'rmol')) + lambda x: get_sfc_value(x, self.plume_rise_pahts['obukhov_length_path'], 'rmol')) catalog['obukhov_len'] = 1. / catalog['obukhov_len'] self.logger.write_log('\t\tGetting layer thickness from {0}'.format( - self.plume_rise_pahts['layer_thickness_dir']), message_level=3) + self.plume_rise_pahts['layer_thickness_path']), message_level=3) catalog['layers'] = catalog.groupby('date_utc')['X', 'Y'].apply( - lambda x: get_layers(x, self.plume_rise_pahts['layer_thickness_dir'], 'layer_thickness')) + lambda x: get_layers(x, self.plume_rise_pahts['layer_thickness_path'], 'layer_thickness')) self.logger.write_log('\t\tGetting temperatue at the top from {0}'.format( - self.plume_rise_pahts['temperature_4d_dir']), message_level=3) + self.plume_rise_pahts['temperature_4d_path']), message_level=3) catalog['temp_top'] = catalog.groupby('date_utc')['X', 'Y', 'Height', 'layers', 'temp_sfc'].apply( - lambda x: get_temp_top(x, self.plume_rise_pahts['temperature_4d_dir'], 't')) + lambda x: get_temp_top(x, self.plume_rise_pahts['temperature_4d_path'], 't')) self.logger.write_log('\t\tGetting wind speed at 10 m', message_level=3) catalog['wSpeed_10'] = catalog.groupby('date_utc')['X', 'Y'].apply( - lambda x: get_wind_speed_10m(x, self.plume_rise_pahts['u10_wind_speed_dir'], - self.plume_rise_pahts['v10_wind_speed_dir'], 'u10', 'v10')) + lambda x: get_wind_speed_10m(x, self.plume_rise_pahts['u10_wind_speed_path'], + self.plume_rise_pahts['v10_wind_speed_path'], 'u10', 'v10')) self.logger.write_log('\t\tGetting wind speed at the top', message_level=3) catalog['wSpeed_top'] = catalog.groupby('date_utc')['X', 'Y', 'Height', 'layers', 'wSpeed_10'].apply( - lambda x: get_wind_speed_top(x, self.plume_rise_pahts['u_wind_speed_4d_dir'], - self.plume_rise_pahts['v_wind_speed_4d_dir'], 'u', 'v')) + lambda x: get_wind_speed_top(x, self.plume_rise_pahts['u_wind_speed_4d_path'], + self.plume_rise_pahts['v_wind_speed_4d_path'], 'u', 'v')) catalog.drop(columns=['wSpeed_10', 'layers', 'X', 'Y'], inplace=True) self.logger.write_time_log('PointSourceSector', 'get_plumerise_meteo', timeit.default_timer() - spent_time) return catalog diff --git a/hermesv3_bu/sectors/sector.py b/hermesv3_bu/sectors/sector.py index bac090a..73ec946 100755 --- a/hermesv3_bu/sectors/sector.py +++ b/hermesv3_bu/sectors/sector.py @@ -638,3 +638,10 @@ class Sector(object): self.comm.Barrier() return True + + @staticmethod + def parse_netcdf_path(path, date_aux=None): + if date_aux is not None: + path = path.replace('', date_aux.strftime("%Y%m%d%H")) + path = path.replace('', date_aux.strftime("%Y%m%d")) + return path diff --git a/hermesv3_bu/sectors/sector_manager.py b/hermesv3_bu/sectors/sector_manager.py index 575730a..c7a53a6 100755 --- a/hermesv3_bu/sectors/sector_manager.py +++ b/hermesv3_bu/sectors/sector_manager.py @@ -162,16 +162,16 @@ class SectorManager(object): arguments.speciation_map, arguments.point_source_speciation_profiles, arguments.point_source_snaps, arguments.point_source_measured_emissions, arguments.molecular_weights, plume_rise=arguments.plume_rise, plume_rise_pahts={ - 'friction_velocity_dir': arguments.friction_velocity_dir, - 'pblh_dir': arguments.pblh_dir, - 'obukhov_length_dir': arguments.obukhov_length_dir, - 'layer_thickness_dir': arguments.layer_thickness_dir, - 'temperature_sfc_dir': arguments.temperature_sfc_dir, - 'temperature_4d_dir': arguments.temperature_4d_dir, - 'u10_wind_speed_dir': arguments.u10_wind_speed_dir, - 'v10_wind_speed_dir': arguments.v10_wind_speed_dir, - 'u_wind_speed_4d_dir': arguments.u_wind_speed_4d_dir, - 'v_wind_speed_4d_dir': arguments.v_wind_speed_4d_dir}) + 'friction_velocity_path': arguments.friction_velocity_path, + 'pblh_path': arguments.pblh_path, + 'obukhov_length_path': arguments.obukhov_length_path, + 'layer_thickness_path': arguments.layer_thickness_path, + 'temperature_sfc_path': arguments.temperature_sfc_path, + 'temperature_4d_path': arguments.temperature_4d_path, + 'u10_wind_speed_path': arguments.u10_wind_speed_path, + 'v10_wind_speed_path': arguments.v10_wind_speed_path, + 'u_wind_speed_4d_path': arguments.u_wind_speed_4d_path, + 'v_wind_speed_4d_path': arguments.v_wind_speed_4d_path}) elif sector == 'traffic' and comm_world.Get_rank() in sector_procs: from hermesv3_bu.sectors.traffic_sector import TrafficSector self.sector = TrafficSector( diff --git a/hermesv3_bu/tools/checker.py b/hermesv3_bu/tools/checker.py index 69d124a..f7dd6b8 100644 --- a/hermesv3_bu/tools/checker.py +++ b/hermesv3_bu/tools/checker.py @@ -24,6 +24,7 @@ def check_files(file_path_list, warning=False): warn(error_message.replace('ERROR', 'WARNING')) return False else: + # raise FileNotFoundError(error_message) error_exit(error_message) return True -- GitLab From 5b834021b356ba9d2597c6219a80d2c260de8e17 Mon Sep 17 00:00:00 2001 From: ctena Date: Tue, 29 Mar 2022 12:37:12 +0200 Subject: [PATCH 4/4] PlumeRise with global meteo: Ready to test --- hermesv3_bu/config/config.py | 2 +- hermesv3_bu/sectors/point_source_sector.py | 99 ++++++++++++++++------ hermesv3_bu/sectors/sector_manager.py | 2 +- 3 files changed, 76 insertions(+), 27 deletions(-) diff --git a/hermesv3_bu/config/config.py b/hermesv3_bu/config/config.py index fc11549..6ea62b7 100755 --- a/hermesv3_bu/config/config.py +++ b/hermesv3_bu/config/config.py @@ -789,7 +789,7 @@ class Config(ArgParser): if arguments.plume_rise: arguments.plume_rise_output = self._parse_bool(arguments.plume_rise_output) else: - arguments.plume_rise_output = False + arguments.plume_rise_output = None # Point Source lists arguments.point_source_pollutants = self._parse_list(arguments.point_source_pollutants) diff --git a/hermesv3_bu/sectors/point_source_sector.py b/hermesv3_bu/sectors/point_source_sector.py index edee48a..bef285e 100755 --- a/hermesv3_bu/sectors/point_source_sector.py +++ b/hermesv3_bu/sectors/point_source_sector.py @@ -52,7 +52,7 @@ class PointSourceSector(Sector): 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_filename=False, plume_rise_pahts=None): + molecular_weights_path, plume_rise=False, plume_rise_filename=None, plume_rise_pahts=None): spent_time = timeit.default_timer() logger.write_log('===== POINT SOURCES SECTOR =====') check_files( @@ -70,6 +70,7 @@ class PointSourceSector(Sector): self.catalog_measured = self.read_catalog_for_measured_emissions(catalog_path, sector_list) self.measured_path = measured_emission_path self.plume_rise_pahts = plume_rise_pahts + self.meteo_xy = self.get_meteo_xy() self.logger.write_time_log('PointSourceSector', '__init__', timeit.default_timer() - spent_time) @@ -423,7 +424,29 @@ class PointSourceSector(Sector): timeit.default_timer() - spent_time) return catalog - def get_meteo_xy(self, dataframe, netcdf_path): + def get_meteo_xy(self): + if self.plume_rise: + meteo_xy_aux_path = os.path.join(self.auxiliary_dir, 'point_sources', 'netcdf_locations.csv') + if not os.path.exists(meteo_xy_aux_path): + meteo_xy = self.calculate_meteo_xy( + self.catalog.groupby('Code').first(), + self.parse_netcdf_path(self.plume_rise_pahts['temperature_sfc_path'], date_aux=self.date_array[0])) + + meteo_xy_aux = self.comm.gather(meteo_xy, root=0) + if self.comm.Get_rank() == 0: + if not os.path.exists(os.path.dirname(meteo_xy_aux_path)): + os.makedirs(os.path.dirname(meteo_xy_aux_path)) + meteo_xy_aux = pd.concat(meteo_xy_aux) + meteo_xy_aux.to_csv(meteo_xy_aux_path) + else: + meteo_xy = pd.read_csv(meteo_xy_aux_path, index_col='Code') + else: + meteo_xy = None + print(meteo_xy) + sys.stdout.flush() + return meteo_xy + + def calculate_meteo_xy(self, dataframe, netcdf_path): spent_time = timeit.default_timer() def nearest(row, geom_union, df1, df2, geom1_col='geometry', geom2_col='geometry', src_column=None): @@ -443,25 +466,48 @@ class PointSourceSector(Sector): import pandas as pd import geopandas as gpd check_files(netcdf_path) - nc = Dataset(netcdf_path, mode='r') - try: - lats = nc.variables['lat'][:] - lons = nc.variables['lon'][:] - except KeyError as e: - error_exit("{0} variable not found in {1} file.".format(str(e), netcdf_path)) - x = np.array([np.arange(lats.shape[1])] * lats.shape[0]) - y = np.array([np.arange(lats.shape[0]).T] * lats.shape[1]).T - - nc_dataframe = pd.DataFrame.from_dict({'X': x.flatten(), 'Y': y.flatten()}) - nc_dataframe = gpd.GeoDataFrame(nc_dataframe, - geometry=[Point(xy) for xy in list(zip(lons.flatten(), lats.flatten()))], - crs={'init': 'epsg:4326'}) - nc_dataframe['index'] = nc_dataframe.index - - union = nc_dataframe.unary_union + if self.comm.Get_rank() == 0: + self.logger.write_log("\t\t\tGetting nearest neighbours", message_level=0) + nc = Dataset(netcdf_path, mode='r') + self.logger.write_log("\t\t\t\tReading coordinates", message_level=3) + try: + lats = nc.variables['lat'][:] + lons = nc.variables['lon'][:] + except KeyError as e: + error_exit("{0} variable not found in {1} file.".format(str(e), netcdf_path)) + + if len(lats.shape) > 1: # Rotated + x = np.array([np.arange(lats.shape[1])] * lats.shape[0]) + y = np.array([np.arange(lats.shape[0]).T] * lats.shape[1]).T + else: # Global + x = np.array([np.arange(lons.shape[0])] * lats.shape[0]) + y = np.array([np.arange(lats.shape[0])] * lons.shape[0]).T + + lats = np.array([lats] * lons.shape[0]).T + lons = np.array([lons] * lats.shape[0]) + # print('X =====', x.shape, x) + # print('Y =====', y.shape, y) + # + # print('Lats =====', lats.shape, lats) + # print('Lons =====', lons.shape, lons) + self.logger.write_log("\t\t\t\tTo GeoDataframe", message_level=3) + nc_dataframe = pd.DataFrame.from_dict({'X': x.flatten(), 'Y': y.flatten()}) + nc_dataframe = gpd.GeoDataFrame(nc_dataframe, + geometry=[Point(xy) for xy in list(zip(lons.flatten(), lats.flatten()))], + crs={'init': 'epsg:4326'}) + nc_dataframe['index'] = nc_dataframe.index + self.logger.write_log("\t\t\t\tUnion", message_level=3) + union = nc_dataframe.unary_union + else: + nc_dataframe = None + union = None + self.logger.write_log("\t\t\t\tBroadcasting", message_level=3) + nc_dataframe = self.comm.bcast(nc_dataframe, root=0) + union = self.comm.bcast(union, root=0) + self.logger.write_log("\t\t\t\tStarting nearest neighbour", message_level=3) dataframe['meteo_index'] = dataframe.apply( nearest, geom_union=union, df1=dataframe, df2=nc_dataframe, src_column='index', axis=1) - + self.logger.write_log("\t\t\t\tAssigning X & Y positions", message_level=3) dataframe['X'] = nc_dataframe.loc[dataframe['meteo_index'], 'X'].values dataframe['Y'] = nc_dataframe.loc[dataframe['meteo_index'], 'Y'].values @@ -469,6 +515,8 @@ class PointSourceSector(Sector): return dataframe[['X', 'Y']] def get_plumerise_meteo(self, catalog): + self.logger.write_log("\t\tPlume Raise selected", message_level=0) + def get_sfc_value(dataframe, filepath, var_name): from netCDF4 import Dataset, num2date nc_path = self.parse_netcdf_path(filepath, date_aux=dataframe.name) @@ -663,12 +711,13 @@ class PointSourceSector(Sector): # TODO Use IoNetCDF spent_time = timeit.default_timer() - print(self.parse_netcdf_path(self.plume_rise_pahts['temperature_sfc_path'], date_aux=self.date_array[0])) - meteo_xy = self.get_meteo_xy(catalog.groupby('Code').first(), - self.parse_netcdf_path(self.plume_rise_pahts['temperature_sfc_path'], - date_aux=self.date_array[0])) - catalog = catalog.merge(meteo_xy, left_index=True, right_index=True) + # meteo_xy = self.calculate_meteo_xy(catalog.groupby('Code').first(), + # self.parse_netcdf_path(self.plume_rise_pahts['temperature_sfc_path'], + # date_aux=self.date_array[0])) + # + # catalog = catalog.merge(meteo_xy, left_index=True, right_index=True) + catalog = catalog.merge(self.meteo_xy, left_index=True, right_index=True, how='left') # ===== 3D Meteo variables ===== # Adding stc_temp @@ -759,7 +808,7 @@ class PointSourceSector(Sector): def write_plume_rise_output(self, catalog): catalog_aux = catalog.drop( - columns=['index', 'step', 'Diameter', 'Speed', 'Temp', 'P_month', 'P_week', 'P_hour', 'P_spec', + columns=['index', 'tstep', 'Diameter', 'Speed', 'Temp', 'P_month', 'P_week', 'P_hour', 'P_spec', 'geometry', 'FID', 'month', 'weekday', 'hour', 'date_as_date', 'temp_sfc', 'friction_v', 'pbl', 'obukhov_len', 'temp_top', 'wSpeed_top', 'Fb', 'S', ]).copy() catalog_aux.rename(columns={'Height': 'h_stck'}, inplace=True) diff --git a/hermesv3_bu/sectors/sector_manager.py b/hermesv3_bu/sectors/sector_manager.py index 6ce18c8..3f35e30 100755 --- a/hermesv3_bu/sectors/sector_manager.py +++ b/hermesv3_bu/sectors/sector_manager.py @@ -165,7 +165,7 @@ class SectorManager(object): arguments.point_source_weekly_profiles, arguments.point_source_hourly_profiles, arguments.speciation_map, arguments.point_source_speciation_profiles, arguments.point_source_snaps, arguments.point_source_measured_emissions, arguments.molecular_weights, - plume_rise=arguments.plume_rise, plume_rise_pahts={ + plume_rise=arguments.plume_rise, plume_rise_filename=plume_rise_out_file, plume_rise_pahts={ 'friction_velocity_path': arguments.friction_velocity_path, 'pblh_path': arguments.pblh_path, 'obukhov_length_path': arguments.obukhov_length_path, -- GitLab