Solar zenith and azimuth angles calculation

[BaptisteVandecrux]

Turns out that there are many different algorithms to calculate them.

JAWS is using sun-position but it is not vectorized and therefore very slow for large files.

I also tried to play with suncalc-py which is very fast, but it does not take the site elevation as input and give solar elevation as output instead of zenith angle.

Eventually, I found solar position which work on vectors. It produces similar zenith angles as JAWS but the azimuth does not seem to pass by 0 nor 360. That is surprising because in polar summer the sun should go through all azimuth angles.

Right now I implement solar position and future investigation will be needed. Feedbacks are welcome!

[BaptisteVandecrux]

Not really matching the one calculated by JAWS

[BaptisteVandecrux]

Koni's code for SZA:
Zenith = 90-SolarElevation(Dataset(i,3,4),Latitude,Longitude); % Gives the Zenith Angle (DEGREES)

function Alpha = SolarElevation(JulianDay,Latitude,Longitude)

% -------------------------------------------------------------------------
% --- Returns the Solar Elevation Angle as a function of the follwing inputs:
% ---
% ---   JulianDay   Decimal Julian Day
% ---               Greenwich Mean Time
% ---   Latitude    Degree Latitude of Grid Point
% ---               Counted positive Northern Hemisphere
% ---   Longitude   Degree Longitude of Grid Point
% ---               Counted positive West of Greenwich
% ---
% ---   Output:     ZenithAngle Zenith Angle, in degrees
% ---
% ---   Reference:  Woolf, H.M. NASA TM X-1646
% ---               On the computational of solar elevation
% ---               angles and the determination of 
% ---               sunrise and sunset times.
% ---
% --- http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19680025707_1968025707.pdf
% ---
% --- Version 2.10: February 15th, 2011 - N. Bayou, Boulder
% ---               * Original Radian/Degrees error corrected...
% --- Version 2.20: No changes
% --- Version 2.21: No changes
% -------------------------------------------------------------------------

% --- SIN(Alpha) = SIN(Phi)*SIN(D)+COS(Phi)*COS(D)*COS(h)

% --- Alpha:    Zenith Angle
% --- Phi:      Latitude

DpY = 365.242; % Number of Days per year

D0 = 23+26/60+37.8/3600; % D0 = 23deg26'37.8"

Hour = (JulianDay - floor(JulianDay))*24; % Hour of the day

d = 360*(JulianDay-1)/DpY;      % Measured in DEGREES
l = 279.9348 + d;               % Measured in DEGREES

Sigma   = l + 0.4087*sind(l) + 1.8724*cosd(l) - 0.0182*sind(2*l) + 0.0083*cosd(2*l);

SinD  = sind(D0)*sind(Sigma);

M = 12 + 0.123570*sind(d) - 0.004289*cosd(d) + 0.153809*sind(2*d) + 0.060783*cosd(2*d);
h = 15*(Hour-M)-Longitude;

SinAlpha = sind(Latitude)*SinD+cosd(Latitude)*cosd(asind(SinD))*cosd(h);
Alpha = asind(SinAlpha);

[BaptisteVandecrux]

Ken's code:

    # Calculating zenith and hour angle of the sun
    doy = ds['time'].to_dataframe().index.dayofyear.values
    hour = ds['time'].to_dataframe().index.hour.values
    minute = ds['time'].to_dataframe().index.minute.values
    # lat = ds['gps_lat']
    # lon = ds['gps_lon']
    lat = ds.attrs['latitude']
    lon = ds.attrs['longitude']
    
    d0_rad = 2 * np.pi * (doy + (hour + minute / 60) / 24 -1) / 365
    
    Declination_rad = np.arcsin(0.006918 - 0.399912
                                * np.cos(d0_rad) + 0.070257
                                * np.sin(d0_rad) - 0.006758
                                * np.cos(2 * d0_rad) + 0.000907
                                * np.sin(2 * d0_rad) - 0.002697
                                * np.cos(3 * d0_rad) + 0.00148
                                * np.sin(3 * d0_rad))
    
    HourAngle_rad = 2 * np.pi * (((hour + minute / 60) / 24 - 0.5) - lon/360)
    # ; - 15.*timezone/360.) ; NB: Make sure time is in UTC and longitude is positive when west! Hour angle should be 0 at noon.
    
    # This is 180 deg at noon (NH), as opposed to HourAngle.
    DirectionSun_deg = HourAngle_rad * 180/np.pi - 180
    
    DirectionSun_deg[DirectionSun_deg < 0] += 360
    DirectionSun_deg[DirectionSun_deg < 0] += 360
    
    ZenithAngle_rad = np.arccos(np.cos(lat * deg2rad)
                                * np.cos(Declination_rad)
                                * np.cos(HourAngle_rad)
                                + np.sin(lat * deg2rad)
                                * np.sin(Declination_rad))
    
    ZenithAngle_deg = ZenithAngle_rad * rad2deg

[BaptisteVandecrux]

Using Ken's code give good result: within +/- 1 degree from JAWS slow output: