Evolution of Brightness and Color of the Night Sky in Madrid

Major schemes to replace other streetlight technologies with Light-Emitting Diode (LED) lamps are being undertaken across much of the world. This is predicted to have important consequences for nighttime sky brightness and color. Here, we report the results of a long-term study of these characteristics focused on the skies above Madrid. The sky brightness and color monitoring station at Universidad Complutense de Madrid (inside the city) collected Johnson B, V, and R sky brightness data, Sky Quality Meter (SQM), and Telescope Encoder Sky Sensor-WiFi (TESS-W) broadband photometry throughout the night, every night between 2010–2020. Our analysis includes a data filtering process that can be used with other similar sky brightness monitoring data. Major changes in sky brightness and color took place during 2015–2016, when a sizable fraction of the streetlamps in Madrid changed from High-Pressure Sodium (HPS) to LEDs. The sky brightness detected in the Johnson B band darkened by 14% from 2011 to 2015 and brightened by 32% from 2015 to 2019.


Introduction
The nighttime environment is being altered across increasingly large areas of the Earth as a consequence of the introduction of artificial lighting. Indeed, the vast majority of the human population now lives under light-polluted skies [1]. Such light pollution is exacerbated by the widespread failure to limit the use of artificial lighting to the places, forms, and times in which it is required to deliver the intended human benefits. Artificial lighting is, for example, often excessively bright, and is allowed to spread into areas (spatial trespass) and times (temporal trespass) beyond those in which it is actually needed. This has a wide array of consequences for living organisms, including measurable effects on human health [2,3] and on wild plants and animals [2][3][4], many of which arise from the impacts of artificial nighttime lighting on the timings of biological activities [5].
Given these environmental impacts, it is important to understand the effects of widespread and ongoing street lighting changes from legacy (mainly High-Pressure Sodium; HPS) to Light-Emitting Diode (LED) technology [6][7][8]. A common method for assessing the extent and impact of light pollution is using remote sensing and ground-based techniques to measure night sky brightness at the zenith as a proxy. Whilst the consequences of street lighting changes for direct light emissions are relatively straightforward to determine, this is more challenging for the brightness and color of the night sky (skyglow). First, although skyglow can be detected from nighttime satellite imagery (e.g., Defense Meteorological Satellite Program/Operational Line-Scan (DMSP/OLS) and Suomi National Polar Orbiting Partnership/Visible Infrared Imaging Radiometer Suite-Day-Night Band (SNPP/VIIRS-DNB); [9,10]), the associated sensors are broadband and do not detect blue light. Thus, nighttime satellite imagery does not provide color information. This is important both because these colors can change with lighting technology and because emissions of different colors can influence skyglow in different ways. Indeed, concerns have arisen that temporal trends in the emissions measured using such imagery can be misleading if the consequences of shifts in the spectrum of artificial light emissions are not accounted for [11]. Second, whilst ground-based sensors are available for measuring skyglow (e.g., Sky Quality Meter (SQM); STARS4ALL Telescope Encoder Sky Sensor-WiFi (TESS-W) photometer [12]), and networks of such devices are being developed, these also tend to be broadband and, again, proper adjustments need to be made for spectral shifts to evaluate temporal trends [6,11,13]. SQMs, broadband photometers used for network monitoring stations, are an example of a device in need of such adjustments. Sánchez de Miguel et al. [11] showed that properly assessing night sky quality with SQMs needed adjustments via color indices.
The contributions of both night sky brightness and color have been recognized by researchers studying light pollution. Kyba et al. [14] used color as an indicator of light pollution and pointed out the need for photometer measuring in several bands to monitor sky color evolution. The authors used standard SQMs with red, green, and blue filters to monitor night sky color evolution, however, the filters were not the standard ones used in astronomy. In particular, the red filter did not properly cover the HPS spectra produced by streetlights. More recently, Barentine et al. [7] used several different data collection instruments, including a Sky Quality Meter (SQM-L) and an unmodified Canon T5i DLSR, to measure night sky brightness changes in Tucson, AZ (USA) following a street lighting retrofit from a mixture of high-and low-pressure sodium (HPS/LPS) streetlamps to fully shielded 3000 K LED. They reported a decrease in skyglow, however, they did not carry out a color analysis.
In this paper, we report how the brightness and color of the night sky changed following a major retrofit of the street lighting in the urban area of Madrid, Spain. Madrid is well-suited for such research, as it has been a focus for multiple previous studies of night sky brightness and light pollution [9,11,[15][16][17][18][19][20][21][22][23]. Here, we used a professional astronomical device (All-Sky Transmission Monitoring Camera; AstMon) that provides information in the Johnson B, V, and R photometric bands along with their associated color, and was operated for a long period (a decade). Finally, we provide a framework that could be used for a collective analysis with data from the International Space Station (ISS) and the Spectrometer for Aerosol Night Detection (SAND) spectrometer to better understand the general problem of the spatial reach of artificial light scattering from cities [10].

Background
Madrid is a densely populated urban center with approximately 3.3 million inhabitants and with about 5.4 million inhabitants within a 27 km radius from downtown [23]. The entire Madrid region includes 6.6 million inhabitants. It has four other province capitals (Avila, Segovia, Toledo, and Guadalajara), which act as satellite cities. The entire Madrid region retrofitted their streetlights. Therefore, the skyglow changes in the Madrid region cannot be explained by the skyglow changes within the city of Madrid alone, although the city is the primary source of light pollution.
It is relatively straightforward to determine the impact of a streetlight retrofit on night sky brightness in a small, isolated town that changes its public street illumination over the course of, say, a week. Using monitoring photometers before, during, and after the streetlight retrofit, one can easily determine the resultant change in brightness. Large cities, however, present a more complicated case, as retrofits have inevitably to be performed in phases. Madrid has been upgrading its street lighting system from HPS to LED technology since 2011, and this process is ongoing. Like many other cities, it has also been implementing other street lighting changes (e.g., changing part of its street lighting system from high power HPS to lower power HPS, testing and implementing regulations to reduce the intensity of street lighting by some percentage after midnight).
implementing other street lighting changes (e.g., changing part of its street lighting system from high power HPS to lower power HPS, testing and implementing regulations to reduce the intensity of street lighting by some percentage after midnight).
We will focus on events before, during, and after the major streetlight retrofit that occurred in Madrid in 2015. This year was notable in the extent of the upgrade. Photographs taken by astronauts aboard the International Space Station (ISS) before and after the main changes in street illumination show visible variation in color from 2011, when street lighting in Madrid was dominated by HPS lamps, to 2017, when street lighting in some areas changed to primarily white LEDs, following the 2015 retrofit ( Figure 1). For example, Rivas Vaciamadrid, at the lower right of Figure 1, completely changed from HPS to LED. Changes to both, manufacturing technology and lamp power ( Figure 2) occurred before and after the major streetlight retrofit in 2015. Changes included replacing HPS lamps bulbs with lower power LED bulbs (from 250W/150W to 150W/70W). Additionally, a portion of 150W HPS globe lamps were changed to 56W demi-globe LED lamps. The diffusers were different between these demi-globe lamps, the 150W had a 0% Upward Light Output Ratio (ULOR), and the 56W LED lamps had a 5% ULOR. The ISS image of 2017 ( Figure 1, right panel) shows an apparent decrease in street brightness due to the change in power of HPS lamps (Right panel, Figure 2). It should be noted that there was not a complete change from HPS to LED lamps. The streetlight retrofits were made in steps throughout the year, and the fraction of technology change was inconsistent from year to year ( Figure  2). According to the Madrid City Hall, there is a reduction in light power after midnight when HPS street lights are dimmed to 80% power and LEDs to 60%. Changes to both, manufacturing technology and lamp power ( Figure 2) occurred before and after the major streetlight retrofit in 2015. Changes included replacing HPS lamps bulbs with lower power LED bulbs (from 250 W/150 W to 150 W/70 W). Additionally, a portion of 150 W HPS globe lamps were changed to 56 W demi-globe LED lamps. The diffusers were different between these demi-globe lamps, the 150 W had a 0% Upward Light Output Ratio (ULOR), and the 56W LED lamps had a 5% ULOR. The ISS image of 2017 ( Figure 1, right panel) shows an apparent decrease in street brightness due to the change in power of HPS lamps (Right panel, Figure 2). It should be noted that there was not a complete change from HPS to LED lamps. The streetlight retrofits were made in steps throughout the year, and the fraction of technology change was inconsistent from year to year ( Figure 2). According to the Madrid City Hall, there is a reduction in light power after midnight when HPS street lights are dimmed to 80% power and LEDs to 60%.

Available Data
The data used for this work were obtained via observations from fixed monitoring instrumentation on top of the Physics building of Universidad Complutense de Madrid (UCM Observatorio UCM, IAU-MPC code I86, 40°27′04″N 03°43′34″W) in Ciudad Universitaria (red circle in Figure 1), northwest of Madrid's city center ( Figure 3). The location of this observatory within a metropolitan area makes it possible to study the high level of light pollution of the surrounding area, thus, making it a unique place dedicated to monitoring night sky brightness with first-class instrumentation.

Available Data
The data used for this work were obtained via observations from fixed monitoring instrumentation on top of the Physics building of Universidad Complutense de Madrid (UCM Observatorio UCM, IAU-MPC code I86, 40 • 27 04 N 03 • 43 34 W) in Ciudad Universitaria (red circle in Figure 1), northwest of Madrid's city center ( Figure 3). The location of this observatory within a metropolitan area makes it possible to study the high level of light pollution of the surrounding area, thus, making it a unique place dedicated to monitoring night sky brightness with first-class instrumentation.

Available Data
The data used for this work were obtained via observations from fixed monitoring instrumentation on top of the Physics building of Universidad Complutense de Madrid (UCM Observatorio UCM, IAU-MPC code I86, 40°27′04″N 03°43′34″W) in Ciudad Universitaria (red circle in Figure 1), northwest of Madrid's city center ( Figure 3). The location of this observatory within a metropolitan area makes it possible to study the high level of light pollution of the surrounding area, thus, making it a unique place dedicated to monitoring night sky brightness with first-class instrumentation.  The UCM astronomical observatory uses two Sky Quality Meters (SQM-LE and SQM-LU), and more recently, two TESS-W photometers, continuously to monitor the night sky brightness of Madrid. The SQM-LE connects to a computer with a Local Area Network, LAN wire, while the SQM-LU connects with a USB connection. TESS-W is a photometer equipped with a dichroic filter that allows an extension of the spectral response to the red region of the spectrum. Thus, the TESS-W photometer is better suited than the SQM for monitoring polluted night sky brightness variation. While these photometers have a broad single photometric band (panchromatic detector), the UCM All-Sky Transmission Monitor (AstMon) provides all-sky maps in the Johnson B, V, and R astronomical photometric bands [25]. The UCM astronomical observatory also hosts a Spectrometer for Aerosol Night Detection (SAND) for monitoring the night sky spectra [21]. Figure 4 demonstrates the detection range of our instrumentation. It shows a comparison of the mean night sky spectra after midnight in January 2014 (before the retrofit) and for the same month in 2019 (after the retrofit). In the top panel, we plot the detection range available by TESS-W (note that the spectral response extends to the red region) and SQM, along with the spectra obtained by the SAND spectrometer. In the bottom panel, we present the Johnson B, V, and R photometric bands, which are within AstMon's spectral range. For reference, we include the prominent emission lines of typical lamps (middle panel). The UCM astronomical observatory uses two Sky Quality Meters (SQM-LE and SQM-LU), and more recently, two TESS-W photometers, continuously to monitor the night sky brightness of Madrid. The SQM-LE connects to a computer with a Local Area Network, LAN wire, while the SQM-LU connects with a USB connection. TESS-W is a photometer equipped with a dichroic filter that allows an extension of the spectral response to the red region of the spectrum. Thus, the TESS-W photometer is better suited than the SQM for monitoring polluted night sky brightness variation. While these photometers have a broad single photometric band (panchromatic detector), the UCM All-Sky Transmission Monitor (AstMon) provides all-sky maps in the Johnson B, V, and R astronomical photometric bands [25]. The UCM astronomical observatory also hosts a Spectrometer for Aerosol Night Detection (SAND) for monitoring the night sky spectra [21]. Figure 4 demonstrates the detection range of our instrumentation. It shows a comparison of the mean night sky spectra after midnight in January 2014 (before the retrofit) and for the same month in 2019 (after the retrofit). In the top panel, we plot the detection range available by TESS-W (note that the spectral response extends to the red region) and SQM, along with the spectra obtained by the SAND spectrometer. In the bottom panel, we present the Johnson B, V, and R photometric bands, which are within AstMon's spectral range. For reference, we include the prominent emission lines of typical lamps (middle panel).  Data collection has been ongoing since the installation of each device (Table 1; [18][19][20] 2013), and data from Aerosol Optical Depth (AOD) from AERONET (since 2012). To our knowledge, there is currently no other such systematic information accounting for the evolution of night sky brightness in Johnson B, V, and R astronomical photometric bands and its associated color for a city with a large population. A detailed analysis of the evolution of the spectra provided by SAND and the AOD data collected by AERONET is beyond the scope of this paper and will be presented elsewhere.

Data Management and Pre-Processing
The data from the SQMs were recorded with the freeware software PySQM [26], while TESS-W photometers use the Internet of Things (IoT) protocols to insert the data in realtime via STARS4ALL repositories. Both photometers point to the zenith. Their data were archived using the standard International Dark-Sky Association (IDA) format [13].
Night sky brightness data from the three photometers were filtered to remove values that were too bright, too dark, or otherwise unreliable, leaving values in the range 14.0 < m < 19.0 mag/arcsec 2 . The resulting files contain records with the date-time and sky brightness in their respective bands for each measure, typically one value per minute.
AstMon proprietary software provides calibrated all-sky images along with sky brightness values for several locations in the sky. In this work, we used sky brightness at the zenith (3.8 arcminutes pixel scale) in order to compare with the SQM (Field of View (FoV) = 20 • ) and TESS-W (FoV = 17 • ) photometer data. The observational procedure of AstMon at UCM consists of a sequence of Johnson B, V, and R measurements (in that order) followed by a small interval, after which the sequence repeats. AstMon at UCM only observes during the astronomical night. AstMon data are stored in simple ASCII files whose records contain date and time of observation (Julian dates) and zenith sky brightness in the corresponding B, V, and R bands. Color indices (differences in magnitudes from two photometric bands) were derived later to determine changes in the color of the sky. In order to determine B-V and V-R color indices, the data have been grouped according to the observation time. The time of the V observation has been assigned to every B, V, and R sequence. An ASCII file containing date-time and color index data was created for the whole series.

Filtering Process
The sky brightness at the zenith measured at the UCM astronomical observatory is a proxy for light pollution in the Madrid urban area. The relative contribution of different light sources depends on their intensity and distance from the photometers. Our purpose is to determine the changes in the number, intensity, and color of pollution sources using night sky brightness values. There are other factors that increase the brightness of the sky, particularly the presence of the moon over the horizon and variation in atmospheric conditions, particularly cloud cover and aerosol content. In this paper, we address the long-term effects of night sky brightness, produced by changes in street lighting systems. Effects produced by factors such as aerosols and atmospheric conditions are often short-term Remote Sens. 2021, 13, 1511 7 of 25 and can be interrelated (such as wind diluting aerosol content) and will be studied in a future paper.
One approach to avoid the influence of clouds would be to select only the photometric nights in the astronomical sense: transparent nights with good conditions during the whole night. This procedure drastically reduces the dataset and throws away a large portion of otherwise useful data, including observations with clear skies during some intervals and non-photometric nights. Another approach is to use the complete dataset and perform an automatic filtering process to eliminate observations taken with cloud cover over the observatory. Since the monitoring station is located inside a highly polluted city, any cloud will increase the brightness. It would be straightforward to write code to filter out unexpected bright values by comparing them with previous and posterior observations. In our case, we elected to filter the data by visual inspection. This procedure is lengthy and tedious but yielded a filtered dataset that we could better trust. Our method offers the opportunity to visualize all the data, observe the variation in sky brightness throughout the night for different sky conditions, and learn about the effects of clouds and the presence of the moon on sky brightness.
Each record of night sky brightness data collected by AstMon and the photometers was categorized into three conditions: clear sky without the moon (no moon, no clouds), clouds during the observation, and data obtained with the moon over the horizon. The procedure is as follows: (a) For each night, a graph is displayed showing the Johnson B, V, and R photometric data throughout the night (x-axis). The magnitude of the broadband provided by the photometers is displayed as a comparison. A second panel of the graph shows the variation of the B-V and V-R color indices. The Julian date when the moon is over the horizon is calculated and the corresponding records flagged. The data with the moon over the horizon are displayed with a different marker. (b) After visual inspection, the data suspected to include clouds were flagged in the dataset. (c) The data are displayed again to ensure that the process was done properly. After the process, all the records of the dataset (both magnitudes and colors) were flagged, and we performed statistics on the observations that were free of external influences (no moon, no clouds). Figure 5 shows a typical graph of the variation throughout one night with observations during a cloudy sky and the presence of moonlight labeled. Figure 6 shows the graph for a photometric night without the moon over the horizon. During the first part of the astronomical night (after twilight), the sky brightness decreases due to the progressive reduction in human activity, including traffic, occupation of office buildings and shops being open, among other anthropogenic factors. A step around 23 h is due to the switching-off of ornamental lights. After midnight, most of the light sources are those from public street illumination. According to Madrid city hall, the intensity of a fraction of the street lamps is reduced in the late evening, but the number has been inconsistent over time. The brightening from 04:45 h to 05:30 h is the result of street illumination returning to normal operational intensity (increase in electric power) before citizens leave for the daily commute. Finally, the end of the plot shows the twilight as recorded by SQM-LE/SQM-LU. The lower panel shows the B-V and V-R color indices. Lower values represent bluer skies. The sky is bluer at the beginning of the night due to the types of lights used in traffic, ornamentation, sports stadiums, and other light pollution sources, which are generally bluer. After midnight, the color reflects street illumination showing low night sky brightness variability. The observations during this part of the night are better suited to determine the nature of the lamps used for street illumination.
After the filtering process, we used 119,296 measurements (no moon, no clouds), accounting for 28% of AstMon total observations. As regards the photometer data, we used 283,922 measurements, accounting for 41% of the SQM measurements, and 234,543 measurements, accounting for 39% of the TESS-W measurements. The percentage of useful data for this research (no moon, no clouds) depends on the experiment starting date, device operability, and the hourly instrument sampling of night sky brightness (Table 1).

Madrid Night Sky Spectra
In Figure 4, we present a graph comparing the mean night sky spectra after midnight in January 2014 and January 2019. Examining this figure, we see that the 2019 night spectrum shows the blue bump of the LED spectrum. The Johnson B, V, and R photometric band responses show how the changes in spectra are related to the change in the sky's color.

AstMon
In Figure 7, we show histograms of the Johnson B, V, and R night sky brightness and color indices at the zenith measured using AstMon. A progressive brightening of the Madrid sky is detected in the Johnson B photometric band. The Johnson B, V, and R band brightness values for nights with both moon and clouds (lighter bars) were found to be to the left of values for moonless and cloudless nights (darker bars), which is expected since clouds brighten the night sky. The B-V color index value increases up to B-V = 2.0 with clouds (more yellow, or less blue, sky). This results in a sky that is much redder than the natural one, which may have a relevant impact on the biological environment.
counting for 28% of AstMon total observations. As regards the photometer data, we used 283,922 measurements, accounting for 41% of the SQM measurements, and 234,543 measurements, accounting for 39% of the TESS-W measurements. The percentage of useful data for this research (no moon, no clouds) depends on the experiment starting date, device operability, and the hourly instrument sampling of night sky brightness (Table 1).

Madrid Night Sky Spectra
In Figure 4, we present a graph comparing the mean night sky spectra after midnight in January 2014 and January 2019. Examining this figure, we see that the 2019 night spectrum shows the blue bump of the LED spectrum. The Johnson B, V, and R photometric band responses show how the changes in spectra are related to the change in the sky's color.

AstMon
In Figure 7, we show histograms of the Johnson B, V, and R night sky brightness and color indices at the zenith measured using AstMon. A progressive brightening of the Madrid sky is detected in the Johnson B photometric band. The Johnson B, V, and R band brightness values for nights with both moon and clouds (lighter bars) were found to be to the left of values for moonless and cloudless nights (darker bars), which is expected since clouds brighten the night sky. The B-V color index value increases up to B-V = 2.0 with clouds (more yellow, or less blue, sky). This results in a sky that is much redder than the natural one, which may have a relevant impact on the biological environment. Comparing Figures 8 and 9 shows that, before midnight, the sky is brighter and bluer than after midnight when it is dimmer and more yellow. The data were, therefore, split into before and after midnight to remove the ornamental light component, headlights Comparing Figures 8 and 9 shows that, before midnight, the sky is brighter and bluer than after midnight when it is dimmer and more yellow. The data were, therefore, split into before and after midnight to remove the ornamental light component, headlights from cars, and light from buildings that contribute to brighter and bluer skies before midnight. This process helped isolate the contribution of the street lighting component. In these figures, we present histograms that only include filtered data. The range of values is reduced, and the changes can be seen more clearly. We used a non-parametric Gaussian Kernel Density Estimation (KDE) to derive a representative value for each histogram (mode) and the underlying distribution. The KDE uses an optimization via the Scott rule of thumb, a databased procedure for choosing an optimal bandwidth [27]. The bandwidth of the Gaussian kernel is chosen to minimize the magnitude mean integrated square error (goodness-of-fit). There are several other options for selecting a Kernel, but we used the Gaussian reference because of the comparison with the optimal bin widths and different densities (within the kernel). The densities can then be rescaled into a standard form, which facilitates the comparison between the optimal bin and the different densities yearly by Johnson B, V, and R photometric bands, and associated color indices.
(mode) and the underlying distribution. The KDE uses an optimization via the Scott rule of thumb, a data-based procedure for choosing an optimal bandwidth [27]. The bandwidth of the Gaussian kernel is chosen to minimize the magnitude mean integrated square error (goodness-of-fit). There are several other options for selecting a Kernel, but we used the Gaussian reference because of the comparison with the optimal bin widths and different densities (within the kernel). The densities can then be rescaled into a standard form, which facilitates the comparison between the optimal bin and the different densities yearly by Johnson B, V, and R photometric bands, and associated color indices.
It should be noted that the Johnson B, V, and R photometric bands used a logarithmic unit of mag/arcsec 2 , which is related inversely to the radiance. The results may differ after conversion to linear units.  The evolution of night sky brightness and color after midnight with reduced ornamental lights, headlights, and light from buildings in Figure 9 shows a progressive evolution toward blue after 2015, when a significant streetlight retrofit replaced HPS lamps with LEDs. The light in the blue band brightened in comparison with less variable changes in the V band and a compensated change in the R band. It is interesting to note that the B-V color index highlights this bluing of the night sky (lower values of the B-V index). The It should be noted that the Johnson B, V, and R photometric bands used a logarithmic unit of mag/arcsec 2 , which is related inversely to the radiance. The results may differ after conversion to linear units.
The evolution of night sky brightness and color after midnight with reduced ornamental lights, headlights, and light from buildings in Figure 9 shows a progressive evolution toward blue after 2015, when a significant streetlight retrofit replaced HPS lamps with LEDs. The light in the blue band brightened in comparison with less variable changes in the V band and a compensated change in the R band. It is interesting to note that the B-V color index highlights this bluing of the night sky (lower values of the B-V index). The change is also clear in the V-R color index.
After midnight, the Johnson B photometric band brightened following the street light retrofit, seen in the movement of the peaks to the left (Figure 9), while the Johnson R photometric band dimmed, as seen in the movement of the peaks to the right. The Johnson V photometric band was less sensitive to the streetlight retrofit. The B-V and V-R color indices showed changes from reddish to bluish, allowing us to identify HPS, streetlight retrofit, and LED data clusters ( Figure 10). For yearly night sky brightness and color statistics, see Tables S1-S5 of the Supplementary Materials. We visualized these changes using a color-color diagram, where the observations are color-coded and binned by observation date (Figure 10). The two panels in Figure 10 show how the change of technology transformed the night sky color in the expected direction: bluer skies as the LEDs are being introduced. These diagrams show the HPS, street light retrofit, and LED data clusters and the shift from reddish in the HPS region to bluish in the LED region. Splitting the data before and after midnight allowed us, again, to better account for changes in the color indices.

Sky Quality Meter SQM
The SQM photometer registered brighter values before the streetlight retrofit in 2015 ( Figure 11); the values were similar to AstMon's Johnson V band measurement in the same period (City Hall decreased bulb power of the HPS lamps (150W to 50W)). After the streetlight retrofit during 2015-2016, SQM registered dimmer values, while Johnson V registered brighter values of night sky brightness. These differences between the Johnson V and SQM measurements are later explained in terms of color indices correction applied to We visualized these changes using a color-color diagram, where the observations are color-coded and binned by observation date (Figure 10). The two panels in Figure 10 show how the change of technology transformed the night sky color in the expected direction: bluer skies as the LEDs are being introduced. These diagrams show the HPS, street light retrofit, and LED data clusters and the shift from reddish in the HPS region to bluish in the LED region. Splitting the data before and after midnight allowed us, again, to better account for changes in the color indices.

Sky Quality Meter SQM
The SQM photometer registered brighter values before the streetlight retrofit in 2015 ( Figure 11); the values were similar to AstMon's Johnson V band measurement in the same period (City Hall decreased bulb power of the HPS lamps (150 W to 50 W)). After the streetlight retrofit during 2015-2016, SQM registered dimmer values, while Johnson V registered brighter values of night sky brightness. These differences between the Johnson V and SQM measurements are later explained in terms of color indices correction applied to the SQM photometer (Section 3.5).
Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 26 Figure 11. Evolution of SQM magnitudes by grouping the data into years. Darker bars represent filtered data (no moon over the horizon and no clouds), and lighter bars are unfiltered data (with moon and clouds). As expected, the unfiltered data were distributed to the left of the filtered data, indicating brighter skies on nights with moon and/or clouds (similar to Figure 7). All measures excluded twilight.
The SQM data recorded a brighter night sky before midnight and dimmer values after midnight ( Figure 12). As with AstMon, the procedure of splitting the data in this manner guarantees a reduction in the ornamental lights and other contributors to light pollution during the first part of the night (See Table S6 of the Supplementary Materials).
The year 2015 registered the highest magnitude (lower radiance) of the series, while numerical trends of the median were variable. Figure 11. Evolution of SQM magnitudes by grouping the data into years. Darker bars represent filtered data (no moon over the horizon and no clouds), and lighter bars are unfiltered data (with moon and clouds). As expected, the unfiltered data were distributed to the left of the filtered data, indicating brighter skies on nights with moon and/or clouds (similar to Figure 7). All measures excluded twilight.
The SQM data recorded a brighter night sky before midnight and dimmer values after midnight ( Figure 12). As with AstMon, the procedure of splitting the data in this manner guarantees a reduction in the ornamental lights and other contributors to light pollution during the first part of the night (See Table S6

TESS-W
For Madrid night skies, the TESS-W registered brighter values compared to the SQM. Consistent with AstMon and SQM, the data with clouds showed a bimodal distribution (Figure 13), while the data with no moon and no clouds showed a left-skewed distribution.

TESS-W
For Madrid night skies, the TESS-W registered brighter values compared to the SQM. Consistent with AstMon and SQM, the data with clouds showed a bimodal distribution (Figure 13), while the data with no moon and no clouds showed a left-skewed distribution.

TESS-W
For Madrid night skies, the TESS-W registered brighter values compared to the SQM. Consistent with AstMon and SQM, the data with clouds showed a bimodal distribution (Figure 13), while the data with no moon and no clouds showed a left-skewed distribution.  TESS-W reported brighter night sky brightness values before midnight and dimmer values after midnight ( Figure 13). As with the SQM, nothing can be said about the evolution of color indices. Table S7 shows a decreasing trend for the median magnitude values (brighter skies) for data after midnight. This result is expected because of the removal of HPS streetlights (which project in the red spectrum), changes in streetlight bulb power, and changes in the diffuser technology. TESS-W is better suited for detecting night sky brightness variation due to streetlight retrofit because it has an extended capacity to measure in red spectra. The progressive brightening of the sky detected by TESS-W ( Figure 14) is not apparent in the statistics of SQM (Figure 13). TESS-W reported brighter night sky brightness values before midnight and dimmer values after midnight ( Figure 13). As with the SQM, nothing can be said about the evolution of color indices. Table S7 shows a decreasing trend for the median magnitude values (brighter skies) for data after midnight. This result is expected because of the removal of HPS streetlights (which project in the red spectrum), changes in streetlight bulb power, and changes in the diffuser technology. TESS-W is better suited for detecting night sky brightness variation due to streetlight retrofit because it has an extended capacity to measure in red spectra. The progressive brightening of the sky detected by TESS-W ( Figure 14) is not apparent in the statistics of SQM (Figure 13).

SQM and Johnson Photometry Comparison
Sánchez de Miguel et al. [11] examined the light response properties of a synthetic SQM using the Johnson V magnitude and the B-V color index. The synthetic formula requires a conversion from the AB magnitude system to radiance, allowing a comparison with the Johnson V band. They showed that irrespective of the lamps used for street illumination, the formula SQMs = V + 0.5 *(B-V), provides a color index correction, (B-V) < 2 (see Figure 8 of that paper). We plot in Figure 15 a graph with the Johnson bands and SQM yearly measurements for Madrid. We add the results for the synthetic SQM (SQMs). Table  2 shows the differences in magnitude for SQM-SQMs. Figure 15 shows that the SQM photometer did not follow the Johnson B and V trends. After performing the color index correction, the SQMs trends, these agreed with the Johnson ones. This discrepancy in trends likely indicates a tendency of SQM optical elements to malfunction as a consequence of the filter becoming less transparent over time. The TESS-W photometer did not need a correction, as it showed the expected brightening of the night sky, supporting the experimental outcome that the SQM presented opposite trends from the ones expected [28].

SQM and Johnson Photometry Comparison
Sánchez de Miguel et al. [11] examined the light response properties of a synthetic SQM using the Johnson V magnitude and the B-V color index. The synthetic formula requires a conversion from the AB magnitude system to radiance, allowing a comparison with the Johnson V band. They showed that irrespective of the lamps used for street illumination, the formula SQM s = V + 0.5 *(B-V), provides a color index correction, (B-V) < 2 (see Figure 8 of that paper). We plot in Figure 15 a graph with the Johnson bands and SQM yearly measurements for Madrid. We add the results for the synthetic SQM (SQM s ). Table 2 shows the differences in magnitude for SQM-SQM s . Figure 15 shows that the SQM photometer did not follow the Johnson B and V trends. After performing the color index correction, the SQM s trends, these agreed with the Johnson ones. This discrepancy in trends likely indicates a tendency of SQM optical elements to malfunction as a consequence of the filter becoming less transparent over time. The TESS-W photometer did not need a correction, as it showed the expected brightening of the night sky, supporting the experimental outcome that the SQM presented opposite trends from the ones expected [28].    Figure 6 shows a typical graph of the variation in sky brightness and color throughout the night for a photometric night. The position and shape of these curves has been changing in recent years. In Figure 16 we present graphs with the composition of the curves for filtered data collected using AstMon in Johnson B separated into seasons. The graphs have been color-coded with the number of observations that represent the average curve for each season. Since winter has the most prolonged nighttime duration in Madrid (more than twelve hours), we use this season for comparison purposes.  Figure 6 shows a typical graph of the variation in sky brightness and color throughout the night for a photometric night. The position and shape of these curves has been changing in recent years. In Figure 16 we present graphs with the composition of the curves for filtered data collected using AstMon in Johnson B separated into seasons. The graphs have been color-coded with the number of observations that represent the average curve for each season. Since winter has the most prolonged nighttime duration in Madrid (more than twelve hours), we use this season for comparison purposes.

Variation throughout the Night
There are four regions of interest related to anthropogenic factors ( Figure 17). After twilight, the sky progressively darkens as human activity decreases from 18:30 to 20:30, as it is the beginning of the night. At 21:30, some ornamental lights are turned off. The second part of the night (after midnight) shows an almost flat curve. As the ornamental lights are off and the traffic is low, it is the best time to measure the streetlight component. Finally, after 4:30, the streetlight power is increased for the morning commute. To show the evolution of night curves in winter from 2010 to 2019, we binned the data from AstMon in hourly intervals and the color index B-V ( Figure 18). We highlight (in insets) on each graph in Figure 18 the time of 03:30 UT, which was optimal for streetlight signal detection. The winter data from SQM and TESS-W are also divided into hourly bins and presented below (Figures 19 and 20). We are reporting the night sky brightness median values and interquartile ranges (IQR)/2sigma for non-normally distributed data. The IQR represents a robust measure for scaling the variability in magnitudes. There are some interesting changes related to the Madrid streetlight retrofit that we will comment on in the discussion. There are four regions of interest related to anthropogenic factors ( Figure 17). After twilight, the sky progressively darkens as human activity decreases from 18:30 to 20:30, as it is the beginning of the night. At 21:30, some ornamental lights are turned off. The second part of the night (after midnight) shows an almost flat curve. As the ornamental lights are off and the traffic is low, it is the best time to measure the streetlight component. Finally, after 4:30, the streetlight power is increased for the morning commute. To show the evolution of night curves in winter from 2010 to 2019, we binned the data from AstMon in hourly intervals and the color index B-V ( Figure 18). We highlight (in insets) on each graph in Figure 18 the time of 03:30 UT, which was optimal for streetlight signal detection. The winter data from SQM and TESS-W are also divided into hourly bins and presented below (Figures 19 and 20). We are reporting the night sky brightness median values and interquartile ranges (IQR)/2sigma for non-normally distributed data. The IQR represents a robust measure for scaling the variability in magnitudes. There are some interesting changes related to the Madrid streetlight retrofit that we will comment on in the discussion.   There are four regions of interest related to anthropogenic factors ( Figure 17). After twilight, the sky progressively darkens as human activity decreases from 18:30 to 20:30, as it is the beginning of the night. At 21:30, some ornamental lights are turned off. The second part of the night (after midnight) shows an almost flat curve. As the ornamental lights are off and the traffic is low, it is the best time to measure the streetlight component. Finally, after 4:30, the streetlight power is increased for the morning commute. To show the evolution of night curves in winter from 2010 to 2019, we binned the data from AstMon in hourly intervals and the color index B-V ( Figure 18). We highlight (in insets) on each graph in Figure 18 the time of 03:30 UT, which was optimal for streetlight signal detection. The winter data from SQM and TESS-W are also divided into hourly bins and presented below (Figures 19 and 20). We are reporting the night sky brightness median values and interquartile ranges (IQR)/2sigma for non-normally distributed data. The IQR represents a robust measure for scaling the variability in magnitudes. There are some interesting changes related to the Madrid streetlight retrofit that we will comment on in the discussion.   We report in Tables 3-5  . The values indicate that most sky brightness variability happens in the first part of the night (before midnight). On average, the night sky brightness variability before midnight is 45%, 27%, and 20% in the Johnson B, V, and R photometric bands, respectively. Before midnight, Johnson B represents a combination of ornamental and commercial lighting, along with the LED streetlights. The Johnson R captures the dimming of HPS streetlights. These findings support previously reported values measured by the ISS [29].  The color index for each winter curve before and after midnight was smoother throughout the night sky from 2010 to 2019 (Bottom of Figure 20). Before midnight, there was no marked change in the B-V color (~0.8). After the streetlight retrofit in 2015, the color index after midnight changed considerably from 1.2 to 1.0, indicating bluer skies. Thus, the color difference throughout the night is now significantly lower. Finally, higher variability in the IQR was due to anthropogenic activities before midnight, and the brightness change at the end of the night when lights return to the nominal power.

Night Sky Brightness Changes over Time
We have presented a long-term time series of observations resulting from continuous, not sporadic, monitoring of the Madrid night sky. The observations from the astronomical observatory of Universidad Complutense de Madrid (inside the city) began in 2010 and was ongoing in 2021. Both broadband panchromatic photometers (SQM and TESS-W) and the All-Sky Transmission Monitoring Camera (AstMon) provided measurements in the Johnson B, V, and R photometric bands throughout the night, every night.
The ongoing streetlight retrofits in Madrid primarily replace HPS lamps with 3000K LEDs. Changes include the bulb, the lamp post model, and a modified lighting regime (programmed change in the power of street illumination) throughout the night. We obtained data before, during, and after a major streetlight retrofit that took place during the year 2015. Our methodology for analysis includes a manual visual filtering process to avoid noise in the data caused by measures taken during cloud coverage. The statistical analysis of the dataset allows us to determine the impact of the Madrid streetlight retrofits on the brightness and color of the night sky. Below, we discuss our findings on the evolution of night sky brightness and color over time in Madrid together with an analysis of the instrumentation used to collect our data. The values indicate that most sky brightness variability happens in the first part of the night (before midnight). On average, the night sky brightness variability before midnight is 45%, 27%, and 20% in the Johnson B, V, and R photometric bands, respectively. Before midnight, Johnson B represents a combination of ornamental and commercial lighting, along with the LED streetlights. The Johnson R captures the dimming of HPS streetlights. These findings support previously reported values measured by the ISS [29]. The color index for each winter curve before and after midnight was smoother throughout the night sky from 2010 to 2019 (Bottom of Figure 20). Before midnight, there was no marked change in the B-V color (~0.8). After the streetlight retrofit in 2015, the color index after midnight changed considerably from 1.2 to 1.0, indicating bluer skies. Thus, the color difference throughout the night is now significantly lower. Finally, higher variability in the IQR was due to anthropogenic activities before midnight, and the brightness change at the end of the night when lights return to the nominal power.

Night Sky Brightness Changes over Time
We have presented a long-term time series of observations resulting from continuous, not sporadic, monitoring of the Madrid night sky. The observations from the astronomical observatory of Universidad Complutense de Madrid (inside the city) began in 2010 and was ongoing in 2021. Both broadband panchromatic photometers (SQM and TESS-W) and the All-Sky Transmission Monitoring Camera (AstMon) provided measurements in the Johnson B, V, and R photometric bands throughout the night, every night.
The ongoing streetlight retrofits in Madrid primarily replace HPS lamps with 3000 K LEDs. Changes include the bulb, the lamp post model, and a modified lighting regime (programmed change in the power of street illumination) throughout the night. We obtained data before, during, and after a major streetlight retrofit that took place during the year 2015. Our methodology for analysis includes a manual visual filtering process to avoid noise in the data caused by measures taken during cloud coverage. The statistical analysis of the dataset allows us to determine the impact of the Madrid streetlight retrofits on the brightness and color of the night sky. Below, we discuss our findings on the evolution of night sky brightness and color over time in Madrid together with an analysis of the instrumentation used to collect our data.

SQM and TESS-W
The SQM is a popular photometer in the light pollution community due to affordability and an easy-to-use user interface. The SQM bandpass was designed to mimic the human eye response and provide values similar to the Johnson V photometric band [30]. The TESS-W photometer bandpass is defined by a dichroic filter [12] that extends the spectral response to the red range, which allows measurement of the wide HPS line range as shown in Figure 4 (upper panel). Sánchez de Miguel et al. [11] demonstrated how SQM leads to misleading results by comparing sky brightness measured with SQM during a retrofit from HPS to LED. They proposed addressing this problem using a color index adjustment for night sky brightness data measured with SQM and encouraged the use of color observations to determine the real evolution of sky brightness.
In this work, we had the opportunity to compare the data collected by TESS-W and SQM, as well as an astronomical camera with color information (AstMon). The average winter evolution curve throughout the night for TESS-W data shows a steady progressive brightening from 2016 to 2019 (Figure 20). On the other hand, the SQM data ( Figure 19) show a brighter sky in 2016 (following 2015, the retrofit's critical year), and a progressive darkening of the night sky brightness from 2016 to 2019. This is a prototypical example of how a streetlight retrofit from HPS to LED increases the nocturnal sky brightness, but a simple analysis of the SQM data yields the opposite conclusion. Figure 21 shows a comparison of the yearly median brightness measured by SQM and TESS-W. In addition to SQM data requiring a color index correction (Figure 15), Bará et al. [28] and Puschnig et al. [31] recently found that SQM photometers do not properly detect light over several years (an aging effect). They found a correction drift of around +0.05(6) magSQM/arcsec 2 per year. This correction likely explains the discrepancy in SQM brightness trends with AstMon and TESS data. The SQM was four years older than the TESS-W. Therefore, a time-dependent sensitivity change might also play a role in these trends. shown in Figure 4 (upper panel). Sánchez de Miguel et al. [11] demonstrated how SQM leads to misleading results by comparing sky brightness measured with SQM during a retrofit from HPS to LED. They proposed addressing this problem using a color index adjustment for night sky brightness data measured with SQM and encouraged the use of color observations to determine the real evolution of sky brightness.
In this work, we had the opportunity to compare the data collected by TESS-W and SQM, as well as an astronomical camera with color information (AstMon). The average winter evolution curve throughout the night for TESS-W data shows a steady progressive brightening from 2016 to 2019 (Figure 20). On the other hand, the SQM data ( Figure 19) show a brighter sky in 2016 (following 2015, the retrofit's critical year), and a progressive darkening of the night sky brightness from 2016 to 2019. This is a prototypical example of how a streetlight retrofit from HPS to LED increases the nocturnal sky brightness, but a simple analysis of the SQM data yields the opposite conclusion. Figure 21 shows a comparison of the yearly median brightness measured by SQM and TESS-W. In addition to SQM data requiring a color index correction (Figure 15), Bará et al. [28] and Puschnig et al. [31] recently found that SQM photometers do not properly detect light over several years (an aging effect). They found a correction drift of around +0.05(6) magSQM/arcsec 2 per year. This correction likely explains the discrepancy in SQM brightness trends with AstMon and TESS data. The SQM was four years older than the TESS-W. Therefore, a time-dependent sensitivity change might also play a role in these trends.

ASTMON Brightness and Color Evolution
To show the evolution in brightness and color of the Madrid sky during the night, in Figure 22 we present changes in the median values of brightness and color indices for each year in the first and second parts of the night using AstMon data. The dispersion of values is higher during the first part of the night when there is variation in human activity, as expected. The second part of the night is consistently darker. The evolution of median night sky brightness shows a compensated difference of magnitude in the Johnson B and

ASTMON Brightness and Color Evolution
To show the evolution in brightness and color of the Madrid sky during the night, in Figure 22 we present changes in the median values of brightness and color indices for each year in the first and second parts of the night using AstMon data. The dispersion of values is higher during the first part of the night when there is variation in human activity, as expected. The second part of the night is consistently darker. The evolution of median night sky brightness shows a compensated difference of magnitude in the Johnson B and Johnson R photometric bands. While the Johnson R band increases in magnitude (towards lower brightness), the Johnson B presents the opposite trend, showing increasing brightness. The B-V and V-R color indices decrease in value over the years, showing a bluer night sky. The Johnson V photometric band was less sensitive to the streetlight retrofit and did not show a marked trend. night sky. The Johnson V photometric band was less sensitive to the streetlight retrofit and did not show a marked trend.
To check the significance of the trends in Figure 22, we performed a model selection using leave-one-out cross-validation [32]. The two competing models are: a simple linear trend and a constant value, both including the variance of the data points. We fit the models before and after midnight, and for each magnitude and color. Before midnight, the best model for the Johnson R band and V-R is a linear trend (increase in magnitude and decrease in V-R color), whereas for the rest of magnitudes and colors the best model is flat. The best model in all cases after midnight is a linear trend, except for Johnson V which was less sensitive and did not show a marked trend.
.  To check the significance of the trends in Figure 22, we performed a model selection using leave-one-out cross-validation [32]. The two competing models are: a simple linear trend and a constant value, both including the variance of the data points. We fit the models before and after midnight, and for each magnitude and color. Before midnight, the best model for the Johnson R band and V-R is a linear trend (increase in magnitude and decrease in V-R color), whereas for the rest of magnitudes and colors the best model is flat. The best model in all cases after midnight is a linear trend, except for Johnson V which was less sensitive and did not show a marked trend.
For the first part of the night, while the increase in Johnson R is significant, the dispersion in the median values of Johnson B is enough to explain the changes in this magnitude. In the second part of the night, the evolution of median night sky brightness shows a compensated difference of magnitude in the Johnson B and Johnson R photometric bands. While the Johnson B band presents the opposite trend, showing increasing brightness. The B-V and V-R color indices decrease in value over the years in the second part of the night, showing bluer night sky.
We observed that the night sky brightness increased primarily in the blue part of the spectrum (Table 6), and most of the detected change took place during 2015 when the major streetlight retrofit from HPS to LEDs occurred. The winter curve is uniquely suited to demonstrate the effect of anthropogenic factors on sky brightness data throughout the night because the night is longer than in other seasons, providing more sampling points. Examining the winter curve, we could identify the part of the night when anthropogenic factors were reduced to an extent that the remaining light was primarily from streetlights. These winter curves (Figures 18-20) allow us to understand the effect of the 2015 streetlight retrofits on the Madrid night sky brightness and color. We emphasize the importance of identifying the streetlight signal with reduced noise. The large sample, available due to continuous nightly measurement, allowed us to identify this signal and is a strength of our design.

Comparison with Other Studies
Using data gathered by satellites, Kyba et al. [33] detected a 2% increase in global night sky brightness in lit outdoor areas. However, they found that measurements by local devices reported greater increases in some regions, suggesting that satellite data sometimes underestimate radiance on a local scale. Our study supports these findings, where TESS-W, a device that measures in the yellowish-red spectra, reported a 10% increase in brightness following a major streetlight retrofit in 2015.
We found that the transition towards solid-state LED technology impacts night spectra, supporting findings by Posch et al. [34]. According to Bará et al. [6], night sky brightness responses behave differently for narrower bands. Therefore, an assessment of increases in sky radiance due to transitions from HPS to LED has to include measures of the brightness in the channels of major changes. We measured in the Johnson B, V, and R photometric bands and found that the bands did indeed behave differently.
Madrid City Hall did not make significant changes to the number or manufacturing type of streetlights in Madrid after 2015. Hence, the night sky brightening reported in the B and V Johnson photometric band and the TESS-W photometer following the 2015 streetlight retrofit have to be produced by other sources. Potential sources for brightening in the night sky at the UCM observatory after 2015 include: (a) inaccurate Madrid City Hall data (we consider this unlikely, but not impossible), (b) changes to street light diffusers in Madrid that increased the upper light ratio, and (c) street lighting changes in Madrid's surrounding areas have impacted the night sky brightness measured at the UCM observatory.

Conclusions
Examining filtered data (no moon, no clouds) from the winter season in Madrid at 03:30 (when other anthropogenic light pollutants are minimized, allowing optimal streetlight signal detection), we found: As regards SQM:

1.
A 38% darkening of the night sky from 2012 to 2019, as measured by SQM. After the streetlight retrofit in 2015, the SQM photometer did not follow the trends detected by AstMon in the Johnson V band and the broadband TESS-W photometer. Although the change in shifting lighting types could play a small role, the main effect is the loss of the SQM photometer sensitivity. This experimental outcome has far-reaching implications for long-term global monitoring of night sky brightness, if confirmed by similar studies.
As regards TESS-W: We conclude that due to the inclusion of the dichroic filter, TESS-W is better suited to tracking changes in night sky brightness compared to SQM, allowing for the detection of night sky brightness in the yellowish-reddish range of spectra (HPS range).
Regarding the AstMon data: The streetlight retrofit from HPS to LED in Madrid:

7.
Madrid's HPS streetlamps remain a strong contributor to the night sky brightness V and R components. The new LED streetlights changed the proportions of the night sky brightness B and V components. 8.
Without careful design, a white LED streetlight retrofit is not optimal to reduce light pollution, even after reducing the electric power because of the high content of blue light. The relevant energy savings are not due to the lamps' energy efficiency but the dimming capabilities and the lower ULOR of the LEDs compared to traditional lighting. An alternative way to reduce light pollution consists of reducing the power of in-place sodium lamps. 9.
There was a small decrease in sky brightness following the streetlight retrofit in 2015, but these changes were later reversed.

Data Availability Statement:
The data presented in this research are openly available in https: //zenodo.org/record/4633001 (accessed on 9 April 2021).

Acknowledgments:
The original manuscript was improved with the comments and suggestions received from four referees We thank The Cities at Night project for its supports.

Conflicts of Interest:
The Authors declare no conflict of interest.

Abbreviations
The