Night-sky imaging as a potential tool for characterization of total lumen output from small and medium-sized cities

In this article, the asymptotic formula developed in past work and applied to predict skyglow due to distant sources was evolved, with the objective of characterizing small and medium-sized cities in the observer’s surroundings. To enable this, a combination of theoretical computations and in situ measurements is needed, aiming to distinguish between dominant and smaller light-emitting sources, with the latter usually being camouﬂaged when measuring the night sky. Furthermore, for numerical modelling of skyglow, few of the most important parameters, speciﬁcally the amount of total lumens installed and radiated to the upward hemisphere, can be derived. Astronomical observatories, in particular, can proﬁt from this concept, since they are usually situated far away from large cities but can still be surrounded by smaller villages and towns. We present a detailed description of how theoretical computations are combined with all-sky photometry in order to obtain the properties mentioned. Results are compared with satellite data, showing that, regarding approximations undertaken for processing, they are comparable, underlining the functionality of our approach. The idea of including in situ observations enables us to quantify the impact of small and medium-sized cities globally and independent of location, as long as measurements were conducted outside light domes. In addition, the presented work may be of major interest to the light-pollution community if conducting long-term observations of cities, since the quality of commonly used satellite data is going to decrease in the future, due to blindness in short wavelengths and upcoming conversions of public lighting systems to blue-enlightened LEDs.

form, to name only a few examples out of the many existing. Moreover, daily observational astronomical work suffers, since light sources on the ground are impairing the visibility of our night skies and thus faint objects become barely visible. The prediction of skyglow variation in different territories is therefore of great importance from a nature conservation and human health perspective.
Among light-source spectral power distributions and their optical distortion due to light scattering and absorption phenomena in the atmospheric environment, the total light output of a city is a key element in all skyglow modelling tools. The amount of light radiated towards the upper hemisphere is critical to evaluate the resulting impacts from skyglow on living organisms, due to providing a variety of sensitivities and threshold levels in their perception of light intensities (see e.g. Perkin et al. 2011). Although satellite data can provide information on upward directed light, the large amount escaping the urban area in other directions and at other angles is largely unknown. Therefore, the light output contributing to skyglow is difficult, or rather impossible, to infer solely from satellite observations. Some estimates are possible using information on the total power consumption. As an example, Joseph et al. (1991) inferred the light output from street lighting of Tel Aviv based on the average power consumption of the luminaires used and evaluated the upward light intensity in photons per second. Light-output appraisal was also performed in a similar way for the city of Vienna, taking into account the total power consumption of the city's street lighting system (Puschnig, Posch & Uttenthaler 2014). However, this approach can only provide approximate values for total lumen outputs, while contributions coming from other light sources, e.g. vehicles, residential lighting, wall-mounted security lighting or the floodlights used in sports stadiums, are excluded from such derivations (even if they are likely equally important). In most cases, the number of total lumens installed is only known for a a small number of cities, due to having a publicly unavailable or incomplete inventory of light sources. Consequently, there is a strong need for the development of an efficient technique to quantify the lumen output, which is applicable to many cities and independent of whether official statistics for public illumination are available or not. Duriscoe, Luginbuhl and Elvidge (2014) used satellite data obtained by the Suomi-NPP Visible Infrared Imaging Radiometer (VIIRS) as inputs to compute the frequently used Garstang (1986) model, in order to predict the total lumens of outdoor lighting from an anthropogenic sky luminance observed at locations on the ground. A new approximate formula simulating the night-sky brightness at intermediate and larger distances from light sources has been derived only recently (Kocifaj, Wallner and Solano-Lamphar 2019) and is used here to characterize light emissions in terms of total lumen outputs. The formula appeared useful, not only to estimate the amount of light escaping from dominant cities but also to save significant computational resources when modelling skyglow over a large territory. A retrieval method based on the aforementioned formula is preferred, primarily because the detection technique is performed rapidly, the experimental equipment is low-priced and there is no requirement for satellite data. However, use of the formula for mid-sized cities remains largely unexamined. The method introduced here requires measurements of the sky brightness (whatever unit is selected, i.e. radiometric for radiance or photometric for luminance) along a vertical plane crossing the hemisphere in the azimuthal position of a city. Such data are usually difficult to interpreted if acquired toward a medium-sized (satellite) city, due to the optical signal potentially being compounded by the superposition of light emission from both the small or medium-sized city and a dominant city in its vicinity. The key issue behind the new technique is whether the contribution from these smaller cities/towns can be displayed individually, with the aim of assessing their total output emission in lumens or Watts. Astronomical observatories, in particular, are usually located far away from large cities to keep their influence on sky brightness low. Still, smaller light-emitting sources, like villages or towns, can be situated around them and possibly evolve to increase their impact on skyglow. Consequently, it is of great importance for observatories to characterize small and medium-sized cities and distinguish their long-term development independent of other dominating or larger light sources. This article addresses the above question and reports that success may be achieved with the proper combination of in situ optical measurements, theoretical analysis and numerical modelling. The theoretical model we have developed in 2019 is designed to simulate both radiometric and photometric variables (see the discussion below and equation 16 in Kocifaj et al. 2019). In the present study, the data recorded comprise night-sky luminance, but we also use the more common term 'brightness' or the radiometric term 'radiance' wherever it is needed.

S K Y L U M I NA N C E AT L A R G E D I S TA N C E S F RO M A L I G H T S O U R C E
It has been shown by Kocifaj et al. (2019) that night-sky luminance due to a ground-based light source (city or town) can be approximated for intermediate and/or long distances D by the formula where z is the observational zenith angle, μ = cos(z), I 0 the luminance of a city normalized to its area, i.e. the city interpreted as a point light source (I 0 is measured in cd or lm sr −1 ). The functions Si(x) and Ci(x) introduced above are sine and cosine integrals, respectively (Press et al. 2007).The argument x is the product of distance D and parameter a, with where μ 0 is the mean effective cosine of the emission angle (see section 2 of Kocifaj et al. 2019) andH is the altitude up to which a homogeneous atmosphere would extend. Here we assume thatH is composed of the weighed contributions of molecular and aerosol constituents, i.e.
whereτ R andτ a are the Rayleigh and aerosol components of the total optical depth of a cloudless atmosphere except for ozone or watervapour absorption bands, i.e.τ =τ R +τ a (Utrillas et al. 2000). Here, H R and H a are the Rayleigh and aerosol components of the atmospheric scaleheights. Aerosol scaleheight typically ranges from 1-3 km (Wu et al. 2011), while H R = 8 km is a commonly accepted value in many studies (see e.g. Waquet et al. 2009). The emission function for a city, i.e. the luminous intensity emitted by that city, is generally unknown, so an isotropic-like radiator appears to be a useful zero-order approximation, resulting in the following formula: where W 0 (lm) is the total amount of lumens directed above the horizontal; consequently, it is obtained as a superposition of both direct emission upwards and isotropic reflection from the ground. If the entirety of light approaching the sky is due only to diffuse reflection from the ground, then W 0 = G · W installed 0 , where G is the diffuse reflectance and W installed 0 is the total lumen output from all light sources in the city. A potential error in J (z) can result from imperfections of the approximate theoretical model, or from uncertainty in some key parameters, such as W 0 and μ 0 , which are notoriously unknown. The optical thicknessτ and/or scaleheightH can be either measured or inferred by indirect methods. Nevertheless, these parameters are only known with limited accuracy, because the aerosol content is highly variable in time and can change significantly with the point of observation. Such variability links to atmospheric processes, including their dynamics, chemistry, humidity, turbulence, etc. The correct procedure to treat the errors in J (z) is to combine the errors of particular variables in quadrature, i.e.
After a bit of manipulation, we obtain Here, is valid for an arbitrary x. For large values of x, the function of η(x) approaches asymptotically to −1, so J V relative to the theoretical luminance J V , as derived equivalently to equation (6), would read The above formula is applicable to sky luminance observed from intermediate or large distances (for aD>5). Normally, the main goal is to estimate the error of luminous flux W 0 , using the experimental data of sky luminance (J V ). The upper threshold for the error of W 0 , analogously to equation (8), is defined as Note that aerosol optical depth is by far the greatest modulator ofH , so we can easily find from equation (3) that The Rayleigh component (τ R ) of the atmospheric optical depth at the ground is known to vary only slightly with air pressure and temperature, so τ R τ a . Therefore, the uncertainty in the total optical depth is due mostly to uncertainty in its aerosol component, i.e. τ ∼ = τ a . Then equation (9) transforms to the new form Taking into account that H a is typically less than half ofH , the relative error W 0 /W 0 never exceeds the value On the one hand, the theoretical derivations in Kocifaj et al. (2019) dictate that the above equations cannot be used at short distances from a light source. On the other hand, it also follows from the second term under the square root of equation (12) that the error in W 0 increases as D approaches very large values. Consequently, the optimum conditions to retrieve W 0 can be achieved if luminance data are captured at an intermediate distance, i.e. at D satisfying the condition (τ D/2H ) < 1. In such a case, the relative error due to aerosol optical depth uncertainty is smaller than that due to potential measurement errors of J V (comparing the first and second terms under the square root of equation 12). The data need to be collected under clear-sky conditions and ideally when the atmospheric turbidity is low. For a clear-sky optical depth of approximatelyτ ≈0.2 at the scotopic peak andH ≈ 5 km, the above condition dictates D < 50 km. Assuming the luminance data are taken at distances of a few tens of kilometres from a city, the aforementioned second term in equation (12) becomes negligible compared with other terms. As an example, for D = 20 km the square ofτ D/2H is 0.16, resulting in a weighting factor of aerosol uncertainty six times smaller than that of the first term on the right-hand side of equation (12). The scaling parameter μ 0 is unknown and has to be obtained from the least-squares method , and the relative error μ 0 /μ 0 is computed from the sample standard deviation by matching modelled and observed data. The inequality μ 0 /(μ + μ 0 ) < 1 is valid for all combinations of μ 0 and μ, because both of these parameters range from 0-1. Therefore, we easily find that The relation can be applied, if the relative error of the measured luminance data is below 10 per cent and the same is valid for the sample standard deviation and optical depth uncertainty, respectively. The error margin of lumen output retrieval would be reduced to less than 10 per cent, with input data showing an uncertainty within 5 per cent. Nevertheless, the optical depth is frequently unknown or inferred from the visibility using Koschmieder theory (Molnár et al. 2016). The following formula that links aerosol optical depth τ a at a specific wavelength λ and visibility V has been introduced by Wu et al. (2014): where h is the station altitude, H 1 = 0.866 + 0.222V , H 2 = 3.77 km, λ = 0.55 μm and ν is the so-called Ångström exponent (Wagner and Silva 2008), being most typically ν = 3. The value of 3.912 is known as the Koschmieder constant for a distant target contrast equal to 0.02, a common reference value for the human detection threshold (Horvath 1995). The estimation of the limiting magnitude of visibility is too subjective for an individual observer, thus provoking potentially large errors in determination of τ a . It is therefore not unlikely that τ a /τ a could be as large as some tens of per cent. Assuming the same applies to μ 0 /μ 0 , the uncertainty in the total lumen output can rise up to 40 per cent or even more. Indeed, such an error margin yields a useful outcome, considering the fact that lumen outputs (W 0 ) for two different cities can differ by several orders of magnitude. Consequently, an estimate of W 0 within several tens of per cent of its correct value can generally be regarded as a suitable result.

D E R I VAT I O N O F T H E I N F L U E N C E C AU S E D B Y S M A L L / M E D I U M -S I Z E D C I T I E S
The great challenge we face is to distinguish between effects from possible existing dominating light-emitting sources and smaller ones. However, to model the night-sky brightness (NSB) distribution in all zenith angles z and azimuth angles A due solely to dominant light sources, we found that can easily be used for N number i of light-emitting sources at all scattering angles θ i . G(z) describes the gradation function, which depends exclusively on zenith angle and atmospheric transmission function, with the latter being related to the optical depth. Furthermore, J i is the brightness (luminance or radiance) of the ith source of light. As derived in Kittler, Kocifaj and Darula (2012), the scattering indicatrix S(θ ) is defined as with c 1 , c 2 and d being scaling constants. Equipped with substantial knowledge of azimuths A i of all i apparent light sources, the scattering angles, as illustrated in Fig. 1, can be computed as MNRAS 494, 5008-5017 (2020) Downloaded from https://academic.oup.com/mnras/article-abstract/494/4/5008/5831326 by guest on 08 May 2020

Field measurements
Field measurements were carried out in mid-June 2018 at locations in Eastern Austria, more precisely in a small area around the National Park Neusiedler See-Seewinkel. Care was taken to keep distances to the occurring dominating light sources, Vienna and Bratislava, high, i.e. 66 and 49 km, respectively. In order to provide an example of our approach, we show the derivations introduced before as calculated for one location, situated in the triangle of Vienna-Bratislava-Sopron. The observation point is illustrated in Fig. 2. Measurements were taken using a digital single-lens reflex camera (Canon EOS 5D Mark II) equipped with a fisheye lens. The analysis of captured images, yielding a matrix indicating the luminance over the whole sky, was performed by the software 'SKY QUALITY CAMERA' (SQC), developed by Andrej Mohar from Euromix Ltd., Slovenia. Fig. 3 shows the input image as captured in cyclic format and the resulting luminance analysis. In the latter, light domes and their impact on the skyglow near the horizon become apparent. Since SQC provides the possibility of collecting data for the whole sky (i.e. a matrix consisting of zenith angles, azimuth angles and luminance values), fixed azimuths only, e.g. city centres, can also be examined.

Combining measurements with a theoretical approach
As can be seen from the radiance maps in Figs 2 and 3, there are two dominating cities, Bratislava and Vienna. Having collected luminance data from the azimuths of both city centres -Bratislava at A B = 18 • , Vienna at A V = 322 • -the first step is to model the NSB at a local meridian that intersects a horizontal circle at the azimuth position of a radiating city. For Bratislava, i.e. θ B = (π/2) − z, we have In analogy, for Vienna, i.e. θ V = (π/2) − z, this results in The scattering angle θ z , as illustrated in Fig. 4, is computed as follows: For all zenith angles, the products G(z)J B and G(z)J V are generally unknown, but can be determined from the values NSB V (z, A V ) and NSB B (z, A B ) as measured at local meridians of Vienna and Bratislava. Specifically, provides the possibility of computing G(z) (23) Figure 4. The same as in Fig. 1, but for scattering angle θ z , which is the angle between the radius vector of Vienna (OC V ) and a sky element on the principal plane of Bratislava (OE B ). The above formulae allow modelling of the theoretical NSB distribution (luminance or radiance map) due to two dominating light-emitting cities, Bratislava and Vienna, as Fig. 5 compares the digitized input image, seen in Fig. 3 (lower), with the result of the computation provided by equation (24), which is now modelling the net contribution from both dominant cities.
The contribution from M small or medium-sized cities to the NSB is now obtained by subtracting the theoretical night-sky brightness caused by dominant cities solely from the measured NSB, i.e.
with M being the total number of small cities. Subsequently, the result now provides the contribution of small/medium-sized cities. For the location of our observation point, we chose three smaller cities for exemplar computation: Mosonmagyaróvár (A center ≈61 • ), Sopron (A center ≈258 • ) and Eisenstadt (A center ≈293 • ), so M = 3. The scaling coefficients used for calculation of the S function are: c 1 = 16, c 2 = 0.3 and d = −3 (Kittler et al. 2012). The result is shown in Fig. 6.
To estimate the amount W V0 of produced lumens radiated upwards by the cities, we use the formula introduced by (Garstang 1986), as follows:  Table 1.  . Satellite data obtained from VIIRS taken on 2018 May 17, for Eisenstadt, Sopron and Mosonmagyaróvár. Only pixels brighter than 10 −8 W cm −2 sr −1 were considered, to focus on city centres and pixels radiating rather than scattering light. Taking into account that the satellite detection limit is 5 × 10 −11 W cm −2 sr −1 (Cao & Bai 2014), there is also the need to distinguish between scattered and radiated light. To ensure that only latter was included in our analysis, we focus exclusively on pixels brighter than 10 −8 W cm −2 sr −1 . In the conversion, we assume a uniform usage of high-pressure sodium (HPS) luminaires. In the case of at least Eisenstadt, the assumption of containing mostly HPS luminaires was certainly correct at the time of observations (Wallner 2019). In order to ensure comparability, satellite data taken on 2018 May 17, just shortly before the measurements took place, were collected. Furthermore, we also propose that cities produce a uniform angular output pattern, resulting from the calculation of emissions to the upward hemisphere, i.e.
L mean is the mean luminance of the total area of artificially lit pixels A tot . The results are illustrated in Fig. 7 and the estimations of total lumens installed are given in Table 2.

D I S C U S S I O N
The computations performed show that night-sky imagery via all-sky photometry, linked to the asymptotic formula introduced, can be used to provide information about the influence of small and medium-sized cities on the night-sky luminance. However, when comparing the simulated results with obtained satellite data, i.e. Tables 1 and 2, there are still some differences, e.g. for Eisenstadt the value of total lumens installed is less by 40 per cent than estimated from satellite. As mentioned before in Section 2, an estimate of this value within several tens of per cent of its correct number can generally be applied as a suitable result. The differences found can be traced back to approximations made when satellite data were processed (see Section 3.3). One is uniform emission in all directions, which is clearly not the case seen realistically, but, since the emission function of cities is for the most part unknown, this is a native zero-order approximation. To go into detail, this means that if the satellite obtains data from low elevation angles, i.e. not from directly above a city, values for the theoretical uniform emission increase, since many cities appear to have stronger emission at low angles. Satellites are capable of detecting 'cold' and 'hot' spots, at which radiance is measurably weaker or stronger due to blockings (Li et al. 2019), but nevertheless cannot cover large viewing zenith angles close to 90 • . Furthermore, limitations to HPS luminaires and minimum pixel brightness could also influence the results. Surprisingly, we found Mosonmagyaróvár to emit more light than Sopron, even if only containing half of the population. This outcome, however, appears correct when inspecting the localities in more detail. On the one hand, the surface of artificially lit areas seems to be smaller in Sopron, hence the density of luminaires appears to be higher and provides issues considering less efficient emissionsto low elevation angles, e.g. due to blockings. On the other hand, small towns close to the observation point, e.g. Jánossomorja (18 per cent population of Mosonmagyaróvár) or Mosonszolnok (5 per cent population of Mosonmagyaróvár), could have contributed to the NSB. Consequently, taking all comparisons into account, one might argue that satellites are not very accurate in providing information on light conditions perceived by an individual observer on the ground. However, considering that many approximations and assumptions had to be implemented to process the VIIRS dataset, both results, measurements and computations, are compatible.
In Eisenstadt, LED conversion started in 2018 July, so measurements were taken even before this process has started. Still, due to finished or upcoming LED conversions, VIIRS data cannot be used for comparison if cities use luminaires strongly emitting in short wavelengths. This provides an aggregated situation as regards long-term development of light pollution on Earth's surface. The majority of data gathered originate either from one-dimensional photometric measurements detecting the luminance of the sky's zenith luminosity, logged by so-called 'Sky Quality Meter' devices, or, as mentioned, from satellite data. Both methods will lead to difficulties in the long run, as their spectral sensitivities show weaknesses, being mostly blind to short wavelengths. Even if cities become brighter, by changing the colour of lightemitting sources, existing devices potentially predict a decrease in skyglow. This issue especially makes our approach even more essential for use as a possible measurement standard in the future.
The simulation performed in this work was conducted for two dominant light sources far away from the observation point. Congruously, the approach could be ineffective near dominant light-emitting sources, due to their strong influence on the night-sky luminance and their light domes camouflaging smaller sources, which are indistinguishable during measurement.

C O N C L U S I O N S
The approach used in this article combines in situ measurements via all-sky photometry and theoretical computations using the asymptotic formula introduced in past work. The aim is to provide information about the impact of small and medium-sized cities on the night-sky luminance, even with other occurring dominating sources, and additionally, an estimation of total lumens installed. This could be of great importance to astronomical observatories, since smaller cities or towns located in the surroundings can evolve and their impact on the NSB can increase. A comparison with performed computations was conducted via satellite data and has shown to be coherent, considering that approximations for processing satellite data were inevitable. We think that our approach can be of great value for quantifying and distinguishing different small or medium-sized light-emitting sources, but also especially for studying their long-term development, since, in future, the quality of data from e.g. VIIRS or Luojia 1-01 will decrease, due to their blindness in short wavelengths. The modelling approach can be used for all cities globally, although it is necessary not to obtain measurements inside light domes, to retain the possibility of these being separated from dominating ones within the course of the calculation. In addition, clear-sky conditions are necessary. However, only one measurement is required for inclusion in the theoretical modelling and therefore results are independent of any possible meteorological variations us, which usually occur if obtaining multiple images throughout the night. Another strength of the method used is that usually unknown atmospheric parameters, e.g. aerosol optical depth, which are necessary to include in models, need not be known, since the implementation of measurements from real conditions contains information about the environmental conditions. Consequently, the expenditure of computational time is lowered significantly by achieving even greater accuracy.