From 797297a0bed8f74643c32a8de67cfe9afa330b86 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Wed, 15 Jan 2020 18:47:32 +0100 Subject: [PATCH 1/3] Adding Rotated nested --- CHANGELOG | 4 + conf/hermes.conf | 16 +- hermesv3_gr/config/config.py | 17 ++ hermesv3_gr/hermes.py | 4 +- hermesv3_gr/modules/grids/grid.py | 174 ++++-------- hermesv3_gr/modules/grids/grid_mercator.py | 4 +- .../modules/grids/grid_rotated_nested.py | 251 ++++++++++++++++++ 7 files changed, 341 insertions(+), 129 deletions(-) create mode 100755 hermesv3_gr/modules/grids/grid_rotated_nested.py diff --git a/CHANGELOG b/CHANGELOG index b1ff499..4f80e7f 100755 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,3 +1,7 @@ +2.0.3 + XXXX/XX/XX + - Rotated nested grid + 2.0.2 2020/01/14 - Corrected bug on GFAS as point emissions. diff --git a/conf/hermes.conf b/conf/hermes.conf index 7a12597..3d57fc5 100755 --- a/conf/hermes.conf +++ b/conf/hermes.conf @@ -22,12 +22,8 @@ output_model = MONARCH # output_model = WRF_CHEM output_attributes = /data/global_attributes.csv -# ***** domain_type=[global, regular, lcc, rotated, mercator] ***** -domain_type = global -# domain_type = regular -# domain_type = lcc -# domain_type = rotated -# domain_type = mercator +# domain_type=[lcc, rotated, mercator, regular, rotated_nested] +domain_type = rotated_nested vertical_description = /data/profiles/vertical/Benchmark_15layers_vertical_description.csv auxiliary_files_path = /data/auxiliar_files/_ @@ -52,6 +48,14 @@ auxiliary_files_path = /data/auxiliar_files/_', '{1}_{2}'.format( item, arguments.inc_rlat, arguments.inc_rlon)) + elif arguments.domain_type == 'rotated_nested': + arguments.__dict__[item] = arguments.__dict__[item].replace('', '{1}_{2}'.format( + item, arguments.n_rlat, arguments.n_rlon)) elif arguments.domain_type == 'lcc' or arguments.domain_type == 'mercator': arguments.__dict__[item] = arguments.__dict__[item].replace('', '{1}_{2}'.format( item, arguments.inc_x, arguments.inc_y)) diff --git a/hermesv3_gr/hermes.py b/hermesv3_gr/hermes.py index 1df2970..a7deb83 100755 --- a/hermesv3_gr/hermes.py +++ b/hermesv3_gr/hermes.py @@ -34,7 +34,7 @@ class HermesGr(object): Interface class for HERMESv3_GR. """ def __init__(self, config, new_date=None, comm=None): - from hermesv3_gr.modules.grids.grid import Grid + from hermesv3_gr.modules.grids.grid import select_grid from hermesv3_gr.modules.temporal.temporal import TemporalDistribution from hermesv3_gr.modules.writing.writer import Writer global full_time @@ -65,7 +65,7 @@ class HermesGr(object): self.levels = VerticalDistribution.get_vertical_output_profile(self.options.vertical_description) - self.grid = Grid.select_grid( + self.grid = select_grid( self.comm, self.options.domain_type, self.options.vertical_description, self.options.output_timestep_num, self.options.auxiliary_files_path, self.options.inc_lat, self.options.inc_lon, self.options.lat_orig, diff --git a/hermesv3_gr/modules/grids/grid.py b/hermesv3_gr/modules/grids/grid.py index 7fe4bd5..5e057c5 100755 --- a/hermesv3_gr/modules/grids/grid.py +++ b/hermesv3_gr/modules/grids/grid.py @@ -26,6 +26,61 @@ import ESMF import hermesv3_gr.config.settings as settings +def select_grid(comm, arguments): + + st_time = timeit.default_timer() + settings.write_log('Selecting grid', level=1) + + # Creating a different object depending on the grid type + if arguments.grid_type == 'global': + from hermesv3_gr.modules.grids.grid_global import GlobalGrid + grid = GlobalGrid(arguments.grid_type, arguments.vertical_description, arguments.output_timestep_num, + arguments.auxiliary_files_path, arguments.inc_lat, arguments.inc_lon, comm=comm) + + elif arguments.domain_type == 'regular': + from hermesv3_gr.modules.grids.grid_latlon import LatLonGrid + grid = LatLonGrid(arguments.grid_type, arguments.vertical_description, arguments.output_timestep_num, + arguments.auxiliary_files_path, arguments.inc_lat, arguments.inc_lon, + arguments.lat_orig, arguments.lon_orig, arguments.n_lat, arguments.n_lon, comm=comm) + + elif arguments.domain_type == 'rotated': + from hermesv3_gr.modules.grids.grid_rotated import RotatedGrid + grid = RotatedGrid(arguments.grid_type, arguments.vertical_description, arguments.output_timestep_num, + arguments.auxiliary_files_path, arguments.centre_lat, arguments.centre_lon, + arguments.west_boundary, arguments.south_boundary, arguments.inc_rlat, arguments.inc_rlon, + comm=comm) + + elif arguments.domain_type == 'rotated_nested': + from hermesv3_gr.modules.grids.grid_rotated_nested import RotatedNestedGrid + grid = RotatedNestedGrid(arguments.grid_type, arguments.vertical_description, arguments.output_timestep_num, + arguments.auxiliary_files_path, arguments.parent_grid_path, arguments.parent_ratio, + arguments.i_parent_start, arguments.j_parent_start, arguments.n_rlat, arguments.n_rlon, + comm=comm) + + elif arguments.domain_type == 'lcc': + from hermesv3_gr.modules.grids.grid_lcc import LccGrid + grid = LccGrid(arguments.grid_type, arguments.vertical_description, arguments.output_timestep_num, + arguments.auxiliary_files_path, arguments.lat_1, arguments.lat_2, arguments.lon_0, + arguments.lat_0, arguments.nx, arguments.ny, arguments.inc_x, arguments.inc_y, arguments.x_0, + arguments.y_0, comm=comm) + + elif arguments.domain_type == 'mercator': + from hermesv3_gr.modules.grids.grid_mercator import MercatorGrid + grid = MercatorGrid(arguments.grid_type, arguments.vertical_description, arguments.output_timestep_num, + arguments.auxiliary_files_path, arguments.lat_ts, arguments.lon_0, arguments.nx, + arguments.ny, arguments.inc_x, arguments.inc_y, arguments.x_0, arguments.y_0, comm=comm) + else: + settings.write_log('ERROR: Check the .err file to get more info.') + if settings.rank == 0: + raise NotImplementedError("The grid type {0} is not implemented.".format(arguments.domain_type) + + " Use 'global', 'regular, 'rotated', 'lcc' or 'mercator.") + sys.exit(1) + + settings.write_time('Grid', 'select_grid', timeit.default_timer() - st_time, level=3) + + return grid + + class Grid(object): """ Grid object that contains the information of the output grid. @@ -110,125 +165,6 @@ class Grid(object): settings.write_time('Grid', 'create_esmf_grid_from_file', timeit.default_timer() - st_time, level=3) return grid - @staticmethod - 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 - """ - Create a Grid object depending on the grid type. - - :param grid_type: type of grid to create [global, rotated, lcc, mercator] - :type grid_type: str - - :param vertical_description_path: Path to the file that contains the vertical description. - :type vertical_description_path: str - - :param timestep_num: Number of timesteps. - :type timestep_num: int - - :param temporal_path: Path to the temporal folder. - :type temporal_path: str - - :param inc_lat: [global] Increment between latitude centroids (degrees). - :type inc_lat: float - - :param inc_lon: [global] Increment between longitude centroids (degrees). - :type inc_lon: float - - :param centre_lat: [rotated] - :type centre_lat: float - - :param centre_lon: [rotated] - :type centre_lon: float - - :param west_boundary: [rotated] - :type west_boundary: float - - :param south_boundary: [rotated] - :type south_boundary: float - - :param inc_rlat: [rotated] Increment between rotated latitude centroids (degrees). - :type inc_rlat: float - - :param inc_rlon: [rotated] Increment between rotated longitude centroids (degrees). - :type inc_rlon: float - - :param lat_ts: [mercator] - :type lat_ts: float - - :param lat_1: [lcc] Value of the Lat1 for the LCC grid type. - :type lat_1: float - - :param lat_2: [lcc] Value of the Lat2 for the LCC grid type. - :type lat_2: float - - :param lon_0: [lcc, mercator] Value of the Lon0 for the LCC grid type. - :type lon_0: float - - :param lat_0: [lcc] Value of the Lat0 for the LCC grid type. - :type lat_0: float - - :param nx: [lcc, mercator] Number of cells on the x dimension. - :type nx: int - - :param ny: [lcc, mercator] Number of cells on the y dimension. - :type ny: int - - :param inc_x: [lcc, mercator] Increment between x dimensions cell centroids (metres). - :type inc_x: int - - :param inc_y: [lcc, mercator] Increment between y dimensions cell centroids (metres). - :type inc_y: int - - :param x_0: [lcc, mercator] Value of the X0 for the LCC grid type. - :type x_0: float - - :param y_0: [lcc, mercator] Value of the Y0 for the LCC grid type. - :type y_0: float - - :return: Grid object. It will return a GlobalGrid, RotatedGrid or LccGrid depending on the type. - :rtype: Grid - """ - - st_time = timeit.default_timer() - settings.write_log('Selecting grid', level=1) - - # 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, - 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, 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, 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, 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, comm=comm) - else: - settings.write_log('ERROR: Check the .err file to get more info.') - if settings.rank == 0: - raise NotImplementedError("The grid type {0} is not implemented.".format(grid_type) - + " Use 'global', 'regular, 'rotated', 'lcc' or 'mercator.") - sys.exit(1) - - settings.write_time('Grid', 'select_grid', timeit.default_timer() - st_time, level=3) - - return grid - @staticmethod def set_vertical_levels(vertical_description_path): """ diff --git a/hermesv3_gr/modules/grids/grid_mercator.py b/hermesv3_gr/modules/grids/grid_mercator.py index e956d17..1714439 100755 --- a/hermesv3_gr/modules/grids/grid_mercator.py +++ b/hermesv3_gr/modules/grids/grid_mercator.py @@ -68,12 +68,12 @@ 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): + 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 Mercator grid.', level=2) # Initialises with parent class - super(MercatorGrid, self).__init__(grid_type, vertical_description_path, temporal_path) + super(MercatorGrid, self).__init__(grid_type, vertical_description_path, temporal_path, comm) # Setting parameters self.lat_ts = lat_ts diff --git a/hermesv3_gr/modules/grids/grid_rotated_nested.py b/hermesv3_gr/modules/grids/grid_rotated_nested.py new file mode 100755 index 0000000..60a9d2e --- /dev/null +++ b/hermesv3_gr/modules/grids/grid_rotated_nested.py @@ -0,0 +1,251 @@ +#!/usr/bin/env python + +# Copyright 2018 Earth Sciences Department, BSC-CNS +# +# This file is part of HERMESv3_GR. +# +# HERMESv3_GR is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# HERMESv3_GR is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with HERMESv3_GR. If not, see . + + +import sys +import os +import timeit +import hermesv3_gr.config.settings as settings +from .grid import Grid + + +class RotatedNestedGrid(Grid): + # TODO Rotated options description + """ + :param grid_type: Type of the output grid [global, rotated, lcc, mercator]. + :type grid_type: str + + :param vertical_description_path: Path to the file that contains the vertical description. + :type vertical_description_path: str + + + :param timestep_num: Number of timesteps. + :type timestep_num: int + """ + + def __init__(self, grid_type, vertical_description_path, timestep_num, temporal_path, parent_grid_path, + parent_ratio, i_parent_start, j_parent_start, n_rlat, n_rlon, comm=None): + st_time = timeit.default_timer() + settings.write_log('\tCreating Rotated grid.', level=2) + + # Initialises with parent class + super(RotatedNestedGrid, self).__init__(grid_type, vertical_description_path, temporal_path, comm) + + # Setting parameters + attributes = {'parent_grid_path': parent_grid_path, 'parent_ratio': parent_ratio, + 'i_parent_start': i_parent_start, 'j_parent_start': j_parent_start, + 'n_rlat': n_rlat, 'n_rlon': n_rlon, 'crs': {'init': 'epsg:4326'}} + attributes = self.get_parent_attributes(attributes) + + self.new_pole_longitude_degrees = self.attributes['parent']['new_pole_longitude_degrees'] + self.new_pole_latitude_degrees = attributes['parent']['new_pole_latitude_degrees'] + self.rlat_1st = attributes['1st_rlat'] + self.rlon_1st = attributes['1st_rlon'] + self.inc_rlat = attributes['inc_rlat'] + self.inc_rlon = attributes['inc_rlon'] + self.n_rlat = attributes['n_rlat'] + self.n_rlon = attributes['n_rlon'] + + # Rotated coordinates + self.rlat = None + self.rlon = None + + # Create coordinates + self.crs = {'init': 'epsg:4326'} + self.create_coords() + + if not os.path.exists(self.coords_netcdf_file): + if settings.rank == 0: + self.write_coords_netcdf() + settings.comm.Barrier() + + 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) + + total_area = self.get_cell_area() + self.full_shape = (timestep_num, len(self.vertical_description), total_area.shape[-2], total_area.shape[-1]) + self.cell_area = total_area[self.x_lower_bound:self.x_upper_bound, self.y_lower_bound:self.y_upper_bound] + + settings.write_time('RotatedNestedGrid', 'Init', timeit.default_timer() - st_time, level=1) + + @staticmethod + def get_parent_attributes(attributes): + from netCDF4 import Dataset + + netcdf = Dataset(attributes['parent_grid_path'], mode='r') + + rlat = netcdf.variables['rlat'][:] + attributes['inc_rlat'] = (rlat[1] - rlat[0]) / attributes['parent_ratio'] + attributes['1st_rlat'] = rlat[int(attributes['j_parent_start'])] + + rlon = netcdf.variables['rlon'][:] + attributes['inc_rlon'] = (rlon[1] - rlon[0]) / attributes['parent_ratio'] + attributes['1st_rlon'] = rlon[attributes['i_parent_start']] + + rotated_pole = netcdf.variables['rotated_pole'] + attributes['parent'] = {'new_pole_longitude_degrees': rotated_pole.grid_north_pole_longitude, + 'new_pole_latitude_degrees': 90 - rotated_pole.grid_north_pole_latitude} + + return attributes + + def create_coords(self): + """ + Create the coordinates for a rotated domain. + """ + from hermesv3_gr.tools.coordinates_tools import create_regular_rotated + import numpy as np + + st_time = timeit.default_timer() + settings.write_log('\t\tCreating rotated coordinates.', level=3) + + # Create rotated coordinates + (self.rlat, self.rlon, br_lats_single, br_lons_single) = create_regular_rotated( + self.rlat_1st, self.rlon_1st, self.inc_rlat, self.inc_rlon, self.n_rlat, self.n_rlon) + if len(self.rlon)//2 < settings.size: + settings.write_log('ERROR: Check the .err file to get more info.') + if settings.rank == 0: + raise AttributeError("ERROR: Maximum number of processors exceeded. " + + "It has to be less or equal than {0}.".format(len(self.rlon)//2)) + sys.exit(1) + # 1D to 2D + c_lats = np.array([self.rlat] * len(self.rlon)).T + c_lons = np.array([self.rlon] * len(self.rlat)) + + # Create rotated boundary coordinates + b_lats = super(RotatedNestedGrid, self).create_bounds(c_lats, self.inc_rlat, number_vertices=4, inverse=True) + b_lons = super(RotatedNestedGrid, self).create_bounds(c_lons, self.inc_rlon, number_vertices=4) + + # Rotated to Lat-Lon + self.boundary_longitudes, self.boundary_latitudes = self.rotated2latlon(b_lons, b_lats) + self.center_longitudes, self.center_latitudes = self.rotated2latlon(c_lons, c_lats) + + settings.write_time('RotatedNestedGrid', 'create_coords', timeit.default_timer() - st_time, level=2) + + def rotated2latlon(self, lon_deg, lat_deg, lon_min=-180): + """ + Calculate the unrotated coordinates using the rotated ones. + + :param lon_deg: Rotated longitude coordinate. + :type lon_deg: numpy.array + + :param lat_deg: Rotated latitude coordinate. + :type lat_deg: numpy.array + + :param lon_min: Minimum value for the longitudes: -180 (-180 to 180) or 0 (0 to 360) + :type lon_min: float + + :return: Unrotated coordinates. Longitudes, Latitudes + :rtype: tuple(numpy.array, numpy.array) + """ + import numpy as np + import math + + st_time = timeit.default_timer() + settings.write_log('\t\t\tTransforming rotated coordinates to latitude, longitude coordinates.', level=3) + + # TODO Document this function + degrees_to_radians = math.pi / 180. + # radians_to_degrees = 180. / math.pi + + # Positive east to negative east + # self.new_pole_longitude_degrees -= 180 + + tph0 = self.new_pole_latitude_degrees * degrees_to_radians + tlm = lon_deg * degrees_to_radians + tph = lat_deg * degrees_to_radians + tlm0d = self.new_pole_longitude_degrees + ctph0 = np.cos(tph0) + stph0 = np.sin(tph0) + + stlm = np.sin(tlm) + ctlm = np.cos(tlm) + stph = np.sin(tph) + ctph = np.cos(tph) + + # Latitude + sph = (ctph0 * stph) + (stph0 * ctph * ctlm) + # if sph > 1.: + # sph = 1. + # if sph < -1.: + # sph = -1. + # print type(sph) + sph[sph > 1.] = 1. + sph[sph < -1.] = -1. + + aph = np.arcsin(sph) + aphd = aph / degrees_to_radians + + # Longitude + anum = ctph * stlm + denom = (ctlm * ctph - stph0 * sph) / ctph0 + relm = np.arctan2(anum, denom) - math.pi + almd = relm / degrees_to_radians + tlm0d + + # if almd < min_lon: + # almd += 360 + # elif almd > max_lon: + # almd -= 360 + almd[almd > (lon_min + 360)] -= 360 + almd[almd < lon_min] += 360 + + settings.write_time('RotatedNestedGrid', 'rotated2latlon', timeit.default_timer() - st_time, level=3) + + return almd, aphd + + def write_coords_netcdf(self): + """ + Writes the temporal file with the coordinates of the output needed to generate the weight matrix. + If it is already well created it will only add the cell_area parameter. + """ + from hermesv3_gr.modules.writing.writer import Writer + + st_time = timeit.default_timer() + settings.write_log('\tWriting {0} file.'.format(self.coords_netcdf_file), level=3) + + if not self.chech_coords_file(): + # Writes an auxiliary empty NetCDF only with the coordinates and an empty variable. + Writer.write_netcdf(self.coords_netcdf_file, self.center_latitudes, self.center_longitudes, + [{'name': 'var_aux', 'units': '', 'data': 0}], + boundary_latitudes=self.boundary_latitudes, + boundary_longitudes=self.boundary_longitudes, + roated=True, rotated_lats=self.rlat, rotated_lons=self.rlon, + north_pole_lat=self.new_pole_latitude_degrees, + north_pole_lon=self.new_pole_longitude_degrees) + + # Calculate the cell area of the auxiliary NetCDF file + self.cell_area = self.get_cell_area() + + # Re-writes the NetCDF adding the cell area + Writer.write_netcdf(self.coords_netcdf_file, self.center_latitudes, self.center_longitudes, + [{'name': 'var_aux', 'units': '', 'data': 0}], + boundary_latitudes=self.boundary_latitudes, + boundary_longitudes=self.boundary_longitudes, cell_area=self.cell_area, + roated=True, rotated_lats=self.rlat, rotated_lons=self.rlon, + north_pole_lat=self.new_pole_latitude_degrees, + north_pole_lon=self.new_pole_longitude_degrees) + else: + self.cell_area = self.get_cell_area() + + settings.write_time('RotatedNestedGrid', 'write_coords_netcdf', timeit.default_timer() - st_time, level=3) + + +if __name__ == '__main__': + pass -- GitLab From 31e75db871e4bb12d0a08c2075a6f43d06144d49 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Thu, 13 Feb 2020 10:46:01 +0100 Subject: [PATCH 2/3] Bug solving --- hermesv3_gr/config/config.py | 2 +- hermesv3_gr/hermes.py | 10 +--------- hermesv3_gr/modules/grids/grid.py | 12 ++++++------ hermesv3_gr/modules/grids/grid_rotated_nested.py | 8 +++++--- hermesv3_gr/modules/writing/writer.py | 2 +- 5 files changed, 14 insertions(+), 20 deletions(-) diff --git a/hermesv3_gr/config/config.py b/hermesv3_gr/config/config.py index 9cb8f9d..c5adae9 100755 --- a/hermesv3_gr/config/config.py +++ b/hermesv3_gr/config/config.py @@ -79,7 +79,7 @@ class Config(ArgParser): help='Path to the file that contains the global attributes.') p.add_argument('--domain_type', required=True, help='Type of domain to simulate.', - choices=['global', 'lcc', 'rotated', 'mercator', 'regular']) + choices=['global', 'lcc', 'rotated', 'mercator', 'regular', 'rotated_nested']) 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.') diff --git a/hermesv3_gr/hermes.py b/hermesv3_gr/hermes.py index a7deb83..5a3344b 100755 --- a/hermesv3_gr/hermes.py +++ b/hermesv3_gr/hermes.py @@ -65,15 +65,7 @@ class HermesGr(object): self.levels = VerticalDistribution.get_vertical_output_profile(self.options.vertical_description) - self.grid = select_grid( - self.comm, - self.options.domain_type, self.options.vertical_description, self.options.output_timestep_num, - 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, - self.options.nx, self.options.ny, self.options.inc_x, self.options.inc_y, self.options.x_0, - self.options.y_0, self.options.lat_ts) + self.grid = select_grid(self.comm, self.options) self.emission_list = EmissionInventory.make_emission_list(self.options, self.grid, self.levels, self.options.start_date) diff --git a/hermesv3_gr/modules/grids/grid.py b/hermesv3_gr/modules/grids/grid.py index 5e057c5..794ad96 100755 --- a/hermesv3_gr/modules/grids/grid.py +++ b/hermesv3_gr/modules/grids/grid.py @@ -32,41 +32,41 @@ def select_grid(comm, arguments): settings.write_log('Selecting grid', level=1) # Creating a different object depending on the grid type - if arguments.grid_type == 'global': + if arguments.domain_type == 'global': from hermesv3_gr.modules.grids.grid_global import GlobalGrid grid = GlobalGrid(arguments.grid_type, arguments.vertical_description, arguments.output_timestep_num, arguments.auxiliary_files_path, arguments.inc_lat, arguments.inc_lon, comm=comm) elif arguments.domain_type == 'regular': from hermesv3_gr.modules.grids.grid_latlon import LatLonGrid - grid = LatLonGrid(arguments.grid_type, arguments.vertical_description, arguments.output_timestep_num, + grid = LatLonGrid(arguments.domain_type, arguments.vertical_description, arguments.output_timestep_num, arguments.auxiliary_files_path, arguments.inc_lat, arguments.inc_lon, arguments.lat_orig, arguments.lon_orig, arguments.n_lat, arguments.n_lon, comm=comm) elif arguments.domain_type == 'rotated': from hermesv3_gr.modules.grids.grid_rotated import RotatedGrid - grid = RotatedGrid(arguments.grid_type, arguments.vertical_description, arguments.output_timestep_num, + grid = RotatedGrid(arguments.domain_type, arguments.vertical_description, arguments.output_timestep_num, arguments.auxiliary_files_path, arguments.centre_lat, arguments.centre_lon, arguments.west_boundary, arguments.south_boundary, arguments.inc_rlat, arguments.inc_rlon, comm=comm) elif arguments.domain_type == 'rotated_nested': from hermesv3_gr.modules.grids.grid_rotated_nested import RotatedNestedGrid - grid = RotatedNestedGrid(arguments.grid_type, arguments.vertical_description, arguments.output_timestep_num, + grid = RotatedNestedGrid(arguments.domain_type, arguments.vertical_description, arguments.output_timestep_num, arguments.auxiliary_files_path, arguments.parent_grid_path, arguments.parent_ratio, arguments.i_parent_start, arguments.j_parent_start, arguments.n_rlat, arguments.n_rlon, comm=comm) elif arguments.domain_type == 'lcc': from hermesv3_gr.modules.grids.grid_lcc import LccGrid - grid = LccGrid(arguments.grid_type, arguments.vertical_description, arguments.output_timestep_num, + grid = LccGrid(arguments.domain_type, arguments.vertical_description, arguments.output_timestep_num, arguments.auxiliary_files_path, arguments.lat_1, arguments.lat_2, arguments.lon_0, arguments.lat_0, arguments.nx, arguments.ny, arguments.inc_x, arguments.inc_y, arguments.x_0, arguments.y_0, comm=comm) elif arguments.domain_type == 'mercator': from hermesv3_gr.modules.grids.grid_mercator import MercatorGrid - grid = MercatorGrid(arguments.grid_type, arguments.vertical_description, arguments.output_timestep_num, + grid = MercatorGrid(arguments.domain_type, arguments.vertical_description, arguments.output_timestep_num, arguments.auxiliary_files_path, arguments.lat_ts, arguments.lon_0, arguments.nx, arguments.ny, arguments.inc_x, arguments.inc_y, arguments.x_0, arguments.y_0, comm=comm) else: diff --git a/hermesv3_gr/modules/grids/grid_rotated_nested.py b/hermesv3_gr/modules/grids/grid_rotated_nested.py index 60a9d2e..1ae474f 100755 --- a/hermesv3_gr/modules/grids/grid_rotated_nested.py +++ b/hermesv3_gr/modules/grids/grid_rotated_nested.py @@ -53,7 +53,7 @@ class RotatedNestedGrid(Grid): 'n_rlat': n_rlat, 'n_rlon': n_rlon, 'crs': {'init': 'epsg:4326'}} attributes = self.get_parent_attributes(attributes) - self.new_pole_longitude_degrees = self.attributes['parent']['new_pole_longitude_degrees'] + self.new_pole_longitude_degrees = attributes['parent']['new_pole_longitude_degrees'] self.new_pole_latitude_degrees = attributes['parent']['new_pole_latitude_degrees'] self.rlat_1st = attributes['1st_rlat'] self.rlon_1st = attributes['1st_rlon'] @@ -94,11 +94,13 @@ class RotatedNestedGrid(Grid): rlat = netcdf.variables['rlat'][:] attributes['inc_rlat'] = (rlat[1] - rlat[0]) / attributes['parent_ratio'] - attributes['1st_rlat'] = rlat[int(attributes['j_parent_start'])] + # j_parent_start start the index at 1 so we must - 1 + attributes['1st_rlat'] = rlat[int(attributes['j_parent_start']) - 1] rlon = netcdf.variables['rlon'][:] attributes['inc_rlon'] = (rlon[1] - rlon[0]) / attributes['parent_ratio'] - attributes['1st_rlon'] = rlon[attributes['i_parent_start']] + # i_parent_start start the index at 1 so we must - 1 + attributes['1st_rlon'] = rlon[attributes['i_parent_start'] - 1] rotated_pole = netcdf.variables['rotated_pole'] attributes['parent'] = {'new_pole_longitude_degrees': rotated_pole.grid_north_pole_longitude, diff --git a/hermesv3_gr/modules/writing/writer.py b/hermesv3_gr/modules/writing/writer.py index 9fb1d1a..c32517b 100755 --- a/hermesv3_gr/modules/writing/writer.py +++ b/hermesv3_gr/modules/writing/writer.py @@ -577,7 +577,7 @@ class Writer(object): # Rotated pole mapping = netcdf.createVariable('rotated_pole', 'c') mapping.grid_mapping_name = 'rotated_latitude_longitude' - mapping.grid_north_pole_latitude = north_pole_lat + mapping.grid_north_pole_latitude = 90 - north_pole_lat mapping.grid_north_pole_longitude = north_pole_lon elif lcc: # CRS -- GitLab From 237581e062f734fd91e7794681bafec2f62202c9 Mon Sep 17 00:00:00 2001 From: Carles Tena Date: Thu, 13 Feb 2020 15:10:40 +0100 Subject: [PATCH 3/3] Added rotated-nested to the monarch writer --- hermesv3_gr/modules/writing/writer_monarch.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hermesv3_gr/modules/writing/writer_monarch.py b/hermesv3_gr/modules/writing/writer_monarch.py index 1bd6758..8c5a551 100755 --- a/hermesv3_gr/modules/writing/writer_monarch.py +++ b/hermesv3_gr/modules/writing/writer_monarch.py @@ -129,7 +129,7 @@ class WriterMonarch(Writer): LambertConformalConic = False if self.grid.grid_type == 'global' or self.grid.grid_type == 'regular': RegularLatLon = True - elif self.grid.grid_type == 'rotated': + elif self.grid.grid_type in ['rotated', 'rotated_nested']: Rotated = True elif self.grid.grid_type == 'lcc': LambertConformalConic = True @@ -430,7 +430,7 @@ class WriterMonarch(Writer): if self.grid.grid_type == 'global' or self.grid.grid_type == 'regular': regular_latlon = True - elif self.grid.grid_type == 'rotated': + elif self.grid.grid_type in ['rotated', 'rotated_nested']: rotated = True elif self.grid.grid_type == 'lcc': lcc = True -- GitLab