Source code for nes.nes_formats.wrf_chem_format

#!/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
[docs] 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
[docs] 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
[docs] 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(overwrite=False) 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
[docs] 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
[docs] 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
[docs] 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
[docs] 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
[docs] 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