# -*- coding: utf-8 -*-
import numpy as np
import math
from snowtools.utils.resources import print_used_memory
# L. Roussel dec. 2024 : modification start slope_aspect_correction when extracting dates
# + write comprehensive commented code
def interp1d(x, y, tab):
"""
Return the linear interpolation of y for the tab values of the x coordinate
Method designed to avoid the use of scipy.interpolate (avoid the dependency to an external module and environment problems on Bull)
.. warning::
This is not relevant anymore
:meta private:
"""
ind_after = np.searchsorted(x, tab)
ind_before = ind_after - 1
x_before = np.take(x, ind_before)
x_after = np.take(x, ind_after)
y_before = np.take(y, ind_before)
y_after = np.take(y, ind_after)
return ((x_after - tab) * y_before + (tab - x_before) * y_after) / ((x_after - x_before) * 1.)
[docs]
class sun():
"""
A class for computation around sun and radiations
"""
printmemory = False
def __init__(self):
self.missingvalue = -9999999.
[docs]
def slope_aspect_correction(self, direct, diffus, time, lat_in, lon_in, aspect_in, slope_in, list_list_azim=None, list_list_mask=None, lnosof_surfex=True, convert_time = True, return_angles = False):
"""
This routine corrects the direct solar radiation because of explicit slope or surrounding masks
:param direct: array of direct solar radiation over an infinite flat surface (time,loc) or (time,x,y)
:param time: time vector
:param latitude: latitude vector
:param logitude: longitude vector
:param aspect: aspect vector
:param slope: slope vector
:param list_list_azim: list of azimuths of masks on each point(radian or degrees ?)
:param list_list_mask: list of heights of masks on each point(radian or degrees ?)
:param lnosof_surfex: (subgrid orography deactivated in Surfex)
:type lnosof_surfex: bool
:param convert_time: take time at the middle of the time step to compute angles more representative of the fluxes
:type convert_time: bool
:param return_angles: if True, also return the sun angle
:type return_angles: bool
:returns: corrected direct solar radiation and optionally, the angular positions of sun
"""
# Docs about solar radiation :
# "Basics in solar radiation at earth surface" (L. Wald) : https://minesparis-psl.hal.science/hal-01676634
# "Wikipedia: solar irradience" : https://en.wikipedia.org/wiki/Solar_irradiance
tab_direct = direct[:]
tab_diffus = diffus[:]
tab_global = tab_direct + tab_diffus
# --------- Upscale tabs ---------
# (all tab will have the same size as input tab_direct)
slope = self.upscale_tab(slope_in, tab_direct.shape)
aspect = self.upscale_tab(aspect_in, tab_direct.shape)
lon = self.upscale_tab(lon_in, tab_direct.shape)
lat = self.upscale_tab(lat_in, tab_direct.shape)
# --------- Extraction of days and time ---------
# extraction of date and computation of fractional Julian day, j_2 (1 january is 0)
if convert_time is True:
# M. Lafaysse : convert date in date at the middle of the time step to compute angles more representative of the fluxes
# L. Roussel 2024 : SAFRAN gives the hourly solar radiation at hour-30 minutes,
# then SURFEX uses this flux between H-1 and H.
deltatime = time[1] - time[0]
tab_time_date = time - deltatime / 2
else:
tab_time_date = time
julian_days = np.ones(tab_time_date.shape, 'f')
decimal_hours = np.ones(tab_time_date.shape, 'f')
for i in range(len(tab_time_date)):
# timetuple list description:
# j[0]: year, j[1]: month, j[2]: day,
# j[3]: hour, j[4]: minute, j[5]: second,
# j[6]: day of the week, j[7]: day of the year,
# j[8]: daylight saving time
time_tuple = tab_time_date[i].timetuple()
julian_days[i] = time_tuple[7] # extract Julian day (integer, 1st january is 1)
# L. Roussel: fix to decimal hour instead of integer hour
decimal_hours[i] = time_tuple[3] + time_tuple[4] / 60 + time_tuple[5] / 3600 # extrac time in hour
j = self.upscale_tab_time(julian_days, tab_direct.shape)
h = self.upscale_tab_time(decimal_hours, tab_direct.shape)
# all tabs now have the same dimension, which is that of tab_direct.
# method from crocus meteo.f90 original file (v2.4)
# --------- Math constants ---------
eps = 0.0001
PI = math.pi
DG2RD = PI / 180. # degree to radian
RD2DG = 180. / PI # radian to degree
# --------- Convert to radian ---------
slope, aspect = slope * DG2RD, aspect * DG2RD
# Later computations are almost all approximations used for engineering
# --------- Equation of Time ---------
# (time shift in minutes, take in account orbit eccentricity, and obliquity)
eot = 9.9 * np.sin(2 * (j - 101.4) / 365 * 2 * PI) \
- 7.7 * np.sin((j - 2.027) / 365 * 2 * PI)
# --------- Sun declination (delta) ---------
# (angular distance of the sun's rays north (or south) of the equator)
# L. Roussel 03/2024, small changes here, but easier to understand:
# 81th days of the year is the spring equinox
# equivalent to : 0.4 * np.sin((0.999694 * j - 81.111) / 365 * 2 * PI) ~ 0.4 * sin((j - 81) / 365 * 2.PI)
sin_delta = 0.4 * np.sin((0.986 * j - 80) * DG2RD)
delta = np.arcsin(sin_delta)
# --------- Solar angular time ---------
# Solar angle at time h of the day
# This do not take in account the shift from year to year (bissextil year and so on)
# h - 12: centered at noontime
# eot / 60: equation of time in hour
# lon / 15: position from Greenwhich (french Alps around 7° East, i.e ~30 min shift)
# / 24: in days
# * 2 * PI: to angle
omega_time = (h - 12 + eot / 60 + lon / 15) / 24 * 2 * PI
# --------- Solar angular height ---------
# cos(theta) = sin(gamma)
#
# zenith, vertical
# |
# | solar zenith angle, sza (theta)
# |---->x
# | x^
# | x |
# | x | solar altitude angle (gamma)
# | x | ~ solar angular height
# |x |
# +--------- horizon
#
lat_radian = lat * DG2RD
sin_gamma = np.sin(lat_radian) * sin_delta + np.cos(lat_radian) * np.cos(delta) * np.cos(omega_time)
# L. Roussel: set as 0 if negative (before 0.001)
sin_gamma = np.where(sin_gamma < eps, 0., sin_gamma)
gamma = np.arcsin(sin_gamma)
# M Lafaysse : remove this threshold because there are controls in SAFRAN and because there are numerical issues at sunset.
# --------- Theoretical maximum radiation ---------
max_theor_radiation = 1370 * (1 - sin_delta / 11.7)
# F. Besson/M. Lafaysse : Threshold on ZSINPS (bug on some cells of SIM-FRANCE)
# solar zenith angle
ZSINPS = np.cos(delta) * np.sin(omega_time) / np.cos(gamma)
ZSINPS = np.where(ZSINPS < -1.0, -1.0, ZSINPS)
ZSINPS = np.where(ZSINPS > 1.0, 1.0, ZSINPS)
# Not used : ML comment this instruction
# ZCOSPS=(np.sin(ZLAT)*ZSINGA-ZSINDL)/(np.cos(ZLAT)*np.cos(ZGAMMA))
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# Before summer 2017, we had this mistake :
# # # # ZPSI = np.arcsin(ZSINPS) # solar azimuth, 0. is South # # # # # NEVER USE THIS WRONG FORMULA
# This gives a wrong trajectory in summer when the azimuth should be <90° or >270°
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# The new computation of azimuth below is extracted from the theoricRadiation method
# Who knows where it comes from ?
# --------- Azimuth ---------
# L. Roussel 2024: precedent computation was fine, just rename variables
# component x of azimuth angle:
x_azimuth_angle = - np.sin(omega_time) * np.cos(delta)
# component y of azimuth angle:
y_azimuth_angle = np.sin(lat_radian) * np.cos(omega_time) * np.cos(delta) - np.cos(lat_radian) * np.sin(delta)
azimuth = np.where((x_azimuth_angle == 0) & (y_azimuth_angle == 0), 0.,
(np.pi) - np.arctan(x_azimuth_angle / y_azimuth_angle))
azimuth = np.where(y_azimuth_angle <= 0., azimuth + np.pi, azimuth)
azimuth = np.where(azimuth >= 2. * np.pi, azimuth - 2. * np.pi, azimuth) # back into [0, 2 * PI]
##########################################
######## Recompute diffuse/direct ########
# diffuse/global theorical ratio for clear sky (SBDART modelling by M. Dumont for Col de Porte)
a = -2.243613
b = 5.199838
c = -4.472389
d = -0.276815
mu = sin_gamma.copy()
ratio_clearsky = np.exp(a * (mu ** 3) + b * (mu ** 2) + c * mu + d)
ratio_clearsky = np.where(ratio_clearsky <= 1, ratio_clearsky, 1)
# Note that it is not sufficient to apply the threshold on the ratio, it is necessary to apply a threshold to the direct radiation
# Compute theorical maximum global radiation (clear sky)
a3 = -0.25070845442
a2 = 0.56643807637
a1 = 0.620060537777
a0 = -0.025767794921
ZTHEOR = max_theor_radiation * (a3 * mu**3 + a2 * mu**2 + a1 * mu + a0)
ZTHEOR = np.where(ZTHEOR < 0., 0., ZTHEOR)
# Compute theorical components
theor_diffus = ratio_clearsky * ZTHEOR
theor_direct = ZTHEOR - theor_diffus
# Apply a threshold to the direct radiation (can not exceed the theorical direct component, otherwise, the projection can provide very strange values
tab_direct = np.where(tab_direct <= theor_direct, tab_direct, theor_direct)
# Conserve energy by a transfer to the diffuse component
tab_diffus = tab_global - tab_direct
# direct incident radiation (maximum value authorized is the theoretical maximum radiation)
direct_incident = np.divide(tab_direct, sin_gamma, out=np.zeros_like(tab_direct), where=sin_gamma!=0)
# M Lafaysse : remove this threshold because there are controls in SAFRAN and because there are numerical issues at sunset.
# L. Roussel 2024: this one looks useless, still keep it to control anyway
direct_incident = np.where(direct_incident >= max_theor_radiation, max_theor_radiation, direct_incident)
# --------- Projection on slope/aspect surface ---------
# (note that aspect is taken from North)
# old code 1: ZRSIP=ZRSI*(np.cos(ZGAMMA)*np.sin(slope)*np.cos(ZPSI - (np.pi + aspect))+ZSINGA*np.cos(slope))
# old code 2: ZRSIP = ZRSI * (np.cos(ZGAMMA) * np.sin(slope) * np.cos(azimuth - aspect) + ZSINGA * np.cos(slope))
direct_plane_projected = direct_incident * (np.cos(gamma) * np.sin(slope) * np.cos(azimuth - aspect) + sin_gamma * np.cos(slope))
direct_plane_projected = np.where(direct_plane_projected <= 0., 0., direct_plane_projected)
# --------- Solar masking ---------
# if mask are available, mask direct when relief between sun and location
# (S. Morin 2014/06/27, taken from meteo.f90 in Crocus)
# Matthieu 2014/09/16 : réécriture parce que les masques changent d'un point à l'autre
direct_plane_projected_masked = direct_plane_projected.copy()
if list_list_mask is not None: # if mask are given
simu_azimuth_degree = azimuth * RD2DG # solar azimuths of the simulation, 0. is North
simu_interp_mask = np.zeros_like(simu_azimuth_degree) # init for computed solar mask (time, loc) or (time, x, y)
for i, list_azim in enumerate(list_list_azim): # for each point
# interp1d(x, y, interp_x)
# return interp_y with x, y and interp_x as arguments
simu_interp_mask[:, i] = interp1d(list_azim, list_list_mask[i], simu_azimuth_degree[:, i])
# set to zero direct radiation values when solar angle is below mask angle
direct_plane_projected_masked = np.where(simu_interp_mask > gamma * RD2DG, 0., direct_plane_projected)
if lnosof_surfex:
# Now this is the normal case
tab_direct = direct_plane_projected_masked
else:
# Not recommended
# put the result back on the horizontal ; surfex will carry out the inverse operation when lnosof=f.
tab_direct = direct_plane_projected_masked / np.cos(slope)
if self.printmemory:
print_used_memory()
if return_angles:
return tab_direct, tab_diffus, gamma, azimuth
else:
return tab_direct, tab_diffus
# The two following routines from JM Willemet should be rewritten without the loops
def upscale_tab(self, var, theshape):
# array var can be considered at the same shape than the other arrays (simplify matricial traitment)
bigvar = np.ma.zeros(theshape) + self.missingvalue
shapein = len(var.shape)
if len(theshape) == 2:
for i in range(0, var.shape[0]):
bigvar[:, i] = var[i]
elif len(theshape) == 3:
if shapein == 2:
for ilat in range(0, var.shape[0]):
for ilon in range(0, var.shape[1]):
bigvar[:, ilat, ilon] = var[ilat, ilon]
elif shapein == 1:
for i in range(0, var.shape[0]):
bigvar[:, :, i] = var[i]
elif len(theshape) == 4:
for i in range(0, var.shape[0]):
bigvar[:, :, :, i] = var[i]
elif len(theshape) == 5:
for ilat in range(0, var.shape[0]):
for ilon in range(0, var.shape[1]):
bigvar[:, :, :, ilat, ilon] = var[ilat, ilon]
else:
print("error on indices in upscale_tab")
return bigvar
def upscale_tab_time(self, var, theshape):
# array slope can be considered at the same shape than the other arrays (simplify matricial traitment)
bigvar = np.ma.zeros(theshape)
if len(theshape) == 2:
for i in range(0, theshape[1]):
bigvar[:, i] = var[:]
elif len(theshape) == 3:
for ilat in range(0, theshape[1]):
for ilon in range(0, theshape[2]):
bigvar[:, ilat, ilon] = var[:]
else:
print("error on indices in upscale_tab_time")
return bigvar
[docs]
def coszenith(self, tab_time_date, lat, lon, slope, aspect):
"""
Cosinus of solar zenith angle
(detailed explaination in function slope_aspect_correction)
Take a certain amount of times (Nt) and points (Np described by lat, lon, slope, aspect)
and return for each time and point the cosinus of zenith angle of sun (shape Nt, Np).
:param tab_time_date: time values, numpy array of size Nt of python datetime.datetime elements
(or other type with timetuple function available)
:param lat: latitude, numpy array (floats), length Np
:param lon: longitude, numpy array (floats), length Np
:param slope: slope (degrees), numpy array (floats), length Np
:param aspect: aspect (degrees), numpy array (floats), length Np
:returns: Cosinus of zenith angle
:rtype: numpy array of shape (Nt, Np)
"""
julian_days = np.ones(tab_time_date.shape, 'f')
decimal_hours = np.ones(tab_time_date.shape, 'f')
for i in range(len(tab_time_date)):
# tab_time_date_2[i] = time_init + datetime.timedelta(seconds = tab_time_date[i])
timetup = tab_time_date[i].timetuple()
julian_days[i] = timetup[7] # extract Julian day (integer)
# L. Roussel: fix to decimal hour instead of integer hour
decimal_hours[i] = timetup[3] + timetup[4] / 60 + timetup[5] / 3600 # extrac time in hour
j = self.upscale_tab_time(julian_days, (tab_time_date.shape[0], 1))
h = self.upscale_tab_time(decimal_hours, (tab_time_date.shape[0], 1))
# all tabs now have the same dimension, which is that of tab_direct.
# method from crocus meteo.f90 original file (v2.4)
# Variables needed for solar computations:
PI = math.pi
DG2RD = PI / 180.
eps = 0.001
# --------- Convert to radian ---------
slope, aspect = slope * DG2RD, aspect * DG2RD
# --------- Equation of Time ---------
eot = 9.9 * np.sin(2 * (j - 101.4) / 365 * 2 * PI) \
- 7.7 * np.sin((j - 2.027) / 365 * 2 * PI)
# --------- Sun declination (delta) ---------
sin_delta = 0.4 * np.sin((0.986 * j - 80) * DG2RD)
delta = np.arcsin(sin_delta)
# --------- Solar angular time ---------
omega_time = (h - 12 + eot / 60 + lon / 15) / 24 * 2 * PI
# --------- Solar angular height ---------
# sza (solar zenith angle) = PI / 2 - gamma
# cos(sza) = sin(gamma)
lat_radian = lat * DG2RD
sin_gamma = np.sin(lat_radian) * sin_delta + np.cos(lat_radian) * np.cos(delta) * np.cos(omega_time)
# L. Roussel: set as 0 if negative (before 0.001)
sin_gamma = np.where(sin_gamma < eps, 0, sin_gamma)
return sin_gamma
[docs]
def directdiffus(self, SWglo, time, lat, lon, slope, aspect, site):
"""
Separation of direct and diffuse short wave radiations
:param SWglo: global short wave
:param time: time
:param lat: latitude
:param lon: longitude
:param slope: The slope angle (degrees)
:param aspect: The aspect (degrees)
:param site: The site code for clear sky decomposition. Must be one of wfj, snb, swa, sap, sod, rme, oas, obs, ojp or cdp)
:returns: direct short wave, diffuse short wave
"""
# Clear sky decomposition by Marie Dumont
if site == "wfj":
a = -2.449042
b = 5.639542
c = -4.745762
d = -0.555133
elif site == "snb":
a = -2.492804
b = 5.743704
c = -4.870253
d = -0.464349
elif site == "swa":
a = -2.481555
b = 5.717005
c = -4.838604
d = -0.484631
elif site == "sap":
a = -1.444924
b = 3.506267
c = -3.342956
d = 0.050087
elif site == "sod":
a = -2.171899
b = 5.040324
c = -4.335702
d = -0.256246
elif site == "rme":
a = -2.375236
b = 5.481626
c = -4.653772
d = -0.421807
elif site == "oas":
a = -2.252006
b = 5.208891
c = -4.443042
d = -0.339813
elif site == "obs":
a = -2.254170
b = 5.213592
c = -4.446855
d = -0.340464
elif site == "ojp":
a = -2.250787
b = 5.205990
c = -4.440420
d = -0.339316
elif site == "cdp":
a = -2.243613
b = 5.199838
c = -4.472389
d = -0.276815
coszenith = self.coszenith(time, lat, lon, slope, aspect)[:, 0]
ratio = np.exp(a * (coszenith**3) + b * (coszenith**2) + c * coszenith + d)
SWdif = np.where(ratio <= 1, ratio * SWglo, SWglo)
SWdir = SWglo - SWdif
# for hour in range(1, 24):
# print hour, coszenith[ndays * 24 + hour - 1], ratio[ndays * 24 + hour - 1], SWdir[ndays * 24 + hour - 1], SWdif[ndays * 24 + hour - 1], SWglo[ndays * 24 + hour - 1]
return SWdir, SWdif