From 7630cc8d97dc0c2b84e6a84a8bcc47581001a540 Mon Sep 17 00:00:00 2001 From: ctena Date: Mon, 17 Apr 2023 14:12:03 +0200 Subject: [PATCH 1/6] WIP: MONARCH out_format --- nes/nc_projections/default_nes.py | 30 +++++++--- nes/nes_formats/__init__.py | 1 + nes/nes_formats/monarch_format.py | 98 +++++++++++++++++++++++++++++++ 3 files changed, 121 insertions(+), 8 deletions(-) create mode 100644 nes/nes_formats/monarch_format.py diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py index 9f9d361..2523b37 100644 --- a/nes/nc_projections/default_nes.py +++ b/nes/nc_projections/default_nes.py @@ -16,7 +16,7 @@ from copy import deepcopy, copy import datetime import pyproj from ..methods import vertical_interpolation, horizontal_interpolation, cell_measures, spatial_join -from ..nes_formats import to_netcdf_cams_ra +from ..nes_formats import to_netcdf_cams_ra, to_netcdf_monarch, to_monarch_units class Nes(object): @@ -2056,7 +2056,7 @@ class Nes(object): Parameters ---------- - data_type : str + data_type : str or Type Data type, by default 'float32' """ @@ -2541,7 +2541,7 @@ class Nes(object): return None - def append_time_step_data(self, i_time): + def append_time_step_data(self, i_time, out_format='DEFAULT'): """ Fill the netCDF data for the indicated index time. @@ -2561,6 +2561,8 @@ class Nes(object): self.serial_nc.append_time_step_data(i_time) self.comm.Barrier() else: + if out_format == 'MONARCH': + self.variables = to_monarch_units(self.variables) for i, (var_name, var_dict) in enumerate(self.variables.items()): for att_name, att_value in var_dict.items(): if att_name == 'data': @@ -2641,6 +2643,8 @@ class Nes(object): Path to the output netCDF file. chunking: bool Indicates if you want to chunk the output netCDF. + keep_open : bool + Indicates if you want to keep open the NetCDH to fill the data by time-step """ # Open NetCDF @@ -2686,6 +2690,7 @@ class Nes(object): def __to_netcdf_cams_ra(self, path): return to_netcdf_cams_ra(self, path) + def to_netcdf(self, path, compression_level=0, serial=False, info=False, chunking=False, type='NES', keep_open=False): """ @@ -2705,6 +2710,8 @@ class Nes(object): Indicates if you want a chunked netCDF output. Only available with non-serial writes. Default: False. type : str Type to NetCDf to write. 'CAMS_RA' or 'NES' + keep_open : bool + Indicates if you want to keep open the NetCDH to fill the data by time-step """ nc_type = type old_info = self.info @@ -2729,23 +2736,30 @@ class Nes(object): new_nc.set_communicator(MPI.COMM_SELF) new_nc.variables = data new_nc.cell_measures = c_measures - if type == 'NES': + if type in ['NES', 'DEFAULT']: new_nc.__to_netcdf_py(path, keep_open=keep_open) elif type == 'CAMS_RA': new_nc.__to_netcdf_cams_ra(path) + elif type == 'MONARCH': + to_netcdf_monarch(new_nc, path, chunking=chunking, keep_open=keep_open) else: - raise ValueError( - "Unknown NetCDF type '{0}'. Use 'CAMS_RA' or 'NES'; default='NES'".format(nc_type)) + msg = "Unknown NetCDF type '{0}'. ".format(nc_type) + msg += "Use CAMS_RA, MONARCH or NES (or DEFAULT)" + raise ValueError(msg) self.serial_nc = new_nc else: self.serial_nc = True else: - if nc_type == 'NES': + if type in ['NES', 'DEFAULT']: self.__to_netcdf_py(path, chunking=chunking, keep_open=keep_open) elif nc_type == 'CAMS_RA': self.__to_netcdf_cams_ra(path) + elif nc_type == 'MONARCH': + to_netcdf_monarch(self, path, chunking=chunking, keep_open=keep_open) else: - raise ValueError("Unknown NetCDF type '{0}'. Use 'CAMS_RA' or 'NES'; default='NES'".format(nc_type)) + msg = "Unknown NetCDF type '{0}'. ".format(nc_type) + msg += "Use CAMS_RA, MONARCH or NES (or DEFAULT)" + raise ValueError(msg) self.info = old_info diff --git a/nes/nes_formats/__init__.py b/nes/nes_formats/__init__.py index 0b0a484..f597c54 100644 --- a/nes/nes_formats/__init__.py +++ b/nes/nes_formats/__init__.py @@ -1 +1,2 @@ from .cams_ra_format import to_netcdf_cams_ra +from .monarch_format import to_netcdf_monarch, to_monarch_units diff --git a/nes/nes_formats/monarch_format.py b/nes/nes_formats/monarch_format.py new file mode 100644 index 0000000..3a77335 --- /dev/null +++ b/nes/nes_formats/monarch_format.py @@ -0,0 +1,98 @@ +#!/usr/bin/env python + +import sys +import warnings +import numpy as np +import os +import nes +from netCDF4 import Dataset +from mpi4py import MPI +from copy import copy + + +def to_netcdf_monarch(self, path, chunking=False, keep_open=False): + """ + Create the NetCDF using netcdf4-python methods. + + Parameters + ---------- + self : nes.Nes + Source projection Nes Object. + path : str + Path to the output netCDF file. + chunking: bool + Indicates if you want to chunk the output netCDF. + keep_open : bool + Indicates if you want to keep open the NetCDH to fill the data by time-step + """ + + self.to_dtype(np.float32) + + # Open NetCDF + if self.info: + print("Rank {0:03d}: Creating {1}".format(self.rank, path)) + if self.size > 1: + netcdf = Dataset(path, format="NETCDF4", mode='w', parallel=True, comm=self.comm, info=MPI.Info()) + else: + netcdf = Dataset(path, format="NETCDF4", mode='w', parallel=False) + if self.info: + print("Rank {0:03d}: NetCDF ready to write".format(self.rank)) + + # Create dimensions + self._create_dimensions(netcdf) + + # Create dimension variables + self._create_dimension_variables(netcdf) + if self.info: + print("Rank {0:03d}: Dimensions done".format(self.rank)) + + # Create cell measures + self._create_cell_measures(netcdf) + + # Create variables + self._create_variables(netcdf, chunking=chunking) + + # Create metadata + self._create_metadata(netcdf) + + # Close NetCDF + if self.global_attrs is not None: + for att_name, att_value in self.global_attrs.items(): + netcdf.setncattr(att_name, att_value) + netcdf.setncattr('Conventions', 'CF-1.7') + + if keep_open: + self.netcdf = netcdf + else: + netcdf.close() + + return None + + +def to_monarch_units(variables): + """ + Change the data values according to the MONARCH conventions + + Parameters + ---------- + variables : dict + Variables data + + Returns + ------- + dict + Variable in the MONARCH units + """ + for var_name in variables.keys(): + if variables['units'] == 'mol.s-1.m-2': + # Kmol to mol + variables['data'] *= 1000 + elif variables['units'] == 'kg.s-1.m-2': + # No action needed + pass + else: + raise TypeError("The unit '{0}' of specie {1} is not defined correctly. ".format( + variables['units'], var_name) + + "Should be 'mol.s-1.m-2' or 'kg.s-1.m-2'") + return variables + -- GitLab From b315764712623d0e38f7e8444161f1c17ecf1bf5 Mon Sep 17 00:00:00 2001 From: ctena Date: Fri, 28 Apr 2023 11:25:21 +0200 Subject: [PATCH 2/6] MONARCH outputs with NES seems to be working. --- nes/nc_projections/default_nes.py | 33 ++++++++++++++--------- nes/nes_formats/monarch_format.py | 45 ++++++++++++++++++++----------- 2 files changed, 50 insertions(+), 28 deletions(-) diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py index 2523b37..cef4939 100644 --- a/nes/nc_projections/default_nes.py +++ b/nes/nc_projections/default_nes.py @@ -378,7 +378,7 @@ class Nes(object): if self.cell_measures[cell_measure]['data'] is not None: del self.cell_measures[cell_measure]['data'] del self.cell_measures - except AttributeError: + except (AttributeError, KeyError): pass del self @@ -2061,10 +2061,10 @@ class Nes(object): """ for var_name, var_info in self.variables.items(): - if var_info['data'] is not None: + if isinstance(var_info['data'], np.ndarray): self.variables[var_name]['data'] = self.variables[var_name]['data'].astype(data_type) - self.variables[var_name]['dtype'] = data_type - + self.variables[var_name]['dtype'] = data_type + return None def concatenate(self, aux_nessy): @@ -2307,7 +2307,8 @@ class Nes(object): time_bnds_var[:] = date2num(self._time_bnds, time_var.units, calendar='standard') # LEVELS - lev = netcdf.createVariable('lev', np.float64, ('lev',), zlib=self.zip_lvl > 0, complevel=self.zip_lvl) + lev = netcdf.createVariable('lev', self._lev['data'].dtype, ('lev',), + zlib=self.zip_lvl > 0, complevel=self.zip_lvl) if 'units' in self._lev.keys(): lev.units = Units(self._lev['units'], formatted=True).units else: @@ -2320,7 +2321,8 @@ class Nes(object): lev[:] = self._lev['data'] # LATITUDES - lat = netcdf.createVariable('lat', np.float64, self._lat_dim, zlib=self.zip_lvl > 0, complevel=self.zip_lvl) + lat = netcdf.createVariable('lat', self._lat['data'].dtype, self._lat_dim, + zlib=self.zip_lvl > 0, complevel=self.zip_lvl) lat.units = 'degrees_north' lat.axis = 'Y' lat.long_name = 'latitude coordinate' @@ -2333,14 +2335,16 @@ class Nes(object): # LATITUDES BOUNDS if self._lat_bnds is not None: - lat_bnds_var = netcdf.createVariable('lat_bnds', np.float64, self._lat_dim + ('spatial_nv',), + lat_bnds_var = netcdf.createVariable('lat_bnds', self._lat_bnds['data'].dtype, + self._lat_dim + ('spatial_nv',), zlib=self.zip_lvl > 0, complevel=self.zip_lvl) if self.size > 1: lat_bnds_var.set_collective(True) lat_bnds_var[:] = self._lat_bnds['data'] # LONGITUDES - lon = netcdf.createVariable('lon', np.float64, self._lon_dim, zlib=self.zip_lvl > 0, complevel=self.zip_lvl) + lon = netcdf.createVariable('lon', self._lon['data'].dtype, self._lon_dim, + zlib=self.zip_lvl > 0, complevel=self.zip_lvl) lon.units = 'degrees_east' lon.axis = 'X' lon.long_name = 'longitude coordinate' @@ -2353,7 +2357,8 @@ class Nes(object): # LONGITUDES BOUNDS if self._lon_bnds is not None: - lon_bnds_var = netcdf.createVariable('lon_bnds', np.float64, self._lon_dim + ('spatial_nv',), + lon_bnds_var = netcdf.createVariable('lon_bnds', self._lon_bnds['data'].dtype, + self._lon_dim + ('spatial_nv',), zlib=self.zip_lvl > 0, complevel=self.zip_lvl) if self.size > 1: lon_bnds_var.set_collective(True) @@ -2365,13 +2370,13 @@ class Nes(object): # CELL AREA if 'cell_area' in self.cell_measures.keys(): - cell_area = netcdf.createVariable('cell_area', np.float64, self._var_dim, + cell_area = netcdf.createVariable('cell_area', self.cell_measures['cell_area']['data'].dtype, self._var_dim, zlib=self.zip_lvl > 0, complevel=self.zip_lvl) if self.size > 1: cell_area.set_collective(True) cell_area[self.read_axis_limits['y_min']:self.read_axis_limits['y_max'], - self.read_axis_limits['x_min']:self.read_axis_limits['x_max']] = self.cell_measures[ - 'cell_area']['data'] + self.read_axis_limits['x_min']:self.read_axis_limits['x_max']] = \ + self.cell_measures['cell_area']['data'] cell_area.long_name = 'area of grid cell' cell_area.standard_name = 'cell_area' @@ -2549,6 +2554,8 @@ class Nes(object): ---------- i_time : int index of the time step to write + out_format : str + Indicates the output format type to change the units (if needed) """ if self.serial_nc is not None: try: @@ -2558,7 +2565,7 @@ class Nes(object): data = self.__gather_data_py_object(self.variables) if self.master: self.serial_nc.variables = data - self.serial_nc.append_time_step_data(i_time) + self.serial_nc.append_time_step_data(i_time, out_format=out_format) self.comm.Barrier() else: if out_format == 'MONARCH': diff --git a/nes/nes_formats/monarch_format.py b/nes/nes_formats/monarch_format.py index 3a77335..778c4ec 100644 --- a/nes/nes_formats/monarch_format.py +++ b/nes/nes_formats/monarch_format.py @@ -1,13 +1,9 @@ #!/usr/bin/env python -import sys -import warnings import numpy as np -import os import nes from netCDF4 import Dataset from mpi4py import MPI -from copy import copy def to_netcdf_monarch(self, path, chunking=False, keep_open=False): @@ -42,11 +38,28 @@ def to_netcdf_monarch(self, path, chunking=False, keep_open=False): self._create_dimensions(netcdf) # Create dimension variables + self._lev['data'] = np.array(self._lev['data'], dtype=np.float32) + self._lat['data'] = np.array(self._lat['data'], dtype=np.float32) + self._lat_bnds['data'] = np.array(self._lat_bnds['data'], dtype=np.float32) + if hasattr(self, '_rlat'): + self._rlat['data'] = np.array(self._rlat['data'], dtype=np.float32) + if hasattr(self, '_y'): + self._y['data'] = np.array(self._y['data'], dtype=np.float32) + + self._lon['data'] = np.array(self._lon['data'], dtype=np.float32) + self._lon_bnds['data'] = np.array(self._lon_bnds['data'], dtype=np.float32) + if hasattr(self, '_rlon'): + self._rlon['data'] = np.array(self._rlon['data'], dtype=np.float32) + if hasattr(self, '_x'): + self._x['data'] = np.array(self._x['data'], dtype=np.float32) + self._create_dimension_variables(netcdf) if self.info: print("Rank {0:03d}: Dimensions done".format(self.rank)) # Create cell measures + if 'cell_area' in self.cell_measures.keys(): + self.cell_measures['cell_area']['data'] = np.array(self.cell_measures['cell_area']['data'], dtype=np.float32) self._create_cell_measures(netcdf) # Create variables @@ -84,15 +97,17 @@ def to_monarch_units(variables): Variable in the MONARCH units """ for var_name in variables.keys(): - if variables['units'] == 'mol.s-1.m-2': - # Kmol to mol - variables['data'] *= 1000 - elif variables['units'] == 'kg.s-1.m-2': - # No action needed - pass - else: - raise TypeError("The unit '{0}' of specie {1} is not defined correctly. ".format( - variables['units'], var_name) + - "Should be 'mol.s-1.m-2' or 'kg.s-1.m-2'") + if isinstance(variables[var_name]['data'], np.ndarray): + if variables[var_name]['units'] == 'mol.s-1.m-2': + # Kmol to mol + variables[var_name]['data'] = np.array(variables[var_name]['data'] * 1000, dtype=np.float32) + elif variables[var_name]['units'] == 'kg.s-1.m-2': + # No unit change needed + variables[var_name]['data'] = np.array(variables[var_name]['data'], dtype=np.float32) + + else: + raise TypeError("The unit '{0}' of specie {1} is not defined correctly. ".format( + variables[var_name]['units'], var_name) + + "Should be 'mol.s-1.m-2' or 'kg.s-1.m-2'") + variables[var_name]['dtype'] = np.float32 return variables - -- GitLab From d00ebe5e61823187ec3f299fbacf9072cc13c25d Mon Sep 17 00:00:00 2001 From: ctena Date: Tue, 2 May 2023 15:21:36 +0200 Subject: [PATCH 3/6] Bugfix when selecting lat/lon min/max and have lat/lon bounds --- nes/nc_projections/default_nes.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py index 2523b37..243476a 100644 --- a/nes/nc_projections/default_nes.py +++ b/nes/nc_projections/default_nes.py @@ -888,18 +888,20 @@ class Nes(object): self._lon['data'] = self._lon['data'][idx['idx_x_min']:idx['idx_x_max']] if self._lat_bnds is not None: - self._lat_bnds = self._lat_bnds[idx['idx_y_min']:idx['idx_y_max'], :] + self._lat_bnds['data'] = self._lat_bnds['data'][idx['idx_y_min']:idx['idx_y_max'], :] if self._lon_bnds is not None: - self._lon_bnds = self._lon_bnds[idx['idx_x_min']:idx['idx_x_max'], :] + self._lon_bnds['data'] = self._lon_bnds['data'][idx['idx_x_min']:idx['idx_x_max'], :] else: # Irregular projections self._lat['data'] = self._lat['data'][idx['idx_y_min']:idx['idx_y_max'], idx['idx_x_min']:idx['idx_x_max']] self._lon['data'] = self._lon['data'][idx['idx_y_min']:idx['idx_y_max'], idx['idx_x_min']:idx['idx_x_max']] if self._lat_bnds is not None: - self._lat_bnds = self._lat_bnds[idx['idx_y_min']:idx['idx_y_max'], idx['idx_x_min']:idx['idx_x_max'], :] + self._lat_bnds['data'] = self._lat_bnds['data'][idx['idx_y_min']:idx['idx_y_max'], + idx['idx_x_min']:idx['idx_x_max'], :] if self._lon_bnds is not None: - self._lon_bnds = self._lon_bnds[idx['idx_y_min']:idx['idx_y_max'], idx['idx_x_min']:idx['idx_x_max'], :] + self._lon_bnds['data'] = self._lon_bnds['data'][idx['idx_y_min']:idx['idx_y_max'], + idx['idx_x_min']:idx['idx_x_max'], :] self.hours_start = 0 self.hours_end = 0 -- GitLab From 1ce4aecee6a113978b868434b1029df3535dc556 Mon Sep 17 00:00:00 2001 From: ctena Date: Tue, 2 May 2023 15:37:58 +0200 Subject: [PATCH 4/6] Bugfix when selecting and concatenating --- nes/nc_projections/default_nes.py | 1 + 1 file changed, 1 insertion(+) diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py index f84c503..0ed13b0 100644 --- a/nes/nc_projections/default_nes.py +++ b/nes/nc_projections/default_nes.py @@ -2094,6 +2094,7 @@ class Nes(object): new = False for var_name, var_info in aux_nessy.variables.items(): if var_info['data'] is None: + aux_nessy.read_axis_limits = self.read_axis_limits aux_nessy.load(var_name) new_vars_added = [] -- GitLab From f7e6450bd97a410dbcee49418c8a671b45294504 Mon Sep 17 00:00:00 2001 From: ctena Date: Fri, 5 May 2023 16:29:38 +0200 Subject: [PATCH 5/6] CMAQ outputs with NES seems to be working. --- nes/nc_projections/default_nes.py | 10 +- nes/nes_formats/__init__.py | 1 + nes/nes_formats/cmaq_format.py | 345 ++++++++++++++++++++++++++++++ nes/nes_formats/monarch_format.py | 35 ++- 4 files changed, 370 insertions(+), 21 deletions(-) create mode 100644 nes/nes_formats/cmaq_format.py diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py index cef4939..da66b99 100644 --- a/nes/nc_projections/default_nes.py +++ b/nes/nc_projections/default_nes.py @@ -16,7 +16,7 @@ from copy import deepcopy, copy import datetime import pyproj from ..methods import vertical_interpolation, horizontal_interpolation, cell_measures, spatial_join -from ..nes_formats import to_netcdf_cams_ra, to_netcdf_monarch, to_monarch_units +from ..nes_formats import to_netcdf_cams_ra, to_netcdf_monarch, to_monarch_units, to_netcdf_cmaq, to_cmaq_units class Nes(object): @@ -2569,7 +2569,9 @@ class Nes(object): self.comm.Barrier() else: if out_format == 'MONARCH': - self.variables = to_monarch_units(self.variables) + self.variables = to_monarch_units(self) + if out_format == 'CMAQ': + self.variables = to_cmaq_units(self) for i, (var_name, var_dict) in enumerate(self.variables.items()): for att_name, att_value in var_dict.items(): if att_name == 'data': @@ -2749,6 +2751,8 @@ class Nes(object): new_nc.__to_netcdf_cams_ra(path) elif type == 'MONARCH': to_netcdf_monarch(new_nc, path, chunking=chunking, keep_open=keep_open) + elif type == 'CMAQ': + to_netcdf_cmaq(new_nc, path, keep_open=keep_open) else: msg = "Unknown NetCDF type '{0}'. ".format(nc_type) msg += "Use CAMS_RA, MONARCH or NES (or DEFAULT)" @@ -2763,6 +2767,8 @@ class Nes(object): self.__to_netcdf_cams_ra(path) elif nc_type == 'MONARCH': to_netcdf_monarch(self, path, chunking=chunking, keep_open=keep_open) + elif nc_type == 'MONARCH': + to_netcdf_cmaq(self, path, keep_open=keep_open) else: msg = "Unknown NetCDF type '{0}'. ".format(nc_type) msg += "Use CAMS_RA, MONARCH or NES (or DEFAULT)" diff --git a/nes/nes_formats/__init__.py b/nes/nes_formats/__init__.py index f597c54..b8083ee 100644 --- a/nes/nes_formats/__init__.py +++ b/nes/nes_formats/__init__.py @@ -1,2 +1,3 @@ from .cams_ra_format import to_netcdf_cams_ra from .monarch_format import to_netcdf_monarch, to_monarch_units +from .cmaq_format import to_netcdf_cmaq, to_cmaq_units diff --git a/nes/nes_formats/cmaq_format.py b/nes/nes_formats/cmaq_format.py new file mode 100644 index 0000000..27a7df7 --- /dev/null +++ b/nes/nes_formats/cmaq_format.py @@ -0,0 +1,345 @@ +#!/usr/bin/env python + +import numpy as np +import nes +from netCDF4 import Dataset +from mpi4py import MPI +from copy import deepcopy +from datetime import datetime + +GLOBAL_ATTRIBUTES_ORDER = [ + 'IOAPI_VERSION', 'EXEC_ID', 'FTYPE', 'CDATE', 'CTIME', 'WDATE', 'WTIME', 'SDATE', 'STIME', 'TSTEP', 'NTHIK', + 'NCOLS', 'NROWS', 'NLAYS', 'NVARS', 'GDTYP', 'P_ALP', 'P_BET', 'P_GAM', 'XCENT', 'YCENT', 'XORIG', 'YORIG', + 'XCELL', 'YCELL', 'VGTYP', 'VGTOP', 'VGLVLS', 'GDNAM', 'UPNAM', 'FILEDESC', 'HISTORY', 'VAR-LIST'] + + +# noinspection DuplicatedCode +def to_netcdf_cmaq(self, path, chunking=False, keep_open=False): + """ + Create the NetCDF using netcdf4-python methods. + + Parameters + ---------- + self : nes.Nes + Source projection Nes Object. + path : str + Path to the output netCDF file. + chunking: bool + Indicates if you want to chunk the output netCDF. + keep_open : bool + Indicates if you want to keep open the NetCDH to fill the data by time-step + """ + self.to_dtype(np.float32) + + set_global_attributes(self) + change_variable_attributes(self) + + # Open NetCDF + if self.info: + print("Rank {0:03d}: Creating {1}".format(self.rank, path)) + if self.size > 1: + netcdf = Dataset(path, format="NETCDF4", mode='w', parallel=True, comm=self.comm, info=MPI.Info()) + else: + netcdf = Dataset(path, format="NETCDF4", mode='w', parallel=False) + if self.info: + print("Rank {0:03d}: NetCDF ready to write".format(self.rank)) + + # Create dimensions + create_dimensions(self, netcdf) + + create_dimension_variables(self, netcdf) + if self.info: + print("Rank {0:03d}: Dimensions done".format(self.rank)) + + # Create variables + create_variables(self, netcdf) + + for att_name in GLOBAL_ATTRIBUTES_ORDER: + netcdf.setncattr(att_name, self.global_attrs[att_name]) + + # Close NetCDF + if keep_open: + self.netcdf = netcdf + else: + netcdf.close() + + return None + + +def change_variable_attributes(self): + """ + + Parameters + ---------- + self : nes.Nes + """ + """ + Modify the emission list to be consistent to use the output as input for CMAQ model. + + :return: Emission list ready for CMAQ + :rtype: dict + """ + for var_name in self.variables.keys(): + + if self.variables[var_name]['units'] == 'mol.s-1': + self.variables[var_name]['units'] = "{:<16}".format('mole/s') + self.variables[var_name]['var_desc'] = "{:<80}".format(self.variables[var_name]['long_name']) + self.variables[var_name]['long_name'] = "{:<16}".format(var_name) + elif self.variables[var_name]['units'] == 'g.s-1': + self.variables[var_name]['units'] = "{:<16}".format('g/s') + self.variables[var_name]['var_desc'] = "{:<80}".format(self.variables[var_name]['long_name']) + self.variables[var_name]['long_name'] = "{:<16}".format(var_name) + + else: + raise TypeError("The unit '{0}' of specie {1} is not defined correctly. ".format( + self.variables[var_name]['units'], var_name) + + "Should be 'mol.s-1.m-2' or 'kg.s-1.m-2'") + return None + + +def to_cmaq_units(self): + """ + Change the data values according to the CMAQ conventions + + Parameters + ---------- + self : nes.Nes + + Returns + ------- + dict + Variable in the MONARCH units + """ + self.calculate_grid_area() + for var_name in self.variables.keys(): + if isinstance(self.variables[var_name]['data'], np.ndarray): + if self.variables[var_name]['units'] == 'mol.s-1': + # Kmol.m-2.s-1 to mol.s-1 + self.variables[var_name]['data'] = np.array( + self.variables[var_name]['data'] * 1000 * self.cell_measures['cell_area']['data'], dtype=np.float32) + elif self.variables[var_name]['units'] == 'g.s-1': + # Kg.m-2.s-1 to g.s-1 + self.variables[var_name]['data'] = np.array( + self.variables[var_name]['data'] * 1000 * self.cell_measures['cell_area']['data'], dtype=np.float32) + + else: + raise TypeError("The unit '{0}' of specie {1} is not defined correctly. ".format( + self.variables[var_name]['units'], var_name) + + "Should be 'mol.s-1.m-2' or 'kg.s-1.m-2'") + self.variables[var_name]['dtype'] = np.float32 + + return self.variables + + +def create_tflag(self): + """ + Create the content of the CMAQ variable TFLAG + + Parameters + ---------- + self : nes.Nes + + Returns + ------- + numpy.ndarray + Array with the content of TFLAG + """ + t_flag = np.empty((len(self.time), len(self.variables), 2)) + + for i_d, aux_date in enumerate(self.time): + y_d = int(aux_date.strftime('%Y%j')) + hms = int(aux_date.strftime('%H%M%S')) + for i_p in range(len(self.variables)): + t_flag[i_d, i_p, 0] = y_d + t_flag[i_d, i_p, 1] = hms + + return t_flag + + +def str_var_list(self): + """ + Transform the list of variable names to a string with the elements with 16 white spaces. + + Parameters + ---------- + self : nes.Nes + + Returns + ------- + str + List of variable names transformed on string. + """ + str_var_list = "" + for var in self.variables.keys(): + str_var_list += "{:<16}".format(var) + + return str_var_list + + +def set_global_attributes(self): + """ + Set the NetCDF global attributes + + Parameters + ---------- + self : nes.Nes + """ + now = datetime.now() + if len(self.time) > 1: + tstep = ((self.time[1] - self.time[0]).seconds // 3600) * 10000 + else: + tstep = 1 * 10000 + + current_attributes = deepcopy(self.global_attrs) + del self.global_attrs + + self.global_attrs = {'IOAPI_VERSION': 'None: made only with NetCDF libraries', + 'EXEC_ID': "{:<80}".format('0.1alpha'), # Editable + 'FTYPE': np.int32(1), # Editable + 'CDATE': np.int32(now.strftime('%Y%j')), + 'CTIME': np.int32(now.strftime('%H%M%S')), + 'WDATE': np.int32(now.strftime('%Y%j')), + 'WTIME': np.int32(now.strftime('%H%M%S')), + 'SDATE': np.int32(self.time[0].strftime('%Y%j')), + 'STIME': np.int32(self.time[0].strftime('%H%M%S')), + 'TSTEP': np.int32(tstep), + 'NTHIK': np.int32(1), # Editable + 'NCOLS': None, # Projection dependent + 'NROWS': None, # Projection dependent + 'NLAYS': np.int32(len(self.lev['data'])), + 'NVARS': None, # Projection dependent + 'GDTYP': None, # Projection dependent + 'P_ALP': None, # Projection dependent + 'P_BET': None, # Projection dependent + 'P_GAM': None, # Projection dependent + 'XCENT': None, # Projection dependent + 'YCENT': None, # Projection dependent + 'XORIG': None, # Projection dependent + 'YORIG': None, # Projection dependent + 'XCELL': None, # Projection dependent + 'YCELL': None, # Projection dependent + 'VGTYP': np.int32(7), # Editable + 'VGTOP': np.float32(5000.), # Editable + 'VGLVLS': np.array([1., 0.], dtype=np.float32), # Editable + 'GDNAM': "{:<16}".format(''), # Editable + 'UPNAM': "{:<16}".format('HERMESv3'), + 'FILEDESC': "", # Editable + 'HISTORY': "", # Editable + 'VAR-LIST': str_var_list(self)} + + # Editable attributes + for att_name, att_value in current_attributes.items(): + if att_name == 'EXEC_ID': + self.global_attrs[att_name] = "{:<80}".format(att_value) # Editable + elif att_name == 'FTYPE': + self.global_attrs[att_name] = np.int32(att_value) # Editable + elif att_name == 'NTHIK': + self.global_attrs[att_name] = np.int32(att_value) # Editable + elif att_name == 'VGTYP': + self.global_attrs[att_name] = np.int32(att_value) # Editable + elif att_name == 'VGTOP': + self.global_attrs[att_name] = np.float32(att_value) # Editable + elif att_name == 'VGLVLS': + self.global_attrs[att_name] = np.array(att_value.split(), dtype=np.float32) # Editable + elif att_name == 'GDNAM': + self.global_attrs[att_name] = "{:<16}".format(att_value) # Editable + elif att_name == 'FILEDESC': + self.global_attrs[att_name] = att_value # Editable + elif att_name == 'HISTORY': + self.global_attrs[att_name] = att_value # Editable + + # Projection dependent attributes + if isinstance(self, nes.LCCNes): + self.global_attrs['NCOLS'] = np.int32(len(self.x['data'])) + self.global_attrs['NROWS'] = np.int32(len(self.y['data'])) + self.global_attrs['NVARS'] = np.int32(len(self.variables)) + self.global_attrs['GDTYP'] = np.int32(2) + + self.global_attrs['P_ALP'] = np.float64(self.projection_data['standard_parallel'][0]) + self.global_attrs['P_BET'] = np.float64(self.projection_data['standard_parallel'][1]) + self.global_attrs['P_GAM'] = np.float64(self.projection_data['longitude_of_central_meridian']) + self.global_attrs['XCENT'] = np.float64(self.projection_data['longitude_of_central_meridian']) + self.global_attrs['YCENT'] = np.float64(self.projection_data['latitude_of_projection_origin']) + self.global_attrs['XORIG'] = np.float64( + self.x['data'][0]) - (np.float64(self.x['data'][1] - self.x['data'][0]) / 2) + self.global_attrs['YORIG'] = np.float64( + self.y['data'][0]) - (np.float64(self.y['data'][1] - self.y['data'][0]) / 2) + self.global_attrs['XCELL'] = np.float64(self.x['data'][1] - self.x['data'][0]) + self.global_attrs['YCELL'] = np.float64(self.y['data'][1] - self.y['data'][0]) + + return None + + +def create_dimensions(self, netcdf): + """ + Create 'time', 'time_bnds', 'lev', 'lon' and 'lat' dimensions. + + Parameters + ---------- + self : nes.Nes + netcdf : Dataset + netcdf4-python open dataset. + """ + netcdf.createDimension('TSTEP', len(self._time)) + netcdf.createDimension('DATE-TIME', 2) + netcdf.createDimension('LAY', len(self._lev['data'])) + netcdf.createDimension('VAR', len(self.variables)) + if isinstance(self, nes.LCCNes): + netcdf.createDimension('COL', len(self._x['data'])) + netcdf.createDimension('ROW', len(self._y['data'])) + + return None + + +def create_dimension_variables(self, netcdf): + """ + Create the 'y' and 'x' variables. + + Parameters + ---------- + self : nes.Nes + netcdf : Dataset + NetCDF object. + """ + + tflag = netcdf.createVariable('TFLAG', 'i', ('TSTEP', 'VAR', 'DATE-TIME',)) + tflag.setncatts({'units': "{:<16}".format(''), 'long_name': "{:<16}".format('TFLAG'), + 'var_desc': "{:<80}".format('Timestep-valid flags: (1) YYYYDDD or (2) HHMMSS')}) + tflag[:] = create_tflag(self) + + return None + + +def create_variables(self, netcdf): + """ + Create the netCDF file variables. + + Parameters + ---------- + self : nes.Nes + netcdf : Dataset + netcdf4-python open dataset. + """ + for var_name, var_info in self.variables.items(): + var = netcdf.createVariable(var_name, 'f', ('TSTEP', 'LAY', 'ROW', 'COL',), + zlib=self.zip_lvl > 0, complevel=self.zip_lvl) + var.units = var_info['units'] + var.long_name = str(var_info['long_name']) + var.var_desc = str(var_info['var_desc']) + if var_info['data'] is not None: + if self.info: + print("Rank {0:03d}: Filling {1})".format(self.rank, var_name)) + + if isinstance(var_info['data'], int) and var_info['data'] == 0: + var[self.write_axis_limits['t_min']:self.write_axis_limits['t_max'], + self.write_axis_limits['z_min']:self.write_axis_limits['z_max'], + self.write_axis_limits['y_min']:self.write_axis_limits['y_max'], + self.write_axis_limits['x_min']:self.write_axis_limits['x_max']] = 0 + + elif len(var_info['data'].shape) == 4: + var[self.write_axis_limits['t_min']:self.write_axis_limits['t_max'], + self.write_axis_limits['z_min']:self.write_axis_limits['z_max'], + self.write_axis_limits['y_min']:self.write_axis_limits['y_max'], + self.write_axis_limits['x_min']:self.write_axis_limits['x_max']] = var_info['data'] + + return None diff --git a/nes/nes_formats/monarch_format.py b/nes/nes_formats/monarch_format.py index 778c4ec..0087472 100644 --- a/nes/nes_formats/monarch_format.py +++ b/nes/nes_formats/monarch_format.py @@ -6,6 +6,7 @@ from netCDF4 import Dataset from mpi4py import MPI +# noinspection DuplicatedCode def to_netcdf_monarch(self, path, chunking=False, keep_open=False): """ Create the NetCDF using netcdf4-python methods. @@ -41,16 +42,13 @@ def to_netcdf_monarch(self, path, chunking=False, keep_open=False): self._lev['data'] = np.array(self._lev['data'], dtype=np.float32) self._lat['data'] = np.array(self._lat['data'], dtype=np.float32) self._lat_bnds['data'] = np.array(self._lat_bnds['data'], dtype=np.float32) - if hasattr(self, '_rlat'): - self._rlat['data'] = np.array(self._rlat['data'], dtype=np.float32) - if hasattr(self, '_y'): - self._y['data'] = np.array(self._y['data'], dtype=np.float32) - self._lon['data'] = np.array(self._lon['data'], dtype=np.float32) self._lon_bnds['data'] = np.array(self._lon_bnds['data'], dtype=np.float32) - if hasattr(self, '_rlon'): + if isinstance(self, nes.RotatedNes): + self._rlat['data'] = np.array(self._rlat['data'], dtype=np.float32) self._rlon['data'] = np.array(self._rlon['data'], dtype=np.float32) - if hasattr(self, '_x'): + if isinstance(self, nes.LCCNes) or isinstance(self, nes.MercatorNes): + self._y['data'] = np.array(self._y['data'], dtype=np.float32) self._x['data'] = np.array(self._x['data'], dtype=np.float32) self._create_dimension_variables(netcdf) @@ -82,32 +80,31 @@ def to_netcdf_monarch(self, path, chunking=False, keep_open=False): return None -def to_monarch_units(variables): +def to_monarch_units(self): """ Change the data values according to the MONARCH conventions Parameters ---------- - variables : dict - Variables data + self : nes.Nes Returns ------- dict Variable in the MONARCH units """ - for var_name in variables.keys(): - if isinstance(variables[var_name]['data'], np.ndarray): - if variables[var_name]['units'] == 'mol.s-1.m-2': + for var_name in self.variables.keys(): + if isinstance(self.variables[var_name]['data'], np.ndarray): + if self.variables[var_name]['units'] == 'mol.s-1.m-2': # Kmol to mol - variables[var_name]['data'] = np.array(variables[var_name]['data'] * 1000, dtype=np.float32) - elif variables[var_name]['units'] == 'kg.s-1.m-2': + self.variables[var_name]['data'] = np.array(self.variables[var_name]['data'] * 1000, dtype=np.float32) + elif self.variables[var_name]['units'] == 'kg.s-1.m-2': # No unit change needed - variables[var_name]['data'] = np.array(variables[var_name]['data'], dtype=np.float32) + self.variables[var_name]['data'] = np.array(self.variables[var_name]['data'], dtype=np.float32) else: raise TypeError("The unit '{0}' of specie {1} is not defined correctly. ".format( - variables[var_name]['units'], var_name) + + self.variables[var_name]['units'], var_name) + "Should be 'mol.s-1.m-2' or 'kg.s-1.m-2'") - variables[var_name]['dtype'] = np.float32 - return variables + self.variables[var_name]['dtype'] = np.float32 + return self.variables -- GitLab From db1d476742cf92ef85450e4769c7b44919eb9f06 Mon Sep 17 00:00:00 2001 From: ctena Date: Thu, 11 May 2023 13:13:20 +0200 Subject: [PATCH 6/6] WRF_CHEM outputs with NES seems to be working. --- nes/nc_projections/default_nes.py | 15 +- nes/nes_formats/__init__.py | 1 + nes/nes_formats/cmaq_format.py | 25 +- nes/nes_formats/wrf_chem_format.py | 385 +++++++++++++++++++++++++++++ 4 files changed, 406 insertions(+), 20 deletions(-) create mode 100644 nes/nes_formats/wrf_chem_format.py diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py index 9961191..ea8afe6 100644 --- a/nes/nc_projections/default_nes.py +++ b/nes/nc_projections/default_nes.py @@ -16,7 +16,8 @@ from copy import deepcopy, copy import datetime import pyproj from ..methods import vertical_interpolation, horizontal_interpolation, cell_measures, spatial_join -from ..nes_formats import to_netcdf_cams_ra, to_netcdf_monarch, to_monarch_units, to_netcdf_cmaq, to_cmaq_units +from ..nes_formats import to_netcdf_cams_ra, to_netcdf_monarch, to_monarch_units, to_netcdf_cmaq, to_cmaq_units, \ + to_netcdf_wrf_chem, to_wrf_chem_units class Nes(object): @@ -341,7 +342,7 @@ class Nes(object): Parameters ---------- - strlen : int + strlen : int or None Max length of the string """ @@ -2573,8 +2574,10 @@ class Nes(object): else: if out_format == 'MONARCH': self.variables = to_monarch_units(self) - if out_format == 'CMAQ': + elif out_format == 'CMAQ': self.variables = to_cmaq_units(self) + elif out_format == 'WRF_CHEM': + self.variables = to_wrf_chem_units(self) for i, (var_name, var_dict) in enumerate(self.variables.items()): for att_name, att_value in var_dict.items(): if att_name == 'data': @@ -2756,6 +2759,8 @@ class Nes(object): to_netcdf_monarch(new_nc, path, chunking=chunking, keep_open=keep_open) elif type == 'CMAQ': to_netcdf_cmaq(new_nc, path, keep_open=keep_open) + elif type == 'WRF_CHEM': + to_netcdf_wrf_chem(new_nc, path, keep_open=keep_open) else: msg = "Unknown NetCDF type '{0}'. ".format(nc_type) msg += "Use CAMS_RA, MONARCH or NES (or DEFAULT)" @@ -2770,8 +2775,10 @@ class Nes(object): self.__to_netcdf_cams_ra(path) elif nc_type == 'MONARCH': to_netcdf_monarch(self, path, chunking=chunking, keep_open=keep_open) - elif nc_type == 'MONARCH': + elif nc_type == 'CMAQ': to_netcdf_cmaq(self, path, keep_open=keep_open) + elif nc_type == 'WRF_CHEM': + to_netcdf_wrf_chem(self, path, keep_open=keep_open) else: msg = "Unknown NetCDF type '{0}'. ".format(nc_type) msg += "Use CAMS_RA, MONARCH or NES (or DEFAULT)" diff --git a/nes/nes_formats/__init__.py b/nes/nes_formats/__init__.py index b8083ee..e09d18e 100644 --- a/nes/nes_formats/__init__.py +++ b/nes/nes_formats/__init__.py @@ -1,3 +1,4 @@ from .cams_ra_format import to_netcdf_cams_ra from .monarch_format import to_netcdf_monarch, to_monarch_units from .cmaq_format import to_netcdf_cmaq, to_cmaq_units +from .wrf_chem_format import to_netcdf_wrf_chem, to_wrf_chem_units diff --git a/nes/nes_formats/cmaq_format.py b/nes/nes_formats/cmaq_format.py index 27a7df7..649d9d5 100644 --- a/nes/nes_formats/cmaq_format.py +++ b/nes/nes_formats/cmaq_format.py @@ -68,17 +68,12 @@ def to_netcdf_cmaq(self, path, chunking=False, keep_open=False): def change_variable_attributes(self): """ + Modify the emission list to be consistent to use the output as input for CMAQ model. Parameters ---------- self : nes.Nes """ - """ - Modify the emission list to be consistent to use the output as input for CMAQ model. - - :return: Emission list ready for CMAQ - :rtype: dict - """ for var_name in self.variables.keys(): if self.variables[var_name]['units'] == 'mol.s-1': @@ -92,8 +87,7 @@ def change_variable_attributes(self): else: raise TypeError("The unit '{0}' of specie {1} is not defined correctly. ".format( - self.variables[var_name]['units'], var_name) + - "Should be 'mol.s-1.m-2' or 'kg.s-1.m-2'") + self.variables[var_name]['units'], var_name) + "Should be 'mol.s-1' or 'g.s-1'") return None @@ -124,8 +118,7 @@ def to_cmaq_units(self): else: raise TypeError("The unit '{0}' of specie {1} is not defined correctly. ".format( - self.variables[var_name]['units'], var_name) + - "Should be 'mol.s-1.m-2' or 'kg.s-1.m-2'") + self.variables[var_name]['units'], var_name) + "Should be 'mol.s-1' or 'g.s-1'") self.variables[var_name]['dtype'] = np.float32 return self.variables @@ -250,8 +243,8 @@ def set_global_attributes(self): # Projection dependent attributes if isinstance(self, nes.LCCNes): - self.global_attrs['NCOLS'] = np.int32(len(self.x['data'])) - self.global_attrs['NROWS'] = np.int32(len(self.y['data'])) + self.global_attrs['NCOLS'] = np.int32(len(self._x['data'])) + self.global_attrs['NROWS'] = np.int32(len(self._y['data'])) self.global_attrs['NVARS'] = np.int32(len(self.variables)) self.global_attrs['GDTYP'] = np.int32(2) @@ -261,11 +254,11 @@ def set_global_attributes(self): self.global_attrs['XCENT'] = np.float64(self.projection_data['longitude_of_central_meridian']) self.global_attrs['YCENT'] = np.float64(self.projection_data['latitude_of_projection_origin']) self.global_attrs['XORIG'] = np.float64( - self.x['data'][0]) - (np.float64(self.x['data'][1] - self.x['data'][0]) / 2) + self._x['data'][0]) - (np.float64(self._x['data'][1] - self._x['data'][0]) / 2) self.global_attrs['YORIG'] = np.float64( - self.y['data'][0]) - (np.float64(self.y['data'][1] - self.y['data'][0]) / 2) - self.global_attrs['XCELL'] = np.float64(self.x['data'][1] - self.x['data'][0]) - self.global_attrs['YCELL'] = np.float64(self.y['data'][1] - self.y['data'][0]) + self._y['data'][0]) - (np.float64(self._y['data'][1] - self._y['data'][0]) / 2) + self.global_attrs['XCELL'] = np.float64(self._x['data'][1] - self._x['data'][0]) + self.global_attrs['YCELL'] = np.float64(self._y['data'][1] - self._y['data'][0]) return None diff --git a/nes/nes_formats/wrf_chem_format.py b/nes/nes_formats/wrf_chem_format.py new file mode 100644 index 0000000..a74da56 --- /dev/null +++ b/nes/nes_formats/wrf_chem_format.py @@ -0,0 +1,385 @@ +#!/usr/bin/env python + +import numpy as np +import nes +from netCDF4 import Dataset +from mpi4py import MPI +from copy import deepcopy +from datetime import datetime + +GLOBAL_ATTRIBUTES_ORDER = [ + 'TITLE', 'START_DATE', 'WEST-EAST_GRID_DIMENSION', 'SOUTH-NORTH_GRID_DIMENSION', 'BOTTOM-TOP_GRID_DIMENSION', 'DX', + 'DY', 'GRIDTYPE', 'DIFF_OPT', 'KM_OPT', 'DAMP_OPT', 'DAMPCOEF', 'KHDIF', 'KVDIF', 'MP_PHYSICS', 'RA_LW_PHYSICS', + 'RA_SW_PHYSICS', 'SF_SFCLAY_PHYSICS', 'SF_SURFACE_PHYSICS', 'BL_PBL_PHYSICS', 'CU_PHYSICS', 'SF_LAKE_PHYSICS', + 'SURFACE_INPUT_SOURCE', 'SST_UPDATE', 'GRID_FDDA', 'GFDDA_INTERVAL_M', 'GFDDA_END_H', 'GRID_SFDDA', + 'SGFDDA_INTERVAL_M', 'SGFDDA_END_H', 'WEST-EAST_PATCH_START_UNSTAG', 'WEST-EAST_PATCH_END_UNSTAG', + 'WEST-EAST_PATCH_START_STAG', 'WEST-EAST_PATCH_END_STAG', 'SOUTH-NORTH_PATCH_START_UNSTAG', + 'SOUTH-NORTH_PATCH_END_UNSTAG', 'SOUTH-NORTH_PATCH_START_STAG', 'SOUTH-NORTH_PATCH_END_STAG', + 'BOTTOM-TOP_PATCH_START_UNSTAG', 'BOTTOM-TOP_PATCH_END_UNSTAG', 'BOTTOM-TOP_PATCH_START_STAG', + 'BOTTOM-TOP_PATCH_END_STAG', 'GRID_ID', 'PARENT_ID', 'I_PARENT_START', 'J_PARENT_START', 'PARENT_GRID_RATIO', 'DT', + 'CEN_LAT', 'CEN_LON', 'TRUELAT1', 'TRUELAT2', 'MOAD_CEN_LAT', 'STAND_LON', 'POLE_LAT', 'POLE_LON', 'GMT', 'JULYR', + 'JULDAY', 'MAP_PROJ', 'MMINLU', 'NUM_LAND_CAT', 'ISWATER', 'ISLAKE', 'ISICE', 'ISURBAN', 'ISOILWATER'] + + +# noinspection DuplicatedCode +def to_netcdf_wrf_chem(self, path, chunking=False, keep_open=False): + """ + Create the NetCDF using netcdf4-python methods. + + Parameters + ---------- + self : nes.Nes + Source projection Nes Object. + path : str + Path to the output netCDF file. + chunking: bool + Indicates if you want to chunk the output netCDF. + keep_open : bool + Indicates if you want to keep open the NetCDH to fill the data by time-step + """ + self.to_dtype(np.float32) + + set_global_attributes(self) + change_variable_attributes(self) + + # Open NetCDF + if self.info: + print("Rank {0:03d}: Creating {1}".format(self.rank, path)) + if self.size > 1: + netcdf = Dataset(path, format="NETCDF4", mode='w', parallel=True, comm=self.comm, info=MPI.Info()) + else: + netcdf = Dataset(path, format="NETCDF4", mode='w', parallel=False) + if self.info: + print("Rank {0:03d}: NetCDF ready to write".format(self.rank)) + + # Create dimensions + create_dimensions(self, netcdf) + + create_dimension_variables(self, netcdf) + if self.info: + print("Rank {0:03d}: Dimensions done".format(self.rank)) + + # Create variables + create_variables(self, netcdf) + + for att_name in GLOBAL_ATTRIBUTES_ORDER: + netcdf.setncattr(att_name, self.global_attrs[att_name]) + + # Close NetCDF + if keep_open: + self.netcdf = netcdf + else: + netcdf.close() + + return None + + +def change_variable_attributes(self): + """ + Modify the emission list to be consistent to use the output as input for WRF-CHEM model. + + Parameters + ---------- + self : nes.Nes + """ + for var_name in self.variables.keys(): + if self.variables[var_name]['units'] == 'mol.h-1.km-2': + self.variables[var_name]['FieldType'] = np.int32(104) + self.variables[var_name]['MemoryOrder'] = "XYZ" + self.variables[var_name]['description'] = "EMISSIONS" + self.variables[var_name]['units'] = "mol km^-2 hr^-1" + self.variables[var_name]['stagger'] = "" + self.variables[var_name]['coordinates'] = "XLONG XLAT" + + elif self.variables[var_name]['units'] == 'ug.s-1.m-2': + self.variables[var_name]['FieldType'] = np.int32(104) + self.variables[var_name]['MemoryOrder'] = "XYZ" + self.variables[var_name]['description'] = "EMISSIONS" + self.variables[var_name]['units'] = "ug/m3 m/s" + self.variables[var_name]['stagger'] = "" + self.variables[var_name]['coordinates'] = "XLONG XLAT" + + else: + raise TypeError("The unit '{0}' of specie {1} is not defined correctly. ".format( + self.variables[var_name]['units'], var_name) + "Should be 'mol.h-1.km-2' or 'ug.s-1.m-2'") + + if 'long_name' in self.variables[var_name].keys(): + del self.variables[var_name]['long_name'] + + return None + + +def to_wrf_chem_units(self): + """ + Change the data values according to the WRF-CHEM conventions + + Parameters + ---------- + self : nes.Nes + + Returns + ------- + dict + Variable in the MONARCH units + """ + self.calculate_grid_area() + for var_name in self.variables.keys(): + if isinstance(self.variables[var_name]['data'], np.ndarray): + if self.variables[var_name]['units'] == 'mol.h-1.km-2': + # 10**6 -> from m2 to km2 + # 10**3 -> from kmol to mol + # 3600 -> from s to h + self.variables[var_name]['data'] = np.array( + self.variables[var_name]['data'] * 10 ** 6 * 10 ** 3 * 3600, dtype=np.float32) + elif self.variables[var_name]['units'] == 'ug.s-1.m-2': + # 10**9 -> from kg to ug + self.variables[var_name]['data'] = np.array( + self.variables[var_name]['data'] * 10 ** 9, dtype=np.float32) + + else: + raise TypeError("The unit '{0}' of specie {1} is not defined correctly. ".format( + self.variables[var_name]['units'], var_name) + "Should be 'mol.h-1.km-2' or 'ug.s-1.m-2'") + self.variables[var_name]['dtype'] = np.float32 + + return self.variables + + +def create_times_var(self): + """ + Create the content of the WRF-CHEM variable times + + Parameters + ---------- + self : nes.Nes + + Returns + ------- + numpy.ndarray + Array with the content of TFLAG + """ + aux_times = np.chararray((len(self.time), 19), itemsize=1) + + for i_d, aux_date in enumerate(self.time): + aux_times[i_d] = list(aux_date.strftime("%Y-%m-%d_%H:%M:%S")) + + return aux_times + + +def set_global_attributes(self): + """ + Set the NetCDF global attributes + + Parameters + ---------- + self : nes.Nes + """ + now = datetime.now() + if len(self.time) > 1: + tstep = ((self.time[1] - self.time[0]).seconds // 3600) * 10000 + else: + tstep = 1 * 10000 + + current_attributes = deepcopy(self.global_attrs) + del self.global_attrs + + self.global_attrs = {'TITLE': None, + 'START_DATE': self.time[0].strftime("%Y-%m-%d_%H:%M:%S"), + 'WEST-EAST_GRID_DIMENSION': None, # Projection dependent attributes + 'SOUTH-NORTH_GRID_DIMENSION': None, # Projection dependent attributes + 'BOTTOM-TOP_GRID_DIMENSION': np.int32(45), + 'DX': None, # Projection dependent attributes + 'DY': None, # Projection dependent attributes + 'GRIDTYPE': 'C', + 'DIFF_OPT': np.int32(1), + 'KM_OPT': np.int32(4), + 'DAMP_OPT': np.int32(3), + 'DAMPCOEF': np.float32(0.2), + 'KHDIF': np.float32(0.), + 'KVDIF': np.float32(0.), + 'MP_PHYSICS': np.int32(6), + 'RA_LW_PHYSICS': np.int32(4), + 'RA_SW_PHYSICS': np.int32(4), + 'SF_SFCLAY_PHYSICS': np.int32(2), + 'SF_SURFACE_PHYSICS': np.int32(2), + 'BL_PBL_PHYSICS': np.int32(8), + 'CU_PHYSICS': np.int32(0), + 'SF_LAKE_PHYSICS': np.int32(0), + 'SURFACE_INPUT_SOURCE': None, # Projection dependent attributes + 'SST_UPDATE': np.int32(0), + 'GRID_FDDA': np.int32(0), + 'GFDDA_INTERVAL_M': np.int32(0), + 'GFDDA_END_H': np.int32(0), + 'GRID_SFDDA': np.int32(0), + 'SGFDDA_INTERVAL_M': np.int32(0), + 'SGFDDA_END_H': np.int32(0), + 'WEST-EAST_PATCH_START_UNSTAG': None, # Projection dependent attributes + 'WEST-EAST_PATCH_END_UNSTAG': None, # Projection dependent attributes + 'WEST-EAST_PATCH_START_STAG': None, # Projection dependent attributes + 'WEST-EAST_PATCH_END_STAG': None, # Projection dependent attributes + 'SOUTH-NORTH_PATCH_START_UNSTAG': None, # Projection dependent attributes + 'SOUTH-NORTH_PATCH_END_UNSTAG': None, # Projection dependent attributes + 'SOUTH-NORTH_PATCH_START_STAG': None, # Projection dependent attributes + 'SOUTH-NORTH_PATCH_END_STAG': None, # Projection dependent attributes + 'BOTTOM-TOP_PATCH_START_UNSTAG': None, + 'BOTTOM-TOP_PATCH_END_UNSTAG': None, + 'BOTTOM-TOP_PATCH_START_STAG': None, + 'BOTTOM-TOP_PATCH_END_STAG': None, + 'GRID_ID': np.int32(1), + 'PARENT_ID': np.int32(0), + 'I_PARENT_START': np.int32(1), + 'J_PARENT_START': np.int32(1), + 'PARENT_GRID_RATIO': np.int32(1), + 'DT': np.float32(18.), + 'CEN_LAT': None, # Projection dependent attributes + 'CEN_LON': None, # Projection dependent attributes + 'TRUELAT1': None, # Projection dependent attributes + 'TRUELAT2': None, # Projection dependent attributes + 'MOAD_CEN_LAT': None, # Projection dependent attributes + 'STAND_LON': None, # Projection dependent attributes + 'POLE_LAT': None, # Projection dependent attributes + 'POLE_LON': None, # Projection dependent attributes + 'GMT': np.float32(self.time[0].hour), + 'JULYR': np.int32(self.time[0].year), + 'JULDAY': np.int32(self.time[0].strftime("%j")), + 'MAP_PROJ': None, # Projection dependent attributes + 'MMINLU': 'MODIFIED_IGBP_MODIS_NOAH', + 'NUM_LAND_CAT': np.int32(41), + 'ISWATER': np.int32(17), + 'ISLAKE': np.int32(-1), + 'ISICE': np.int32(15), + 'ISURBAN': np.int32(13), + 'ISOILWATER': np.int32(14), + 'HISTORY': "", # Editable + } + + # Editable attributes + float_atts = ['DAMPCOEF', 'KHDIF', 'KVDIF', 'CEN_LAT', 'CEN_LON', 'DT'] + int_atts = ['BOTTOM-TOP_GRID_DIMENSION', 'DIFF_OPT', 'KM_OPT', 'DAMP_OPT', + 'MP_PHYSICS', 'RA_LW_PHYSICS', 'RA_SW_PHYSICS', 'SF_SFCLAY_PHYSICS', 'SF_SURFACE_PHYSICS', + 'BL_PBL_PHYSICS', 'CU_PHYSICS', 'SF_LAKE_PHYSICS', 'SURFACE_INPUT_SOURCE', 'SST_UPDATE', + 'GRID_FDDA', 'GFDDA_INTERVAL_M', 'GFDDA_END_H', 'GRID_SFDDA', 'SGFDDA_INTERVAL_M', 'SGFDDA_END_H', + 'BOTTOM-TOP_PATCH_START_UNSTAG', 'BOTTOM-TOP_PATCH_END_UNSTAG', 'BOTTOM-TOP_PATCH_START_STAG', + 'BOTTOM-TOP_PATCH_END_STAG', 'GRID_ID', 'PARENT_ID', 'I_PARENT_START', 'J_PARENT_START', + 'PARENT_GRID_RATIO', 'NUM_LAND_CAT', 'ISWATER', 'ISLAKE', 'ISICE', 'ISURBAN', 'ISOILWATER'] + str_atts = ['GRIDTYPE', 'MMINLU', 'HISTORY'] + for att_name, att_value in current_attributes.items(): + if att_name in int_atts: + self.global_attrs[att_name] = np.int32(att_value) + elif att_name in float_atts: + self.global_attrs[att_name] = np.float32(att_value) + elif att_name in str_atts: + self.global_attrs[att_name] = str(att_value) + + # Projection dependent attributes + if isinstance(self, nes.LCCNes) or isinstance(self, nes.MercatorNes): + self.global_attrs['WEST-EAST_GRID_DIMENSION'] = np.int32(len(self._x['data']) + 1) + self.global_attrs['SOUTH-NORTH_GRID_DIMENSION'] = np.int32(len(self._y['data']) + 1) + self.global_attrs['DX'] = np.float32(self._x['data'][1] - self._x['data'][0]) + self.global_attrs['DY'] = np.float32(self._y['data'][1] - self._y['data'][0]) + self.global_attrs['SURFACE_INPUT_SOURCE'] = np.int32(1) + self.global_attrs['WEST-EAST_PATCH_START_UNSTAG'] = np.int32(1) + self.global_attrs['WEST-EAST_PATCH_END_UNSTAG'] = np.int32(len(self._x['data'])) + self.global_attrs['WEST-EAST_PATCH_START_STAG'] = np.int32(1) + self.global_attrs['WEST-EAST_PATCH_END_STAG'] = np.int32(len(self._x['data']) + 1) + self.global_attrs['SOUTH-NORTH_PATCH_START_UNSTAG'] = np.int32(1) + self.global_attrs['SOUTH-NORTH_PATCH_END_UNSTAG'] = np.int32(len(self._y['data'])) + self.global_attrs['SOUTH-NORTH_PATCH_START_STAG'] = np.int32(1) + self.global_attrs['SOUTH-NORTH_PATCH_END_STAG'] = np.int32(len(self._y['data']) + 1) + + self.global_attrs['POLE_LAT'] = np.float32(90) + self.global_attrs['POLE_LON'] = np.float32(0) + + if isinstance(self, nes.LCCNes): + self.global_attrs['MAP_PROJ'] = np.int32(1) + self.global_attrs['TRUELAT1'] = np.float32(self.projection_data['standard_parallel'][0]) + self.global_attrs['TRUELAT2'] = np.float32(self.projection_data['standard_parallel'][1]) + self.global_attrs['MOAD_CEN_LAT'] = np.float32(self.projection_data['latitude_of_projection_origin']) + self.global_attrs['STAND_LON'] = np.float32(self.projection_data['longitude_of_central_meridian']) + self.global_attrs['CEN_LAT'] = np.float32(self.projection_data['latitude_of_projection_origin']) + self.global_attrs['CEN_LON'] = np.float32(self.projection_data['longitude_of_central_meridian']) + elif isinstance(self, nes.MercatorNes): + self.global_attrs['MAP_PROJ'] = np.int32(3) + self.global_attrs['TRUELAT1'] = np.float32(self.projection_data['standard_parallel']) + self.global_attrs['TRUELAT2'] = np.float32(0) + self.global_attrs['MOAD_CEN_LAT'] = np.float32(self.projection_data['standard_parallel']) + self.global_attrs['STAND_LON'] = np.float32(self.projection_data['longitude_of_projection_origin']) + self.global_attrs['CEN_LAT'] = np.float32(self.projection_data['standard_parallel']) + self.global_attrs['CEN_LON'] = np.float32(self.projection_data['longitude_of_projection_origin']) + + return None + + +def create_dimensions(self, netcdf): + """ + Create 'time', 'time_bnds', 'lev', 'lon' and 'lat' dimensions. + + Parameters + ---------- + self : nes.Nes + netcdf : Dataset + netcdf4-python open dataset. + """ + netcdf.createDimension('Time', len(self._time)) + netcdf.createDimension('DateStrLen', 19) + netcdf.createDimension('emissions_zdim', len(self._lev['data'])) + if isinstance(self, nes.LCCNes): + netcdf.createDimension('west_east', len(self._x['data'])) + netcdf.createDimension('south_north', len(self._y['data'])) + + return None + + +def create_dimension_variables(self, netcdf): + """ + Create the 'y' and 'x' variables. + + Parameters + ---------- + self : nes.Nes + netcdf : Dataset + NetCDF object. + """ + + times = netcdf.createVariable('Times', 'S1', ('Time', 'DateStrLen', )) + times[:] = create_times_var(self) + + return None + + +def create_variables(self, netcdf): + """ + Create the netCDF file variables. + + Parameters + ---------- + self : nes.Nes + netcdf : Dataset + netcdf4-python open dataset. + """ + for var_name, var_info in self.variables.items(): + var = netcdf.createVariable(var_name, 'f', ('Time', 'emissions_zdim', 'south_north', 'west_east',), + zlib=self.zip_lvl > 0, complevel=self.zip_lvl) + var.FieldType = var_info['FieldType'] + var.MemoryOrder = var_info['MemoryOrder'] + var.description = var_info['description'] + var.units = var_info['units'] + var.stagger = var_info['stagger'] + var.coordinates = var_info['coordinates'] + + if var_info['data'] is not None: + if self.info: + print("Rank {0:03d}: Filling {1})".format(self.rank, var_name)) + + if isinstance(var_info['data'], int) and var_info['data'] == 0: + var[self.write_axis_limits['t_min']:self.write_axis_limits['t_max'], + self.write_axis_limits['z_min']:self.write_axis_limits['z_max'], + self.write_axis_limits['y_min']:self.write_axis_limits['y_max'], + self.write_axis_limits['x_min']:self.write_axis_limits['x_max']] = 0 + + elif len(var_info['data'].shape) == 4: + var[self.write_axis_limits['t_min']:self.write_axis_limits['t_max'], + self.write_axis_limits['z_min']:self.write_axis_limits['z_max'], + self.write_axis_limits['y_min']:self.write_axis_limits['y_max'], + self.write_axis_limits['x_min']:self.write_axis_limits['x_max']] = var_info['data'] + + return None -- GitLab