diff --git a/.codacy.yml b/.codacy.yml index bf12fc2d77140a852ff1b9197370766a0ecec7c3..0bc0d39b0f2175fbc21f61847dace53b65460c48 100755 --- a/.codacy.yml +++ b/.codacy.yml @@ -13,7 +13,7 @@ engines: enabled: true pylint: enabled: true - python_version: 2 + python_version: 3 exclude_paths: [ 'doc/**', diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 08c1355a02d158efce8b3c8438df86fdc3275b6e..d4676c915a2d3b2e503a4a7e293c64e664d1875b 100755 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -20,11 +20,11 @@ prepare: test_python2: stage: test script: - - conda env update -f environment.yml -n hermesv3_gr python=2.7 + - conda env update -f environment.yml -n hermesv3_gr python=3.7 - source activate hermesv3_gr - python run_test.py - pip install codacy-coverage --upgrade - - python-codacy-coverage -r tests/report/python2/coverage.xml + - python-codacy-coverage -r tests/report/python3/coverage.xml #test_python3: # stage: test diff --git a/CHANGELOG b/CHANGELOG index 7d4bbafa65b02f3ce2e21662713287f5f90fbc05..947c8d87748f019855defcc52447189e43c1317a 100755 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,3 +1,8 @@ +2.0.0 + - Python 3 + - first_time option + - erase_auxiliary_files option + 1.0.4 2019/07/19 - Specified version on timezonefinder to be less than 4.0.0 (not further support to python 2.7.X) diff --git a/conf/hermes.conf b/conf/hermes.conf index a9d139395600ea7e6bdb7c016f6b4d649469b192..7a125976b8999b572c72dcc012ebf68037b17ee4 100755 --- a/conf/hermes.conf +++ b/conf/hermes.conf @@ -11,6 +11,8 @@ start_date = 2017/01/01 00:00:00 output_timestep_type = hourly output_timestep_num = 24 output_timestep_freq = 1 +first_time = 0 +erase_auxiliary_files = 1 [DOMAIN] @@ -27,7 +29,7 @@ domain_type = global # domain_type = rotated # domain_type = mercator vertical_description = /data/profiles/vertical/Benchmark_15layers_vertical_description.csv -auxiliar_files_path = /data/auxiliar_files/_ +auxiliary_files_path = /data/auxiliar_files/_ # if domain_type == global: inc_lat = 1. diff --git a/environment.yml b/environment.yml index bebcce37b2e46490fde22e27f8e5d1a8de3ebaf5..dc46ab7a228c1d91821ffda6d4e2e27ca5f5abc3 100755 --- a/environment.yml +++ b/environment.yml @@ -6,7 +6,7 @@ channels: - conda-forge dependencies: - - python = 2 + - python = 3 - numpy - netcdf4 >= 1.3.1 - python-cdo >= 1.3.6 diff --git a/hermesv3_gr/config/config.py b/hermesv3_gr/config/config.py index f8d1d4acdde44858f55501290fcc387a763110eb..7775770f7c6802c6534e0a658724e0f76b17b193 100755 --- a/hermesv3_gr/config/config.py +++ b/hermesv3_gr/config/config.py @@ -20,13 +20,20 @@ from configargparse import ArgParser import sys +import os +from shutil import rmtree +from mpi4py import MPI class Config(ArgParser): """ Initialization of the arguments that the parser can handle. """ - def __init__(self): + def __init__(self, comm=None): + if comm is None: + comm = MPI.COMM_WORLD + self.comm = comm + super(Config, self).__init__() self.options = self.read_options() @@ -61,6 +68,10 @@ class Config(ArgParser): type=str, choices=['hourly', 'daily', 'monthly', 'yearly']) p.add_argument('--output_timestep_num', required=True, help='Number of timesteps to simulate.', type=int) p.add_argument('--output_timestep_freq', required=True, help='Frequency between timesteps.', type=int) + p.add_argument('--first_time', required=False, default='False', type=str, + help='Indicates if you want to run it for first time (only create auxiliary files).') + p.add_argument('--erase_auxiliary_files', required=False, default='False', type=str, + help='Indicates if you want to start from scratch removing the auxiliary files already created.') p.add_argument('--output_model', required=True, help='Name of the output model.', choices=['MONARCH', 'CMAQ', 'WRF_CHEM']) @@ -69,7 +80,7 @@ class Config(ArgParser): p.add_argument('--domain_type', required=True, help='Type of domain to simulate.', choices=['global', 'lcc', 'rotated', 'mercator', 'regular']) - p.add_argument('--auxiliar_files_path', required=True, + p.add_argument('--auxiliary_files_path', required=True, help='Path to the directory where the necessary auxiliary files will be created if them are ' + 'not created yet.') @@ -177,8 +188,16 @@ class Config(ArgParser): arguments.start_date = self._parse_start_date(arguments.start_date) arguments.end_date = self._parse_end_date(arguments.end_date, arguments.start_date) + arguments.first_time = self._parse_bool(arguments.first_time) + arguments.erase_auxiliary_files = self._parse_bool(arguments.erase_auxiliary_files) + self.create_dir(arguments.output_dir) - self.create_dir(arguments.auxiliar_files_path) + if arguments.erase_auxiliary_files: + if os.path.exists(arguments.auxiliary_files_path): + if self.comm.Get_rank() == 0: + rmtree(arguments.auxiliary_files_path) + self.comm.Barrier() + self.create_dir(arguments.auxiliary_files_path) return arguments @@ -207,8 +226,7 @@ class Config(ArgParser): full_path = os.path.join(self.options.output_dir, file_name) return full_path - @staticmethod - def create_dir(path): + def create_dir(self, path): """ Create the given folder if it is not created yet. @@ -217,8 +235,7 @@ class Config(ArgParser): """ import os from mpi4py import MPI - icomm = MPI.COMM_WORLD - comm = icomm.Split(color=0, key=0) + comm = self.comm.Split(color=0, key=0) rank = comm.Get_rank() if rank == 0: @@ -304,12 +321,12 @@ class Config(ArgParser): else: return self._parse_start_date(end_date) - def set_log_level(self): + def set_log_level(self, comm): """ Defines the log_level using the common script settings. """ from . import settings - settings.define_global_vars(self.options.log_level) + settings.define_global_vars(self.options.log_level, comm) if __name__ == '__main__': diff --git a/hermesv3_gr/config/settings.py b/hermesv3_gr/config/settings.py index 18b6a6a65a2bfd481fcb7c1d8c48389e6f094cc7..038f004e3febd6d15bd8bb3ccc08601e1021a05a 100755 --- a/hermesv3_gr/config/settings.py +++ b/hermesv3_gr/config/settings.py @@ -46,7 +46,7 @@ global log_file global df_times -def define_global_vars(in_log_level): +def define_global_vars(in_log_level, mycomm): # TODO Documentation from mpi4py import MPI @@ -55,7 +55,11 @@ def define_global_vars(in_log_level): global rank global size - icomm = MPI.COMM_WORLD + if mycomm is None: + icomm = MPI.COMM_WORLD + else: + icomm = mycomm + comm = icomm.Split(color=0, key=0) rank = comm.Get_rank() size = comm.Get_size() diff --git a/hermesv3_gr/hermes.py b/hermesv3_gr/hermes.py index 52d46cba127a6672cd2439c69a37ee53b23d0e22..bb047b9e7f8dd0f0ddc82cba836ea3248f09b8e7 100755 --- a/hermesv3_gr/hermes.py +++ b/hermesv3_gr/hermes.py @@ -34,13 +34,16 @@ class Hermes(object): """ Interface class for HERMESv3. """ - def __init__(self, config, new_date=None): + def __init__(self, config, new_date=None, comm=None): from hermesv3_gr.modules.grids.grid import Grid from hermesv3_gr.modules.temporal.temporal import TemporalDistribution from hermesv3_gr.modules.writing.writer import Writer global full_time st_time = full_time = timeit.default_timer() + if comm is None: + comm = MPI.COMM_WORLD + self.comm = comm self.config = config self.options = config.options @@ -48,7 +51,7 @@ class Hermes(object): if new_date is not None: self.options.start_date = new_date - config.set_log_level() + config.set_log_level(self.comm) settings.define_log_file(self.options.output_dir, self.options.start_date) settings.define_times_file() @@ -64,8 +67,9 @@ class Hermes(object): self.levels = VerticalDistribution.get_vertical_output_profile(self.options.vertical_description) self.grid = Grid.select_grid( + self.comm, self.options.domain_type, self.options.vertical_description, self.options.output_timestep_num, - self.options.auxiliar_files_path, self.options.inc_lat, self.options.inc_lon, self.options.lat_orig, + self.options.auxiliary_files_path, self.options.inc_lat, self.options.inc_lon, self.options.lat_orig, self.options.lon_orig, self.options.n_lat, self.options.n_lon, self.options.centre_lat, self.options.centre_lon, self.options.west_boundary, self.options.south_boundary, self.options.inc_rlat, self.options.inc_rlon, self.options.lat_1, self.options.lat_2, self.options.lon_0, self.options.lat_0, @@ -89,47 +93,53 @@ class Hermes(object): settings.write_time('HERMES', 'Init', timeit.default_timer() - st_time, level=1) # @profile - def main(self): + def main(self, return_emis=False): """ Main functionality of the model. """ from datetime import timedelta - st_time = timeit.default_timer() - settings.write_log('') - settings.write_log('***** Starting HERMESv3 *****') - num = 1 - for ei in self.emission_list: - settings.write_log('Processing emission inventory {0} for the sector {1} ({2}/{3}):'.format( - ei.inventory_name, ei.sector, num, len(self.emission_list))) - num += 1 - - ei.do_regrid() - - if ei.vertical is not None: - settings.write_log("\tCalculating vertical distribution.", level=2) - if ei.source_type == 'area': - ei.vertical_factors = ei.vertical.calculate_weights() - elif ei.source_type == 'point': - ei.calculate_altitudes(self.options.vertical_description) - ei.point_source_by_cell() - # To avoid use point source as area source when is going to apply vertical factors while writing. - ei.vertical = None - else: - settings.write_log('ERROR: Check the .err file to get more info.') - if settings.rank == 0: - raise AttributeError('Unrecognized emission source type {0}'.format(ei.source_type)) - sys.exit(1) - - if ei.temporal is not None: - ei.temporal_factors = ei.temporal.calculate_3d_temporal_factors() - if ei.speciation is not None: - ei.emissions = ei.speciation.do_speciation(ei.emissions) - - self.writer.write(self.emission_list) - - settings.write_log("***** HERMESv3 simulation finished successfully *****") - settings.write_time('HERMES', 'main', timeit.default_timer() - st_time) + if self.options.first_time: + # Stop run + settings.write_log('***** HERMESv3_BU First Time finished successfully *****') + else: + st_time = timeit.default_timer() + settings.write_log('') + settings.write_log('***** Starting HERMESv3_GR *****') + num = 1 + for ei in self.emission_list: + settings.write_log('Processing emission inventory {0} for the sector {1} ({2}/{3}):'.format( + ei.inventory_name, ei.sector, num, len(self.emission_list))) + num += 1 + + ei.do_regrid() + + if ei.vertical is not None: + settings.write_log("\tCalculating vertical distribution.", level=2) + if ei.source_type == 'area': + ei.vertical_factors = ei.vertical.calculate_weights() + elif ei.source_type == 'point': + ei.calculate_altitudes(self.options.vertical_description) + ei.point_source_by_cell() + # To avoid use point source as area source when is going to apply vertical factors while writing + ei.vertical = None + else: + settings.write_log('ERROR: Check the .err file to get more info.') + if settings.rank == 0: + raise AttributeError('Unrecognized emission source type {0}'.format(ei.source_type)) + sys.exit(1) + + if ei.temporal is not None: + ei.temporal_factors = ei.temporal.calculate_3d_temporal_factors() + if ei.speciation is not None: + ei.emissions = ei.speciation.do_speciation(ei.emissions) + if return_emis: + return self.writer.get_emis(self.emission_list) + else: + self.writer.write(self.emission_list) + + settings.write_log("***** HERMESv3_GR simulation finished successfully *****") + settings.write_time('HERMES', 'main', timeit.default_timer() - st_time) settings.write_time('HERMES', 'TOTAL', timeit.default_timer() - full_time) settings.finish_logs(self.options.output_dir, self.options.start_date) diff --git a/hermesv3_gr/modules/emision_inventories/emission_inventory.py b/hermesv3_gr/modules/emision_inventories/emission_inventory.py index 43e0db7b102637b11e6d43791960a948e182c47f..167042479a4a6f79db5b1c727b89eec59dfd8e35 100755 --- a/hermesv3_gr/modules/emision_inventories/emission_inventory.py +++ b/hermesv3_gr/modules/emision_inventories/emission_inventory.py @@ -103,7 +103,7 @@ class EmissionInventory(object): # It will also create the WoldMasks necessaries self.masking = Masking( options.world_info, factors, regrid_mask, grid, - world_mask_file=os.path.join(os.path.dirname(options.auxiliar_files_path), + world_mask_file=os.path.join(os.path.dirname(options.auxiliary_files_path), '{0}_WorldMask.nc'.format(inventory_name))) self.input_pollutants = pollutants @@ -116,19 +116,19 @@ class EmissionInventory(object): if self.source_type == 'area': self.regrid = ConservativeRegrid( self.pollutant_dicts, - os.path.join(options.auxiliar_files_path, + os.path.join(options.auxiliary_files_path, "Weight_Matrix_{0}_{1}.nc".format(self.inventory_name, settings.size)), grid, masking=self.masking) - - # Creating Vertical Object - if p_vertical is not None: - self.vertical = VerticalDistribution( - self.get_profile(p_vertical), vertical_profile_path=options.p_vertical, - vertical_output_profile=vertical_output_profile) - else: - self.vertical = None - settings.write_log('\t\tNone vertical profile set.', level=2) - self.vertical_factors = None + if not options.first_time: + # Creating Vertical Object + if p_vertical is not None: + self.vertical = VerticalDistribution( + self.get_profile(p_vertical), vertical_profile_path=options.p_vertical, + vertical_output_profile=vertical_output_profile) + else: + self.vertical = None + settings.write_log('\t\tNone vertical profile set.', level=2) + self.vertical_factors = None # Creating Temporal Object # It will also create the necessaries timezone files @@ -136,18 +136,18 @@ class EmissionInventory(object): self.temporal = TemporalDistribution( current_date, options.output_timestep_type, options.output_timestep_num, options.output_timestep_freq, options.p_month, p_month, options.p_week, p_week, options.p_day, p_day, options.p_hour, p_hour, - options.world_info, options.auxiliar_files_path, grid) + options.world_info, options.auxiliary_files_path, grid, options.first_time) else: self.temporal = None settings.write_log('\t\tNone temporal profile set.', level=2) self.temporal_factors = None - - # Creating Speciation Object - if p_speciation is not None: - self.speciation = Speciation(p_speciation, options.p_speciation, options.molecular_weights) - else: - self.speciation = None - settings.write_log('\t\tNone speciation profile set.', level=2) + if not options.first_time: + # Creating Speciation Object + if p_speciation is not None: + self.speciation = Speciation(p_speciation, options.p_speciation, options.molecular_weights) + else: + self.speciation = None + settings.write_log('\t\tNone speciation profile set.', level=2) self.vertical_weights = None diff --git a/hermesv3_gr/modules/emision_inventories/gfas_emission_inventory.py b/hermesv3_gr/modules/emision_inventories/gfas_emission_inventory.py index 33a532844a47aafc823a39759ebf0474763a6ad1..c06261561fa526b69de98473bb325f2aaec0914d 100755 --- a/hermesv3_gr/modules/emision_inventories/gfas_emission_inventory.py +++ b/hermesv3_gr/modules/emision_inventories/gfas_emission_inventory.py @@ -82,13 +82,13 @@ class GfasEmissionInventory(EmissionInventory): vertical_output_profile, reference_year=reference_year, factors=factors, regrid_mask=regrid_mask, p_vertical=None, p_month=p_month, p_week=p_week, p_day=p_day, p_hour=p_hour, p_speciation=p_speciation) + if not options.first_time: + self.approach = self.get_approach(p_vertical) + self.method = self.get_method(p_vertical) - self.approach = self.get_approach(p_vertical) - self.method = self.get_method(p_vertical) + self.altitude = self.get_altitude() - self.altitude = self.get_altitude() - - self.vertical = GfasVerticalDistribution(vertical_output_profile, self.approach, self.get_altitude()) + self.vertical = GfasVerticalDistribution(vertical_output_profile, self.approach, self.get_altitude()) settings.write_time('GFAS_EmissionInventory', 'Init', timeit.default_timer() - st_time, level=3) diff --git a/hermesv3_gr/modules/grids/grid.py b/hermesv3_gr/modules/grids/grid.py index dca201366752607cc310f379c699d3129dea1273..520b56b800ddbd25c60e7c3cbc78dbdc35474df0 100755 --- a/hermesv3_gr/modules/grids/grid.py +++ b/hermesv3_gr/modules/grids/grid.py @@ -40,10 +40,10 @@ class Grid(object): :type temporal_path: str """ - def __init__(self, grid_type, vertical_description_path, temporal_path): + def __init__(self, grid_type, vertical_description_path, temporal_path, comm): st_time = timeit.default_timer() # settings.write_log('Creating Grid...', level=1) - + self.comm = comm # Defining class atributes self.procs_array = None self.nrows = 0 @@ -66,7 +66,7 @@ class Grid(object): self.temporal_path = temporal_path self.shapefile_path = None - self.esmf_grid = None + # self.esmf_grid = None self.x_lower_bound = None self.x_upper_bound = None self.y_lower_bound = None @@ -78,6 +78,27 @@ class Grid(object): settings.write_time('Grid', 'Init', timeit.default_timer() - st_time, level=1) + def calculate_bounds(self): + x_min = 0 + y_min = 0 + if len(self.center_latitudes.shape) == 1: + x_max = self.center_latitudes.shape[0] + y_max = self.center_longitudes.shape[0] + else: + x_max = self.center_latitudes.shape[0] + y_max = self.center_longitudes.shape[1] + y_len = y_max // self.comm.Get_size() + exceeded_ranks = y_max % self.comm.Get_size() + self.x_lower_bound = x_min + self.x_upper_bound = x_max + if self.comm.Get_rank() < exceeded_ranks: + y_len = y_len + 1 + self.y_lower_bound = self.comm.Get_rank() * y_len + self.y_upper_bound = (self.comm.Get_rank() + 1) * y_len + else: + self.y_lower_bound = (self.comm.Get_rank() * y_len) + exceeded_ranks + self.y_upper_bound = ((self.comm.Get_rank() + 1) * y_len) + exceeded_ranks + @staticmethod def create_esmf_grid_from_file(file_name, sphere=True): import ESMF @@ -94,7 +115,7 @@ class Grid(object): return grid @staticmethod - def select_grid(grid_type, vertical_description_path, timestep_num, temporal_path, inc_lat, inc_lon, lat_orig, + def select_grid(comm, grid_type, vertical_description_path, timestep_num, temporal_path, inc_lat, inc_lon, lat_orig, lon_orig, n_lat, n_lon, centre_lat, centre_lon, west_boundary, south_boundary, inc_rlat, inc_rlon, lat_1, lat_2, lon_0, lat_0, nx, ny, inc_x, inc_y, x_0, y_0, lat_ts): # TODO describe better the rotated parameters @@ -180,26 +201,27 @@ class Grid(object): # Creating a different object depending on the grid type if grid_type == 'global': from hermesv3_gr.modules.grids.grid_global import GlobalGrid - grid = GlobalGrid(grid_type, vertical_description_path, timestep_num, temporal_path, inc_lat, inc_lon) + grid = GlobalGrid(grid_type, vertical_description_path, timestep_num, temporal_path, inc_lat, inc_lon, + comm=comm) elif grid_type == 'regular': from hermesv3_gr.modules.grids.grid_latlon import LatLonGrid grid = LatLonGrid(grid_type, vertical_description_path, timestep_num, temporal_path, inc_lat, inc_lon, - lat_orig, lon_orig, n_lat, n_lon) + lat_orig, lon_orig, n_lat, n_lon, comm=comm) elif grid_type == 'rotated': from hermesv3_gr.modules.grids.grid_rotated import RotatedGrid grid = RotatedGrid(grid_type, vertical_description_path, timestep_num, temporal_path, - centre_lat, centre_lon, west_boundary, south_boundary, inc_rlat, inc_rlon) + centre_lat, centre_lon, west_boundary, south_boundary, inc_rlat, inc_rlon, comm=comm) elif grid_type == 'lcc': from hermesv3_gr.modules.grids.grid_lcc import LccGrid grid = LccGrid(grid_type, vertical_description_path, timestep_num, temporal_path, lat_1, lat_2, lon_0, - lat_0, nx, ny, inc_x, inc_y, x_0, y_0) + lat_0, nx, ny, inc_x, inc_y, x_0, y_0, comm=comm) elif grid_type == 'mercator': from hermesv3_gr.modules.grids.grid_mercator import MercatorGrid grid = MercatorGrid(grid_type, vertical_description_path, timestep_num, temporal_path, lat_ts, lon_0, - nx, ny, inc_x, inc_y, x_0, y_0) + nx, ny, inc_x, inc_y, x_0, y_0, comm=comm) else: settings.write_log('ERROR: Check the .err file to get more info.') if settings.rank == 0: @@ -389,8 +411,9 @@ class Grid(object): st_time = timeit.default_timer() settings.write_log('\t\tGetting 2D coordinates from ESMPy Grid', level=3) - lat = self.esmf_grid.get_coords(1, ESMF.StaggerLoc.CENTER).T - lon = self.esmf_grid.get_coords(0, ESMF.StaggerLoc.CENTER).T + esmf_grid = self.create_esmf_grid_from_file(self.coords_netcdf_file) + lat = esmf_grid.get_coords(1, ESMF.StaggerLoc.CENTER).T + lon = esmf_grid.get_coords(0, ESMF.StaggerLoc.CENTER).T settings.write_time('Grid', 'get_coordinates_2d', timeit.default_timer() - st_time, level=3) @@ -427,21 +450,13 @@ class Grid(object): settings.write_log('\t\tGrid shapefile not done. Lets try to create it.', level=3) # Create Shapefile - # Use the meters coordiantes to create the shapefile - y = self.boundary_latitudes x = self.boundary_longitudes - # sys.exit() if self.grid_type in ['global', 'regular']: x = x.reshape((x.shape[1], x.shape[2])) y = y.reshape((y.shape[1], y.shape[2])) - # x_aux = np.empty((x.shape[0], y.shape[0], 4)) - # x_aux[:, :, 0] = x[:, np.newaxis, 0] - # x_aux[:, :, 1] = x[:, np.newaxis, 1] - # x_aux[:, :, 2] = x[:, np.newaxis, 1] - # x_aux[:, :, 3] = x[:, np.newaxis, 0] aux_shape = (y.shape[0], x.shape[0], 4) x_aux = np.empty(aux_shape) x_aux[:, :, 0] = x[np.newaxis, :, 0] @@ -452,12 +467,6 @@ class Grid(object): x = x_aux del x_aux - # y_aux = np.empty((x.shape[0], y.shape[0], 4)) - # y_aux[:, :, 0] = y[np.newaxis, :, 0] - # y_aux[:, :, 1] = y[np.newaxis, :, 0] - # y_aux[:, :, 2] = y[np.newaxis, :, 1] - # y_aux[:, :, 3] = y[np.newaxis, :, 1] - y_aux = np.empty(aux_shape) y_aux[:, :, 0] = y[:, np.newaxis, 0] y_aux[:, :, 1] = y[:, np.newaxis, 0] @@ -467,8 +476,6 @@ class Grid(object): y = y_aux del y_aux - # exit() - if not full_grid: y = y[self.x_lower_bound:self.x_upper_bound, self.y_lower_bound:self.y_upper_bound, :] x = x[self.x_lower_bound:self.x_upper_bound, self.y_lower_bound:self.y_upper_bound, :] @@ -476,29 +483,12 @@ class Grid(object): aux_b_lats = y.reshape((y.shape[0] * y.shape[1], y.shape[2])) aux_b_lons = x.reshape((x.shape[0] * x.shape[1], x.shape[2])) - # The regular lat-lon projection has only 2 (laterals) points for each cell instead of 4 (corners) - # if aux_b_lats.shape[1] == 2: - # aux_b = np.empty((aux_b_lats.shape[0], 4)) - # aux_b[:, 0] = aux_b_lats[:, 0] - # aux_b[:, 1] = aux_b_lats[:, 0] - # aux_b[:, 2] = aux_b_lats[:, 1] - # aux_b[:, 3] = aux_b_lats[:, 1] - # aux_b_lats = aux_b - # - # if aux_b_lons.shape[1] == 2: - # aux_b = np.empty((aux_b_lons.shape[0], 4)) - # aux_b[:, 0] = aux_b_lons[:, 0] - # aux_b[:, 1] = aux_b_lons[:, 1] - # aux_b[:, 2] = aux_b_lons[:, 1] - # aux_b[:, 3] = aux_b_lons[:, 0] - # aux_b_lons = aux_b - # Create one dataframe with 8 columns, 4 points with two coordinates each one df_lats = pd.DataFrame(aux_b_lats, columns=['b_lat_1', 'b_lat_2', 'b_lat_3', 'b_lat_4']) df_lons = pd.DataFrame(aux_b_lons, columns=['b_lon_1', 'b_lon_2', 'b_lon_3', 'b_lon_4']) df = pd.concat([df_lats, df_lons], axis=1) - # Substituate 8 columns by 4 with the two coordinates + # Substitute 8 columns by 4 with the two coordinates df['p1'] = list(zip(df.b_lon_1, df.b_lat_1)) del df['b_lat_1'], df['b_lon_1'] df['p2'] = list(zip(df.b_lon_2, df.b_lat_2)) @@ -511,18 +501,11 @@ class Grid(object): # Make a list of list of tuples # [[(point_1.1), (point_1.2), (point_1.3), (point_1.4)], # [(point_2.1), (point_2.2), (point_2.3), (point_2.4)], ...] - list_points = df.as_matrix() + list_points = df.values del df['p1'], df['p2'], df['p3'], df['p4'] # List of polygons from the list of points geometry = [Polygon(list(points)) for points in list_points] - # geometry = [] - # for point in list_points: - # print point - # geometry.append(Polygon(list(point))) - # print geometry[0] - # sys.exit() - # print len(geometry), len(df), gdf = gpd.GeoDataFrame(df, crs={'init': 'epsg:4326'}, geometry=geometry) gdf = gdf.to_crs(self.crs) diff --git a/hermesv3_gr/modules/grids/grid_global.py b/hermesv3_gr/modules/grids/grid_global.py index a3ca2368e8134a9bf3690c460597fc2a803b6c6b..2f63c3cb213105ac43e778f5f1c04c6162e7af26 100755 --- a/hermesv3_gr/modules/grids/grid_global.py +++ b/hermesv3_gr/modules/grids/grid_global.py @@ -54,14 +54,13 @@ class GlobalGrid(Grid): """ def __init__(self, grid_type, vertical_description_path, timestep_num, temporal_path, inc_lat, inc_lon, - center_longitude=float(0)): - import ESMF + center_longitude=float(0), comm=None): st_time = timeit.default_timer() settings.write_log('\tCreating Global grid.', level=2) # Initialize the class using parent - super(GlobalGrid, self).__init__(grid_type, vertical_description_path, temporal_path) + super(GlobalGrid, self).__init__(grid_type, vertical_description_path, temporal_path, comm) self.center_lat = float(0) self.center_lon = center_longitude @@ -76,12 +75,7 @@ class GlobalGrid(Grid): super(GlobalGrid, self).write_coords_netcdf() settings.comm.Barrier() - self.esmf_grid = super(GlobalGrid, self).create_esmf_grid_from_file(self.coords_netcdf_file) - - self.x_lower_bound = self.esmf_grid.lower_bounds[ESMF.StaggerLoc.CENTER][1] - self.x_upper_bound = self.esmf_grid.upper_bounds[ESMF.StaggerLoc.CENTER][1] - self.y_lower_bound = self.esmf_grid.lower_bounds[ESMF.StaggerLoc.CENTER][0] - self.y_upper_bound = self.esmf_grid.upper_bounds[ESMF.StaggerLoc.CENTER][0] + self.calculate_bounds() self.shape = (timestep_num, len(self.vertical_description), self.x_upper_bound-self.x_lower_bound, self.y_upper_bound-self.y_lower_bound) diff --git a/hermesv3_gr/modules/grids/grid_latlon.py b/hermesv3_gr/modules/grids/grid_latlon.py index 0409d7bc63d90cbc3c99235f9d02c115593257a4..df0a676794c2367748905522dd2f817e0334cfa1 100755 --- a/hermesv3_gr/modules/grids/grid_latlon.py +++ b/hermesv3_gr/modules/grids/grid_latlon.py @@ -63,14 +63,13 @@ class LatLonGrid(Grid): """ def __init__(self, grid_type, vertical_description_path, timestep_num, temporal_path, inc_lat, inc_lon, lat_orig, - lon_orig, n_lat, n_lon): - import ESMF + lon_orig, n_lat, n_lon, comm=None): st_time = timeit.default_timer() settings.write_log('\tCreating Regional regular lat-lon grid.', level=2) # Initialize the class using parent - super(LatLonGrid, self).__init__(grid_type, vertical_description_path, temporal_path) + super(LatLonGrid, self).__init__(grid_type, vertical_description_path, temporal_path, comm) self.inc_lat = inc_lat self.inc_lon = inc_lon @@ -87,12 +86,7 @@ class LatLonGrid(Grid): self.write_coords_netcdf() settings.comm.Barrier() - self.esmf_grid = self.create_esmf_grid_from_file(self.coords_netcdf_file, sphere=False) - - self.x_lower_bound = self.esmf_grid.lower_bounds[ESMF.StaggerLoc.CENTER][1] - self.x_upper_bound = self.esmf_grid.upper_bounds[ESMF.StaggerLoc.CENTER][1] - self.y_lower_bound = self.esmf_grid.lower_bounds[ESMF.StaggerLoc.CENTER][0] - self.y_upper_bound = self.esmf_grid.upper_bounds[ESMF.StaggerLoc.CENTER][0] + self.calculate_bounds() self.shape = (timestep_num, len(self.vertical_description), self.x_upper_bound-self.x_lower_bound, self.y_upper_bound-self.y_lower_bound) diff --git a/hermesv3_gr/modules/grids/grid_lcc.py b/hermesv3_gr/modules/grids/grid_lcc.py index d360d74087dbc7c312ae0b8d887d747e8647b1df..42aaebcc32e3829f60ad72d693232624e0cf7e99 100755 --- a/hermesv3_gr/modules/grids/grid_lcc.py +++ b/hermesv3_gr/modules/grids/grid_lcc.py @@ -77,13 +77,12 @@ class LccGrid(Grid): """ def __init__(self, grid_type, vertical_description_path, timestep_num, temporal_path, lat_1, lat_2, lon_0, lat_0, - nx, ny, inc_x, inc_y, x_0, y_0, earth_radius=6370000.000): - import ESMF + nx, ny, inc_x, inc_y, x_0, y_0, earth_radius=6370000.000, comm=None): st_time = timeit.default_timer() settings.write_log('\tCreating Lambert Conformal Conic (LCC) grid.', level=2) # Initialises with parent class - super(LccGrid, self).__init__(grid_type, vertical_description_path, temporal_path) + super(LccGrid, self).__init__(grid_type, vertical_description_path, temporal_path, comm) # Setting parameters self.lat_1 = lat_1 @@ -113,12 +112,7 @@ class LccGrid(Grid): self.write_coords_netcdf() settings.comm.Barrier() - self.esmf_grid = super(LccGrid, self).create_esmf_grid_from_file(self.coords_netcdf_file, sphere=False) - # - self.x_lower_bound = self.esmf_grid.lower_bounds[ESMF.StaggerLoc.CENTER][1] - self.x_upper_bound = self.esmf_grid.upper_bounds[ESMF.StaggerLoc.CENTER][1] - self.y_lower_bound = self.esmf_grid.lower_bounds[ESMF.StaggerLoc.CENTER][0] - self.y_upper_bound = self.esmf_grid.upper_bounds[ESMF.StaggerLoc.CENTER][0] + self.calculate_bounds() self.shape = (timestep_num, len(self.vertical_description), self.x_upper_bound-self.x_lower_bound, self.y_upper_bound-self.y_lower_bound) diff --git a/hermesv3_gr/modules/grids/grid_mercator.py b/hermesv3_gr/modules/grids/grid_mercator.py index 38bb722d883faf92ebebcc4a64cf32d2feeadf7b..d18c32de25fe21f36737c24a6285149216e185ac 100755 --- a/hermesv3_gr/modules/grids/grid_mercator.py +++ b/hermesv3_gr/modules/grids/grid_mercator.py @@ -69,7 +69,6 @@ class MercatorGrid(Grid): def __init__(self, grid_type, vertical_description_path, timestep_num, temporal_path, lat_ts, lon_0, nx, ny, inc_x, inc_y, x_0, y_0, earth_radius=6370000.000): - import ESMF st_time = timeit.default_timer() settings.write_log('\tCreating Mercator grid.', level=2) @@ -101,12 +100,7 @@ class MercatorGrid(Grid): self.write_coords_netcdf() settings.comm.Barrier() - self.esmf_grid = super(MercatorGrid, self).create_esmf_grid_from_file(self.coords_netcdf_file, sphere=False) - # - self.x_lower_bound = self.esmf_grid.lower_bounds[ESMF.StaggerLoc.CENTER][1] - self.x_upper_bound = self.esmf_grid.upper_bounds[ESMF.StaggerLoc.CENTER][1] - self.y_lower_bound = self.esmf_grid.lower_bounds[ESMF.StaggerLoc.CENTER][0] - self.y_upper_bound = self.esmf_grid.upper_bounds[ESMF.StaggerLoc.CENTER][0] + self.calculate_bounds() self.shape = (timestep_num, len(self.vertical_description), self.x_upper_bound-self.x_lower_bound, self.y_upper_bound-self.y_lower_bound) diff --git a/hermesv3_gr/modules/grids/grid_rotated.py b/hermesv3_gr/modules/grids/grid_rotated.py index 24c0f9a83ffc580e1dd34ae4deb8050dcade7d16..8c9903612a1ef51efbff669a46d1f9edb0c159b4 100755 --- a/hermesv3_gr/modules/grids/grid_rotated.py +++ b/hermesv3_gr/modules/grids/grid_rotated.py @@ -40,14 +40,12 @@ class RotatedGrid(Grid): """ def __init__(self, grid_type, vertical_description_path, timestep_num, temporal_path, centre_lat, centre_lon, - west_boundary, south_boundary, inc_rlat, inc_rlon): - import ESMF - + west_boundary, south_boundary, inc_rlat, inc_rlon, comm=None): st_time = timeit.default_timer() settings.write_log('\tCreating Rotated grid.', level=2) # Initialises with parent class - super(RotatedGrid, self).__init__(grid_type, vertical_description_path, temporal_path) + super(RotatedGrid, self).__init__(grid_type, vertical_description_path, temporal_path, comm) # Setting parameters self.new_pole_longitude_degrees = -180 + centre_lon @@ -71,18 +69,10 @@ class RotatedGrid(Grid): if not os.path.exists(self.coords_netcdf_file): if settings.rank == 0: - # super(RotatedGrid, self).write_coords_netcdf() self.write_coords_netcdf() settings.comm.Barrier() - # self.write_coords_netcdf() - - self.esmf_grid = super(RotatedGrid, self).create_esmf_grid_from_file(self.coords_netcdf_file, sphere=False) - - self.x_lower_bound = self.esmf_grid.lower_bounds[ESMF.StaggerLoc.CENTER][1] - self.x_upper_bound = self.esmf_grid.upper_bounds[ESMF.StaggerLoc.CENTER][1] - self.y_lower_bound = self.esmf_grid.lower_bounds[ESMF.StaggerLoc.CENTER][0] - self.y_upper_bound = self.esmf_grid.upper_bounds[ESMF.StaggerLoc.CENTER][0] + self.calculate_bounds() self.shape = (timestep_num, len(self.vertical_description), self.x_upper_bound-self.x_lower_bound, self.y_upper_bound-self.y_lower_bound) diff --git a/hermesv3_gr/modules/regrid/regrid_conservative.py b/hermesv3_gr/modules/regrid/regrid_conservative.py index 35c718f0cad174da3de337dbc07b0f0541f40ef9..74c219e8ecbb9ec7a7f883e07ba1d69e2b25af55 100755 --- a/hermesv3_gr/modules/regrid/regrid_conservative.py +++ b/hermesv3_gr/modules/regrid/regrid_conservative.py @@ -48,7 +48,7 @@ class ConservativeRegrid(Regrid): src_field.read(filename=self.pollutant_dicts[0]['path'], variable=self.pollutant_dicts[0]['name'], timeslice=0) - dst_grid = self.grid.esmf_grid + dst_grid = self.grid.create_esmf_grid_from_file(self.grid.coords_netcdf_file) dst_field = ESMF.Field(dst_grid, name='my outut field') ESMF.Regrid(src_field, dst_field, filename=self.weight_matrix_file, regrid_method=ESMF.RegridMethod.CONSERVE,) diff --git a/hermesv3_gr/modules/speciation/speciation.py b/hermesv3_gr/modules/speciation/speciation.py index 9d0c1b8d63363971a46c41eb8c9807d5b6fc3e71..713a4848bda6b6ca9ef4d7d5f4bf5aaf4cd0ffbe 100755 --- a/hermesv3_gr/modules/speciation/speciation.py +++ b/hermesv3_gr/modules/speciation/speciation.py @@ -146,7 +146,7 @@ class Speciation(object): raise KeyError('ERROR: {0} pollutant is not in the molecular weights file {1} .'.format( emission['name'], self.molecular_weights_path)) sys.exit(1) - exec ("{0} = np.array(emission['data'], dtype=settings.precision)".format(emission['name'])) + exec("{0} = np.array(emission['data'], dtype=settings.precision)".format(emission['name'])) emission['units'] = '' input_pollutants.append(emission['name']) @@ -182,7 +182,7 @@ class Speciation(object): raise AttributeError( "Error in speciation profile {0}: ".format(self.id) + "The output specie {0} cannot be calculated ".format(pollutant['name']) + - "with the expression {0} because{1}".format(formula, e.message)) + "with the expression {0} because{1}".format(formula, str(e))) else: sys.exit(1) speciated_emissions.append(dict_aux) diff --git a/hermesv3_gr/modules/temporal/temporal.py b/hermesv3_gr/modules/temporal/temporal.py index 457d0c87e06adb06bf15fee4963886dad2b9883f..6f2c457973921cff91ec458d2175cb17c5c8173f 100755 --- a/hermesv3_gr/modules/temporal/temporal.py +++ b/hermesv3_gr/modules/temporal/temporal.py @@ -70,7 +70,7 @@ class TemporalDistribution(object): """ def __init__(self, starting_date, timestep_type, timestep_num, timestep_freq, monthly_profile_path, monthly_profile_id, weekly_profile_path, weekly_profile_id, daily_profile_path, daily_profile_id, - hourly_profile_path, hourly_profile_id, world_info_path, auxiliar_files_dir, grid): + hourly_profile_path, hourly_profile_id, world_info_path, auxiliar_files_dir, grid, first_time=False): from timezonefinder import TimezoneFinder import pandas as pd @@ -80,40 +80,34 @@ class TemporalDistribution(object): self.grid = grid - # self.starting_date = starting_date + if not first_time: + self.date_array = self.caclulate_date_array(starting_date, timestep_type, timestep_num, timestep_freq) - self.date_array = self.caclulate_date_array(starting_date, timestep_type, timestep_num, timestep_freq) + if timestep_type in ['daily', 'hourly']: + self.daily_profile = self.get_daily_profile(daily_profile_id, daily_profile_path) + else: + self.daily_profile = None - # self.timestep_type = timestep_type - # self.timestep_num = timestep_num - # self.timestep_freq = timestep_freq - # - # self.ending_date = self.calculate_ending_date() - if timestep_type in ['daily', 'hourly']: - self.daily_profile = self.get_daily_profile(daily_profile_id, daily_profile_path) - else: - self.daily_profile = None + if self.daily_profile is None: + if timestep_type in ['monthly', 'daily', 'hourly']: + self.monthly_profile = self.get_monthly_profile(monthly_profile_id, monthly_profile_path) + else: + self.monthly_profile = None - if self.daily_profile is None: - if timestep_type in ['monthly', 'daily', 'hourly']: - self.monthly_profile = self.get_monthly_profile(monthly_profile_id, monthly_profile_path) + if timestep_type in ['daily', 'hourly']: + self.weekly_profile = self.get_weekly_profile(weekly_profile_id, weekly_profile_path) + else: + self.weekly_profile = None else: self.monthly_profile = None - - if timestep_type in ['daily', 'hourly']: - self.weekly_profile = self.get_weekly_profile(weekly_profile_id, weekly_profile_path) - else: self.weekly_profile = None - else: - self.monthly_profile = None - self.weekly_profile = None - if monthly_profile_id is not None or weekly_profile_id is not None: - warn('WARNING: Daily profile is set. Monthly or Weekly profiles will be ignored.') + if monthly_profile_id is not None or weekly_profile_id is not None: + warn('WARNING: Daily profile is set. Monthly or Weekly profiles will be ignored.') - if timestep_type in ['hourly']: - self.hourly_profile = self.get_hourly_profile(hourly_profile_id, hourly_profile_path) - else: - self.hourly_profile = None + if timestep_type in ['hourly']: + self.hourly_profile = self.get_hourly_profile(hourly_profile_id, hourly_profile_path) + else: + self.hourly_profile = None self.world_info = world_info_path self.netcdf_timezones = os.path.join(auxiliar_files_dir, 'timezones.nc') @@ -688,7 +682,7 @@ class TemporalDistribution(object): nc_in = Dataset(self.netcdf_timezones) timezones = nc_in.variables['timezone_id'][:, self.grid.x_lower_bound:self.grid.x_upper_bound, - self.grid.y_lower_bound:self.grid.y_upper_bound].astype(int) + self.grid.y_lower_bound:self.grid.y_upper_bound].astype(int) nc_in.close() tz_list = np.chararray(timezones.shape, itemsize=32) diff --git a/hermesv3_gr/modules/writing/writer.py b/hermesv3_gr/modules/writing/writer.py index 8cb46f5e95356b65340544f079ce4e82eace43ee..535d50e7bc1fe2fed76c84b51619e0264fc6ccbe 100755 --- a/hermesv3_gr/modules/writing/writer.py +++ b/hermesv3_gr/modules/writing/writer.py @@ -25,6 +25,8 @@ from mpi4py import MPI from netCDF4 import Dataset from hermesv3_gr.config import settings from functools import reduce +import pandas as pd +from pandas import DataFrame, MultiIndex class Writer(object): @@ -186,7 +188,7 @@ class Writer(object): return True - def calculate_data_by_var(self, variable, inventory_list, shape): + def calculate_data_by_var(self, variable, inventory_list, shape, change_units=True): """ Calculate the date of the given variable throw the inventory list. @@ -243,7 +245,8 @@ class Writer(object): settings.write_time('TemporalDistribution', 'calculate_data_by_var', timeit.default_timer() - temporal_time, level=2) # Unit changes - data = self.unit_change(variable, data) + if change_units: + data = self.unit_change(variable, data) if data is not None: data[data < 0] = 0 settings.write_time('Writer', 'calculate_data_by_var', timeit.default_timer() - st_time, level=3) @@ -253,7 +256,7 @@ class Writer(object): """ Implement on inner class """ - return np.array([0]) + return data @staticmethod def calculate_displacements(counts): @@ -603,3 +606,39 @@ class Writer(object): netcdf.setncatts(global_attributes) netcdf.close() + + def get_emis(self, emission_list): + self.set_variable_attributes(emission_list) + self.change_variable_attributes() + + data_list = [] + + fid_array = np.array(range(self.grid.full_shape[-2] * self.grid.full_shape[-1])) + fid_array = fid_array.reshape((self.grid.full_shape[-2], self.grid.full_shape[-1])) + fid_array = fid_array[self.grid.x_lower_bound:self.grid.x_upper_bound, + self.grid.y_lower_bound:self.grid.y_upper_bound].flatten() + + zero_vars = [] + for var_name in self.variables_attributes.keys(): + var_data = self.calculate_data_by_var(var_name, emission_list, self.grid.shape, change_units=False) + if var_data is None: + zero_vars.append(var_name) + else: + var_data_list = [] + for tstep in range(var_data.shape[0]): + for layer in range(var_data.shape[1]): + aux_data = DataFrame(data=var_data[tstep, layer, :, :].flatten(), columns=[var_name], + index=MultiIndex.from_product([fid_array, [layer], [tstep]], + names=['FID', 'layer', 'tstep'])) + aux_data = aux_data.replace(0, np.nan) + aux_data.dropna(how='all', inplace=True) + var_data_list.append(aux_data) + data_list.append(pd.concat(var_data_list, sort=False)) + + data = pd.concat(data_list, axis=1, sort=False) + data = data.replace(np.nan, 0) + + for var_name in zero_vars: + data[var_name] = 0 + + return data diff --git a/hermesv3_gr/tools/netcdf_tools.py b/hermesv3_gr/tools/netcdf_tools.py index 48869a1d4c86b80da769416b995a9627038cc37f..8678e9f8c130fd50c8ed4fe91f96e7aae9b57655 100755 --- a/hermesv3_gr/tools/netcdf_tools.py +++ b/hermesv3_gr/tools/netcdf_tools.py @@ -23,11 +23,6 @@ from netCDF4 import Dataset from mpi4py import MPI from functools import reduce -ICOMM = MPI.COMM_WORLD -COMM = ICOMM.Split(color=0, key=0) -RANK = COMM.Get_rank() -SIZE = COMM.Get_size() - def open_netcdf(netcdf_path): """ @@ -606,9 +601,8 @@ def create_netcdf(netcdf_path, center_latitudes, center_longitudes, data_list, # if variable['data'] is not 0: # print var[:].shape, variable['data'].shape, variable['data'].max() shape = tuple() - exec ("shape = (len(hours), {0}.size, {1}.size, {2}.size)".format(var_dim[0], var_dim[1], var_dim[2])) - # exit() - print(shape) + exec("shape = (len(hours), {0}.size, {1}.size, {2}.size)".format(var_dim[0], var_dim[1], var_dim[2])) + var[:] = np.zeros(shape) # Grid mapping diff --git a/run_test.py b/run_test.py index 6eb2926793f63a9c62eee79c06273d12f95fab82..655465f0f26d182068b5ffe928597a593a781365 100755 --- a/run_test.py +++ b/run_test.py @@ -1,5 +1,4 @@ # coding=utf-8 -"""Script to run the tests for EarthDiagnostics and generate the code coverage report""" import os import sys