solar_time.py

Contenu du fichier

#!/usr/bin/env python3

""""This script reads a NetCDF file containing the diurnal cycle as a
function of Greenwich time and creates a NetCDF file containing the
diurnal cycle as a function of solar time.

 Author: Lionel GUEZ

"""

import netCDF4
import numpy as np
import shutil
from scipy import interpolate
import sys
from numpy import ma

if len(sys.argv) != 3:
    sys.exit("Required arguments: input-netCDF-file output-netCDF-file")
    
shutil.copyfile(sys.argv[1], sys.argv[2])

with netCDF4.Dataset(sys.argv[1]) as f_in, \
     netCDF4.Dataset(sys.argv[2], "r+") as f_out:
    f_out.renameDimension("time_counter", "solar_time")
    f_out.renameVariable("time_counter", "solar_time")
    f_out.variables["solar_time"].units = "hours"
    

    for name in ["bounds", "calendar", "cell_methods"]:
        if name in f_out.variables["solar_time"].ncattrs():
            f_out.variables["solar_time"].delncattr(name)

    t_solar = np.arange(1800, 86400, 3600) # in s
    f_out.variables["solar_time"][:] = t_solar / 3600

    if isinstance(f_in.variables["lon"][:], ma.MaskedArray):
        lon = f_in.variables["lon"][:].filled()
        # So that t_UTC is not a masked array. interpolate.interp1d
        # does not accept a masked array.
    else:
        lon = f_in.variables["lon"][:]    

    t_UTC = t_solar[:, np.newaxis] - lon * 240

    # Add a multiple of 86400 s to t_UTC to get into the interval
    # [f_in.variables["time_counter"][0] - 1800;
    # f_in.variables["time_counter"][23] + 1800]:
    m = np.ceil((f_in.variables["time_counter"][0] - t_UTC) / 86400 - 1 / 48)
    t_UTC += m * 86400

    # Check that the reference date of NetCDF variable time_counter is
    # at 00:00. Then t_UTC can be used as a time_counter value.
    d = netCDF4.num2date(0, f_in["time_counter"].units,
                         f_in["time_counter"].calendar)
    assert np.all(np.array((d.hour, d.minute, d.second, d.microsecond)) == 0)
    
    # Since t_UTC can be half an hour before
    # f_in.variables["time_counter"][0] or after
    # f_in.variables["time_counter"][23], extend the time series:
    x = np.block([f_in.variables["time_counter"][0] - 3600,
                  f_in.variables["time_counter"],
                  f_in.variables["time_counter"][23] + 3600])

    for name in f_out.variables:
        if set(('solar_time', 'lat', 'lon')).issubset(f_out.variables[name].
                                                      dimensions):
            print("solar_time.py: name =", name)
       
            for i in range(f_in.dimensions["lon"].size):
                # Extend f_in.variables[name][..., i] in time taking
                # into account periodicity. "y" values correspond to
                # "x" values. We assume that time is the first
                # dimension, longitude is the last dimension and
                # create "y" as having all dimensions except
                # longitude, and an extended time dimension:
                my_shape = f_in.variables[name].shape
                y = np.empty((my_shape[0] + 2, *my_shape[1:-1]))
                y[0] = f_in.variables[name][23, ..., i]
                y[1:- 1] = f_in.variables[name][..., i]
                y[-1] = f_in.variables[name][0, ..., i]

                interpol = interpolate.interp1d(x, y, axis = 0, copy = False,
                                                assume_sorted = True)
                f_out.variables[name][..., i] = interpol(t_UTC[:, i])