diff --git a/CHANGELOG b/CHANGELOG index 9a76d15812d2cd918101180c2eb5f194f640d161..27c4272d8f98b3966b3234bc609beea895a80d5a 100755 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,3 +1,13 @@ +1.0.8 + 2023/06/02 + + - Shipping port emissions hardcoded to 2nd layer (layer = 1) + +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 37fd0691eb2301f81b19bb86554c602ee3dfce1e..8d2ac1969d10fa9b2c788d4becabaed80ebe203f 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] @@ -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 382021f30ff6ce35872aa0de4ebe2a43f50c21d7..e13bd590c858eab2ba3506af867008063b736c2e 100755 --- a/hermesv3_bu/__init__.py +++ b/hermesv3_bu/__init__.py @@ -1 +1 @@ -__version__ = "1.0.6" +__version__ = "1.0.8" diff --git a/hermesv3_bu/config/config.py b/hermesv3_bu/config/config.py index 9e6e07a08831c81fdc8f80d91316a912a4ab8350..6ea62b7e13386d33e947f54c335cc69c695801a9 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 ***** @@ -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 = 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 6a404ba6284b9efb09602a068f7f495dab8b5afc..bef285e9329438c3ea0c40933f056b9006c4cb5a 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=None, plume_rise_pahts=None): spent_time = timeit.default_timer() logger.write_log('===== POINT SOURCES SECTOR =====') check_files( @@ -64,11 +64,13 @@ 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) 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) @@ -422,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): @@ -442,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 @@ -468,10 +515,11 @@ class PointSourceSector(Sector): return dataframe[['X', 'Y']] def get_plumerise_meteo(self, catalog): - def get_sfc_value(dataframe, dir_path, var_name): + 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 = 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 +540,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 +575,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 +618,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 +638,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 +658,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 +682,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') @@ -670,48 +712,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")))) - - 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 - 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 @@ -755,12 +798,27 @@ 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', '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) + 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, index=False) + return True + def set_layer(self, catalog): spent_time = timeit.default_timer() diff --git a/hermesv3_bu/sectors/sector.py b/hermesv3_bu/sectors/sector.py index bac090adb382ad76e84f6b7f94ebf6eb29375cfb..73ec9461e906d91efb1cb33925f25b85de3b0180 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 575730a6a150eb3c9f30855a7ad815235f3d8d6d..3f35e30835c87ee94a1f94f686140493c8a0e5fb 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,17 +165,17 @@ 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={ - '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}) + 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, + '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/sectors/shipping_port_sector.py b/hermesv3_bu/sectors/shipping_port_sector.py index 88759a2bc2ae2cde8aa395f63acb6ab6d5a0c464..09d8a1266cd8126b27af517e3080119d6acb974e 100755 --- a/hermesv3_bu/sectors/shipping_port_sector.py +++ b/hermesv3_bu/sectors/shipping_port_sector.py @@ -594,7 +594,7 @@ class ShippingPortSector(Sector): dataframe.rename(columns={'idx2': 'FID'}, inplace=True) dataframe.drop(columns=['src_inter_fraction', 'idx1', 'geometry'], inplace=True) - dataframe['layer'] = 0 + dataframe['layer'] = 1 dataframe = dataframe.loc[:, ~dataframe.columns.duplicated()] dataframe = dataframe.groupby(['FID', 'layer', 'tstep']).sum() self.logger.write_time_log('ShippingPortSector', 'to_grid_geometry', timeit.default_timer() - spent_time) diff --git a/hermesv3_bu/tools/checker.py b/hermesv3_bu/tools/checker.py index 69d124a7f5af62044977cb7cc8393bb1e9b654e1..f7dd6b811844260cd33177fd6cb467d6891fb9da 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