Heure locale
Script pour calculer le cycle diurne en fonction de l'heure locale
solar_time.py
—
Python Source,
3 ko (3466 bytes)
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])