Star photometry with all-sky cameras to retrieve aerosol optical depth at nighttime

The lack of aerosol optical depth (AOD) data at night can be partially addressed through Moon photometer measurements or fully covered with star photometer observations. However, the limited availability and complexity of star photometers have motivated this study to use all-sky cameras to extract starlight signals and derive AOD at night using star photometry. For this purpose, eight all-sky cameras were configured and deployed in nine different locations to capture raw images with varying exposure times every 2 min during the night. This work proposes a novel methodology to extract the starlight signal from the raw data from all-sky cameras and convert it into AOD values. This process consists of the following steps: removing the background image, selecting the pixels, and extracting the signal for each star from a predefined list of 56 stars; performing in situ Langley calibration of the instruments and retrieving the total optical depth (TOD); calculating the effective wavelength for each camera channel; deriving the AOD by subtracting the gas contribution to TOD; and averaging, cloud-screening, and quality-assuring the AOD time series. The AOD time series obtained through this methodology are compared with independent AOD measurements from collocated Moon photometers in the nine locations. The obtained results show that the AOD values derived with the proposed method generally correlate with reference values, often achieving correlation coefficients (r) above 0.90. The AOD values retrieved using the cameras tend to overestimate the reference values by approximately 0.02 and exhibit a precision of around 0.03–0.04. The agreement between the two datasets varies with wavelength and decreases at high-latitude locations, likely due to the poorer performance of Langley calibration in these regions. AOD values align well with day-to-night transitions obtained by solar photometers, demonstrating their reliability. Despite the slight overestimation, the AOD values derived by this new method approximate the real values and provide coverage throughout the entire night, without requiring the presence of the Moon. Therefore, they serve to study and monitor the nocturnal evolution of AOD.

1 Introduction

The study of aerosol optical properties and their temporal variability is essential to better understand their role in the Earth's climate system. Aerosols influence radiative forcing both directly, by scattering and absorbing solar radiation, and indirectly, by modifying cloud properties and lifetimes (Boucher et al., 2013). This is particularly critical in polar regions, where aerosols can have a significant influence due to the unique atmospheric and surface conditions (Law and Stohl, 2007). In these regions, the high albedo of ice and snow amplifies the impact of aerosols on the radiative balance, and the deposition of light-absorbing particles, such as black carbon, can accelerate ice melt and snowmelt (Hansen and Nazarenko, 2004; Bond et al., 2013). It is therefore essential to continuously monitor aerosol properties and generate long-term data series to study potential changes and trends, especially in polar regions, which remain underrepresented in global aerosol monitoring due to logistical challenges and harsh environmental conditions.

One of the most important aerosol properties is the aerosol optical depth (AOD), which is related to the extinction (scattering + absorption) of light caused by aerosols as it crosses the atmosphere. AOD serves as a proxy for the quantity of aerosol particles present in the atmosphere. AOD spectral variation, quantified through the empirical Ångström law (Angström, 1961), also provides information about particle size predominance. AOD is typically calculated using measurements from photometers that, by directly pointing at the Sun, measure the direct solar irradiance reaching the Earth's surface at one wavelength or generally multiple wavelengths separately. If the photometer is properly calibrated, these direct irradiance measurements can be converted to AOD values through the Beer–Bouguer–Lambert law. One method to calibrate photometers for obtaining AOD values is a technique named “Langley” or a “Langley plot” (Shaw, 1983). This method involves performing direct irradiance measurements at different solar zenith angle (SZA) values, fast enough to assume that the direct extraterrestrial solar irradiance and the atmospheric optical depth remain constant (Toledano et al., 2018). It allows calculating AOD without knowing the true extraterrestrial solar irradiance in physical units. This methodology is used by AERONET (AErosol RObotic NETwork; Holben et al., 1998, https://aeronet.gsfc.nasa.gov, last access: 28 June 2025), a network which includes more than 500 photometer stations distributed worldwide. AERONET calibrates its field photometers by intercomparison with reference photometers (masters), which are themselves calibrated using the Langley technique (Giles et al., 2019).

Thanks to this method and direct solar irradiance measurements, AOD can be obtained during the daytime but not at nighttime. However, knowing the AOD at night is crucial, among other issues, to study complete aerosol cycles (Perrone et al., 2022; Herrero del Barrio et al., 2023) and the aerosol interaction with longwave radiation, especially in polar regions, where continuous night periods occur for several consecutive months during the polar night. Fortunately, some photometer models, such as those used in AERONET, were updated a few years ago to measure direct lunar irradiance at the ground from the first Moon quarter to the third one (Berkoff et al., 2011; Barreto et al., 2013). The main challenge with these measurements is that the direct extraterrestrial lunar irradiance changes significantly, even within a single night, making the Langley technique insufficient; it is then necessary to know the direct extraterrestrial lunar irradiance each time. In this regard, the ROLO model (Robotic Lunar Observatory; Kieffer and Stone, 2005), in the case of AERONET, and its implementation RIMO (ROLO Implementation for Moon photometry Observation; Barreto et al., 2019) have been used to estimate the direct extraterrestrial lunar irradiance. These models did not initially provide accurate AOD values (Barreto et al., 2016, 2017), but this issue was resolved after applying corrections to them (Román et al., 2020; González et al., 2020; Uchiyama et al., 2019). However, this does not address another issue with lunar photometry: there are no direct lunar irradiance measurements during half of the lunar cycle (from the third quarter to the first one), and even during the other half, the Moon is not visible throughout the entire night, especially near the quarter phases. This means that lunar photometry only partially fills the gap in our knowledge of nighttime AOD.

This gap in nighttime data series can be filled by using star photometers. These instruments essentially consist of a telescope on a tracking mount equipped with a sensor that records direct star irradiance at the ground. AOD can be derived from these measurements by applying star photometry (see Ivănescu et al., 2021, and references therein), which works similarly to solar or lunar photometry but covers the entire night, as there are always multiple stars visible in the sky. Moreover, the direct extraterrestrial irradiance of stars is highly stable over time, allowing the Langley technique to be applied without requiring prior knowledge of these extraterrestrial values. In addition to Langley calibration, the capability of star photometers to measure direct irradiance from multiple stars enables the application of alternative calibration methods, such as the two-star method, which uses two different stars at varying zenith angles (Leiterer et al., 1995; Pérez-Ramírez et al., 2008a), or the multi-star Langley method proposed by Ivănescu and O'Neill (2023). Although star photometer measurements are still used to study AOD at night (e.g., Graßl et al., 2024), the availability of star photometers is extremely limited due to, among other issues, their high cost and the fact that they are not fully automated. There are currently only about five star photometers in operation worldwide for aerosol monitoring (Barreto et al., 2019; Román et al., 2020).

This scarcity of star photometers has led to the search for alternatives that can be more widely used to obtain AOD at night. The fact that the sensor of one of the existing star photometers is a CCD camera that captures images where the star occupies only a small fraction of the CCD image (Pérez-Ramírez et al., 2008b) suggests that other instruments capable of capturing images with visible stars could also be used for star photometry. All-sky cameras can be used for this purpose. These instruments typically consist of a CMOS sensor with a fisheye lens that captures images of the entire sky dome, as opposed to focusing on a smaller magnified area like star photometers. Most modern all-sky cameras can operate both day and night, and their wide field of view allows them to capture stars, which occupy a very small area in the image but include a large number of stars as they observe the entire sky at once. Therefore, if the starlight can be extracted from these images, the nighttime AOD can be derived in a similar way to that of a star photometer. In fact, it has already been demonstrated that AOD can be estimated during daytime using all-sky cameras with different methods (Cazorla et al., 2008; Román et al., 2022; Scarlatti et al., 2023). Among its advantages, all-sky cameras are fully automated instruments, with no moving parts, and are much more cost-effective.

In this framework, our starting hypothesis is that the starlight signal, proportional to star irradiance, can be extracted from all-sky camera images, and thus, AOD can be obtained using star photometry. The main idea is that while the signal from each individual star – and consequently the AOD derived from it – may be very noisy due to technical limitations of all-sky cameras, the ability to capture the signal from many stars simultaneously and to take continuous images allows a large number of AOD values to be obtained. Although each individual value may have high noise, averaging them reduces the overall noise. With this in mind, the main objective of this work is to develop a new methodology capable of extracting the starlight signal from images captured by all-sky cameras and calculating AOD from these values. Additionally, we aim to test and validate this new method by comparing the results with independently obtained measurements using lunar photometry across various locations.

This paper is structured as follows. Section 2 introduces the technical characteristics and configuration of the all-sky cameras used in this study, as well as the locations where they are installed. The proposed methodology for extracting star direct irradiances and converting them into AOD values is detailed in Sect. 3. Section 4 presents a general comparison between the AOD retrieved from the all-sky cameras and the values independently obtained by Moon photometers, along with two case studies. Finally, the main conclusions are summarized in Sect. 5.

2 Instrumentation and sites

The Group of Atmospheric Optics of the University of Valladolid (Grupo de Óptica Atmosférica; GOA-UVa) has created and developed GOA-SCAN (GOA all-Sky CAmeras Network), a network of all-sky cameras deployed around the world. GOA-SCAN is managed by GOA-UVa, which establishes protocols for measurements, data storage, calibration, and quality assurance to maximize the information recorded by the cameras and to standardize and process the data consistently. This network comprises all-sky cameras from GOA-UVa as well as from other research groups.

One of the most commonly used camera models in GOA-SCAN is the OMEA-3C all-sky camera, manufactured by Alcor System. This camera features a fisheye lens (180° × 180°, equidistant projection) mounted over a CMOS sensor (SONY IMX178) with 3096 × 2080 px (6.44 MP allowing a resolution of 5.4 arcmin px⁻¹) and 14-bit resolution. This sensor has three distinct RGB filters (red, green, and blue; see Fig. 1a) arranged in an RGGB Bayer pattern mosaic, along with an additional IR-cut filter (see Fig. 1b) placed over the RGB filters to block infrared skylight. The spectral responses of the RGB filters, the IR-cut filter, and the combined effect of both (the total spectral response of each of the three RGB channels) are shown in Fig. 1d. In GOA-SCAN, there is a modified version of the OMEA-3C camera that is identical to the original but replaces the infrared filter with a tri-band filter (see Fig. 1c), which narrows the spectral bandwidth of the three color channels (see Fig. 1e). We will refer to this modified camera version as OMEA-3C-TF (TF for tri-band filter). In this work, we have also used grayscale (Gr) images, which are obtained as a weighted combination of the three color channels, as shown in Eq. (1):

$$\begin{array}{}\text{(1)}& \mathbf{Gr}=\mathrm{0.299}\mathbf{R}+\mathrm{0.587}\mathbf{G}+\mathrm{0.114}\mathbf{B}.\end{array}$$

This conversion is directly applied by the COLOR_RGB2GRAY option of the cvtColor function of the OpenCV Python library (Bradski, 2000, https://opencv.org, last access: 28 June 2025). The spectral response of this grayscale channel is presented in Fig. 1 for the OMEA-3C (Fig. 1d) and the OMEA-3C-TF (Fig. 1e).

Figure 1 Normalized spectral transmittance (NST) of different filters: (a) red (R), green (G), and blue (B) filters of the Bayer mosaic over the SONY IMX178 CMOS sensor in the OMEA cameras, normalized to the maximum value among the three channels. (b) IR-cut filter included in the OMEA-3C model. (c) Tri-band filter included in the OMEA-3C-TF model. (d) Combination of the RGB filters and the IR-cut filter in the OMEA-3C model. (e) Combination of the RGB filters and the tri-band filter in the OMEA-3C-TF model. Panels (d) and (e) also show the spectral response of the gray signal, which is computed as a weighted combination of the three color channels using the following weights: 0.299 for R, 0.587 for G, and 0.114 for B.

The camera sensor and lens are enclosed in a fully anodized aluminum casing system, which is completely watertight (IP67) and protected by a BK7 glass dome on the top. This model features an automatic heating system to defrost the dome and to avoid rainwater droplets on its surface. The control of these cameras is carried out from an independent computer directly connected to them via cable. It is typically housed either in a weatherproof enclosure outdoors or inside a nearby building. The computers of the OMEA cameras belonging to GOA-SCAN are equipped with the GOA-OMEA-Capture software, developed by GOA-UVa, to configure settings and manage image capture, downloading, processing, and data transmission (Antuña Sánchez, 2022).

All OMEA cameras in GOA-SCAN are configured to capture multi-exposure image sequences every 5 min during the daytime (with some exceptions set to every 2 min) and every 2 min at nighttime (SZA above 97°). These sequences consist of consecutive image captures taken at different exposure times. For nighttime, there are two distinct multi-exposure configurations: Moon, used when the Moon is between the first and third quarters and its zenith angle is below 80°, and Moonless, used the rest of the time. The exposure times for the OMEA-3C camera in Moon mode are 0.1, 0.5, 2, 5, 10, 15, and 20 s, while for Moonless conditions, the times are 5, 10, 15, 20, and 30 s. For the OMEA-3C-TF camera, the exposure times differ: 1, 5, 10, 20, and 30 s in Moon mode and 3, 6, 12, 24, and 48 s in Moonless conditions. The shortest exposures are used to detect the starlight of the brightest stars without saturation, while the longest exposures ensure sufficient signal from the fainter stars. Nighttime images are configured to be recorded with an ISO amplification set to a gain of 10 dB (0.29 $e-/\text{ADU}$) for OMEA-3C and 15 dB (0.16 $e-/\text{ADU}$) for OMEA-3C-TF. All images are cropped to 2000 × 2000 px to reduce memory usage while still covering the full sky. The image data are stored in raw format (14 bits in 16-bit format, without demosaicing or white-balance correction) in a single “.h5” file that includes all images of the same sequence along with additional metadata. GOA-OMEA-Capture is programmed to frequently transmit these files to the GOA-SCAN servers, where the data are processed to produce products such as high dynamic range (HDR) images (Antuña Sánchez, 2022). These HDR images are made available on the GOA-UVa web page, updated to reflect the latest data from each GOA-SCAN station (https://goa.uva.es, last access: 28 June 2025; “Sites and Measurements” tab).

Image data from eight different GOA-SCAN OMEA cameras installed at nine different locations have been utilized in this work. Figure 2 presents a world map highlighting the mentioned stations. As shown, these nine stations cover a wide range of latitudes, including three polar stations – Ny-Ålesund and Andøya (Norway) as well as the Marambio station (Antarctica, Argentina) – and two subtropical stations in the Canary Islands: Fuencaliente (La Palma, Spain) and Izaña (Tenerife, Spain). Table 1 provides detailed information about these stations and the cameras installed there, each represented by its ID. All cameras are OMEA-3C, except for the C013 camera installed in Valladolid (Spain), which is an OMEA-3C-TF model. Some cameras were relocated, such as the C005 camera, which was moved from Valladolid to Fuencaliente in September 2021 to monitor the Tajogaite volcanic eruption and later installed in Izaña and in Andøya. Although some of the mentioned cameras have been installed since 2018, only data from July 2020 onward have been used in this work, as the GOA-OMEA-Capture software was not available before.

Figure 2 Geographical locations of the nine GOA-SCAN stations utilized in this work.

Table 1 Detailed information about the location and measurement period of each GOA-SCAN all-sky camera used in this work. The acronyms in the table are defined as follows: SMN (Servicio Meteorológico Nacional; Argentinian Meteorological Service), AWI (Alfred Wegener Institute for Polar and Marine Research), GOA-UVa (Grupo de Óptica Atmosférica de la Universidad de Valladolid; Atmospheric Optics Group of the University of Valladolid), AEMet (Agencia Estatal de Meteorología; Spanish Meteorological Agency), ASC (Andøya Space Center), DWD (Deutscher Wetterdienst; German Meteorological Service), GFAT-UGR (Grupo de Física de la Atmósfera de la Universidad de Granada; Atmospheric Physics Group of the University of Granada), and PMOD-WRC (Physikalisch-Meteorologisches Observatorium Davos and World Radiation Center). The present corresponds to 22 January 2025.

All GOA-SCAN cameras are routinely inspected and cleaned. Each camera, horizontally leveled to point directly at the zenith, undergoes geometric calibration whenever it is relocated or experiences even slight changes in its position. This calibration is carried out using the ORION software (Antuña-Sánchez et al., 2022), which determines the azimuth and zenith angles of the celestial vault viewed by each camera pixel, enabling the precise localization of any body within the images if its coordinates are known.

In addition, each of the stations listed in Fig. 2 and Table 1 is equipped with at least one CE318-T Sun–sky–Moon photometer (manufactured by Cimel Electronique). All these photometers are part of AERONET, with periodic calibration, maintenance, and data management overseen by GOA-UVa following the AERONET scheme (Giles et al., 2019). These management tasks are supported by CAELIS (Fuertes et al., 2018; https://www.caelis.uva.es, last access: 28 June 2025), a software tool designed to manage and process data from photometers controlled by GOA-UVa and to facilitate data visualization. The cloud-screened AOD data collected during daytime (Sun) and nighttime (Moon) by these photometers have been directly obtained in this work from CAELIS, with the processing methodology described in detail by González et al. (2020) and Román et al. (2020).

3 Methodology

In this section it is explained how the starlight signal is extracted from the camera raw data and converted into AOD values. Stars are punctual sources of light; however, their image in all-sky photos spans several pixels due to long time exposures (as stars move during the integration time) and primarily due to the camera's point spread function (Piotrowski et al., 2013). Due to the RGGB Bayer pattern, one-fourth of the pixels record the blue signal, another one-fourth record the red signal, and half of the pixels record the green signal. This means that the star signal is not fully captured by a single color channel but is distributed across them. Consequently, it is necessary to estimate the signal for each color channel in the remaining pixels of the camera's sensor with different color filters. To achieve this, each 14-bit 2000 × 2000 raw image is converted into a 2000 × 2000 × 3 RGB image using a demosaicing algorithm: the demosaicing_CFA_Bayer_bilinear function from the Python colour_demosaicing package (Losson et al., 2010, https://github.com/colour-science/colour-demosaicing, last access: 28 June 2025). This method produces three 2000 × 2000 px matrices, one for each color channel, containing the interpolated signal. Additionally, a fourth grayscale matrix is calculated from the R, G, and B matrices using Eq. (1). This grayscale matrix represents a broader spectral response (see Fig. 1), but all its pixels are sensitive to the brightest pixels of the observed stars (this was the ultimate reason to calculate the gray channel).

3.1 Background

To extract the starlight signal from the obtained images, it is necessary to remove the background signal caused by sky glow, scattered moonlight, and other factors, such as auroras in polar regions or dark and readout noise frames. The Background2D class, from the photutils (version 1.13.0) Python library (Bradley et al., 2024; https://photutils.readthedocs.io/en/stable/, last access: 28 June 2025), has been used for each channel since it is generally applied to estimate the 2D background noise in astronomical images. Background2D has been configured to use a 10 × 10 px box in which to estimate the background as the median of the pixels within the box (bkg_estimator = MedianBackground). This method iterates over, each time rejecting values that are less or more than a specified number of standard deviations (SigmaClip), which has been chosen as 3 in this work, from a center value. The obtained result has been then smoothed using a 3 × 3 median filter (filter_size = (3,3)) to suppress local underestimations or overestimations.

This proposed methodology appears to be effective for estimating the background in our case, since the stars in the images typically occupy a small area and exhibit abrupt intensity changes compared to the background. Some parameters, such as the box size and SigmaClip, have been manually selected in this work due to their good performance with the cameras used. However, the optimal values may vary for other cameras with different characteristics, such as pixel resolution.

Some examples of the background calculation are shown in Fig. 3 for different cameras, locations, color channels, and conditions. The left panels of Fig. 3 correspond to the true RGB color composition of each case, after applying a white-balance correction. Note that although the raw camera signal is recorded in 14 bits, it is stored in 16 bits; consequently, the digital counts (DCs) shown in the color bars of Fig. 3 are scaled to 16 bits.

The upper panels of Fig. 3 depict a summer night photo taken at Izaña under the presence of the Moon (Moon phase angle approximately −40°) and Saharan desert dust conditions (AOD at 440 nm about 0.26). The gray image (Fig. 3b) and its calculated background (Fig. 3c) look very similar; however, after subtracting this background from the original signal (Fig. 3d), the presence of stars becomes evident, indicating that the background correction is apparently not removing the star signals. Additionally, the differences between the original image and its background reveal that some parts of the horizon, such as buildings with non-smooth variations, are not perfectly modeled in the background calculation. These horizon parts can be easily masked to avoid their use. However, Fig. 3d also shows that the calculated background underestimates the true one for pixels close to the Moon, and hot pixels cannot be removed with this method (see some pixels in the corners of Fig. 3d), as they exhibit a shape similar to stars. Moreover, some lunar reflections on the camera dome, which appear as circular or linear patterns, are not captured in the background and could be mistaken for stars.

Figure 3e to h present a Moonless case for the green channel at Lindenberg, where the Milky Way is accurately estimated in the background (Fig. 3g), and a large number of stars become visible when the background is removed (Fig. 3h). Another case with the Moon (Moon phase angle of −68°) is shown in Fig. 3i to l for Valladolid using the OMEA-3C-TF camera. In this case, a cloud can be observed at the top of the image, which is also captured in the background but underestimated, similar to the Moon and the horizontal city skyline, as shown in Fig. 3l. On the other hand, the Orion and Ursa Major constellations are clearly visible in the upper left and bottom right parts, respectively, of Fig. 3l, where Jupiter is also visible to the left of the Moon. Finally, a case in Davos featuring a red aurora is shown in Fig. 3m to p. In this case, the background accurately captures the aurora, which is completely removed in Fig. 3p. However, a satellite streak can be observed in the left part of Fig. 3p, highlighting that line-shaped features are not well captured by the background estimation method used.

Figure 3 RGB white-balanced color images for (a) Izaña on 10 July 2022 at 21:16 UTC (camera C005 with an exposure time of 5 s), (e) Lindenberg on 25 May 2024 at 00:00 UTC (camera C009 with an exposure time of 30 s), (i) Valladolid on 19 January 2024 at 23:40 UTC (camera C013 with an exposure time of 10 s), and (m) Davos on 10 May 2024 at 21:30 UTC (camera C038 with an exposure time of 20 s). The second column presents the grayscale raw signal for the gray (b), green (f), blue (j), and red (n) channels, corresponding to the images in the first column. Panels (c), (g), (k), and (o) show the background signal estimated from the respective raw images in the second column. Panels (d), (h), (l), and (p) display the difference between the raw signal (b, f, j, n) and the estimated background (c, g, k, o). The color bars represent pixel digital counts scaled to 16 bits. The maximum values of the color bars in the second (b, f, j, n) and third (c, g, k, o) columns correspond to the 99.5th percentile of the pixel signals in the respective raw images (second column). For the last column, the maximum value is set as the product of the exposure time and 50 digital counts per second. This figure is compressed, so the signal from the stars does not appear as it does in the original images. For this reason, it is recommended to access the high-resolution version of the figure available in the online version of this article.

3.2 Star signal (SS)

Once the background is estimated, it is subtracted from the raw demosaiced image for each channel. This new signal is referred to as the background-corrected pixel signal (BCPS). The light signal of a star should correspond to the sum of the BCPS values of the pixels that capture its light. Therefore, it is necessary to identify these pixels.

To achieve this, the coordinates of the star in the celestial vault for the image must first be determined. This task is performed using the SkyCoord function, based on ephemeris, from the Astropy Python package (version 6.1.1; The Astropy Collaboration et al., 2022; https://www.astropy.org, last access: 28 June 2025), taking as input the camera's location coordinates and the image's timestamp, considering the midpoint along the exposure time. If the star's zenith angle exceeds 80°, it is discarded. Then, the closest pixel to the obtained star coordinates is identified using the previously mentioned camera calibration (see Sect. 2). For each channel, an 80 × 80 px square box centered on this pixel is selected, assuming that the entire star image is contained within this box. An example of this box, centered on the Alnilam star, is shown in Fig. 4a for the red channel of a raw image captured with the C013 camera in Valladolid on 19 January 2024 at 23:40 UTC, with an exposure time of 20 s (assumed timestamp: 19 January 2024 at 23:40:45).

The BCPS subimage within the chosen box is expected to contain the image of the selected star, but occasionally other stars may also appear in the same region. This is the case in Fig. 4a, where the stars Alnitak and Mintaka also appear in the image, showing the Orion Belt together with Alnilam. A segmentation of this subimage is performed using the detect_sources function from photutils to identify the stars within the box. This function identifies sources as regions with at least N connected pixels, each exceeding a specified threshold value. In this work, N is set to 10, mainly to avoid hot pixels, and the threshold for each channel is defined as the BCPS median value of all pixels within the box, plus 3 times the standard deviation of all these pixels. In the example shown in Fig. 4, the calculated threshold is 689.6 DC. Figure 4b displays the pixels where the BCPS exceeds this threshold. While the three observed stars, as well as some isolated pixels, surpass this value, only the three stars are identified as such (segments marked in color in Fig. 4c) because they have a sufficient number of connected pixels (at least 10 px).

Figure 4 (a) Background-corrected pixel signal (BCPS), in digital counts (DCs), for a subimage within a selected 80 × 80 px box, centered on the Alnilam star, derived from the red channel of the raw image recorded in Valladolid on 19 January 2024 at 23:40 UTC (camera C013 with an exposure time of 20 s). Pixels where the BCPS exceeds the threshold of 689.6 DC are highlighted in white in panel (b). Segments identified as stars are marked in color in panel (c), where the yellow segment corresponds to Alnilam, the red cross indicates its centroid, and the white circle denotes the region used to extract the star's BCPS.

Additionally, a deblending process is applied using the deblend_sources function from photutils to separate segments corresponding to overlapping stars, which is a very uncommon occurrence. Subsequently, the centroid coordinates of each segment are calculated using the SourceCatalog function from photutils (see the red cross in Fig. 4). The segment with the shortest distance between its centroid coordinates and the analyzed star's coordinates is identified as the corresponding star. If this distance exceeds 1°, the star is considered not identified. In the case of Fig. 4c, the yellow segment is identified as the Alnilam star.

Once the star is identified in the image, the star signal (SS) can be calculated as the sum of the BCPS values of the segmented pixels identified as star pixels – that is, the sum of the BCPS values of all yellow pixels shown in the example of Fig. 4c. However, this method tends to underestimate the actual star signal, as some pixels surrounding the star segment may still record a fraction of the total signal without exceeding the defined threshold. To address this, a circular aperture with a radius of 5 px, centered on the centroid of the star (see the white circle in Fig. 4c), is used to encompass the full star signal. Finally, the SS is calculated as the integral of the BCPS within this aperture using the aperture_photometry function from photutils.

3.3 Star direct irradiance screening and smoothing

In this work, the SS has been calculated using the described method for all available raw images from the cameras listed in Table 1 and for the following 56 selected stars: Achernar, Acrux, Adhara, Aldebaran, Algol, Alhena, Alioth, Alkaid, Alnair, Alnilam, Alphecca, Alpheratz, Altair, Antares, Arcturus, Atria, Avior, Bellatrix, Betelgeuse, Canopus, Capella, Caph, Castor, Deneb, Diphda, Dubhe, Elnath, Eta Herculis, Fomalhaut, Gacrux, Gamma Herculis, Hadar, Kaus Australis, Kochab, Markab, Merak, Miaplacidus, Mimosa, Mintaka, Mirfak, Mizar, Peacock, Polaris, Pollux, Procyon, Regor, Regulus, Rigel, Rigil Kentaurus, Sargas, Shaula, Sheratan, Sirius, Spica, Vega, and Wezen. This list includes the most prominent stars visible from both hemispheres.

Subsequently, the SS of each star is calculated for all four channels and every available exposure time. However, SS values corresponding to exposure times shorter than 2 s are excluded due to the low signal-to-noise ratio under such conditions. The SS values obtained from at least one saturated pixel are also excluded. In this work, a pixel is considered saturated when its raw signal (without background correction) exceeds 58 982 DC, a threshold corresponding to 90 % of the maximum possible recorded signal in a 16-bit format (65 535 DC). In addition, all SS values obtained for SZA values below 110° are discarded to avoid contamination from sunlight.

The SS values are proportional to the exposure time, representing uncalibrated irradiation (energy per unit area). To standardize the measurements, the obtained SS values are normalized by dividing them by the exposure time, resulting in a metric proportional to the star's irradiance (power per unit area). This new parameter is referred to as the uncalibrated star irradiance (USI).

The left panels of Fig. 5 show the USI values obtained during a winter night at Andøya for various stars. For Procyon (Fig. 5a), the signal exhibits a pattern similar to that of direct solar irradiance during daytime under cloud-free conditions, primarily varying with the relative optical air mass (m). The recorded m values for Procyon ranged from 5.57 (at star set and star rise) to 2.28 that night, with the USI decreasing as m increased. The relative optical air mass in this work is calculated using the equation of Kasten and Young (1989).

Figure 5 Uncalibrated star signal (USI) at Andøya (camera C005) from 3 to 4 January 2024 for four color channels (red, green, blue, and gray). Panels (b), (d), (f), and (h) show the corresponding USI values after applying the time-smoothing process. The USI values are presented for the following stars: Procyon (a, b), Vega (c, d), Kochab (e, f), and Polaris (g, h).

Stars may rise at the beginning of the night, similar to the Sun during the day, but it is also possible for stars to rise at any time during the night or even to set at the beginning of the night and rise again later. The latter case applies to Vega (Fig. 5c), which did not even set in the example shown in Fig. 5 before starting to rise again when the optical air mass reached 3.19.

Additionally, due to the position of the stars in the sky and the observer's coordinates, the optical air mass of the stars varies to a greater or lesser extent depending on the star. This is illustrated in Fig. 5e for Kochab, whose optical air mass ranged only between 1.25 and 1.00 that night, resulting in very stable USI values throughout the night. An extreme case is Polaris, which is closely aligned with the Earth's north rotational axis, causing its optical air mass to vary only between 1.07 and 1.06 during the night shown in Fig. 5g.

In general, the USI values shown in the left panels of Fig. 5 appear noisy, with a high deviation, particularly in the red channel. Noise in the USI values is expected due to shot noise (inherent to light), star scintillation caused by atmospheric turbulence (Klaus et al., 2004), and the uncertainties introduced by the proposed method for extracting these values. However, Kochab and Polaris, which both exhibit similar USI values (2000–5000 DC s⁻¹) and minimal variation in optical air mass, display distinct noise patterns, suggesting an additional noise source.

As noted, the position of Polaris in the sky changes very slightly but is not static; for instance, its azimuth and zenith varied by only 3 and 1°, respectively, during the night shown in Fig. 5g. This limited movement ensures that the position of the pixels capturing Polaris does not change rapidly. Consequently, as stars are point light sources, the maximum star signal is consistently recorded by the same pixel (corresponding to the star's center) for a significant period. This pixel is associated with only one of the three color channels in the Bayer pattern. For instance, if the maximum star signal falls on a green pixel, it will not be detected by the blue or red pixels, and as a result, the demosaiced images for these two channels will not reflect this maximum signal. This explains the behavior of the USI time evolution for Polaris in Fig. 5, where all the channels exhibit a cyclic pattern. The maximum of each channel corresponds to the maximum star signal reaching a pixel of that channel, ensuring that a maximum never appears simultaneously in two different channels. This pronounced oscillation, caused by the use of a single color filter (R, G, or B) per pixel, occurs much faster for all the stars that exhibit more movement across the celestial vault. Consequently, it is also, at least in part, responsible for the high deviation observed in the USI values of these stars.

To reduce this deviation, the USI measurements are smoothed using the following method for each star and channel. For a target USI value, a time window of ±4 min is defined around it. All available USI data within this window are normalized to the optical air mass of the target USI value. This normalization is performed by multiplying each value by its own optical air mass and dividing by the optical air mass of the target USI value. The median of all normalized USI values within the time window, including the target USI value, is then calculated. Any USI value that deviates by more than 20 % (in absolute terms) from this median is considered an outlier and excluded from the time window. A new median is then computed using the remaining normalized USI data, which is taken as the smoothed USI value. Additionally, the standard deviation of the remaining normalized USI values is calculated and used as the uncertainty of the smoothed USI. The number of data points used in the final median and standard deviation calculations is also recorded (N_smooth). If N_smooth is less than 2 or if the final USI value is not above 1 DC s⁻¹, the USI measurement is discarded.

The USI values displayed in the left panels of Fig. 5 have been smoothed and are represented in the right panels of the same figure (panels b, d, f, and h). The smoothing process effectively reduces the noise in the data, as can be observed. Additionally, an abrupt decrease in the non-smoothed USI values, caused by the presence of clouds, is noticeable around 16:00 for Vega, Kochab, and Polaris. This decrease is also reflected in the smoothed values; however, some cloud-contaminated data are filtered out due to the presence of outliers and the limited availability of data under these conditions. Consequently, N_smooth could serve as a proxy for cloud screening.

Finally, it is remarkable that most of the USI values shown in Fig. 5 were derived from images captured under intense auroral activity. However, the effect of these northern lights is not appreciated in the final USI data, highlighting the effectiveness of the proposed background characterization and removal method of Sect. 3.1. Hereafter, the term USI values will refer exclusively to the smoothed USI values.

3.4 Langley calibration and TOD calculation

The direct solar, lunar, or stellar irradiance at one single wavelength reaching the Earth is related to the atmospheric total optical depth (TOD) through the Beer–Lambert–Bouguer law. This law is monochromatic, and there are modifications based on approximations for irradiance measurements performed with broadband filters (Rufener, 1986); however, in this work, the measured irradiance has been approximated as the irradiance at an effective wavelength (see Sect. 3.5) that satisfies the Beer–Lambert–Bouguer law. For stars, this relationship can be expressed in terms of the natural logarithm, as shown in Eq. (2):

$$\begin{array}{}\text{(2)}& \ln \left(\text{USI}\right)=\ln \left({\text{USI}}^{\mathrm{0}}\right)-m\cdot \mathit{\tau },\end{array}$$

where m is the relative optical air mass, τ is the TOD, and USI⁰ is the uncalibrated extraterrestrial star irradiance (i.e., the star irradiance at the top of the atmosphere). The calibration constant that converts USI into well-calibrated physical units is the same constant that converts USI⁰ into the actual extraterrestrial irradiance. This allows Eq. (2) to enable the use of uncalibrated irradiances, thus eliminating the need for tedious laboratory calibrations. If USI is measured for a star, the TOD can then be estimated using Eq. (2); only the constant value of ln (USI⁰) needs to be known for this star.

To this end, the Langley plot calibration method is applied. This method assumes that the TOD remains constant during a period in which the optical air mass varies. Under this assumption, ln (USI⁰) can be estimated as the y intercept of a least-squares linear fit of ln (USI) as a function of the optical air mass. The slope of the fit corresponds to the total optical depth with a negative sign (see Eq. 2).

This method is applied in this study to determine ln (USI⁰) for each star and channel on a nightly basis. First, rising and setting periods for each star are identified within a single night, as the Langley plot calibration is performed separately for each period. Only USI data recorded under optical air masses between 2 and 5 are used for this purpose, similar to the air mass range employed by AERONET for Langley calibrations (Toledano et al., 2018). The selected ln (USI) data are fitted to the optical air mass using a weighted least-squares linear fit, where the weight of each data point is the inverse square of its uncertainty (as determined in Sect. 3.3). This ensures that data points with higher uncertainty contribute less to the fit.

Once the initial fit is obtained, measured data with an absolute difference greater than 1 % compared to the modeled values from the fit are removed. The weighted least-squares fit is then recalculated using the remaining data. This outlier rejection process is iteratively applied until no additional data points are excluded in a given iteration. For the final iteration, the ln (USI⁰) (y intercept) and TOD (negative slope) values and their uncertainty (from the fit) are stored, along with the optical air mass range (maximum m minus minimum m), the weighted correlation coefficient (r_w), and the number of data points used (N_Langley).

From all the obtained ln (USI⁰) values, only those meeting the following criteria are considered: (1) an optical air mass range greater than 2, to ensure sufficient sensitivity to variations in m during the fit; (2) a weighted correlation coefficient lower than −0.9, to guarantee a strong inverse correlation (the threshold is negative due to the inverse relationship between ln (USI) and m); (3) N_Langley greater than 40, to ensure an adequate number of data points; and (4) a TOD value (negative slope) below 0.8, to exclude conditions with high aerosol loads.

As an example, Fig. 6 shows the ln (USI⁰) satisfying the required criteria for different stars and the four channels of camera C011 during more than 3 years. These ln (USI⁰) values exhibit some fluctuations but tend to remain consistent over time for each star and channel. It is noteworthy that this camera was relocated to Granada in 2024, yet no abrupt changes were observed in the obtained ln (USI⁰) values at this new station. Similar behavior was also observed for camera C005 (not shown), which was installed at various locations. These results suggest that the calibration of the cameras may indeed be independent of the observation site.

Figure 6 Natural logarithm of the extraterrestrial uncalibrated star irradiance (ln (USI⁰)) obtained using the Langley plot method for camera C011, shown for four stars (Altair, Castor, Markab, and Rigel) and the four color channels: red (a), green (b), blue (c), and gray (d). The straight lines represent the median of ln (USI⁰) over a 1-year time window centered on the date shown on the x axis. Dashed vertical lines separate the measurement periods in Valladolid and Granada.

As expected, the ln (USI⁰) values depend on the star, with higher values observed for brighter stars. This trend can also be seen in Table 2, which presents the median and standard deviation of all available ln (USI⁰) values for camera C011. Additionally, this table reveals that the red channel has the lowest availability of ln (USI⁰) values.

Table 2 Median (Md) and standard deviation (σ) of the ln (USI⁰) values for different stars and channels for the entire measurement period of camera C011. The quantity of averaged data (N) is included.

The ln (USI⁰) value of a star for a given date is calculated as the median of all available ln (USI⁰) values within a 1-year time window centered on that date. This approach helps to reduce the uncertainty in the ln (USI⁰) values. The median is chosen over the mean to minimize the influence of outliers. The uncertainty of these smoothed ln (USI⁰) values is calculated as the root sum of the squares of two terms: the propagated uncertainty of each individual ln (USI⁰) value included in the median calculation and the standard deviation of all averaged ln (USI⁰) values.

These smoothed ln (USI⁰) values are also represented as straight lines in Fig. 6. Overall, the values remain largely constant, with minor and gradual time variations, capturing the average behavior of the ln (USI⁰) data. Caution should be exercised near the edges of the time interval, as the absence of half the averaging period may result in unrealistic values (e.g., the end of the red channel for Markab in Fig. 6a). Such anomalies may resolve themselves as additional data become available in the future. This dynamic method of determining the calibration each day is useful for accounting for possible camera degradations or changes in light intensity from variable stars such as Betelgeuse.

Once ln (USI⁰) is determined for each available USI measurement, the TOD is calculated using Eq. (2). The uncertainty in these TOD values is estimated by propagating the uncertainties of both ln (USI⁰) and USI through Eq. (2).

3.5 Effective wavelengths

The TOD values are generally associated with an effective wavelength (λ_eff). In photometry, this effective wavelength is usually calculated by considering the spectral response function of the instrument capturing the direct irradiance. In this work, the effective wavelength of each camera channel is determined as the wavelength expected value assuming as a probability density function the incoming direct spectral star irradiance weighted by the channel's spectral response, as explained by Román et al. (2017) (see Eq. 4 in Román et al., 2012):

$$\begin{array}{}\text{(3)}& {\mathit{\lambda }}_{\mathrm{eff}}={\displaystyle \frac{{\int }_{\mathit{\lambda }}\mathit{\lambda }\cdot S\left(\mathit{\lambda }\right)\cdot I\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}{{\int }_{\mathit{\lambda }}S\left(\mathit{\lambda }\right)\cdot I\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}},\end{array}$$

where λ, is the wavelength, I(λ) is the direct star irradiance reaching the instrument, and S(λ) is the spectral response of the camera channel given by the NST of Fig. 1.

To calculate λ_eff and to ensure realistic values of direct irradiances, the spectral direct transmittance of the Earth's atmosphere is simulated under different conditions using the DISORT radiative transfer model (Stamnes et al., 1988; Buras et al., 2011) included in the libRadtran 2.0.4 radiative transfer package (Mayer and Kylling, 2005; Emde et al., 2016). These conditions are the same as those described in Antuña-Sánchez et al. (2021) and include seven star zenith angle bins (ranging from 10 to 70° in 10° steps), five Ångström exponent (α) values (from 0.2 to 1.8 in steps of 0.4), and five Ångström turbidity (β) values (AOD at 1000 nm, ranging from 0.01 to 0.21 in steps of 0.05). These simulations yield a total of 200 spectral transmittance scenarios.

A realistic incoming star irradiance is obtained by multiplying the extraterrestrial star irradiance spectrum by these transmittance values. However, the extraterrestrial spectrum varies with the star. Therefore, the incoming star irradiance is simulated for 30 different stars, whose extraterrestrial spectra are obtained from the Pulkovo spectrophotometric catalog (Alekseeva et al., 1996). These stars, listed in Table 3, were selected based on their brightness and availability to cover various star types. For each selected star, 200 effective wavelengths (corresponding to the 200 atmospheric conditions) are calculated. The median and standard deviation of these effective wavelengths are presented in Table 3.

Table 3 Median of the calculated effective wavelengths (± standard deviation) for the four color channels of the OMEA-3C and OMEA-3C-TF camera models, calculated for different stars.

Differences in the effective wavelengths among stars are evident, such as Adhara, a blue supergiant with shorter effective wavelengths compared to Gacrux, a red giant. The green and gray channels exhibit similar effective wavelengths, although the gray channel shows a higher deviation due to its broader spectral response. The spectral response differences between the two analyzed camera models (see Sect. 2) also contribute to variations in effective wavelengths, with longer wavelengths observed for the OMEA-3C-TF model, particularly in the red channel. The standard deviation is comparable between the two camera models, except for the red channel, where it is lower for the OMEA-3C-TF model. The median and standard deviation of the 30 effective wavelengths are calculated and presented in the last row of Table 3 (“All”). These values are assumed to represent the effective wavelengths of the camera channels and, consequently, the wavelengths used for the TOD values. This assumption is made regardless of the star being observed, even though variations between stars are acknowledged. The deviation of the effective wavelength is lower when the bandwidth of the camera channels is reduced with a tri-band filter like in the OMEA-3C-TF.

3.6 Contribution of Rayleigh and gas absorption

To derive AOD, the optical depths (ODs) of Rayleigh scattering and gaseous absorption must be subtracted from the TOD. In this work, the gases considered to significantly contribute to the TOD at the camera wavelengths are ozone (O₃), nitrogen dioxide (NO₂), and water vapor (H₂O). The Rayleigh optical depth (ROD) is assumed to be proportional to ground atmospheric pressure.

The spectral optical depth of each of the aforementioned components can be obtained using the Beer–Bouguer–Lambert law. This is achieved by simulating the direct spectral star irradiance at the ground under specific conditions, both with and without the presence of the component, and applying Eq. (4):

$$\begin{array}{}\text{(4)}& {\text{OD}}_x\left(\mathit{\lambda }\right)=m^{-\mathrm{1}}\cdot \ln \left({\displaystyle \frac{I_{\cancel{x}}\left(\mathit{\lambda }\right)}{I\left(\mathit{\lambda }\right)}}\right),\end{array}$$

where x represents the compound for which the optical depth is being calculated, m is the optical air mass, I(λ) is the direct star irradiance at the ground, and $I_{\cancel{x}}\left(\mathit{\lambda }\right)$ is the same I irradiance simulated while neglecting the x compound. The OD values are independent of the selected star. In this work, these spectral OD values are derived using simulations performed with the libRadtran 2.0.4 package using the DISORT model and the REPTRAN parameterization of absorption cross sections of various gases (Gasteiger et al., 2014). The OD is calculated using this method for six different values of each component: total ozone column (TOC) from 200 to 400 DU in 40 DU increments, NO₂ column from 0.05 to 0.425 DU in 0.075 DU increments, total water vapor column (TWVC) from 2.5 to 40 mm in 7.5 mm increments, and the ground atmospheric pressure from 750 to 1000 hPa in 50 hPa increments. When the OD is not explicitly calculated for these parameters, the default values are 300 DU for TOC, 0.25 DU for the NO₂ column, 20 mm for TWVC, and 940 hPa for ground atmospheric pressure.

The calculated OD is valid for extracting the OD at single wavelengths; however, in this case, the camera records the star signal using wider bandwidth filters. Therefore, the OD of each x compound is determined as the expected value of the OD using the ground direct star spectral irradiance weighted by the spectral response of the camera as the probability density function (Eq. 5):

$$\begin{array}{}\text{(5)}& {\mathit{\tau }}_x={\displaystyle \frac{{\int }_{\mathit{\lambda }}{\text{OD}}_x\left(\mathit{\lambda }\right)\cdot S\left(\mathit{\lambda }\right)\cdot I\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}{{\int }_{\mathit{\lambda }}S\left(\mathit{\lambda }\right)\cdot I\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}},\end{array}$$

where τ_x represents the optical depth obtained for the x compound, OD_x is the simulated spectral OD for the x compound, S is the spectral response of the camera channel defined by the NST shown in Fig. 1, and I(λ) is the simulated ground direct star irradiance. These spectral I(λ) values are calculated for the 30 stars considered in Sect. 3 (see Table 3) across seven star zenith angle bins (10 to 70°) and the same five bins for the Ångström exponent and turbidity described in Sect. 3.5. This results in a total of 5250 optical depths for each component and load.

Figure 7 illustrates the calculated optical depth of each compound as a function of its amount for a specific case and as a function of surface atmospheric pressure in the case of Rayleigh scattering. The obtained optical depths exhibit a linear behavior. A least-squares fit is performed for each component and channel, revealing high correlation and y-intercept values close to zero, as expected. This indicates that the slope represents an effective absorption coefficient (or a scattering coefficient in the case of ROD), which can be used to derive the optical depth by multiplying it by the total amount of the compound or by the ground atmospheric pressure in the case of Rayleigh scattering.

Figure 7 Optical depth of ozone (${\mathit{\tau }}_{{\mathrm{O}}_{\mathrm{3}}}$), nitrogen dioxide (${\mathit{\tau }}_{{\mathrm{NO}}_{\mathrm{2}}}$), water vapor (${\mathit{\tau }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$), and Rayleigh scattering (τ_R) as a function of total ozone column (TOC; panel a), NO₂ column (b), total water vapor column (TWVC; panel c), and ground atmospheric pressure (d), respectively. These values are calculated for the OMEA-3C camera model for the star Regulus with a zenith angle of 50° and an Ångström exponent and turbidity of 1 and 0.11, respectively. The data are shown separately for the four color channels, with solid lines indicating the least-squares linear fits.

This slope, referred to as the optical depth factor (ODF) in this study, is calculated for the 5250 available cases with varying zenith angles, stars, and aerosol conditions. The mean ODF values for ozone, nitrogen dioxide, water vapor, and Rayleigh scattering are presented in Table 4 for each channel and both camera models. The standard deviation of the ODF values is also included in Table 4, representing the uncertainty associated with the ODF. As expected, the ODF for water vapor is higher in the red channel, while it is higher in the blue channel for Rayleigh scattering. The uncertainty in the ODF is lower for the OMEA-3C-TF due to its narrower spectral response.

Table 4 Mean values (± standard deviation) of the calculated optical depth factors (ODFs) of ozone (O₃), nitrogen dioxide (NO₂), water vapor (H₂O) and Rayleigh scattering for the four color channels of the OMEA-3C and OMEA-3C-TF camera models.

In this work, the four τ values are calculated for each available TOD value by multiplying each ODF by its corresponding amount component value. The uncertainty of these retrieved values is also propagated. The hourly TOC and TWVC values are directly obtained from the ERA5 reanalysis (Hersbach et al., 2020), specifically from the dataset of ERA5 hourly data on single levels from 1940 to the present (reanalysis-era5-single-levels; https://cds.climate.copernicus.eu/datasets/reanalysis-era5-single-levels, last access: 28 June 2025) as the total_column_ozone and total_column_water_vapour products.

The ground atmospheric pressure is adjusted to the altitude of each station using an air temperature correction applied to the sea level pressure, with both values obtained from the same ERA5 reanalysis collection: hourly 2m_temperature and mean_sea_level_pressure products. All ERA5 data are retrieved from the Atmosphere Data Store (ADS; https://ads.atmosphere.copernicus.eu/, last access: 28 June 2025) using the Climate Data Store (CDS) Application Program Interface.

For nitrogen dioxide, data are sourced from the monthly climatological table described by González et al. (2020), derived from the OMI version 3 dataset (OMNO2d gridded, level 3; Krotkov et al., 2017). All downloaded values (TOC, NO₂ column, TWVC and ground pressure) are temporally interpolated to estimate the gas and Rayleigh optical depths and their uncertainties at the times corresponding to the available TOD values. The uncertainty in these downloaded values is not considered and is therefore not propagated in the final uncertainty of the mentioned optical depths.

3.7 AOD and cloud screening

Once the TOD and its corresponding gas and Rayleigh optical depths are available, the AOD is calculated for each star as the TOD minus each of the four mentioned optical depths (O₃, NO₂, H₂O, and Rayleigh). These AOD values are represented in the left panels of Fig. 8 for several stars and the four effective wavelengths during one night at Izaña.

Figure 8 Aerosol optical depth (AOD) at Izaña (camera C005) from 3 to 4 November 2022 at different wavelengths for the following stars: Adhara, Alhena, Dubhe, Mirfak, Pollux, Regulus, and Vega. Panels (b), (d), (f), and (h) show the corresponding AOD values after filtering for clouds and spurious data. The AOD values are presented for the following wavelengths: 591 nm (AOD₅₉₁; panels a and b), 534 nm (AOD₅₃₄; panels c and d), 466 nm (AOD₄₆₆; panels e and f), and 537 nm (AOD₅₃₇; panels g and h).

The AOD is available for 37 different stars that night, but only 7 randomly chosen stars are shown to improve visualization in the plots; error bars are omitted for the same reason. The AOD values generally overlap in time for different stars. However, the AOD occasionally shows significant deviations for some stars compared to others. For instance, this is well observed in Fig. 8 in the AOD from Mirfak in the middle of the night, especially at 591 nm (red channel) and 534 nm (green channel). These deviations of the AOD from other stars can at least partially be attributed to inaccuracies in the ln (USI⁰) calibration constant. Additionally, the presence of clouds in the star's position can also cause AOD differences between stars; however, such effects should be noticeable across all channels.

Moreover, the AOD for some stars exhibits pronounced oscillations, such as Mirfak in the middle of the night or Pollux at the end. These oscillations primarily occur under conditions of low N_smooth values and low optical air mass. This suggests that the oscillations may be caused by the insufficient number of USI values in the smoothing process described in Sect. 3.3. The lack of data could be partially related to low optical air masses, where the field of view of these pixels is larger due to the fisheye lens projection. This causes starlight to be concentrated in fewer pixels at low zenith angles, leading to the star segment not occupying enough pixels to be considered. Additionally, the higher concentration of starlight makes the signal in these pixels more prone to saturation under these conditions of low optical air mass.

A cloud screening and quality assurance (CS/QA) method is proposed to filter the AOD data series for each star and channel to exclude cloud-contaminated, inaccurate, and spurious data. First, any AOD value whose sum with its uncertainty falls below zero is removed, as it lacks physical meaning. Additionally, only AOD data with N_smooth > 10 are considered. To reduce temporal fluctuations, a running mean filter, similar to that proposed by Pérez-Ramírez et al. (2012), is applied: any AOD value showing an absolute difference greater than 0.02 from the mean AOD calculated within a ±5 min window (avoiding the AOD value itself) is discarded. If no other AOD data exist within the time window, the value is also rejected. This running mean method is applied iteratively until no further data are removed. Finally, to eliminate outliers and spurious data, a 3σ criterion, similar to that described by Giles et al. (2019), is employed: the median and standard deviation of the remaining AOD values are calculated, and any AOD value falling outside the range defined by the median ±3σ (σ being the standard deviation) is rejected. This process is repeated iteratively until no further data are discarded or the standard deviation falls below 0.02.

The right panels of Fig. 8 present the same AOD values as those in the left panels, but filtered using the method explained above. As observed, the filtering method successfully removes AOD values with higher deviations and fluctuations, leaving most of the remaining AOD values close to low values. The data reduction is more pronounced in the red channel. These filtered data are shown for only seven stars, but filtered AOD values from a total of 36 different stars are available for that night, one fewer than the unfiltered data series. This indicates that the AOD data for one star are entirely removed after filtering. That star is Sirius, which is one of the brightest stars and frequently appears saturated, resulting in N_smooth values that are consistently below the required threshold.

The differences in the filtered AOD obtained from different stars are likely to be primarily due to calibration characteristics, although some of these differences could also arise from inhomogeneities in the aerosol spatial distribution and certain assumptions, such as using the same effective wavelengths as well as gas and Rayleigh optical depths for all stars. Assuming spatial homogeneity, the AOD values from different stars are averaged to produce a single AOD value per channel at each time. This averaging process is explained next.

First, for each channel, the filtered AOD data from all stars are grouped into 5 min intervals centered on specific times: on the hour, 5 min, 10 min, and so on. For example, one group is formed with data in the interval 19:57:30 < time ≤ 20:02:30, which corresponds to 20:00:00; the next interval is 20:02:30 < time ≤ 20:07:30, assigned to 20:05:00, and so forth. Then, a 3σ criterion is applied to each AOD group to remove outliers in a similar way as Giles et al. (2019): the AOD median and standard deviation are calculated, and any data point outside the interval median ±3σ is removed. This process is iteratively repeated until no further data points are removed or the standard deviation falls below 0.02. Once outliers are rejected, the weighted mean of the remaining data is calculated, with the weight of each data point being the inverse square of its estimated uncertainty. This weighted mean is then assumed to represent the AOD for the analyzed channel and time interval. The uncertainty on this AOD is calculated by propagating the uncertainty of the averaged AOD data. The standard deviation of the averaged data (σ_5 min), the number of averaged data (N_5 min), and the number of different stars involved in the averaged data (N_{star−5 min}) are stored for quality control purposes.

An AOD value is then obtained for each channel every 5 min. Figure 9 presents these averaged AOD values corresponding to the case depicted in Fig. 8. The AOD remained relatively stable throughout the night, with values around 0.04 across all channels. As expected, the red channel generally exhibits lower AOD values compared to the other channels; however, it also shows larger temporal fluctuations and higher uncertainties.

These averaged data, which are analogous to AERONET level 1.0 data, can still be contaminated by clouds or may lack sufficient confidence. Therefore, an additional CS/QA method is applied to filter these data. First, for quality assurance, AOD values below 0 are rejected due to their nonphysical nature. Additionally, all AOD data with N_5 min ≤ 10, N_{star−5 min} < 2, or σ_5 min ≥ 0.08 are removed, as they indicate an insufficient number of observations or a high level of deviation.

For cloud screening, a methodology similar to AERONET version 3 is applied (Giles et al., 2019). This process considers all the remaining AOD data from a full night and is applied independently to each channel. The procedure consists of the following steps.

• 1.

  Smoothness criterion. If two consecutive AOD points show an AOD difference greater than 0.01 min⁻¹, the higher AOD value of the two is removed. This step is applied iteratively until no further data are removed.

• 2.

  Stand-alone criterion. If no other AOD data remain within a ±1 h time window around a given AOD data point, that data point is removed.

• 3.

  3σ criterion. The standard deviation (σ) of all remaining AOD data in the night is calculated. If σ > 0.02, any AOD values falling outside the interval defined by the median ±3σ are discarded. This step is applied only once (no iterative).

• 4.

  Repeat stand-alone criterion. The stand-alone criterion is reapplied to the remaining data.

• 5.

  Minimum data requirement. If the number of remaining AOD values is below three or less than 10 % of the initial available averaged values, all the data for that night are removed.

This methodology ensures that only high-quality AOD data remain after the screening process, analogous to AERONET level 1.5 data.

This CS/QA method is applied to all available data. In the case of the AOD shown in Fig. 9, the data passing this CS/QA method are displayed in blue, while the removed data are shown in red. Green and gray channels do not exhibit any removed data, and only one data point is removed for the blue channel due to a low N_5 min value of 7. However, the red channel presents several data points removed by the quality assurance criteria of the CS/QA method: some points at the beginning and end of the night exhibit AOD values below zero; the data removed at the beginning of the night, although above zero, and the five consecutive points around 20:20 have N_5 min values ≤ 10, consistent with the observed in Fig. 8b; and the data points at 04:45 and 04:55 are removed because N_{star−5 min} = 1 and σ_5 min = 0.09, respectively. The several subsequent screen algorithms to remove outliers may remove legitimate data points and underestimate the measurement uncertainty.

Figure 9 Aerosol optical depth (AOD) at Izaña (camera C005) from 3 to 4 November 2022 at 591 nm (AOD₅₉₁; panel a), 534 nm (AOD₅₃₄; panel b), 466 nm (AOD₄₆₆; panel c), and 537 nm (AOD₅₃₇; panel d). Error bars indicate plus or minus the AOD uncertainty. AOD values that do not meet the cloud screening and quality assurance (CS/QA) criteria are shown in red and are removed from the final data series represented in blue.

4 Results

The methodology explained in the previous section has been applied to all available camera data (Table 1), yielding cloud-screened and QA-filtered AOD values at four effective wavelengths corresponding to the red, green, blue, and gray camera channels. To evaluate the performance of the proposed methodology, the retrieved AOD values are compared against independent AOD measurements taken at nighttime at the same stations using a CE318-T Sun–sky–Moon photometer.

4.1 AOD comparison with Moon photometer data

The AOD data from the photometer, obtained through CAELIS, do not correspond to the same wavelengths as those retrieved from the camera. To match them, the photometer AOD values at 440, 500, and 675 nm are selected, as these are the closest to the camera wavelengths. Using these three AOD values, the Ångström exponent and the turbidity parameter are calculated. The photometer AOD is then interpolated to the camera wavelengths using these two Ångström parameters.

To temporally match the two datasets (camera and photometer), for each AOD data point from the camera, photometer data within a ±2.5 min window for the same wavelength are selected. These data within this window are averaged, and this mean value is assumed as the reference for comparison with the camera AOD.

Once the two datasets are matched in time and wavelength, they are directly compared for all the different cameras and wavelengths in Fig. 10. The color scale in this figure indicates the number of camera–photometer data pairs found within AOD pixels of size 0.02 × 0.02.

Figure 10 Density plots representing the AOD values retrieved from the cameras against those obtained by CAELIS from the Sun–sky–Moon photometer measurements at nighttime. Each row corresponds to a specific camera, and the columns represent the red, green, blue, and gray channels, respectively. The red, green, blue, and gray channels correspond to effective wavelengths of 614, 538, 472, and 548 nm, respectively, for the C013 camera; the wavelengths are 591, 534, 466, and 537 nm for the remaining cameras. The density is represented by the color bar, which indicates the number of available data points (N_p) in a square bin with a pixel size of 0.02 AOD units (2D histogram). The number of available data points (N), the correlation coefficient (r), the mean bias error (MBE), and the standard deviation of the differences between the two AODs (SD) are included in each panel. The solid magenta line represents the least-squares linear fit.

In general, the AOD from the cameras shows values similar to those from the photometer, with most of the data points lying close to the 1:1 line and exhibiting correlation coefficient values above 0.90 in most cases. However, this correlation is significantly lower for cameras C003 (Marambio) and C004 (Ny-Ålesund), especially in Marambio, where a considerable number of outliers can be observed. These outliers may be due to camera data contaminated by clouds or failing to meet quality criteria that were not properly filtered by the developed algorithms. They may also result from photometer data not properly filtered for clouds, as thin clouds are harder to detect at nighttime without lunar aureole measurements. This issue is evident in panels (a), (b), (c), and (d) of Fig. 10, where high photometer values appear in the bottom-right corner for low camera AOD values. Although this problem with outliers is common across all cameras (e.g., panel i), it is more pronounced in Marambio due to the high cloud presence and generally low AOD values at this station. In the case of the camera C004 in Ny-Ålesund, the differences between the camera and the photometer may partly be due to the Langley calibration method, which is ideal for low-latitude regions but performs worse at high latitudes due to slower changes in the optical air mass. This means that the time required to perform a Langley plot is longer at high latitudes, which increases the probability of experiencing some atmospheric variation during this extended period, potentially reducing the accuracy of the calibration.

For a more detailed quantification of the differences between the two datasets, Table 5 presents several statistical indices comparing the two datasets, using the photometer as the reference. Generally, fewer AOD data points are available from the camera for the blue and especially the red channels, a pattern also observed in Fig. 9. The correlation coefficient is also generally lower for these channels. Table 5 presents the mean (mean bias error; MBE), median (Md), and standard deviation (SD) of the differences (Δ) between the AOD from the cameras and the photometer. The MBE, which represents the accuracy of the camera AOD, ranges between 0 and 0.05, with the highest values observed for cameras C038 (Davos) and C003 (Marambio) and the lowest values for cameras C011 and C009 (Lindenberg). In general, the camera AOD overestimates the photometer AOD for all cameras, with this overestimation being approximately 0.02 for all wavelengths when considering data from all cameras.

Table 5 Statistical estimators of the comparison between the AOD obtained from cameras and photometers for different cameras and channels. N represents the number of available data pairs, r is the correlation coefficient, and b and y₀ are the slope and y intercept of the least-squares linear fit between the camera and photometer AODs. Regarding the differences between the AOD from the cameras and photometers (Δ), the table includes the mean (mean bias error, MBE), the standard deviation (SD), and the percentage of Δ absolute values that fall within 1 time the sum of the uncertainties of both AODs (Δ < σ) and within 2 times this sum (Δ < 2σ).

The SD values, which represent precision ranges from 0.02 to 0.06, except for camera C003 at Marambio, which shows values of 0.09 for blue channel, likely due to the aforementioned outliers. The lowest SD values are found for camera C013, which may be linked to its quasi-monochromatic filters. The blue channel generally shows the worst precision, with higher SD values of approximately 0.05 when considering all cameras, while this SD is close to 0.03 for the other wavelengths.

The slope (b) and y intercept (y₀) of the linear fit shown in Fig. 10 are also presented in Table 5. The y₀ values are low but positive, and most b values are below 1, indicating a linear fit that tends to overestimate low AOD values while underestimating high ones.

Table 5 also shows the percentage of the absolute values of Δ that fall below the standard uncertainty (Δ < σ) and the expanded uncertainty (Δ < 2σ). The standard uncertainty is defined as the square root of the sum of the squares of the errors of the camera and photometer AODs, while the expanded uncertainty is twice the standard uncertainty. In this case, the uncertainty associated with the photometer AOD is assumed to be 0.03 based on the results of Román et al. (2020), but it could actually be lower. If the uncertainty associated with the camera AOD properly represents its error as a Gaussian probability density function with the standard deviation equal to the uncertainty, then the values of Δ < σ and Δ < 2σ should be approximately 68 % and 95 %, respectively. The values of Δ < σ vary depending on the camera and wavelength, as shown in Table 5. However, in most cases, except for camera C003, the Δ < 2σ values are close to the mentioned 95 %. In fact, when considering data from all cameras combined, the Δ < 2σ values range between 91.3 % and 96.2 %, indicating that the estimated AOD uncertainty accurately reproduces the real uncertainty in the camera AOD.

Cameras C005 and C011 were installed in various locations; however, the results in Table 5 do not account for this detail. To evaluate the performance of these cameras in each location, Table 6 presents the same statistical indices as Table 5, but distinguishing between the locations for cameras C005 and C011.

Table 6 Same values as in Table 5, but for the different locations of cameras C005 and C011 separately.

The r values are generally lower for Andøya compared to the other locations of camera C005. Similarly, for camera C011, r values for Valladolid are lower than those for Granada. Differences between locations also appear in the MBE, with the most significant difference observed in the blue channel of camera C011, where the MBE in Granada reaches the lowest value across all cameras (−0.01). SD values for Izaña are the lowest among all cameras and locations. This improved precision may be related to Izaña being an ideal site for Langley calibrations due to its low latitude and high altitude, which provide a clean and stable atmosphere (Toledano et al., 2018). On the other hand, the SD values in Andøya are higher than in other locations, likely because Andøya is a less suitable site for Langley calibration due to its higher latitude. This variation in SD between locations observed for camera C005 contrasts with the similar SD values obtained at the two locations (Valladolid and Granada) for camera C011. The Δ < 2σ values remain close to the expected 95 % for all locations.

Therefore, the performance of the cameras in estimating AOD appears to depend on the location where they are installed. On the other hand, the differences in the results from Table 5 between cameras C006 and C009, both located in Lindenberg, are not as pronounced and may be caused by the insufficient number of data available for C006. Overall, it seems that the main differences observed between the various cameras are attributable to the location where they are installed rather than the camera itself, at least when considering OMEA-3C cameras (all except C013). Although some of the observed differences may be due to differences in the available data, the differences between locations are most likely related to how suitable each site is for obtaining highly accurate values of ln (USI⁰) through the Langley technique. This depends on whether the site exhibits a stable atmosphere and is geographically located in a region where the star optical air masses change rapidly enough for the atmospheric conditions to remain nearly constant during the Langley calibration (typically at low latitudes and midlatitudes).

4.2 Case study

In the previous comparison (Sect. 4.1), the agreement between the AOD from the cameras and the photometer at nighttime is discussed. However, the photometer AOD at nighttime is only available when the Moon is visible and between the first and last quarters (less than half of the time), while the camera AOD is available throughout the entire night since stars are always visible in the sky. To evaluate the goodness of the AOD obtained from the cameras during Moonless nighttime, Fig. 11 shows the time series of the camera AOD over 20 d in Izaña, along with the daytime AOD obtained with the photometer, to assess whether the day-to-night transition in the AOD is consistent. During this period, the nighttime AOD values from the photometer, when the Moon was visible, are also included.

Figure 11 AOD at (a) 591 nm (AOD₅₉₁), (b) 534 nm (AOD₅₃₄), (c) 466 nm (AOD₄₆₆), and (d) 537 nm (AOD₅₃₇) at Izaña for 20 d in July 2022. Blue error bars represent the nighttime camera AOD ± its associated error. The red crosses and green dots correspond to the AOD data from the Sun–sky–Moon photometer during daytime and nighttime, respectively.

During the period shown, it can be observed how the daytime AOD varied throughout the days, ranging from values close to zero to values above 0.4. The increases in AOD above the typical low values for this station were caused by the arrival of Saharan dust from North Africa. As can be observed for all wavelengths, the nighttime AOD evolution aligns well with the daytime values, even when Moon data are not available. For instance, the increase in daytime AOD from the end of 7 July to the early morning of 8 July is well captured by the camera AOD, which shows a clear increase starting just after midnight. A similar behavior is observed between 21 and 22 July. In general, AOD from camera correlates with the expected values even under Moonless conditions.

One of the poorest performances of the camera AODs against the photometer occurs in Ny-Ålesund (C004), where half of the year is continuously dark. As an example, Fig. 12 shows the AOD over 5 consecutive days in December 2023 at this location. First, it can be observed that there are no daytime data, while nighttime data are continuously available due to the time of year. On 22, 23, and 24 December, AOD values around 0.1 are observed for all wavelengths, showing a smooth and fairly constant temporal evolution. A similar trend can be seen in the photometer data. Although the camera AOD shows some differences compared to the reference AOD obtained from the photometer, it is still useful for observing the temporal evolution and magnitude of the AOD. Notably, the camera AOD at 466 nm shows lower values compared to the photometer and the other wavelengths. This could be caused by inaccuracies in the Langley calibrations, which, as mentioned above, are less accurate when conducted at high latitudes.

Figure 12 AOD at (a) 591 nm (AOD₅₉₁), (b) 534 nm (AOD₅₃₄), (c) 466 nm (AOD₄₆₆), and (d) 537 nm (AOD₅₃₇) at Ny-Ålesund for 5 d in December 2023. Blue error bars represent the nighttime camera AOD ± its associated error. The green dots correspond to the AOD data from the Sun–sky–Moon photometer during nighttime.

The data from the lidar (not shown) installed in Ny-Ålesund reveal the presence of polar stratospheric clouds (PSCs) between 21 and 23 km altitude from 07:00 UTC on 23 December to 19:00 UTC on 24 December. These clouds have a very low optical depth, which prevents the cloud screening algorithm from filtering out AOD data in their presence. Therefore, part of the AOD observed in Fig. 12 is influenced by these PSCs. These clouds are also not filtered out by the photometer, as high clouds are typically screened during the day using solar aureole radiance measurements (Giles et al., 2019). However, at night, lunar aureole measurements are not available for this cloud screening. Additionally, lidar data indicate that the last AOD measurements on 25 December are likely contaminated by high ice clouds, such as cirrus, which are known to be more challenging to filter. The same issue may also affect the first AOD measurements on 21 December, although no lidar data are available to confirm this. These findings highlight that conditions during the polar night can indeed be variable, and cloud screening algorithms could be further improved to better detect high and thin clouds.

5 Conclusions

This study has demonstrated the potential of using all-sky cameras to estimate atmospheric properties such as the total optical depth (TOD) and the aerosol optical depth (AOD) under nighttime conditions. This paper has proposed a novel method for extracting the starlight reaching the Earth for different stars from nighttime sky images captured by an all-sky camera. This recorded starlight is equivalent to an uncalibrated irradiance, and then it can be used for star photometry, transforming this starlight into the AOD.

Langley calibration has been used to calculate the TOD and, after properly accounting for the effects of gas scattering and absorption, to retrieve the AOD at four different effective wavelengths (red, green, blue, and gray camera channels). A comparison against independent and accurate photometer measurements at nine different locations has revealed that the AOD values retrieved from the cameras correlate with the photometers ones, with correlation coefficients exceeding 0.9 in many cases. Additionally, when considering data from all stations combined, the mean bias error (accuracy) and the standard deviation (precision) of the camera–photometer differences are approximately 0.02 and between 0.03 and 0.04, respectively.

However, the worst agreement appears for higher-latitude locations, likely related to the limitations of Langley calibrations in such regions. This suggests that alternative methods or cross-calibration with other stations may be necessary for such areas. Additionally, the Langley calibration is calculated as an average over a ±6-month window, meaning it is updated when new data become available after at least 6 months. This highlights the importance of exercising extra caution when interpreting data from the most recent 6 months.

The proposed methodology to extract starlight, identify valid Langleys, and filter cloudy and low-quality data is based on a series of approaches that use several threshold values manually established. The method and thresholds are assumed to be the same and independent of the measurement station. However, these thresholds may need to be adjusted differently depending on the measurement station for optimal performance. Additionally, the list of stars could be refined and tailored to each specific location. Moreover, for all the measurements an all-sky image is available, where the presence of clouds can be observed globally, and it could be used for cloud screening in the future.

The use of RGB cameras allows obtaining starlight in different spectral bands and, therefore, at various effective wavelengths. However, since stars are point sources and the camera sensor has a mosaic of color filters, such as the RGGB Bayer pattern, this influences the recorded starlight signal. Additionally, the standard camera color filters are usually spectrally wide, which introduces greater inaccuracies in photometry purposes, such as calculating an effective gaseous absorption. Reducing this spectral width improves the precision of the AOD, as observed when comparing the performance of different camera models in Valladolid. With all this in mind, an improvement to the method could involve using cameras without a filter mosaic and with narrower spectral responses, such as incorporating an interference filter like those used in photometers. While this would address the highlighted issues, it is important to consider that the amount of starlight captured by the camera would be significantly lower. This would require longer exposure times and/or higher amplification gains, leading to fewer images and/or an increase in noise.

Although the use of a camera with a color filter mosaic has enabled the calculation of AOD at different effective wavelengths, the small spectral range (less than 150 nm) and the inherent uncertainty in the calculated AOD do not allow for a reliable estimation of the Ångström exponent.

Even with these limitations, the proposed methodology provides AOD values similar to actual ones, serving at least as a good proxy. This method has proven useful for detecting and monitoring temporal changes in AOD, such as those caused by Saharan dust intrusions, highlighting its potential for studying aerosol dynamics. This, combined with the ability to continuously monitor AOD throughout the night, regardless of moonlight, offers an unprecedented opportunity to fill the gap in nighttime AOD measurements without being restricted to lunar visibility.

Future research should focus on improving data filtering algorithms and calibration methods for high-latitude stations, testing other camera models, and further validating the proposed methodology under different atmospheric conditions. In this regard, the proposed methodology should also be compared with measurements from star photometers, such as those installed in Lindenberg or Ny-Ålesund, in order to assess the performance of the AOD retrieved with the all-sky cameras during Moonless periods. This comparison could also benefit from the cloud screening algorithms for ice clouds that have been developed for these instruments (Baibakov et al., 2015; O'Neill et al., 2016). We encourage scientific researchers involved in atmospheric studies to use the data generated from this work and to adopt the proposed methodology (either by joining GOA-SCAN or implementing it themselves) to calculate AOD with their own all-sky cameras. Additionally, we invite other researchers to compare the methodology presented in this work with measurements performed in other locations and to work on improving it so the entire community can benefit from advancements in this promising research line. Finally, we recommend that all-sky camera manufacturers consider these new proposed ideas to adapt their technical developments to the specific needs of this new application.

Appendix A

Acronyms
AEMet	Agencia Estatal de Meteorología: Meteorological State Agency of Spain
AERONET	AErosol RObotic NETwork
AOD	Aerosol optical depth
ASC	Andøya Space Center, Andøya, Norway
AWI	Alfred Wegener Institute: Helmholtz Centre for Polar and Marine Research (German Polar Institute)
AWIPEV	Alfred Wegener Institute Paul Emile Victor: AWI and IPEV German–French Arctic Research Station (also called “Koldewey station”) in Ny-Å, Svalbard
BCPS	Background-corrected pixel signal
CS	Cloud screening
DC	Digital counts
DU	Dobson units
DWD	Deutscher Wetterdienst: German Meteorological Service
GFAT	Grupo de Física de la Atmósfera de la Universidad de Granada (UGR): Atmospheric Physics Group of the University of Granada (UGR), Granada, Spain
GOA	Grupo de Óptica Atmosférica: Atmospheric Optics Group of UVa, Valladolid, Spain
GOA-SCAN	GOA all-Sky CAmeras Network
Gr	Grayscale
GRASP	Generalized Retrieval of Atmosphere and Surface Properties, Lille, France
HDR	High dynamic range
IARC	Izaña Atmospheric Research Center: AEMET Observatory of Izaña, Tenerife, Spain
IISTA-CEAMA	Instituto Interuniversitario de Investigación del Sistema Tierra – Centro Andaluz de Medio Ambiente: Andalusian Institute for Earth System Research – Andalusian Center for the Environment, UGR, Granada, Spain
IPEV	Institut polaire français Paul Émile Victor: French Polar Institute
IR	Infrared
MBE	Mean bias error
Md	Median
MOL-RAO	Meteorologisches Observatorium Lindenberg – Richard-Aßmann-Observatorium: DWD Observatory Lindenberg, Lindenberg (Tauche), Germany
NST	Normalized Spectral Transmittance
Ny-Å	Ny-Ålesund, Spitzbergen, Svalbard
OD	Optical depth
ODF	Optical depth factor
PMOD/WRC	Physikalisch-Meteorologisches Observatorium Davos/World Radiation Center, Davos, Switzerland
QA	Quality assurance
RGB	Red, green, blue
RIMO	ROLO Implementation for Moon photometry Observation: open lunar reflectance model developed by GOA and AEMET/IARC
ROD	Rayleigh optical depth
ROLO	Robotic Lunar Observatory: lunar reflectance model of the United States Geological Survey (USGS)
SD	Standard deviation
SMN	Servicio Meteorológico Nacional: National Meteorological Service of Argentina
SS	Star signal
SZA	Solar zenith angle
TF	Tri-band filter
TOC	Total ozone column (in DU)
TOD	Total optical depth
TWVC	Total water vapor column
UGR	Universidad de Granada, Granada, Spain
USI	Uncalibrated star irradiance
USI0	Uncalibrated extraterrestrial star irradiance
UTC	Coordinated Universal Time
UVa	Universidad de Valladolid, Valladolid, Spain

Paramater	Description
α	Ångström exponent
β	Ångström turbidity
τ	Total optical depth (TOD)
λ	Wavelength
σ	Standard deviation (SD)
Δ	AOD differences
λeff	Effective wavelength
I(λ)	Spectral irradiance (in W m−2 nm−1)
m	Optical air mass
Md	Median
rw	Weighted correlation coefficient
r	Correlation coefficient
S(λ)	Spectral response