diff --git a/preproc/cams_glob_ship_v31_preproc.py b/preproc/cams_glob_ship_v31_preproc.py new file mode 100755 index 0000000000000000000000000000000000000000..b475aa465690cec1303f21b9fd0171593813360a --- /dev/null +++ b/preproc/cams_glob_ship_v31_preproc.py @@ -0,0 +1,225 @@ +#!/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 os +import sys +import numpy as np + +# ============== README ====================== +""" +downloading website: +reference: +Besides citing HERMESv3_GR, users must also acknowledge the use of the corresponding emission inventories in their works +""" + +# ============== CONFIGURATION PARAMETERS ====================== +INPUT_PATH = '/esarchive/recon/ecmwf/cams_glob_shipv31/original_files/' +OUTPUT_PATH = '/esarchive/recon/ecmwf/cams_glob_shipv31' +POLLUTANT_INFO = {'Ash': {'input_var_name': ['glseas', 'inland'], 'output_var_name': 'ash', 'unit_factor': None}, + 'CH4': {'input_var_name': ['glseas', 'inland'], 'output_var_name': 'ch4', 'unit_factor': None}, + 'CO2': {'input_var_name': ['glseas', 'inland'], 'output_var_name': 'co2', 'unit_factor': None}, + 'CO': {'input_var_name': ['glseas', 'inland'], 'output_var_name': 'co', 'unit_factor': None}, + 'EC': {'input_var_name': ['glseas', 'inland'], 'output_var_name': 'ec', 'unit_factor': None}, + 'NOx': {'input_var_name': ['glseas', 'inland'], 'output_var_name': 'nox_no', 'unit_factor': None}, + 'OC': {'input_var_name': ['glseas', 'inland'], 'output_var_name': 'oc', 'unit_factor': None}, + 'SO4': {'input_var_name': ['glseas', 'inland'], 'output_var_name': 'so4', 'unit_factor': None}, + 'SOx': {'input_var_name': ['glseas', 'inland'], 'output_var_name': 'so2', 'unit_factor': None}, + 'VOC_SUM': {'input_var_name': ['glseas', 'inland'], 'output_var_name': 'nmvoc', 'unit_factor': None}} +# LIST_YEARS = 2019 +LIST_YEARS = [2019] +INPUT_NAME = 'CAMS-GLOB-SHIP_Glb_0.1x0.1_anthro__v3.1_monthly_.nc' +# ============================================================== + + +def get_input_name(filename, pollutant, year): + """ + Gets the path for the input file name + + :param pollutant: Name of the pollutant + :type pollutant: str + + :return: Path to the input file + :rtype: str + """ + + filename = filename.replace('', pollutant) + filename = filename.replace('', str(year)) + + return filename + + +def get_full_year_data(file_name, year, var_name): + """ + Gets the needed date in the input format. + + :param file_name: path to the input file. + :type file_name: str + + :param pollutant: Name of the pollutant. + :type pollutant: str + + :param sector: Name of the sector. + :type sector: str + + :param year: Year to calculate. + :type year: int + + :return: Data of the selected emission. + :rtype: numpy.array + """ + from netCDF4 import Dataset + from datetime import datetime + import cf_units + import numpy as np + + nc = Dataset(file_name, mode='r') + + time = nc.variables['time'] + + time_array = cf_units.num2date(time[:], time.units, time.calendar) + time_array = np.array([datetime(year=x.year, month=x.month, day=1) for x in time_array]) + + i_time = np.where(time_array == datetime(year=year, month=1, day=1))[0][0] + + # if sector in LIST_SECTORS: + if isinstance(var_name, list): + for i, var_name_aux in enumerate(var_name): + if i == 0: + data = nc.variables[var_name_aux][i_time:i_time + 12, :, :] + else: + data += nc.variables[var_name_aux][i_time:i_time + 12, :, :] + else: + data = nc.variables[var_name][i_time:i_time+12, :, :] + # else: + # data = None + nc.close() + + return data + + +def get_global_attributes(file_name): + """ + Gets the global attributes of the input file. + + :param file_name: Path to the NetCDF file + :type file_name: str + + :return: Global attributes + :rtype: dict + """ + from netCDF4 import Dataset + + nc = Dataset(file_name, mode='r') + + atts_dict = {} + for name in nc.ncattrs(): + atts_dict[name] = nc.getncattr(name) + + nc.close() + return atts_dict + + +def create_bounds(coordinates, number_vertices=2): + """ + Calculate the vertices coordinates. + + :param coordinates: Coordinates in degrees (latitude or longitude) + :type coordinates: numpy.array + + :param number_vertices: Non mandatory parameter that informs the number of vertices that must have the boundaries. + (by default 2) + :type number_vertices: int + + :return: Array with as many elements as vertices for each value of coords. + :rtype: numpy.array + """ + interval = coordinates[1] - coordinates[0] + + coords_left = coordinates - interval / 2 + coords_right = coordinates + interval / 2 + if number_vertices == 2: + bound_coords = np.dstack((coords_left, coords_right)) + elif number_vertices == 4: + bound_coords = np.dstack((coords_left, coords_right, coords_right, coords_left)) + else: + raise ValueError('The number of vertices of the boudaries must be 2 or 4') + + return bound_coords + + +def do_transformation(year): + """ + Does the transformation for the selected year + + :param year: Year to calculate + :type year: int + """ + from datetime import datetime + from hermesv3_gr.tools.netcdf_tools import extract_vars, get_grid_area, write_netcdf + for pollutant, pollutant_info in POLLUTANT_INFO.items(): + file_name = get_input_name(os.path.join(INPUT_PATH, INPUT_NAME), pollutant, year) + if os.path.exists(file_name): + c_lats, c_lons = extract_vars(file_name, ['lat', 'lon']) + lat_inc = -0.1 # c_lats['data'][1] - c_lats['data'][0] + c_lats['data'] += lat_inc / 2 + lon_inc = 0.1 # c_lons['data'][1] - c_lons['data'][0] + c_lons['data'] += lon_inc / 2 + + b_lats = create_bounds(c_lats['data']) + b_lons = create_bounds(c_lons['data']) + + # cell_area = get_grid_area(file_name) + + global_attributes = get_global_attributes(file_name) + + data = get_full_year_data(file_name, year, pollutant_info['input_var_name']) + + pollutant_name = pollutant_info['output_var_name'] + + file_path = os.path.join(OUTPUT_PATH, 'monthly_mean', '{0}'.format(pollutant_name)) + if not os.path.exists(file_path): + os.makedirs(file_path) + + for month in range(1, 12 + 1, 1): + emission = { + 'name': pollutant_name, + 'units': 'kg.m-2.s-1', + # 'data': data[month - 1, :, :].reshape((1,) + cell_area.shape) + 'data': data[month - 1, :, :].reshape((1, data.shape[-2], data.shape[-1])) + + } + # write_netcdf( + # os.path.join(file_path, '{0}_{1}{2}.nc'.format(pollutant_name, year, str(month).zfill(2))), + # c_lats['data'], c_lons['data'], [emission], date=datetime(year=year, month=month, day=1), + # boundary_latitudes=b_lats, boundary_longitudes=b_lons, cell_area=cell_area, + # global_attributes=global_attributes) + + write_netcdf( + os.path.join(file_path, '{0}_{1}{2}.nc'.format(pollutant_name, year, str(month).zfill(2))), + c_lats['data'], c_lons['data'], [emission], date=datetime(year=year, month=month, day=1), + boundary_latitudes=b_lats, boundary_longitudes=b_lons, global_attributes=global_attributes) + else: + raise IOError('File not found {0}'.format(file_name)) + return True + + +if __name__ == '__main__': + for y in LIST_YEARS: + do_transformation(y) \ No newline at end of file