 |
0.1.39
|
|
0.1.39
|
| Dependencies | None |
|---|---|
| Python Import | from pyhelios import SolarPosition |
| Main Class | SolarPosition |
| Dependencies | None |
|---|---|
| Platforms | Windows, Linux, macOS |
| GPU | Not required |
This plugin calculates the position of the sun, and also implements other models for solar fluxes as well as longwave fluxes from the sky. Model theory and equations are given in the sections below.
| Constructors |
|---|
| SolarPosition(context) |
| SolarPosition(context, utc_offset, latitude, longitude) |
The SolarPosition class can be initialized by simply passing a Helios context as an argument to the constructor. This gives the class access to the time and date currently set in the Context. The model must also know certain parameters about the simulated location, in particular the offset from UTC time, latitude, and longitude. A description of these parameters are given in the table below. These can be supplied using the second constructor listed in the table above. If only the Context is supplied to the constructor, the plugin uses the Context's currently configured location.
The Context's location can be set ahead of time with context.setLocation(latitude, longitude, utc_offset) or with the immutable Location dataclass from pyhelios.types. See Geographic Location in the user guide for details. With this approach, location is configured once and shared across plugins:
| Input Parameter | Description | Convention | Default Behavior |
|---|---|---|---|
| UTC | Difference in hours between Coordinated Universal Time (UTC) for a particular location. See the figure below to determine a particular UTC offset. | UTC offset value is positive moving West. | +8:00 |
| latitude | Geographic coordinate that specifies the north–south position of a point on the Earth's surface in degrees. | Latitude is positive in the northern hemisphere. | +38.55 |
| longitude | Geographic coordinate that specifies the east-west position of a point on the Earth's surface in degrees. | Longitude is positive in the western hemisphere. | +121.76 |
The solar position model was implemented following the description in Chapter 1 of Iqbal (1983).
The day angle \(\Gamma\) given as the polar angle of the earth relative to the sun ( \(\Gamma=0\) on Jan. 1) is calculated as
where DOY is the Julian Day of the year.
The solar declination angle is then calculated as
The equation of time is calculated as
The hour angle is given by
with
and
Finally, the solar elevation angle is given by
and the solar azimuthal angle is given by
Note that \(\mathrm{cos}^{-1}\) gives angles between 0 and \(\pi\), so to get a \(\phi_s\) between 0 and \(2\pi\), we take \(\phi_s=2\pi-\phi_s\) if \(LST>12\).
Clear-sky solar fluxes are calculated using the 'REST-2' (Reference Evaluation of Solar Transmittance, 2 bands) model of Gueymard (2008). REST-2 is a high-performance broadband radiative transfer model derived from parameterizations of the SMARTS spectral code, and is widely recognized as one of the most accurate clear-sky models available.
The model uses a two-band spectral scheme that separately treats the visible band (290-700 nm) and near-infrared band (700-4000 nm). In the visible band, attenuation is dominated by Rayleigh scattering and aerosol extinction, while the NIR band is primarily affected by water vapor absorption. For each band, the model calculates independent broadband transmittances for:
Direct beam irradiance is computed from the product of these transmittances, while diffuse irradiance uses a two-layer scattering scheme that accounts for aerosol forward scattering and backscattering from the atmosphere-ground system. The model partitions the total radiative flux into direct and diffuse components suitable for agricultural, solar energy, and climate applications.
For applications requiring high spectral resolution (e.g., photosynthesis, remote sensing), calculateDirectSolarSpectrum(), calculateDiffuseSolarSpectrum() and calculateGlobalSolarSpectrum() compute the direct normal, diffuse horizontal and global horizontal spectral irradiance at the ground from 300 to 2600 nm at 1 nm resolution. The model follows the structure of SSolar-GOA (Cachorro et al. 2022): the atmosphere is a single homogeneous layer of molecules and aerosol whose scattering transmittance is given by an analytical two-stream solution, gas absorption multiplies it, and multiple reflection between the ground and the atmosphere is included through the atmosphere's spherical albedo.
Parameter Derivation: The model uses the same atmospheric inputs as the REST-2 model (pressure, temperature, humidity, turbidity). The precipitable water is derived from the air temperature and relative humidity, which must therefore be greater than zero. The total column ozone comes from getOzoneColumn() (see Total Column Ozone) and the ground albedo from setGroundAlbedo() (default 0.2).
Output Format: Results are stored in Context global data as vectors of (wavelength, irradiance) pairs with user-defined labels. Three spectral components are computed: global irradiance on horizontal surface, direct irradiance normal to sun direction, and diffuse irradiance on horizontal surface (all in W/m²/nm).
The gas transmittance tables cover a water vapor path (precipitable water times air mass) up to 320 g/cm², an ozone path (ozone column times air mass) up to 25 atm-cm, and an air mass times the ratio of surface pressure to 1013.25 hPa up to 37.99. A longer path raises an error; with the sun close to the horizon this happens for precipitable water above about 8.4 cm, an ozone column above about 660 DU, or a surface pressure above about 1015 hPa. A sun at or below the horizon also raises an error.
To simulate imagery from a sensor above the atmosphere (e.g., a satellite), the radiance leaving the scene must be converted to the radiance reaching the sensor. The atmosphere attenuates the scene radiance and adds light of its own. Following the formulation of the 6S radiative transfer code (Vermote et al. 1997), the radiance at the sensor is
where \(L_{path}\) is the path radiance scattered into the sensor by the atmosphere without reaching the ground, \(T_{dir}^{\uparrow}\) and \(T_{dif}^{\uparrow}\) are the direct and diffuse transmittances from the ground to the sensor, and the last term (the adjacency effect) is light reflected by the surrounding ground of albedo \(\rho_g\), lit by the global irradiance \(E_g\), and scattered into the sensor's line of sight.
The scattering quantities are taken from a look-up table computed with the vector version of 6S for its continental aerosol model. The table spans surface pressures from 701.2 to 1013 hPa, aerosol optical depths at 550 nm from 0 to 1.2, solar zenith angles up to 70 degrees, view zenith angles up to 60 degrees and all relative azimuths. Gas transmittance, also from 6S, is tabulated separately at 1 nm spacing. In a comparison with direct 6S calculations, top-of-atmosphere reflectances differed by 0.0009 on average and by at most 0.011.
Limitations: a single (continental) aerosol type; the sensor is above the atmosphere; the adjacency effect assumes surroundings of uniform albedo; and clear sky.
In a thermal infrared band the atmosphere absorbs part of the radiance leaving the scene and emits radiance of its own, so the spectral radiance at a sensor above the atmosphere is \(L_{sensor} = \tau L_{surface} + L^{\uparrow}\), where \(\tau\) is the atmospheric transmittance along the view path and \(L^{\uparrow}\) the radiance emitted by the atmosphere along it. The atmosphere also emits the sky radiance that the scene reflects; it is described by \(L^{\downarrow}\), the downwelling irradiance at the ground divided by \(\pi\), so that a Lambertian surface of emissivity \(\varepsilon\) reflects \((1-\varepsilon)L^{\downarrow}\).
The three quantities are interpolated from a look-up table computed with the libRadtran radiative transfer code (REPTRAN absorption parameterization) for a clear sky without aerosol, as averages over 5 cm⁻¹ wavenumber bins stored at the bin centers, 5502-15326 nm. The temperature and humidity profiles are not known from the surface conditions alone, so the table assumes climatological ones, interpolated between standard atmospheres according to the surface air temperature.
The table's axes are the column water vapor (0.1-7.5 g/cm²), the surface air temperature (at sea level, 246.2-321.7 K; the range is lower at higher elevations, e.g. 242.4-304.6 K at 701.2 hPa), the surface pressure (701.2-1013 hPa, extrapolated linearly up to 1050 hPa), the view zenith angle (0-60 degrees) and the ozone column (150-450 DU). The table agrees with the line-by-line model LBLRTM to within 1 K in brightness temperature for the satellite bands tested.
The broadband longwave radiation flux emitted by a clear sky (getAmbientLongwaveFlux()) is modeled following Prata (1996).
The model surmounts to calculating the effective emissivity of the sky as a function of precipitable water in the atmosphere
where \(u\) is the precipitable water in cm. It is estimated from the air temperature \(T_a\) (K) and relative humidity \(h\) set with setAtmosphericConditions() as \(u = 0.465\,e_0/T_a\), with the vapor pressure \(e_0 = 611\,h\,\mathrm{exp}\left(17.502\,(T_a-273)/(T_a-273+240.9)\right)\) in Pa. (The precipitable water used by the spectral solar models and the thermal atmosphere look-up table is instead derived from the dew point.)
The downwelling longwave radiation flux on a horizontal surface is given by
where \(\sigma=5.67\times10^{-8}\) W/m2-K4, and \(T_a\) is the air temperature in Kelvin measured near the ground (say 2 m height).
This flux covers all longwave wavelengths. It is not the right diffuse flux for an emission band with wavelength bounds, which should receive only the sky radiation within the band: getThermalSkyFlux() computes that from the thermal atmosphere look-up table (see Thermal Infrared Atmosphere Between the Scene and a Satellite Sensor), for bands within 5502-15326 nm. The table cannot replace the broadband flux, since much of the sky's longwave emission lies outside that range (e.g., the water vapor rotation band beyond 15 µm).
The direction of the sun can be queried in one of several ways: a Cartesian unit vector pointing in the direction of the sun, a spherical coordinate describing the direction of the sun, the elevation angle of the sun, the zenithal angle of the sun, and the azimuthal angle of the sun. The functions to query these quantities are given in the table below. Each of these functions calculates the solar direction based on the current time and date set in the Context (see setTime() "setTime()" and setDate() "setDate()"), and the UTC, latitude, and longitude specified in the SolarPosition constructor.
| Direction Quantity | Function |
|---|---|
| Unit vector pointing toward the sun. | getSunDirectionVector() |
| Spherical coordinate vector pointing toward the sun. | getSunDirectionSpherical() |
| Elevation angle of the sun (radians). | getSunElevation() |
| Zenithal angle of the sun (radians). | getSunZenith() |
| Azimuthal angle of the sun (radians). | getSunAzimuth() |
Below is an example of how to use the SolarPosition plugin to calculate the sun angle.
The SolarPosition plugin requires atmospheric parameters to calculate solar flux and related quantities. The Python API uses setAtmosphericConditions() to set atmospheric parameters once, then calls parameter-free flux methods.
The turbidity parameter used in the SolarPosition plugin is Ångström's aerosol turbidity coefficient (β), the aerosol optical depth at 1 µm. It quantifies the amount of aerosols (dust, pollution, haze) in the atmosphere that scatter and absorb solar radiation.
Important: This turbidity definition is NOT the same as "Linke turbidity" (TL), which is commonly used in some other solar radiation models and uses a different scale.
The aerosol optical depth at wavelength λ follows the Ångström turbidity formula:
\[ \tau_{aerosol}(\lambda) = \beta \lambda^{-\alpha} \]
where λ is in µm and α is the Ångström exponent, taken as 1.3 by the REST-2 and spectral models and the sensor atmosphere model. The aerosol optical depth at 550 nm is therefore about 2.18β, and that at 500 nm about 2.46β; an aerosol optical depth measured at 500 nm (e.g., by a sun photometer) corresponds to β ≈ 0.41 times that value. Passing an aerosol optical depth at 500 nm directly as the turbidity over-states the aerosol load by about 2.5 times. The default turbidity is 0.02.
The Prague sky model (updatePragueSkyModel()) is parameterized by ground-level visibility rather than aerosol optical depth. The aerosol optical depth at 550 nm is converted to visibility by log-log interpolation through the model's continental polluted, average and clean atmospheres, whose visibilities are 27.6, 59.4 and 131.8 km. The visibility is limited to the model's range of 20-131.8 km, so air cleaner than an aerosol optical depth at 550 nm of 0.064 (turbidity about 0.029, which includes the default of 0.02) is represented by the cleanest Prague atmosphere, and air hazier than about 0.45 (turbidity about 0.21) by the haziest. The visibility used is stored in Context global data as prague_sky_visibility_km.
Higher turbidity values result in:
| Parameter | Description | Validation | Example Value |
|---|---|---|---|
| pressure_Pa | Atmospheric pressure in Pascals (near the ground) | Must be > 0 | 101,325 Pa (1 atm) |
| temperature_K | Air temperature in Kelvin (near the ground) | Must be > 0 | 300 K (27°C) |
| humidity_rel | Air relative humidity (near the ground) | Must be 0-1; must be greater than 0 for the spectral irradiance, sensor atmosphere and thermal atmosphere methods | 0.5 (50%) |
| turbidity | Ångström's aerosol turbidity coefficient (β), the aerosol optical depth at 1 µm (see Understanding the Turbidity Parameter). Note: This is NOT Linke turbidity, which uses a different scale. Higher values indicate more aerosols in the atmosphere, which reduces direct solar flux and increases diffuse fraction. | Must be ≥ 0 | 0.02 (default) |
The REST-2 and spectral irradiance models and the sensor atmosphere models need the total column ozone, which absorbs in the ultraviolet and in the Chappuis band in the visible. It can be set, in Dobson units, with setOzoneColumn() (stored in Context global data as atmosphere_ozone_DU). Otherwise it is taken from a monthly zonal-mean climatology: the monthly mean total ozone of the NASA SBUV Merged Ozone Data Set in 36 latitude bands 5 degrees wide, averaged over 2005-2024. The value at the current date and latitude is interpolated linearly in latitude between the centers of the bands and in the day of the year between the middles of the months. It can be retrieved with getOzoneColumn().
The climatology ranges from about 180 DU (the Antarctic ozone hole, 75-80 degrees south in October) to about 420 DU (the Arctic in spring). It has no data where the sun is too low for the satellite measurements, in high-latitude winter months, nor poleward of 80 degrees; for those dates and latitudes, any method that needs the ozone column (including getSolarFlux()) raises an error, and the ozone column must be set with setOzoneColumn(). Next to a latitude band or month without data, the value of the band or month containing the latitude or date is used rather than interpolating toward the missing one.
The solar flux can be calculated using the REST-2 model of Gueymard (2008) using the getSolarFlux() function. IT IS CRITICAL TO NOTE THAT THE CALCULATED FLUX IS FOR A SURFACE PERPENDICULAR TO THE SUN DIRECTION. To get the flux on a horizontal surface, multiply by the cosine of the solar zenith angle.
Methods are available to get the incoming solar radiation flux perpendicular to the direction of the sun 1) for the entire solar spectrum (getSolarFlux()), 2) for the PAR band (getSolarFluxPAR()), and 3) for the NIR band (getSolarFluxNIR()).
The very similar function getDiffuseFraction() calculates the fraction of the total flux that is diffuse. The fraction that is direct is simply one minus the diffuse fraction.
Example code for using these solar flux functions is given below.
The predicted solar flux may not perfectly match local predicted solar fluxes due to uncertainty in the local turbidity value. There is a built-in routine to calibrate the turbidity based on measured radiative fluxes.
For the calibration, you must load radiation flux data into a timeseries within the Context. There must be at least one clear-sky day in the timeseries data, and the radiative fluxes must be for the entire solar spectrum in units of W/m2. You can then use the calibrateTurbidityFromTimeseries() method. This method takes one argument, which is a string corresponding to the timeseries variable name containing the radiation flux data.
The REST2 model for solar fluxes was developed for clear-sky conditions and cannot directly be used when clouds are present. If incident solar radiation data is available (e.g., from a weather station), this can be used to calibrate the model to account for the possible presence of clouds. A simple model is described below for doing so.
Consider \(R_{meas,h}\) to be the measured all-wave incoming solar radiation flux on a horizontal plane (clear or cloudy conditions), and \(R_{clear}\) to be the predicted all-wave incoming solar radiation flux predicted by the REST2 model for clear-sky conditions perpendicular to the direction of the sun. This flux can be projected onto the horizontal plane according to
\[ R_{clear,h} = R_{clear}\mathrm{cos}\,\theta_s. \]
The diffuse fraction can be approximated as
\[ f_{diff} = 1-\frac{R_{meas,h} - R_{clear,h}}{R_{clear,h}}, \]
where it is enforced that \(0\leq f_{diff} \leq 1\). The resulting flux that is output from the model is (flux perpendicular to the sun)
\[ R_{model} = R_{clear}\frac{R_{meas,h}}{R_{clear,h}}. \]
In order to enable flux calibration for cloudy conditions, you must 1) Load timeseries data containing the measured all-wave solar radiation flux. This data must cover the entire period of the simulation. 2) Call enableCloudCalibration(), which requires a string corresponding to the timeseries data label.
Below is a Python example showing cloud calibration using Context timeseries:
It should be emphasized that the above model is a relatively simple approximation that produces reasonable fluxes, but more accurate predictions are possible and require much more complicated models.
For applications requiring wavelength-resolved irradiance (e.g., photosynthesis models with wavelength-dependent quantum yield, remote sensing, hyperspectral image simulation), the calculateGlobalSolarSpectrum(), calculateDirectSolarSpectrum() and calculateDiffuseSolarSpectrum() methods compute high-resolution spectral irradiance with the model described in Spectral Solar Irradiance (Model Based on SSolar-GOA).
The spectral irradiance methods use the atmospheric conditions set with setAtmosphericConditions() (the relative humidity must be greater than zero), the ozone column (see Total Column Ozone) and the ground albedo. Results are stored in Context global data for use by other plugins (e.g., the radiation plugin for ray tracing with spectral sources), which refer to a spectrum by its label.
The spectral model also depends on the albedo of the ground surrounding the scene. Light reflected by the ground is partly scattered back down by the atmosphere, which raises the global and diffuse irradiance, most strongly at blue and ultraviolet wavelengths where atmospheric scattering is strongest; without cloud calibration the direct beam is unaffected. The ground albedo defaults to 0.2 and can be set with setGroundAlbedo() (stored in Context global data as atmosphere_ground_albedo):
Example code for calculating spectral irradiance:
Three separate methods are available for the different spectral components:
Each method accepts an optional resolution parameter (default 1 nm) allowing wavelength downsampling. For example, resolution_nm=10.0 produces 231 wavelengths instead of the native 2301.
calculateSensorAtmosphereSpectra() computes the spectra needed to convert radiance leaving the scene into radiance at a sensor above the atmosphere (see Atmosphere Between the Scene and a Satellite Sensor). Given the direction from the scene toward the sensor, it stores in Context global data the path radiance, the adjacency radiance, the upward direct and diffuse transmittances and the atmospheric spherical albedo, together with the direct, diffuse and global solar irradiance at the ground computed with the same atmosphere. To keep the illumination of the scene consistent with the atmosphere between it and the sensor, use these irradiance spectra, rather than those of calculateDirectSolarSpectrum() and calculateDiffuseSolarSpectrum(), for the scene's radiation sources.
| Suffix | Description | Units |
|---|---|---|
| _path_radiance | Radiance scattered by the atmosphere into the sensor without reaching the ground | W/m²/sr/nm |
| _adjacency_radiance | Radiance reflected by the surrounding ground and scattered by the atmosphere into the sensor | W/m²/sr/nm |
| _upward_direct_transmittance | Transmittance of the surface-leaving radiance straight to the sensor, including gas absorption | - |
| _upward_diffuse_transmittance | Transmittance of surface-leaving radiance scattered into the sensor direction, including gas absorption | - |
| _spherical_albedo | Spherical albedo of the atmosphere | - |
| _direct_irradiance | Direct solar irradiance at the ground, normal to the sun direction | W/m²/nm |
| _diffuse_irradiance | Diffuse solar irradiance at the ground on a horizontal surface | W/m²/nm |
| _global_irradiance | Global solar irradiance at the ground on a horizontal surface | W/m²/nm |
| _direction_to_sensor | Unit vector toward the sensor (a single vec3) | - |
The turbidity set with setAtmosphericConditions() (the Ångström coefficient β, see Understanding the Turbidity Parameter) is converted to an aerosol optical depth at 550 nm using an Ångström exponent of 1.3. Conditions outside the range of the look-up table (solar zenith above 70 degrees, view zenith above 60 degrees, aerosol optical depth at 550 nm above 1.2, surface pressure below 701.2 hPa or above 1050 hPa, or a water vapor path, the column water vapor times \(1/\mu_s + 1/\mu_v\), above 320 g/cm²) and a sun below the horizon raise an error; pressures between 1013 and 1050 hPa are extrapolated linearly. The model assumes clear sky, so it cannot be used together with cloud calibration.
calculateSensorThermalAtmosphere() computes the thermal infrared transmittance and upwelling and downwelling atmospheric radiance for a sensor viewing the scene from the given direction (see Thermal Infrared Atmosphere Between the Scene and a Satellite Sensor). They are stored in Context global data as spectra of 234 points from 5502 to 15326 nm, with the suffixes "_thermal_transmittance", "_thermal_upwelling_radiance" and "_thermal_downwelling_radiance" (radiances in W/m²/sr/nm). The atmospheric profile is built from the air temperature, humidity and pressure set with setAtmosphericConditions() and the ozone column (Total Column Ozone); the column water vapor derived from the temperature and humidity is stored with the suffix "_thermal_water_vapor_cm", and the unit vector toward the sensor with the suffix "_direction_to_sensor". The model does not depend on the sun's position, so night-time imagery can be simulated.
The column water vapor must lie within 0.1-7.5 g/cm², the pressure within 701.2-1050 hPa, the ozone column within 150-450 DU, and the view zenith angle must not exceed 60 degrees. The supported air temperatures depend on the pressure: 246.2-321.7 K at 1013 hPa. Conditions outside these ranges raise an error.
The sky longwave radiance that the scene reflects in a thermal band should come from the same atmosphere. getThermalSkyFlux() gives the clear-sky longwave irradiance on a horizontal surface within a band, the band-specific counterpart of the broadband getAmbientLongwaveFlux(); it is the diffuse flux to give an emission band with the same wavelength bounds:
The downwelling longwave radiation flux from the sky can be calculated using the getAmbientLongwaveFlux() function. This function is based on the Prata (1996) model and returns the clear-sky downwelling longwave radiation flux on a horizontal surface in W/m2, integrated over all longwave wavelengths (see Ambient Longwave Flux). It is the diffuse flux to give an emission band without wavelength bounds. For an emission band with wavelength bounds, use getThermalSkyFlux() instead, which gives the clear-sky irradiance within the band (see Simulating Satellite Thermal Imagery).
1.13.2