Abstract
This paper presents a sky luminance model for the tropics. The model is mathematically expressed as a multiplication of two functions. These are , which is a function of the zenith angle of a sky element and solar zenith angle, and , which is a function of the angle between the sky element and the sun. To obtain the analytical forms of these functions, the sky luminance data collected at Nakhon Pathom (13.82°N, 100.04°E), Thailand, during a four-year period were analyzed. Additionally, satellite-derived cloud index at the position of the sky luminance measurements during the same period was estimated. Based on values of the cloud index, the skies were classified into 10 types, from clear to overcast skies. By using appropriate grouping and mathematical operation of the sky luminance data, the values of the two functions were obtained and then fitted with empirical equations. The multiplication of these equations gives the final form of the sky luminance model. To validate the model, it was used to calculate the relative sky luminance at other three sites in the tropics. It was found that values of relative sky luminance calculated from the model and those obtained from the measurements were in reasonable agreement.
1. Introduction
A utilization of daylight to illuminate building interior helps to save electricity for lighting [1]. For this reason, daylight-integrated buildings and daylight devices have been actively developed in many countries [2–4]. To design a daylight-integrated building, information on daylight environment of the building is generally required. One of the most important information required is the daylight sky luminance.
Sky luminance can be measured by using a sky scanner. Ideally, information on sky luminance should be obtained from a dense network of sky scanners, but due to costs, only a few instruments are routinely deployed in most parts of the world, and limited works on sky luminance measurements have been reported [5–8]. As a result, sky luminance has to be obtained from model calculations. Therefore, a number of sky luminance models have been proposed to estimate sky luminance. These can be summarized as follows.
Pokrowski [9] was the first researcher who proposed that the angular distribution of clear sky luminance depended on four parameters: a solar zenith angle, a zenith angle of a sky element, an angle between the sun and the sky element, and zenith luminance. Then a sky luminance model was formulated as a function of these parameters. Kittler [10] further developed a sky luminance model for the clear sky by expressing the relative sky luminance as a product of two functions, namely, indicatrix and gradation functions. For the case of the overcast sky, Moon and Spencer [11] proposed a model for calculating the relative sky luminance for this sky condition. Later on, a number of sky luminance models have been proposed [12–23]. As sky luminance depends on sky conditions, a number of researchers have addressed this problem. These include the modeling work of Brunger [13], Harrison and Coombes [15], Perez et al. [17, 18], Matsuura and Iwata [16], and Igawa and Nakamura [20]. Many investigators have tested some of these models against measured data, and a wide range of discrepancies have been reported [24–28].
As sky luminance depends strongly on sky types, Kittler et al. [19] classified sky conditions into 15 types and proposed a general model for estimating sky luminance for all sky types. The values of the coefficients of the model depend on the sky types. Afterwards, Commission Internationale d’Eclairage (CIE) has adopted this model as the CIE standard general sky [21]. To use the CIE model, it is necessary to identify the sky types in order to select the values of the model coefficients. The selection is based on the description of the luminance distribution given by CIE [21], which practically depends on the skill of observers.
To avoid the use of solar radiation and other meteorological data and the problem of subjectivity to identify sky types, we proposed a novel method to identify sky conditions by using satellite-derived cloud index, and a new formula for estimating sky luminance was established using this identification method.
The model was aimed to use in the tropics. The tropical environment offers a different perspective to sky luminance climatology when compared to other environments. In the tropics, the solar altitude angle at local noon is high or at overhead, and cloud structures are different from those of mid and high latitudes, with high frequency of cumulonimbus clouds. In developing the model, a 4-year period of luminance data obtained from the measurement at one station in a tropical environment of Thailand and concurrent satellite data were used to formulate the model. To ensure the generality of the model for the tropics, sky luminance independent data obtained from other 3 stations situated in the tropical environment were used to validate the model. The distance between these stations varies in the range of 600–1,500 km.
2. Measurements and Data Processing
2.1. Sky Luminance Measurements
To obtain necessary data for this work, sky luminance was measured at our four existing solar radiation monitoring stations located in a tropical environment of Thailand. These stations are at Chiang Mai (18.78°N, 98.98°E) in the northern region of the country, Ubon Ratchathani (15.25°N, 104.87°E) in the northeastern region, Nakhon Pathom (13.82°N, 100.04°E) in the central region, and Songkhla (7.20°N, 100.60°E) in the southern region (see Figure 1). For each station, a sky scanner was used to measure sky luminance. It scanned a total luminance of 145 sky elements at the view angle of 5° across the sky hemisphere. These 145 elements of sky luminance were made according to the decomposition of Tregenza sky [29]. The luminance measurements were taken every 30 minutes, which each scan was completed in 3 minutes. A scanning sequence carried out at the zenith angles of 6°, 18°, 30°, 42°, 54°, 66°, 78°, and 90°. The periods of data from each station used in this work are shown in Table 1. As Nakhon Pathom station is located at the central region of the country and it has the longest period of data (4 years), the luminance data from this station were used to formulate the sky luminance model, whereas the luminance data from Chiang Mai, Ubon Ratchathani, and Songkhla were employed for the model validation. Only the sky luminance data taken at the middle of the hour (i.e., 8:30–16:30) were used in the analysis. All sky luminance data were subjected to a quality control. To identify the abnormal sky luminance data, threshold values of sky luminance were set for clear, partly cloudy, and overcast skies. These values were estimated by using the luminance values calculated from the CIE standard general sky [21], as a guideline. The sky luminance data from the four stations were compared with the threshold values, and the luminance data whose values are abnormally deviated from the threshold values were considered as suspected abnormal data. As we have also sky images from sky cameras at the four stations (Figure 1), these images were also used to recheck these abnormal data. The confirmed abnormal data were discarded from the data sets.
As the model is aimed to estimate the sky luminance, not the direct luminance from the sun, the luminance data at the position of the sun and the sky elements close to the sun disk were also discarded from the data sets.
2.2. Processing of Satellite Data
The satellite data used in this work are 8-bit digital data from a visible channel of three geostationary meteorological satellites: GMS5 (January–May, 2003), GOES9 (June, 2003–July, 2005), and MTSAT-1R (August, 2005–December, 2006 and September–December, 2011). In total, the satellite data for the period of 4.3 years were used in this work. The sensor of the satellites measured solar radiance reflected from the earth-atmospheric system. The radiance was converted into irradiance and used for the calculation of reflectivity by assuming that the earth-atmospheric system is a Lambertian surface. This assumption is true only for long-term average basis [30]. As our study uses the satellite data for the period of 4 years in the modeling process to identify the sky condition, on average basis as described in the next section, this assumption is considered to be valid for this work. In addition, the satellite images taken in the early morning and late afternoon were not used to avoid the non-Lambertian effect which causes error in the calculation of the earth-atmospheric reflectivity. Only the satellite images taken at 8:30, 9:30, 10:30, 11:30, 12:30, 13:30, 14:30, 15:30, and 16:30 local time were used in this work. These satellite images have a spatial resolution of 3 km × 3 km. The original forms of the images are in a satellite projection showing a curvature of the earth’s surface. They were remapped into a cylindrical projection, linear in latitude and longitude lines. To control the quality of these data, a visual inspection of all satellite image data was undertaken, and the image data with distortion were discarded from the satellite data set. Each satellite image comprises a matrix of pixels () covering the entire area of Thailand. Each pixel contains information of a gray level varying in the range from 0 to 255 digital counts, which is related to the reflectivity of the earth-atmospheric system. After the map transformation and quality control, they were cut into matrices of pixels centered at the positions of Nakhon Pathom, Chiang Mai, Ubon Ratchathani, and Songkhla measuring stations. Each matrix of pixels covers most parts of the sky over these target stations. A calibration tables supplied by a satellite data agency were used to convert the gray levels into earth-atmospheric reflectivity ().
2.3. Processing of Sky Luminance Data
Although, the satellite-derived earth-atmospheric reflectivity () can be used to identify sky conditions [31], the values of vary only in a small range (from 0.1 to 0.4), resulting in coarse sky classifications. In addition, the effect of the ground reflectivity is also included in . As a results, not only depends on the sky conditions but also on the ground reflectivity. Therefore, to avoid the effect of the ground reflectivity and to increase the precision of the sky classification, we proposed to use the satellite-derived cloud index () to classify the sky conditions. This cloud index (), which was first introduced by Cano et al. [32], was defined as where , , and are the earth-atmospheric reflectivity, ground reflectivity, and maximum cloud reflectivity, respectively. was derived from satellite data, as described in the preceding section. The ground reflectivity was derived from cloud-free satellite image using the method explained by Janjai et al. [33]. For the maximum cloud albedo, , it was calculated from the maximum value of the gray level of the satellite data. As the study area is situated in the tropics with the surface covered by green vegetation year round and the frequent occurrence of tropical clouds, the typical values of ground and maximum cloud reflectivities are 0.12 and 0.8, respectively. In addition, most skies in this region are under partly cloudy conditions, resulting in the typical values of the cloud index and the earth-atmospheric reflectivity of 0.5 and 0.46, respectively.
In this work, a matrix of pixels of the satellite data centered at a target luminance measuring stations was used to indicate the sky type over the target stations. From the size of each pixel, this matrix of pixels covers most parts of the sky at the stations. In the first step, the cloud index at each pixel was calculated by using (1). From (1), in the pixel without a cloud, is equal to and is equal to zero. On the contrary, in a pixel which is totally covered by cloud, is equal to and is equal one. The pixel which is partly covered by cloud, the value of the cloud index will be between these two extreme cases. Therefore, the cloud index indicates the contribution of cloud to the earth-atmospheric reflectivity through its reflective property. In the second step, the values of cloud index of all pixels in the matrix were averaged, and the average value was used to represent the cloud index of the sky. As there are sky cameras at the four stations which routinely take photos of the sky, the images of the sky were also used to visually check the sky types indicated by the cloud index. It was found that when the cloud index is equal to zero, the image obtained from the sky camera shows clear sky condition, and when the cloud index is equal to one, the image shows overcast sky condition. In addition, for the case of , the sky is clearly seen from the image as partly cloudy condition. This means that the sky types indicated by the cloud index agree well with those observed from the images of the sky.
Concurrently, the corresponding sky luminance data obtained from the sky scanner at Nakhon Pathom station were also gathered. The data of each scan consist of values of sky luminance from 145 sky elements. An example of these data and satellite data is shown in Figure 2. In the next step, these 145-sky element luminance data from all scans were classified into 10 groups according to the cloud index , (range from 0 to 1, with the step of 0.1), prevailing at the time of each scan.
The luminance data in each cloud index group were separated into 15 subgroups according to the solar zenith angle () prevailing at time of the scan. The intervals of solar zenith angle for these subgroups are , , , , , , , , , , , , , , and . The luminance data with were not used due to the lack of the corresponding satellite data. For each subgroup, the value of the solar zenith angle at the middle of the interval was used to represent the solar zenith angle of that subgroup. In total, there were 150 subgroups of sky luminance data sets (10 sky types × 15 solar zenith angle intervals).
One sky luminance data set was obtained from one scan. In each subgroup, there are several sky luminance data sets which have the same intervals of the cloud index and solar zenith angle. For a given interval of the cloud index, the data sets which have the same interval of solar zenith angle were averaged. As there are gaps between the sky elements, an interpolation was used to find the values of the sky luminance in the gaps. The sky luminance values obtained from this process were employed in the modeling described in the next section.
3. Modeling
As it is difficult to find a function to represent mathematically the absolute sky luminance (cd·m−2), most investigators [10, 18, 19] preferred to formulate a model for a relative sky luminance () which is defined as a ratio of the absolute sky luminance () to the zenith luminance (). In this study, we followed the approach of Kittler et al. [19] to model the relative sky luminance by expressing it as a multiplication of two normalized functions as where and are gradation and indicatrix functions, respectively. is the zenith angle of the sky element, is the indicatrix angle of the sky elements, and is the zenith angle of the sun (see Figure 3). The indicatrix angle can be calculated by using the following equation [21]: where and are the azimuths of the sky element and the sun, respectively.
Although the relative sky luminance is dependent of functions and , by using an appropriate grouping and mathematical operation of the sky luminance data, mathematical expression of and can be separately determined. To obtain these mathematical expressions, the sky luminance data processed in Section 2 were analysed as follows.
For the given values of a solar zenith angle and a cloud index, the sky was partitioned into 30 concentric zones around the sun or “concentric zones of .” These zones have the indicatrix angle of , , , , , , , , , , , , , , , , , , , , , , , , , , , and .
As mentioned earlier, the luminance data corresponding to the value of were not included in the data set because these data come from the area of the sun disk and the nearest area of the sky around the sun, and the model to be constructed is not for use for these areas. Therefore, the sky luminance starts from.
In addition, the sky was also partitioned into 17 concentric zones around the sky zenith or “concentric zones of .” These zones have the zenith angle , , , , , , , , , , , , , , , and .
The partition step of 5.0° for and was chosen because it corresponds to the angular diameter of each sky element observed by the sky scanners. The interception between the concentric zones of and the concentric zones of gives subzones of the sky (see Figure 4).
The sky luminance of the subzone was denoted as . An example of the values of the sky luminance of some subzones is shown in Table 2.
The determination of an empirical formula for started with the evaluation of relative sky luminance of all subzones in the concentric zone . In this concentric zone, there are 17 subzones, namely, . The subzone was arbitrarily selected as a reference subzone, and it was renamed to be . The sky luminance of this reference subzone was denoted as . By applying (2) to subzones , the relative sky luminance of these subzones can be written as follows:
Note that (7) is the equation for the reference subzone .
In the second step, (4)–(6) and (8)-(9) were divided by (7) of the reference subzone to obtain the following:
Note that there is no equation for sky luminance of the subzone because this subzone was selected as a reference subzone.
As the values of in the left hand side of (10) were obtained from the measurements, the values of were consequently known. Then the process was repeated for the concentric zones . In the third step, the values of obtained from these concentric zones were plotted against , as an example shown in Figure 5.
It is observed that the plots of versus obtained from all concentric zones of are similar, and they can be represented with one graph. The -intercept of this graph gives the value of . In order to eliminate , all values of in the graph was divided by the value of to obtain the values of . Then the values of were plotted against , as shown in Figure 6. The process was repeated for all values of and ; therefore, a set of graphs versus is obtained. An example of these graphs is shown in Figure 7.
In the final step, the graphs versus for all values of and were fitted with an empirical function, and the best-fitted equation is where and are expressed in radian. The values of the coefficient , , , and depend on the sky types indicated by values of the cloud index, as shown in Table 3. Note that the fitting was carried out for each interval value of and the coefficient of determination () of the fitting is also shown in Table 3.
To find the mathematical expression for , the sky luminance in the concentric zones of was analyzed. The analysis started with the concentric zone (Figure 4). The subzone () in the concentric zone was assigned as a reference subzone with the sky luminance . Then, the relative sky luminance in each subzone of the concentric zone was calculated by using (2) as follows: In order to eliminate , (12), (14), (15) were divided by (13) to obtain
The values of terms in the left hand side of (16) were obtained from the sky luminance measurements. The process was repeated for the concentric zone and for all solar zenith angles (). Then the values of were plotted against the values of , and the results are shown Figure 8.
It is noticed that the plots versus are independent of the solar zenith angle (), and the plots can be represented by one graph. The process was repeated for all cloud indices. The results are graphically shown in Figure 9.
The graphs in Figure 9 can be best fitted with the following empirical equation: where is expressed in radian. The values of the coefficients , and obtained from the fitting are shown in Table 4.
For the case of , (17) becomes By dividing (17) with (18), was obtained as
The last step of the modeling is the combination of in (11) with in (19) to obtain the final expression for the relative sky luminance as follows:
For a given set of the coefficient determined by the cloud index, (20) will give a value of . An example of values of across the sky for clear, intermediate, and overcast skies is shown in Figure 10. As expected, the graph for the clear sky is above the graph for the intermediate sky and the graph for the overcast sky is at the bottom.
4. Validation of the Model
The model was validated against the sky luminance data obtained from the sky scanners at three measuring stations in the tropical environments of Thailand. These are Chiang Mai (18.78°N, 98.98°E) in the northern region, Ubon Ratchathani (15.25°N, 104.87°E) in the northeastern region, and Songkhla (7.2°N, 100.60°E) in the southern region (see Figure 1). The distances between these stations are 600–1,500 km. The data from these stations are independent data sets because they were not involved in the model formulation. Values of cloud index over the stations were derived from the satellite data and used to select the coefficient values (Tables 3 and 4). As the model was derived from the spatial average basis and the sky luminance data from sky zones are practically used, the sky luminance data are partitioned into 25 sky zones uniformly distributed in the sky hemisphere and all data within each zone were averaged. The position of the zones is shown in Figure 11. The number of zones was chosen so that it can represent sufficiently the luminance from different parts of the sky. This partition follows the method used in the SATEL-LIGHT project which produced sky luminance of Europe [34]. In the next step, the model (20) was used to calculate the relative sky luminance in each zone. Then the values of the relative sky luminance calculated from the model were plotted against those obtained from the measurements. The results are shown in Figures 12, 13, and 14.
Although the comparison graphs are scattered, the data points exhibit a linear trend around the 1 : 1 line for all stations. The discrepancies between the predicted and the measured relative sky luminances in terms of root mean square difference (RMSD) and mean bias difference (MBD) are 15.1% and 5.1% for Chiang Mai, 15.1% and 4.9% for Songkhla, and 16.6% and 6.7% for Ubon Ratchathani, respectively. Although, the validation was made for three stations in Thailand, the distances between these stations are 600–1,500 km. From this validation result, the model is valid for a wide range of distance in the tropics. It is expected that the model will be valid for other tropical sites. However, it will be interesting to test the model against the measurements in other tropical locations.
Ineichen et al. [25] have tested the performance of seven sky luminance models against four sky luminance data banks obtained from measurements. These are the models proposed by Perez et al. [18], Brunger and Hooper [35], Perez et al. [24], Matsuura and Iwata [16], Harrison [36], Perraudeau [37], and Kittler [38]. They found that values of RMSD varied from 27% for the model of Perez et al. [18] to 65% for the model of Perraudeau [37]. Additionally, values of MBD also varied in a wide range (−4 to 14%), depending on the models, data banks, and conditions for the elimination of the data. Although, the model developed in this study gives relatively lower RMSD and MBD as compared to those of the above-mentioned models for most cases, the comparison cannot be made directly due to differences in the partition of the sky data sets and the conditions for the elimination of the data. In addition to its reasonable accuracy, our model has the advantage that it needs only satellite-derived cloud index as an input. This cloud index can be derived from meteorological satellite data which are available for most parts of the world.
5. Conclusion
A model for estimating a sky luminance in the tropics using a satellite-derived cloud index to identify sky conditions has been developed. Satellite data from GMS5, GOES9, and MTSAT-1R satellites were used in this study. A four-year period of concurrent satellite and sky luminance data collected at Nakhon Pathom, Thailand were used to formulate the model. The model is expressed as a multiplication of two empirical functions. The zenith angle of the sky element, the angle between the sky element and the sun, and the zenith angle of the sun were used as independent variables of the functions. The values of the model coefficients are given according to sky types indicated by the values of the cloud index. To validate its performance, the model was used to calculate the relative sky luminance at Chiang Mai, Ubon Ratchathani and Songkhla. The calculated and the measured relative sky luminances are in reasonable agreement with RMSD and MBD of 15.1% and 5.1% for Chiang Mai, 15.1% and 4.9% for Songkhla, and 16.6% and 6.7% for Ubon Ratchathani, respectively. With its reasonable accuracy and simplicity, this proposed model is expected to be useful for designing a daylight-integrated building in a tropical environment.
Acknowledgments
The author would like to thank Silpakorn University Research and Development Institute for the financial support to this research. The author would also like to thank the regional offices of Thai Meteorological Department for supporting the measurements of the sky luminance. Dr. Manuel Nunez is gratefully acknowledged for his valuable advice. Finally, the author thanks Dr. Somjet Pattarapanitchai for his assistance in the processing of the data.