Large‐scale spatiotemporal calculation of photovoltaic capacity factors using ray tracing: A case study in urban environments

Photovoltaics (PVs) on building facades, either building‐integrated or building‐attached, offer a large energy yield potential especially in densely populated urban areas. Targeting this potential requires the availability of planning tools such as insolation forecasts. However, calculating the PV potential of facade surfaces in an urban environment is challenging. Complex time‐dependent shadowing and light reflections must be considered. In this contribution, we present fast ray tracing calculations for insolation forecasts in large urban environments using clustering of Sun positions into typical days. We use our approach to determine time resolved PV capacity factors for rooftops and facades in a wide variety of environments, which is particularly useful for energy system analyses. The advantage of our approach is that the determined capacity factors for one geographic location can be easily extended to larger geographic regions. In this contribution, we perform calculations in three exemplary environments and extend the results globally. Especially for facade surfaces, we find that there is a pronounced intra‐day and also seasonal distribution of PV potentials that strongly depends on the degree of latitude. The consideration of light reflections in our ray tracing approach causes an increase in calculated full load hours for facade surfaces between 10% and 25% for most geographical locations.


| INTRODUCTION
For achieving the goals of the Paris Agreement, a significant reduction in greenhouse gas emissions is necessary.This requires a transition of the energy system towards renewable energy sources.Low levelized costs of electricity (LCOE) make photovoltaic (PV) energy generation a key technology in the energy system transition. 1 PVs in cities will become of particular importance, since two thirds of the world's population are expected to live in urban areas by 2050. 2 However, with the ongoing densification of urban areas and increasing energy consumption, suitable space for utility-scale PV installations will become increasingly scarce.Therefore, the market share of buildingattached PV (BAPV) and building-integrated PV (BIPV) installations are expected to increase in the future. 3Facade PV installations will likely make an important contribution to the future energy system in addition to the rooftop PV systems that are predominantly used today.The reason is that the ratio of available roof area to living area is comparatively small in densely populated urban areas.In addition, facades offer a significantly larger area potential compared to rooftops. 4nergy system models are an important tool for advancing the energy transition.They are widely used to evaluate the technical, social, and economic feasibility of decarbonized systems with a high share of PV. [5][6][7][8] However, since PV generation is fluctuating in time, spatially resolved and time-resolved PV potential data are required.
The focus of our simulations is not on the PV potential of an individual building but rather the PV potential of a large set of buildings in a large urban area for supporting the energy transition of districts and other larger administrative levels.However, as a by-product, the most attractive sites and regions can also be extracted from the results.
Calculating the yield of PV installations in urban environment is challenging.Such installations are often exposed to complex and timedependent shadowing of insolation.This holds true especially for installations on facades.4][15] More advanced techniques use a model of the surroundings to take objects causing shadowing of both direct and diffuse insolation into account. 16,17A widely used approach is the calculation of a skyline profile that shows in which direction seen from the PV module an object shades the sunlight. 18These models work well as long as light reflections in the surroundings can be neglected for the most part, which is usually fulfilled for rooftop installations.However, for the calculation of the insolation onto a facade installation, the light reflected in the surroundings plays a major role. 19Ray tracing calculations provide an accurate solution to this problem, as they can account for both shadowing and multiple light reflections even in very complex environments.However, a practical disadvantage of the ray tracing approach is the comparatively high computational demand.
Recently, Ernst et al. 20 published a ray tracing approach that uses binning of Sun positions for the reduction of computational demand.We extend this approach by using clustering of Sun positions and show that the introduction of relative irradiances can be used for global PV potential calculations.This is because, in addition to reducing the computational time required, the use of relative irradiances offers the advantage that the results become independent of weather data and can be timeshifted along latitudes.
In general, there are two approaches for the calculation of PV potentials at larger scales.There is a bottom-up approach, in which PV feed-in is calculated for each PV-installation individually and subsequently aggregated across districts, counties, or even countries. 21,22is approach is hampered on the one hand by 3D building data not being consistently available and on the other hand by the large computing time required.4][25] The difficulty of this approach lies in the identification of suitable indicators.Shading can only be considered using a timeindependent shadowing factor here.Recently, Joshi et al. 26  In this work, we pay special attention to facade PV and present a hybrid approach using an advanced ray tracing technique that allows to generalize the results for locations around the globe.Our approach combines the accuracy of ray tracing for the calculation of time resolved PV potentials on a microscopic scale with the possibility of large-scale modeling for use in energy system simulations.This hybrid method allows for a top-down approach, however without the disadvantage of neglecting individual effects like shading that can have a large impact on PV power generation, in particular on the time dependence.

| METHODOLOGY
We define a set of exemplary urban environments with different building densities and building structures and calculate the time resolved PV potentials for each surface of each building while accounting for shading and multiple light reflections.We then extend the results for each environment globally, which results in temporally resolved PV capacity factors for all surfaces in each global region.On top of that, if we assume the building structure in a specific region can be described by a set of representative environments, the PV potentials in that region can be described by a weighted average over the results for each environment in a top-down manner.

| Detailed calculation of PV yield potential using ray tracing
We implement a ray tracing toolbox (CityPV) that calculates diffuse and direct insolation on arbitrary surfaces in complex urban environments based on publicly available geographic data and meteorological timeseries.At first, we build the simulation scene from geographic data.
Information about building geometries is widely available in different level of details (LoD).In this study, we use LoD-2 data, which consist of building models with standard roof structures.Subsequently, the simulation scene is turned into a three-dimensional mesh where each surface is described by a set of triangles.Since a plane is uniquely defined by three points that do not lie on a single line, the ray tracer algorithm works on triangles as primitive structures of the mesh.Each triangle is assigned with an individual albedo value and a scattering characteristic that determines the intensity and direction of the reflected light.The ray tracing algorithm uses a Monte-Carlo integration of the render equation 27,28 ; thus, a large number of light rays are traced, and the interactions between light rays and the environment are recorded.In a basic ray tracing approach, each generated ray must be checked against each triangle in the simulation scene.Thus, the required computation time increases linearly with both the number of rays and the number of triangles.At the same time, the accuracy of the results increases with the number of rays, 29 and the spatial resolution increases with the number of triangles.Hence, accurate high-resolution calculations on large scales are time-consuming and can thereby become impossible in practice.In order to solve this problem and at the same time obtain both accurate and spatially finely resolved results, we implement a bounding volume hierarchy (BVH), which holds the data of all triangles, 30 and use graphic processors (GPU) to benefit from a high degree of parallelization. 31With this approach, the required computation time increases only logarithmically with the number of triangles, which is a great advantage especially for large geometries.For example, on our workstation (NVIDIA RTX 6000), our implementation traces 1.1 Â 10 9 primary rays of direct sunlight per second in a geometry that consists of 8,700 triangles.Increasing the number of triangles in the same geometry to about 35,000, our implementation still traces 0.94 Â 10 9 primary rays per second.
Throughout this study, we use a total of 1 Â 10 9 rays for each Sun position for direct irradiance and 40 Â 10 9 rays for diffuse irradiance.
The actual ray tracing consists of three steps: (i) ray generation, (ii) ray casting and hit detection, and (iii) ray reflection.During ray generation, the starting point and direction of each ray are determined.In the second step, the rays are cast into the simulation scene, and the closest hit of each ray with the triangles of the scene is recorded.
Afterwards, the properties of the triangle hit by the ray determine the reflection probability and the direction into which the light is reflected.Light reflection is handled according to a given bidirectional reflection distribution function (BRDF). 32In the simplest case, the BRDF is constant for all incoming and outgoing ray directions.This models ideal Lambertian scattering, where the scattering object appears equally bright regardless of the viewing direction.The process of ray casting, hit detection, and ray reflection is repeated with all rays until the maximum number of light reflections is reached.Within this study, we include up to three light reflections for each ray and assume each surface with an albedo of α = 0.25.This albedo value is selected because it is in the range of typical albedo values for materials used in urban areas. 33Considering more than three light reflections has a minor impact on the results, because the ratio of non-absorbed rays drops to α 3 ≈ 1.6% after the third reflection.
In CityPV, ray tracing is performed independently for diffuse and direct insolation.For the calculation of diffuse insolation, we assume an isotropic sky hemisphere. 34Recently, Paul et al. 35 found that isotropic models tend to underestimate the global vertical incident irradiance by 6% to 12% of the peak value of the day.On the other hand, there is the advantage of reduced computing time when using the isotropic sky hemisphere, which makes it favorable for the scope of large-scale energy system simulations.Additionally, the accuracy of anisotropic models often depends on the geographic location that makes comparisons between locations difficult.An implementation of an annual light source as proposed by Ernst et al. 36 could allow the use of anisotropic sky models while maintaining high computational speed.We sample the starting points of the rays uniformly on a sky hemisphere of radius r sphere .Subsequently, the ray directions are randomly sampled according to the BRDF of Lambertian scatter.For each triangle, we normalize the hits per triangle area ρ diff.triangle to the density of hits ρ diff.unshaded a horizontal surface would receive if there were no objects causing shading (cf. the derivation in the appendix), with N diff.raysbeing the total number of generated diffuse rays and I diff.
relative being the normalized diffuse insolation in the plane of the triangle.Because we assume an isotropic sky hemisphere, there is only one calculation of diffuse irradiation that is valid for all Sun positions throughout the year.
In contrast, the calculation of direct insolation depends on the position of the Sun.For each Sun position, we uniformly sample the ray starting points on a disk with radius r disk that is perpendicular to the direction of the Sun's rays.We choose r disk sufficiently large to safely cover the whole scene.The direction of each ray is given by the normal vector of the disk.The hits per triangle area ρ direct.trianglecaused by direct irradiation are normalized to the density of rays created on the disk ρ direct.disk with N direct.raysbeing the total number of generated direct rays on the disk and I direct.relativebeing the normalized direct insolation in the plane of the triangle.The absolute values for diffuse and direct insolation are then calculated by multiplying the normalized insolation values with measured meteorological values for diffuse insolation in the horizontal plane and for direct insolation normal to the Sun's direction, respectively.

| Incorporating meteorological conditions
The irradiation available at the surface of the earth depends on meteorological conditions.Measured values for the global horizontal irradiance (GHI), diffuse horizontal irradiance (DHI), and direct normal irradiance (DNI) are available with high temporal and spatial resolution.For the scope of this study, we use meteorological data from the ERA5 reanalysis weather database, which provides worldwide hourly estimates for the solar radiation among other values from 1959 to the present. 37The ERA5 dataset contains the accumulated surface solar radiation downwards (SSRD) and the accumulated direct solar radiation at the surface (FDIR).The SSRD per accumulation period corresponds to the GHI, and FDIR per accumulation period corresponds to the direct solar radiation on a horizontal plane at the surface.We calculate the DNI and DHI from SSRD and FDIR using the implementation in the software package pvlib 38 according to the following equations: with θ the zenith angle of the Sun.However, if the Sun is close to the horizon, the calculation of DNI is prone to errors because of small values for cos(θ).Therefore, we neglect all Sun positions where the Sun is less than 2 above the horizon.Each result in this study is calcu- There is also diffuse insolation onto the ground plane that is partly reflected onto the PV-module.The ground view factor (GVF) describes the exchange of radiation between the PV-module with tilt ϑ and an infinitely extended ground plane Thus, the total diffuse insolation onto the PV-module can be expressed as In the last equation, we express the total diffuse insolation relative to the DHI so the time dependence vanishes.We calculate the direct insolation onto the PV-module in a similar manner.The direct insolation hitting the PV-module, that is, without prior reflection, depends on the DNI and the cosine of the angle ω between the direction of the Sun and the normal vector of the PV module's surface.Further, the ground reflected direct insolation depends on the albedo, the GVF, and the direct horizontal irradiation Schematic of a PV-module with tilt ϑ on an infinite plane with no objects causing shading.The Sun's zenith angle is θ, and the azimuth angle is ϕ.The angle between the direction of the Sun and the normal vector n of the plane of the PV-module is ω.
with θ being the Sun's zenith angle.By using the relative direct insolation, the time dependence on the meteorological data vanishes, and only the time dependence of ω and θ remains.This dependence comes solely from the movement of the Sun's azimuth angle ϕ and the Sun's zenith angle θ over the course of the year.

| Clustering of Sun positions
In the following, we describe our approach for the clustering of Sun positions to reduce the required computation time.The time dependence of the relative direct insolation requires one ray tracing calculation for each Sun position.For instance, an hourly resolution results in about 4,430 calculations, that is, hours per year where the Sun is above the horizon.We use time series aggregation to accelerate the calculation of this computationally intensive task.We apply a k-means cluster algorithm to the Sun positions, that is, to pairs of azimuth angle and zenith angle using the software package tsam. 40The k-means We use the timeseries of the insolation in the plane of the PV module for the estimation of the DC power output timeseries.For this conversion, we assume a commercially available monocrystalline silicon PV module with a nominal power of P nominal = 300 W p .The conversion is based on the Sandia array performance model. 38The solar cell temperature is estimated assuming an open rack mounting with a glass/polymer module. 41Spatially and temporally resolved air temperature and windspeed data for the solar cell temperature estimation are fetched from the ERA5 reanalysis dataset. 37Please note that the model for the estimation of solar cell temperature is based upon empirical data and is likely to be less reliable for PV module mountings with high tilt.Nevertheless, we use this model as it is easy to use and accounts for the impact of temperature and windspeed on a global scale.Subsequently, we calculate the effective irradiance that takes Sun position dependent spectral loss and angle of incidence loss into account. 41Finally, the DC output power P module is calculated based on the module specifications, the effective irradiance, and the solar cell temperature. 41The capacity factor f is then defined as and thus, the time integral of the capacity factor f over 1 year corresponds to the full load hours (FLH) on the DC side of the PV installation.The full load hours are the ratio between annual energy yield of the PV installation and the nominal installed PV power.

| Exemplary urban environments
In this study, we use a set of three exemplary urban environments depicted in Figure 4  an unshaded facade with the respective orientation on a virtually infinite ground plane.The first environment is an excerpt from a mixed district in the city of Berlin (Germany) with partly single-family houses but also some high-density apartment buildings (Environment A). 42 The second environment (Environment B) is an excerpt from a more densely populated district in the city of Hanover (Germany). 43The last environment (Environment C) is a very densely built-up area with many high-rise buildings in the city of New York City (United States). 44,45All environments are roughly square with an edge length of approximately 1 km.Please note that vegetation is not part of the dataset and thus is neglected in our calculation.The consideration of vegetation would typically result in an overall reduction of PV potentials, in particular for facade PV installations near the ground.
We classify each surface with a tilt angle greater than 89 as a facade and each surface with a tilt angle between 15 and 60 as a roof surface.The roof surfaces serve as a reference for the evaluation of the facade potentials.The mean tilt value for the roof surfaces ranges between 28 and 42 for the investigated orientations and environments.All roof surfaces with a tilt angle less than 15 are disregarded, because our ray tracer calculates the insolation in the plane of the roof surface and PV modules on flat roof surfaces would be tilted in most cases.Additionally, we disregard all surfaces with an area of less than 6.25m 2 .This is because these areas are considered too small to accommodate PV modules.We divide all surfaces into east, south, west, and north orientation according to their surface azimuth.For this classification, we use all surfaces with a surface azimuth angle in the range between À30 and +30 of the respective cardinal direction.

| RESULTS
The top part of Figure 5  The bottom part of Figure 5 shows the difference between the hourly capacity factors using typical days and the reference values calculated without clustering.We find that the mean absolute error (MAE) strongly decreases from 38.6 Â 10 À4 using 4 typical days to 2.9 Â 10 À4 using 52 typical days.For this reason, a low number of typical days are not sufficient to accurately determine the timeresolved yield of a PV installation.However, as can be seen in Figure 5, we find that the errors between the values from the clustering approach and the reference values are mostly symmetrical around the cluster centers.Thus, most errors cancel out over the course of 1 year when calculating the annual yield of a PV installation.Hence, we underestimate the full load hours by 0.47% using 4 typical days and by less than 0.04% using 12, 24, or 52 typical days.We validate F I G U R E 5 (Top) Averaged capacity factor for south facing facades at 30 N latitude and 0 longitude in the high-rise building Environment C calculated without clustering of Sun positions with an hourly resolution.(Bottom) Absolute error between the capacity factor shown in the top part of the figure and the capacity factor calculated using clustering with 12 (red line), 24 (green line), and 52 (blue line) typical days.The mean absolute error (MAE) decreases from 12.2 Â 10 À4 using 12 typical days to 2.9 Â 10 À4 using 52 typical days.
our clustering approach for latitudes between 60 S and 90 N using the hourly resolved reference calculations and find that the described behavior of the clustering approach is consistent for all investigated degrees of latitudes.Based on these results, we decided to use 12 typical days for further calculations in this study.The ray tracing of a whole year with hourly resolution without clustering for Environment C (104,611 triangles) takes about 2,070 min on our workstation.In contrast, using 12 typical days, the same calculation only takes about 70 min, thus ≈30 times faster.

| Worldwide full load hours of facade PV installations
Using the clustering approach, we calculate worldwide full load hours of non-shaded facades with different orientations.A non-shaded facade is considered the best case and serves as a reference for the calculations in more complex Environments A, B and C. Figure 6 shows the worldwide full load hours for a south facing facade.In addition to the influence of latitude on the FLH of the facades, the influence of different weather conditions is clearly visible.For example, the number of full load hours of south facades varies by about 36% between the west coast and the east coast of the United States albeit being on the same degree of latitude.This emphasizes the strength of our approach, which clusters only the Sun positions and not the associated weather data, so that the full spatial and temporal resolution of the weather data is preserved.
The full load hours shown in Figure 6 result from the insolation values calculated by our ray tracing approach that considers up to three light reflections in the environment.Figure 7 shows the relative increase in full load hours for a non-shaded south facing facade when up to three reflections are considered instead of none.As can be seen in Figure 7, light reflections in the environment significantly increase the full load hours of facade surfaces.We find that the mean increase of full load hours along the degrees of latitude ranges from about 10% to about 25% for south facing facades in the non-shaded case as well as in Environments A, B, and C. In the case of east facing facades, the mean increase ranges from about 10% to about 20%.This underlines the advantage of ray tracing calculations with multiple light reflections compared to transposition models with simplified light reflection models only.

| Seasonal distribution of monthly full load hours
Please note that the following discussions of the results focus on the northern hemisphere for the sake of clarity.Figure 8

| Intra-day distribution of facade PV power generation
To match the PV power generation in a specific area to the shape of power consumption, information about the seasonal distribution of full load hours and information about the intra-day distribution of PV power generation is necessary.For the calculation of diffuse insolation, we construct a sky hemisphere around the simulation scene as shown in Figure A1.Since we assume an isotropic sky hemisphere, the starting points of rays on the hemisphere are uniformly sampled.The direction for each ray is sampled according to Lambert's cosine law, resulting in an ideal diffuse sky.The ray tracing algorithm yields hit counts for each surface in the simulation scene, which we subsequently normalize to the diffuse horizontal irradiance (DHI).For this normalization, we divide the density of hits on each surface by the density of hits on the horizontal plane if there were no objects that cause shading.Thus, we have to calculate the latter density in dependence of the number of rays N diff.raysand the radius r sphere of the sky hemisphere.The radiant flux Φ 12 that the surface dA 2 exchanges with the surface dA 1 is sphere dA 1 dA 2 , with β 1 the angle between the normal of dA 1 and the connecting line to dA 2 , β 2 the angle between the normal of dA 2 and the connecting line to dA 1 , and L 2 the radiance that is emitted from dA 2 .The ray tracing results yield hits per area; thus, Since we assume the sky as ideal diffuse, we know the relationship between the emitted radiance L 2 , the radiant flux dΦ 2 , and the surface element dA 2 : and thus, dΦ 12 dA 1 ¼ dΦ 2 dA 2 : Hence, the density of radiant flux the surface element dA 1 receives from dA 2 equals the density of radiant flux emitted on the sky hemisphere.In total, we create N diff.rays on the sky hemisphere with a surface area A sphere ¼ 2πr 2 sphere ; thus, we normalize the hits per area on each surface with presented a hybrid approach between bottom-up and top-down for a highresolution global spatiotemporal assessment of rooftop solar PV based on big data, machine learning, and geospatial analysis.Their analysis assumes roof-mounted PV installations without shading and with optimum tilt angle.The result is a global map of the technical PV rooftop potential in terms of the annual sum of generated electricity and associated LCOE values.
Figure1, contains PV-modules with tilts ϑ ranging from 0 (horizontal) to 90 (vertical).The ground plane extents virtually infinite in all directions with an albedo α.In case of the isotropic sky model, the diffuse insolation is directly proportional to the sky view factor (SVF) of each PV-module.The SVF of an unshaded plane with tilt ϑ is given by39

2 . 4 |
Figure 3A) and shift the results either to the west or east (Exemplary Point B in Figure 3A).We interpolate the shifted values to the values at the full hour using a quadratic interpolation.The example depicted in Figure 3A,B shows a shift of 105 to the west, which corresponds to a time shift of 7 h.Since the relative irradiation values are and pay special attention to facade PV potentials.Additionally, we compare our results with the case of F I G U R E 3 (A) The relative direct insolation and the relative diffuse insolation can be shifted along a latitude.The above example shows a shift of results from A at 0 longitude to B at À105 longitude.(B) Relative direct insolation for one exemplary facade shifted from A to B. The shift of results along a latitude equals a time shift of the relative insolation values calculated by our ray tracing approach.The shift from A at 0 longitude to B at À105 longitude equals a time shift of 7 h.F I G U R E 4 Three urban environments we used in this study for the ray tracing of populated areas with different building densities and geometries.The left image shows a mixed district in the city of Berlin (Germany) with single-family houses but also some high-density apartment buildings (Environment A).The center image shows a densely populated district in the city of Hannover (Germany) with largely apartment buildings (Environment B).The right image shows a densely built-up area with high-rise buildings (Environment C) in the city of New York City (United States).All excerpts have an edge length of approximately 1 km.
exemplarily shows the hourly resolved area-weighted average capacity factor f of south facades in densely built-up Environment C at 30 N latitude and 0 longitude.The capacity factor shows a maximum in winter with significantly lower values in summer.This is due to the 90 tilt of the facades and the lower Sun elevation in winter.At the exemplarily chosen geographic location, the Sun reaches a maximum elevation angle of ≈83.4 in summer compared to a maximum elevation of ≈36.6 in winter.Thus, the angle between the facade normal vector and the Sun direction is considerably smaller during winter yielding a higher capacity factor.Subsequently, we compare the hourly capacity factors without clustering of Sun positions with a reproduction of the whole year using 4, 12, 24, and 52 typical days.
shows the seasonal distribution of the monthly full load hours of east and south facing facades (left) and roofs (right) from the equator to the north pole along 10 E of longitude.The results shown are taken from environment A. In general, the annual full load hours of facade installations are lower compared to rooftop installations.This difference between facade and rooftop installations is most pronounced between 5 N and 40 N latitude.Moving further towards the north, the differences in annual full load hours between facade and rooftop installations become smaller.Further, we find pronounced differences in the seasonal distribution of full load hours between facade and rooftop installations depending on their orientation and location.East facing facades show a uniform PV power generation throughout the year from 10 N to 30 N latitude.For the same region, south facing facades show high full load hours during winter but very small PV power generation in summer.This is the reason that east and west facing facades show higher annual full load hours compared to south facing facades up to 22 N latitude.In this region, south facing rooftop installations show a similar seasonal distribution of the full load hours; however, the decrease of PV power generation during summer is much less pronounced.In the case of east facing rooftop installations, we find an increase in PV power generation during summer, making the seasonal distribution less uniform when compared to east facing facades.For locations north of 40 N, the PV power generation of both facade and rooftop installations starts to concentrate during summer.Importantly, we find that south facing facade and rooftop installations show a considerably broader seasonal distribution compared to east facing facade and rooftop installation in this region.This effect is most pronounced when comparing the south facing facade with the east facing facade.Moving further to the north, the full load hours of both facade and rooftop installations during both winter and summer decrease.Interestingly, the full load hours during summer start to increase north of 80 N. The reason is the arctic summer during which the Sun never sets.

Figure 9 Figure 9 A. 1 |
Figure 9 are because the facades in the respective environment are predominantly facing in a slightly more north or south direction.

ρ
diff:unshaded ¼ N diff:rays 2πr 2: sphere F I G U R E A 1 Illustration of the sky hemisphere used in our ray tracing approach for the simulation of an ideal diffuse sky.The surface element dA 2 exchanges radiation with the surface element dA 1 .