Evaluating the potential health impact of obtrusive light: Application to Montréal, Canada

: This paper describes the use of a new obtrusive light module of the Illumina v2 model to estimate the light that may enter bedroom windows. We used as input to the model, 1- the sources’ ﬂux and spectrum derived from the color images taken by astronauts from the international space station, 2- an association between source spectrum and angular emission, and 3- a per zone inventory of obstacles properties and lamp height. The model calculate the spectral irradiance incident to buildings’ windows taking into account for the orientation of the street. By using the color information from an ISS image, we can classify pixels as a function of their spectra. With the same image, it is also possible to determine the upward photopic radiance for each pixel. Both serve as inputs to the model to calculate the spectral irradiance on any window. By having the spectral irradiance, it is possible to determine the Melatonin Suppression Index and the photopic irradiance on the window. Such information can later be used to perform epidemiological studies. The new methodology is applied to the case of Montréal in Canada for a set of houses’ locations. The computations are made for 2013 (pre-LED era).


Introduction
Innovations in the field of lighting have shaped the history of the last century. Due to their numerous uses, our society has become light dependent. As any other innovation, lighting technologies come with their advantages and disadvantages. The negative aspects of these technologies seem to be unknown by most people. The unexpected effects of what is now called Artificial Light At Night (ALAN) have first been mentioned by astronomers. Within the past decade, multidisciplinary researchers became aware of the actual issues. Indeed, abnormal behaviours have been observed among several nocturnal species [1] due to ALAN. Negative impacts on bats [2], turtles [3], fish [4] and songbirds [5] have already been identified. Moreover, it is now known that ALAN does have an effect on plants [6][7][8][9][10]. On the other hand, what seems even more concerning is the possible impacts of ALAN on human health. For humans, the first suspected impacts are related to breast and prostate cancers [11][12][13]. ALAN plays a role in disturbing circadian rhythm [14], even in low intensity lighting [15]. The disturbance in the circadian cycle also plays a role in the increase of obesity [16]. In order to provide a reference to establish the potential impact of a given spectrum of ALAN on the melatonin suppression, Aubé et al. 2013 introduced the Melatonin Suppression Index (MSI). Since its introduction, some studies verified if associations exist between the MSI and a variety of diseases. Among them, Garcia-Saenz et al. [18] and Garcia-Saenz et al. [19] found that, for Spanish cities, regions with high MSI, are associated with increased risk of breast, prostate and colorectal cancers.
The goal of this paper is to present a new methodology combining color remote sensing of the lighting technologies with a radiative transfer model to determine the spectral irradiance falling on a for reflection on buildings. However, when the ground is covered with fresh snow, the reflection on the ground can rise up to the same order of magnitude as the obtrusive direct light.

Modeling domain and period
Throughout this modeling experiment, we used a resolution of 8.5 m which is more or less the resolution of the ISS image used to determine the lamp inventory. For that experiment, we consider the position of each house's window where we want to find the obtrusive light as a new virtual observer. Illumina v2 considers each virtual observer as a different model run and always locates the virtual observer in the center of the modeling domain. In this experiment, we used two layers having a size of 25 x 25 pixels. The first layer at a resolution of 8.5 m (201 x 210 m) and the second at a resolution of 21.25 m (531 x 531 m). This is large enough because we do not want to calculate the contribution of the light scattered by the sky that is considered negligible. One advantage of such a configuration is the rapidity of execution of Illumina v2.
Montréal is the second-biggest city in Canada, and an urban agglomeration of approximately 4 million citizens in 2014 [31]. The city has a population density of 4 438.7 citizens per square kilometre [32]. Overall, the agglomeration has a surface of 499.1 square kilometers of which 365.2 square kilometers is representing the city itself [33]. As for the topography, it is relatively flat, except for a small hill named Mont Royal. The city contains more than forty skyscrapers higher than one hundred meters [34] and the overall building density is high. It is as well frequent to see trees positioned near street lamps, blocking a part of the emitted light during summer.
Before the massive conversion to LED starting in 2017, Montréal has a street lighting mainly composed of High-Pressure Sodium (HPS) and Metal Halides (MH) lamps. The street lights photometry of HPS was mainly assimilated to a ≈50-50 mix of cobrahead-like and helios-like light fixtures with respectively ≈5% and ≈1% of the Upward Light Output Ratio (ULOR). Such a mix results in a 2% ULOR on average. Most MH lamps were cobrahead-like ≈5% ULOR. Some private lights are assumed to be mainly composed of Compact Fluorescent lamps (CFL) with a 15% ULOR. This information is summarized in Table 1. The modeling experiment presented hereby only considers the light sources extracted from an ISS image (see Figure 1). This image covers a surface of about 850 km 2 . All sources falling outside this surface are not considered in our modeling experiment. This is not a problem since we are concerned about the obtrusive light coming from inside a radius of a few hundreds of meters. If one is concerned about the obtrusive light coming from the sky, Simoneau et al. [35] showed that about 90% comes from within 10 km but that situation is not considered in this study.

Modeling parameters
In order to achieve the modeling, a few parameters need to be determined. Those parameters are used as the base information for the model to be executed. We chose year 2013 as a reference period for the pre-conversion to LED lighting. The datasets can be divided into two main types, namely the uniform parameters and the gridded parameters. Those last ones are the parameters that may vary throughout the modeling domain.

Uniform parameters
The first parameter used in the model is the relative humidity. The relevance of this parameter is related to the fact that many aerosols are hydrophilic and consequently their properties, such as their size and their refractive index, are influenced by it. To determine the value that we used for the modeling, we took the average humidity that has been recorded during early 2010s for Montréal City, which is approximately 79.1%, and rounded it to 80 % [36].
Using the same methodology as above, the average atmospheric pressure for Montréal City has been set for the model as 101. 1 kPa [36]. This parameter has an influence over the density of molecules present in the modeling domain, thus affecting their interactions with photons.
Other parameters considered to be uniform over the whole domain are the aerosol optical depth and the Ångström coefficient which are linked to the aerosol density and the size distribution. For this modeling experiment, we used the data available on the aerosol robotic network (AERONET) website [37]. More precisely, the data used, both for the Ångström coefficient and the aerosol optical depth (AOD), is the level 2 daily average values recorded on April 8, 2013 at the closest AERONET site (CARTEL). AOD average values are indexed according to the wavelength. We selected the AOD relative to a wavelength of 500 nm, while the Ångström coefficient is the one recorded for wavelengths 440 nm and 870 nm. A primordial parameter is the ground reflectance. We assumed a mix of 90% asphalt and 10% vegetation to be representative. Table 2 shows the values of each uniform all over the domain parameter used for the modeling.

Gridded parameters
The first non-uniform over the whole domain parameter that is required is the digital elevation model (DEM). This parameter, which is basically the topography of the modeling domain, is obtained from the Shuttle Radar Topography Mission (SRTM) [38]. The absolute accuracy concerning the vertical altitude is of 16 meters, while the accuracy of the relative height between two ground elements is of 10 meters. The horizontal resolution is about 30 meters.

Lamp flux and spectral type from ISS images
The upward ISS photopic radiance and the lamp spectrum classification come from the ISS night images [39]. We used the ISS035-E-17088 image taken during the ISS035 mission with the Nikon D3S Camera (400mm lens) on April 6, 2013 ( Figure 1). We used pre-processed images provided by NOKTOsat [40], from which we developed a filtering process to remove the natural background radiance and the noisy signal over unlit surfaces. This allowed to limit the number of light-emitting pixels to the significant ones. Figure 4 summarizes the steps undertaken in order to obtain filtered ISS images. As a starting point, we had access to images conveying information on the photopic radiance (image ImpactVlG_GR from NOKTOsat) as well as the lighting technologies used throughout Montréal (image Composite from NOKTOsat). As for the second image, the technologies were separated into 5 different classes: number 1 relating to cool white sources (MH or LED higher than 4000 K), number 2 to warm white sources (CFL or MH between 3000 K and 4000 K), number 3 to a mix of 3000 K fluorescent, HPS, MH or white sources between 2400 K and 3000 K, number 4 to phosphor converted (PC) Amber LED and warm HPS, and number 5 to normal HPS. A sixth category representing Low Pressure Sodium (LPS) and pure amber LED is normally included in the image provided by NOKTOsat, but this category has not been detected in Montréal. There are very few pixels belonging to class number 3. Table 1 shows the choice of technology that we have made to meet the above classes in accordance to what was available in the city in 2013.
When looking at the photopic radiance image, it became apparent that there was some background lighting present in the image. This can be explained by the reflection of starlight, moonlight and artificial sky glow on the surfaces along with the upward scattering of ground based lights. Indeed, the images revealed the presence of light in areas where there should not be anthropogenic light emitted at all (such as empty fields or bodies of water). In order to consider only artificial light emissions, we estimated the value of this background light and subtracted it from the value of every pixel in the image. This background was determined by a statistical analysis of the image. At first, a temporary image containing only the pixels of lesser value was created. As points of higher emission intensity had been eliminated, only light emissions in dark areas remained. We then analyzed this remaining data and estimated that the addition of the mode and the standard deviation of this new image would be a good value for the background removal. Indeed, the background light is quite similar throughout the territory. It is therefore a logical assumption that the value of the background light is the most common value in the image. The standard deviation was also added to the background in order to completely eliminate the background light that could be somewhat higher than the mode in certain areas. The use of the standard deviation in the determination of the background was decided by trial and error. This empirical value fits well for the territory of study but might have to be changed when studying a different one with the same method.
Within this filtered image, isolated non-null pixels could still be found within areas that should not emit light. As to not consider these false light sources, this noise had to be eliminated from the data. A binary image was created in which the pixels were given a new value of 1 for non-null pixels and 0 for null pixels. A convolution was then applied to this binary image : a sliding window of 3 pixels by 3 pixels would consider the surrounding pixels of each individual pixel in order to determine if it needed to be eliminated. If a non-null pixel had 4 or more pixels of value 1 surrounding it (meaning that this pixel was not isolated within a dark area), it would remain of the value 1. If it had 3 or less pixels of value 1 in the sliding window, it would be given a new value of 0. This process resulted in a new binary image in which only the pixels that should realistically have a significant value were kept (the others having been converted to zeros). Finally, this binary image was multiplied to both the photopic radiance and the lighting technology images in order to obtain clean images (see  Afterwards, we created an image of the MSI throughout Montréal from the information contained in the filtered lighting technology image ( Figure 2). We estimated the MSI for each of the 5 categories of technology present in Montréal in 2013 (Table 1). We used the data on Lamp Spectral Power Distribution Database (LSPDD) [41] to find the relevant lamp and determined the associated MSI for each spectral class. We created a lookup table associating the MSI to its class. Finally, we substituted the value relating to the class for the value of the MSI in the image ( Table 1). The resulting image therefore presents the MSI for each pixel.

Obstacle properties and lamps heights
The average height of street lamps and properties of the obstacles for each region are determined using Google Maps StreetView [42] and Google Earth [43]. A total of 309 circular zones with various radius and specific central coordinates have been defined with Map Developers [44] in order to cover the entire surface of Montréal. These zones were established according to the relative uniformity of the buildings' size, the distance between them, the lamps' height over the area and the obstacle filling factor. The procedure established to collect the required information is the same for each zone.
The first parameter defined for each zone is the average building height (h o ). This obstacle height was measured using the Path3D function in Google Earth. The average street lamp height (h L ) was also estimated for each zone using the same method. Then, the distance between the obstacles (d o ) (buildings, trees) was determined by using the ruler tool in Google Earth at an angle of ∼45 o from the street orientation. Although the closest obstacles would usually be in front or behind the street lamps, we considered that street lamps are engineered to direct light slightly in front of them but mostly in the direction of the street. Thus, we estimated that the average distance a photon could travel before colliding with an obstacle would be approximately equivalent to the distance of the first obstacle met at an angle of 45 degrees. We finally determined the obstacle filling factor (F). This obstacle property is a ratio out of one corresponding to the peripheral angle of a position where the path of the photons is blocked. For instance, an obstacle filling factor of 0.5 is attributed to a position if ∼180 o out of the total 360 o are blocking light trespass. Table A1 summarizes the various zones defined to characterize the obstacles and lamp heights. In Illumina, the order in which the zones are provided to the model is related to their priority (a new zone having a higher priority over the previous). In case of an intersection between zones, the latest zone overwrites the previous one. In Table A1,   values, is valid for regions where no other circular zone is defined (since any zone superimposed on the first predominates).

Generic light output patterns and spectra
We considered 3 different ULORs which are 2%, 5%, and 15%. The cobrahead-like LOP (5% ULOR) was taken from the Illuminating Engineering Society of North America (IESNA) file for the Cooper's cobrahead model. This LOP was systematically associated to MH lamps. The HPS LOPs are typically in relatively equal number of cobrahead (5% ULOR) and of Helios models (1% ULOR). The average of their IESNA files gives a new LOP with 2% ULOR. For the CFL, we assumed a LOP having 15% ULOR. In the latter case, the LOP was measured by ourselves with a farm lantern (model Globe PL-8120).
We used five different lamp spectra which were taken from the LSPDD. Figure 6 shows the spectra of the lamps used for this experiment. The associations between the spectra and the LOPs are given in Table 1.
For each pixel, we determine the spectrum and its radiant flux Φ e . This is done thanks to the ISS derived spectral classes and upward ISS photopic radiance images. Knowing the typical LOP and spectrum associated to each pixel according to its spectral class (see Table 1), one can convert the ISS derived photopic radiance into spectral flux of the pixel (in Watt). This is done with equation 1.  Table A1.
L e,Ω (θ) is the photopic radiance given by the ISS image Impact_VlG_GR (in units of nw sr -1 cm -2 Å -1 ) multiplied by the nominal photopic bandwidth (∼1000 Å). S is the area of the pixel in cm 2 (which is 722 500 cm 2 for a pixel of 8.5 m). V(λ) is the spectral sensitivity of the photopic band. θ is the zenith angle between the image center and the ISS (31 o in our case). ρ(λ) is the ground reflectance. F 90−180 (λ) (see equation 2), andḠ(θ, λ) are respectively the amount of light going down per unit of wavelength, and the LOP value at any wavelength for the zenith angle z = θ. T(λ) is the atmospheric transmittance but that factor has been already corrected by NOKTOsat so that we set its value to 1.
According to [26], the spectral intensity I e,Ω (z, λ), in [W/sr/nm], leaving the light source pixel at any zenith angle z is given by equation 3.
The radiant flux of a pixel Φ e , in [W], is determined from the ISS photopic radiance image and equation 1 combined with: 1-information from the pixel net LOP and spectrum, 2-the photopic spectral response and 3-the ground spectral reflectance [26].Ḡ(z, λ), in [Wsr −1 nm −1 ], is the light radiation pattern of a pixel. This function gives the relative amount of light emitted at any zenith angle and any wavelength per unit of solid angle. In Illumina's usual configuration, where we get the radiance from VIIRS-DNB,Ḡ(z, λ) of a given pixel is a combination of different street lamps having different LOPs and different spectra because of the large size of the pixel footprint. In the present adaptation to the higher-resolution data from ISS images, the LOP of a pixel is the same for any wavelength since we assume that only one light fixture falls inside a pixel (i.e.Ḡ(z, λ) =Ḡ(z)).

Determining the window orientation
Determining the angle and orientation between the window and the street lamps is crucial to model the irradiance entering in a home. We assume a window facing the closest street so that the front of the house is parallel to the street. With the objective of making our modeling approach as general as possible, we developed a numerical method to determine from any address in the world, the closest driving street. This is done by using OpenStreetMap data and an r-tree integrated in the OSMnx python library to calculate the minimal euclidean distances between points [45]. The library is also used as a tool to fragment into straight lines each small street segments. This allows obtaining the bearing angle of a segment even in curved streets. We then add the appropriate angle to get a right angle between the street and the window. This angle is then used as an input value to the Illumina model (the azimuth viewing angle). For the elevation viewing angles, we aim to determine the irradiance on the vertical surface of a window. So that the elevation angle is always equal to zero. Figure 7 shows the resulting angles for a small part of Montréal.

Results for some Montréal sites in 2013
In order to test and compare our modeling and remote sensing methodologies, we identified five different sites across the modeling domain. They are differing by their urban zoning categories. We actually selected a site in the city center, another in an urban residential area, one in a suburb residential area, plus one in a commercial area and finally another in a rural area. Each site differ by their distances between buildings and buildings' heights. All precise locations are on a building facade facing the street. We defined up to three floors but for some sites only one or two were kept with respect to the actual height of the buildings. We used 2, 5, 8 m above ground as the default set of heights. Detailed information about the sites is shown in Table 3.  Table 4 shows the comparative results obtained with the two methods to estimate both the MSI and the amount of photopic light for the five sites. One can notice that for MSI, no significant variation occur according to the window height (Illumina MSI). On the other hand, the Illumina photopic irradiance decreases with the window height. There is no obvious link between the Illumina MSI and the upward ISS MSI. The same observation arises with the comparison of upward ISS photopic radiance and the Illumina photopic irradiance.
In order to better explore if the ISS derived parameters are somehow correlated to the equivalent parameters determined with Illumina v2, we decided to select 500 sampling points across the city. The 500 points locations are shown in Figure 8. If no significant new information is provided by Illumina v2 compared to the one derived from the ISS images, we should observe a relatively good correlation between the values obtained directly from the ISS images and the Illumina v2 derived values. Figure 9 shows v2 provides significant new information compared to the ISS derived parameters. For that reason such modeling results can represent an advantage over direct remote sensing techniques when it comes to verify if there is any associations between the ALAN and health issues. Figure 9 also shows that the correlation coefficient between the ISS MSI and Illumina v2 do not significantly change according to the window height. This result is consistent with the fact that for our five sites, no significant changes were observed with various window height. One can notice that, for the MSI correlations, the density of points is higher around the 1:1 relationship and that Illumina generally predict higher MSI for low ISS MSI values and lower MSI for high ISS MSI values (see panels (a), (b) and (c) of Figure 9). This is normal because the Illumina MSI does not only involve the pixel's lamp but also the nearby ones so that Illumina provides a sort of complex averaging effect on the MSI. However, the relationship between the upward ISS photopic radiance and the Illumina photopic irradiance shows a decreasing correlation with higher window height. Figure 10 shows the irradiance spectra obtained with Illumina v2 for the 5 sites for windows 2 m above ground. We only provided them for 2 m because almost no changes were observed in the spectra shapes with higher window height. It is also clear that the spectrum is highly sensitive to the site location and the respective mix of lamp spectra located around each location. One nice aspect of using Illumina v2 is that it provides spectra. Having spectra allows a verification of any other lighting parameter, not only the MSI. As an example, one can calculate the CCT. Another interesting possibility provided by the model is to determine the photoreceptor response to the given irradiance by humans or any other animal specie. Illumina v2 also calculates the radiance spectra in a specific direction so that it can be useful when the direction of the light is of some importance. We did not use that feature in the present study.

Conclusion
The goal of this paper was to present a new methodology combining color remote sensing of the lighting technologies with the radiative transfer model Illumina v2 to determine the spectral irradiance falling on a window surface facing a street at any location. By using the numerical model Illumina v2, we were able to establish that no clear correlation exists between the model results and the MSI and photopic radiances derived with the ISS images. The low correlation may indicate that Illumina v2 provides new information that is not available from the ISS images. We think that this result highlights that the use of Illumina v2 is promising to perform epidemiological studies related to the color and amount of obtrusive light. For that reason, we aim to provide such information to researchers in the field of epidemiology in order to verify if any association can be found between the Illumina v2 derived obtrusive light parameters and human health issues. One other interesting result is that the Melatonin Suppression Index is almost insensitive to the window height while the spectral irradiance decrease with the window height. Moreover, the correlation between the ISS derived radiances and the Illumina v2 irradiances is worse for higher window heights. The reduction of the irradiance with height combined with the constant MSI is interesting to disentangle the effect of the spectral content versus the amount of light.
In the future, we plan to perform an in situ validation experiment in Montréal in order to determine the accuracy of the modelled approach presented here. That experiment will require the measurement of the irradiance on a set of sites for which we computed the Illumina v2 irradiance spectra. We are expecting to perform that experiment with the LANcube radiometer [46] in the coming year.
Author Contributions: We applied the sequence-determines-credit approach [47]

Conflicts of Interest:
The authors declare no conflicts of interest.     Figure 10. Spectral irradiance on windows 2m above ground for the five test sites identified in Table 3.