From a03812c00e0a1644d794f44b2be5c68ebdf92533 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Wed, 10 Jun 2020 13:20:25 +0200 Subject: [PATCH 1/4] Corrected bug on Clip shapefile not found error message --- hermesv3_bu/clipping/shapefile_clip.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hermesv3_bu/clipping/shapefile_clip.py b/hermesv3_bu/clipping/shapefile_clip.py index 7d4eb46..4a87412 100755 --- a/hermesv3_bu/clipping/shapefile_clip.py +++ b/hermesv3_bu/clipping/shapefile_clip.py @@ -51,7 +51,7 @@ class ShapefileClip(Clip): clip = gpd.GeoDataFrame(geometry=[geom], crs=clip.crs) clip.to_file(self.shapefile_path) else: - error_exit(" Clip shapefile {0} not found.") + error_exit(" Clip shapefile {0} not found.".format(clip_path)) else: clip = gpd.read_file(self.shapefile_path) self.logger.write_log("\tClip created at '{0}'".format(self.shapefile_path), 3) -- GitLab From 2df3b4f3e0c67c42a6a76f6a8089da0ba8f1b14d Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Thu, 11 Jun 2020 11:21:32 +0200 Subject: [PATCH 2/4] Climatology to test --- conf/hermes.conf | 53 +++++++------- hermesv3_bu/hermes.py | 4 +- hermesv3_bu/io_server/io_netcdf.py | 86 ++++++++++++++++++++--- hermesv3_bu/sectors/residential_sector.py | 4 +- hermesv3_bu/sectors/traffic_sector.py | 29 +++++--- 5 files changed, 128 insertions(+), 48 deletions(-) diff --git a/conf/hermes.conf b/conf/hermes.conf index 70e57f6..4be432e 100755 --- a/conf/hermes.conf +++ b/conf/hermes.conf @@ -2,14 +2,14 @@ log_level = 3 input_dir = /home/Earth/ctena/Models/hermesv3_bu_data data_path = /esarchive/recon -output_dir = /scratch/Earth/ctena/HERMESv3_BU_OUT/ -output_name = HERMESv3__p4_w2.nc +output_dir = /scratch/Earth/ctena/HERMESv3_BU_OUT/climatology +output_name = HERMESv3__climatology_crop_fertilizers.nc emission_summary = 0 -start_date = 2016/11/29 00:00:00 +start_date = 2020/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/ctena/HERMESv3_BU_aux/__test +auxiliary_files_path = /scratch/Earth/ctena/HERMESv3_BU_aux/_ first_time = 0 erase_auxiliary_files = 0 @@ -17,7 +17,7 @@ erase_auxiliary_files = 0 [DOMAIN] # domain_type=[lcc, rotated, mercator, regular, rotated_nested] -domain_type = rotated_nested +domain_type = lcc # output_type=[MONARCH, CMAQ, WRF_CHEM, DEFAULT] output_model = DEFAULT compression_level = 0 @@ -67,22 +67,20 @@ vertical_description = /profiles/vertical/MONARCH_Global_48layers_ver #y_0 = 43862.90625 # CATALUNYA test - nx = 28 - ny = 30 - #nx = 4 - #ny = 4 - inc_x = 10000 - inc_y = 10000 - x_0 = 253151.59375 - y_0 = 43862.90625 + # nx = 28 + # ny = 30 + # inc_x = 10000 + # inc_y = 10000 + # x_0 = 253151.59375 + # y_0 = 43862.90625 # IP - #nx = 397 - #ny = 397 - #inc_x = 4000 - #inc_y = 4000 - #x_0 = -807847.688 - #y_0 = -797137.125 + nx = 397 + ny = 397 + inc_x = 4000 + inc_y = 4000 + x_0 = -807847.688 + y_0 = -797137.125 # EUROPA #nx = 478 @@ -98,6 +96,7 @@ vertical_description = /profiles/vertical/MONARCH_Global_48layers_ver #ny = 158 #inc_x = 1000 #inc_y = 1000 + #x_0 = -142848.422 #y_0 = -20137.891 @@ -121,7 +120,7 @@ vertical_description = /profiles/vertical/MONARCH_Global_48layers_ver [CLIPPING] -# clipping = /Shapefiles/barcelona/barcelona_munic.shp +clipping = /shapefiles/barcelona/barcelona_munic.shp # clipping = 2.2 41.41, 2.15 41.41, 2.2 41.42, 2.15 41.42 # clipping = 2.2 41.41, 2.19 41.41 @@ -129,19 +128,19 @@ vertical_description = /profiles/vertical/MONARCH_Global_48layers_ver # SECTORS #################################################################### [SECTOR MANAGEMENT] -writing_processors = 2 +writing_processors = 1 # # aviation_processors = 1 # shipping_port_processors = 1 -# livestock_processors = 12 +# livestock_processors = 4 # crop_operations_processors = 1 # crop_fertilizers_processors = 4 # agricultural_machinery_processors = 1 # residential_processors = 4 -# recreational_boats_processors = 4 -# point_sources_processors = 16 -# traffic_processors = 256 -traffic_area_processors = 4 +# recreational_boats_processors = 1 +# point_sources_processors = 1 +# traffic_processors = 1 +# traffic_area_processors = 1 [SHAPEFILES] nuts3_shapefile = /shapefiles/nuts3/nuts3.shp @@ -293,7 +292,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 = True +plume_rise = False # 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/hermes.py b/hermesv3_bu/hermes.py index d44fba5..fb7ab65 100755 --- a/hermesv3_bu/hermes.py +++ b/hermesv3_bu/hermes.py @@ -70,7 +70,9 @@ class HermesBu(object): self.writer.write(emis) self.comm.Barrier() - + if self.comm.Get_rank() == 0: + print('***** HERMESv3_BU simulation finished successfully *****') + sys.stdout.flush() self.logger.write_log('***** HERMESv3_BU simulation finished successfully *****') self.logger.write_time_log('HERMES', 'TOTAL', timeit.default_timer() - self.initial_time) self.logger.finish_logs() diff --git a/hermesv3_bu/io_server/io_netcdf.py b/hermesv3_bu/io_server/io_netcdf.py index caa9558..c2bd7e5 100755 --- a/hermesv3_bu/io_server/io_netcdf.py +++ b/hermesv3_bu/io_server/io_netcdf.py @@ -3,7 +3,7 @@ import sys import os from mpi4py import MPI -from datetime import timedelta +from warnings import warn import numpy as np import geopandas as gpd from netCDF4 import Dataset @@ -43,7 +43,27 @@ class IoNetcdf(IoServer): :return: GeoDataframe with the data in the desired points. :rtype: geopandas.GeoDataframe """ - check_files(netcdf_path) + # check_files(netcdf_path) + climatology = False + if not os.path.exists(netcdf_path): + if date_type == 'daily': + path_climatology = netcdf_path[:-10] + '_climatology{0}.nc'.format(str(date.month).zfill(2)) + + elif date_type == 'yearly': + path_climatology = netcdf_path[:-8] + '_climatology.nc' + if not os.path.exists(path_climatology): + error_message = "*ERROR* (Rank {0}) File/s not found:".format(MPI.COMM_WORLD.Get_rank()) + error_message += "\n\t{0}\n\tneither {1}".format(netcdf_path, path_climatology) + error_exit(error_message) + else: + error_message = "*WARNING* (Rank {0}) File/s not found:".format(MPI.COMM_WORLD.Get_rank()) + error_message += "\n\t{0}\n\t using {1}".format(netcdf_path, path_climatology) + print(error_message) + sys.stdout.flush() + warn(error_message) + + climatology = True + netcdf_path = path_climatology nc = Dataset(netcdf_path, mode='r') try: lat_o = nc.variables['latitude'][:] @@ -62,7 +82,11 @@ class IoNetcdf(IoServer): error_exit("{0} variable not found in {1} file.".format(str(e), netcdf_path)) # From time array to list of dates. time_array = num2date(time[:], time.units, CALENDAR_STANDARD) - time_array = np.array([aux.date() for aux in time_array]) + + if climatology: + time_array = np.array([aux.date().replace(year=date.year) for aux in time_array]) + else: + time_array = np.array([aux.date() for aux in time_array]) i_time = np.where(time_array == date)[0][0] elif date_type == 'yearly': i_time = 0 @@ -114,30 +138,57 @@ class IoNetcdf(IoServer): :param lat_max: Maximum latitude of the centroid of the road links. :type lat_max: float - :return: Temperature, centroid of the cell and cell identificator (REC). + :return: Temperature, centroid of the cell and cell id (REC). Each time step is each column with the name t_. :rtype: GeoDataFrame """ + climatology = False path = os.path.join(netcdf_dir, '{0}_{1}{2}.nc'.format(var_name, date_array[0].year, str(date_array[0].month).zfill(2))) # self.logger.write_log('Getting temperature from {0}'.format(path), message_level=2) - check_files(path) + if not os.path.exists(path): + path_climatology = os.path.join(netcdf_dir, '{0}_{1}{2}.nc'.format( + var_name, 'climatology', str(date_array[0].month).zfill(2))) + if not os.path.exists(path_climatology): + error_message = "*ERROR* (Rank {0}) File/s not found:".format(MPI.COMM_WORLD.Get_rank()) + error_message += "\n\t{0}\n\tneither {1}".format(path, path_climatology) + if var_name == 'prlr': + raise FileNotFoundError(path) + else: + error_exit(error_message) + else: + error_message = "*WARNING* (Rank {0}) File/s not found:".format(MPI.COMM_WORLD.Get_rank()) + error_message += "\n\t{0}\n\t using {1}".format(path, path_climatology) + print(error_message) + sys.stdout.flush() + warn(error_message) + + climatology = True + path = path_climatology + nc = Dataset(path, mode='r') try: lat_o = nc.variables['latitude'][:] lon_o = nc.variables['longitude'][:] n_lat = len(lat_o) - time = nc.variables['time'] except KeyError as e: try: lat_o = nc.variables['lat'][:] lon_o = nc.variables['lon'][:] n_lat = len(lat_o) - time = nc.variables['time'] except KeyError as e: error_exit("{0} variable not found in {1} file.".format(str(e), path)) + + try: + time = nc.variables['time'] + except KeyError as e: + error_exit("{0} variable not found in {1} file.".format(str(e), path)) # From time array to list of dates. time_array = num2date(time[:], time.units, CALENDAR_STANDARD) + + if climatology: + time_array = np.array([time_aux.replace(year=date_array[0].year) for time_aux in time_array]) + i_time = np.where(time_array == date_array[0])[0][0] # Correction to set the longitudes from -180 to 180 instead of from 0 to 360. @@ -156,7 +207,7 @@ class IoNetcdf(IoServer): lon = np.array([lon_o[:]] * len(lat_o[:])).flatten() # del lat_o, lon_o - # Reads the var variable of the xone and the times needed. + # Reads the var variable of the one and the times needed. try: var = nc.variables[var_name][i_time:i_time + (len(date_array)), j_min:j_max, i_min:i_max] except KeyError as e: @@ -166,10 +217,25 @@ class IoNetcdf(IoServer): # That condition is fot the cases that the needed temperature is in a different NetCDF. while len(var) < len(date_array): aux_date = date_array[len(var)] + climatology = False path = os.path.join(netcdf_dir, '{0}_{1}{2}.nc'.format(var_name, aux_date.year, str(aux_date.month).zfill(2))) - # self.logger.write_log('Getting {0} from {1}'.format(var_name, path), message_level=2) - check_files(path) + if not os.path.exists(path): + path_climatology = os.path.join(netcdf_dir, '{0}_{1}{2}.nc'.format(var_name, 'climatology', + str(date_array[0].month).zfill(2))) + if not os.path.exists(path_climatology): + error_message = "*ERROR* (Rank {0}) File/s not found:".format(MPI.COMM_WORLD.Get_rank()) + error_message += "\n\t{0}\n\tneither {1}".format(path, path_climatology) + error_exit(error_message) + else: + error_message = "*WARNING* (Rank {0}) File/s not found:".format(MPI.COMM_WORLD.Get_rank()) + error_message += "\n\t{0}\n\t using {1}".format(path, path_climatology) + print(error_message) + sys.stdout.flush() + warn(error_message) + + climatology = True + path = path_climatology nc = Dataset(path, mode='r') i_time = 0 try: diff --git a/hermesv3_bu/sectors/residential_sector.py b/hermesv3_bu/sectors/residential_sector.py index 4f34981..8ccad2f 100755 --- a/hermesv3_bu/sectors/residential_sector.py +++ b/hermesv3_bu/sectors/residential_sector.py @@ -318,7 +318,8 @@ class ResidentialSector(Sector): daily_distribution.drop(columns=['centroid', 'REC', 'geometry_y'], axis=1, inplace=True) daily_distribution.rename(columns={'geometry_x': 'geometry'}, inplace=True) - for fuel in self.fuel_list: + for i, fuel in enumerate(self.fuel_list): + self.logger.write_log("\t\t{0} {1}/{2}".format(fuel, i+1, len(self.fuel_list)), message_level=3) # Selection of factor for HDD as a function of fuel type if fuel.startswith('B_'): hdd_f = self.hdd_f_biomass @@ -341,6 +342,7 @@ class ResidentialSector(Sector): daily_distribution = {} for day in self.day_dict.keys(): + self.logger.write_log("\t{0}".format(day), message_level=3) daily_distribution[day] = self.calculate_daily_distribution(day) self.logger.write_time_log('ResidentialSector', 'get_fuel_distribution_by_day', diff --git a/hermesv3_bu/sectors/traffic_sector.py b/hermesv3_bu/sectors/traffic_sector.py index 610518a..35e1e30 100755 --- a/hermesv3_bu/sectors/traffic_sector.py +++ b/hermesv3_bu/sectors/traffic_sector.py @@ -9,7 +9,7 @@ import geopandas as gpd from geopandas import GeoDataFrame import numpy as np from datetime import timedelta -import warnings +from warnings import warn from hermesv3_bu.logger.log import Log from hermesv3_bu.sectors.sector import Sector from hermesv3_bu.io_server.io_netcdf import IoNetcdf @@ -556,7 +556,7 @@ class TrafficSector(Sector): df.drop(columns=['Copert_V_name'], inplace=True) except IOError: self.logger.write_log('WARNING! No mileage correction applied to {0}'.format(pollutant_name)) - warnings.warn('No mileage correction applied to {0}'.format(pollutant_name)) + warn('No mileage correction applied to {0}'.format(pollutant_name)) df = None if downcasting: @@ -841,7 +841,7 @@ class TrafficSector(Sector): if len(resta_1) > 0: self.logger.write_log('WARNING! Exists some fleet codes that not appear on the EF file: {0}'.format( resta_1)) - warnings.warn('Exists some fleet codes that not appear on the EF file: {0}'.format(resta_1), Warning) + warn('Exists some fleet codes that not appear on the EF file: {0}'.format(resta_1), Warning) if len(resta_2) > 0: error_exit('Exists some fleet codes duplicated on the EF file: {0}'.format(resta_2)) @@ -1062,7 +1062,7 @@ class TrafficSector(Sector): else: self.logger.write_log("nmvoc emissions cannot be estimated because voc or ch4 are not selected in " + "the pollutant list.") - warnings.warn( + warn( "nmvoc emissions cannot be estimated because voc or ch4 are not selected in the pollutant list.") compacted = self.speciate_traffic(expanded, self.hot_cold_speciation) @@ -1178,9 +1178,16 @@ class TrafficSector(Sector): road_link_aux['centroid'] = road_link_aux['geometry'].centroid link_lons = road_link_aux['geometry'].centroid.x link_lats = road_link_aux['geometry'].centroid.y - - p_factor = self.calculate_precipitation_factor(link_lons.min(), link_lons.max(), link_lats.min(), - link_lats.max(), self.precipitation_path) + try: + p_factor = self.calculate_precipitation_factor(link_lons.min(), link_lons.max(), link_lats.min(), + link_lats.max(), self.precipitation_path) + except FileNotFoundError as e: + self.resuspension_correction = False + message = "*WARNING* The the rain correction was not applied due to missing time steps in the " + \ + "precipitation dataset {0}".format(e) + self.logger.write_log(message) + warn(message) + if self.resuspension_correction: unary_union = p_factor.unary_union road_link_aux['REC'] = road_link_aux.apply(self.nearest, geom_union=unary_union, df1=road_link_aux, @@ -1239,16 +1246,20 @@ class TrafficSector(Sector): df.rename(columns={p_name: p_name_new}, inplace=True) pollutants_renamed.append(p_name_new) - df_aux = df[['Link_ID', 'Fleet_Code'] + pollutants_renamed] + df_aux = df[['Link_ID', 'Fleet_Code'] + pollutants_renamed].copy() df_aux['tstep'] = tstep df_list.append(df_aux) df.drop(columns=pollutants_renamed, inplace=True) - del df_aux + # del df_aux libc.malloc_trim(0) gc.collect() df = pd.concat(df_list, ignore_index=True) + for to_delete in df_list: + del to_delete + libc.malloc_trim(0) + gc.collect() self.logger.write_time_log('TrafficSector', 'transform_df', timeit.default_timer() - spent_time) return df -- GitLab From 364985b7c45094b7e3b13471051a4ddd0a71fe66 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Thu, 11 Jun 2020 16:09:00 +0200 Subject: [PATCH 3/4] Corrected bug on February Climatology --- conf/hermes.conf | 6 +++--- hermesv3_bu/io_server/io_netcdf.py | 5 +++++ 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/conf/hermes.conf b/conf/hermes.conf index 4be432e..189868e 100755 --- a/conf/hermes.conf +++ b/conf/hermes.conf @@ -3,9 +3,9 @@ log_level = 3 input_dir = /home/Earth/ctena/Models/hermesv3_bu_data data_path = /esarchive/recon output_dir = /scratch/Earth/ctena/HERMESv3_BU_OUT/climatology -output_name = HERMESv3__climatology_crop_fertilizers.nc +output_name = HERMESv3__climatology_livestock.nc emission_summary = 0 -start_date = 2020/11/29 00:00:00 +start_date = 2022/02/02 00:00:00 # ----- end_date = start_date [DEFAULT] ----- # end_date = 2010/01/01 00:00:00 output_timestep_num = 24 @@ -132,7 +132,7 @@ writing_processors = 1 # # aviation_processors = 1 # shipping_port_processors = 1 -# livestock_processors = 4 +livestock_processors = 4 # crop_operations_processors = 1 # crop_fertilizers_processors = 4 # agricultural_machinery_processors = 1 diff --git a/hermesv3_bu/io_server/io_netcdf.py b/hermesv3_bu/io_server/io_netcdf.py index c2bd7e5..a615729 100755 --- a/hermesv3_bu/io_server/io_netcdf.py +++ b/hermesv3_bu/io_server/io_netcdf.py @@ -10,6 +10,7 @@ from netCDF4 import Dataset from shapely.geometry import Point from cf_units import num2date, CALENDAR_STANDARD from geopandas import GeoDataFrame +from calendar import isleap from hermesv3_bu.io_server.io_server import IoServer from hermesv3_bu.tools.checker import check_files, error_exit @@ -84,6 +85,8 @@ class IoNetcdf(IoServer): time_array = num2date(time[:], time.units, CALENDAR_STANDARD) if climatology: + if not isleap(year=date.year) and date.month == 2: + time_array = time_array[:28] time_array = np.array([aux.date().replace(year=date.year) for aux in time_array]) else: time_array = np.array([aux.date() for aux in time_array]) @@ -187,6 +190,8 @@ class IoNetcdf(IoServer): time_array = num2date(time[:], time.units, CALENDAR_STANDARD) if climatology: + if not isleap(year=date_array[0].year) and date_array[0].month == 2: + time_array = time_array[:28*24] time_array = np.array([time_aux.replace(year=date_array[0].year) for time_aux in time_array]) i_time = np.where(time_array == date_array[0])[0][0] -- GitLab From c8a76cd9a0d5ead9b2624ac0e999cb4e68ec9d6d Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Wed, 9 Dec 2020 15:33:11 +0100 Subject: [PATCH 4/4] Added climatology --- CHANGELOG | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG b/CHANGELOG index de2b671..81ff954 100755 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,3 +1,8 @@ +1.0.5 + 2020/12/09 + + - Added Climatic meteorology + 1.0.4 2020/07/17 -- GitLab