Validation of a 3D Local-Scale Adaptive Solar Radiation Model by Using Pyranometer Measurements and High-Resolution DEM

1. Introduction

Knowledge of local solar radiation is essential for designing optimal solar energy systems. Ideally, this knowledge should come from long-term measured data. However, the limited coverage of radiation measurement networks imposes the need to develop solar radiation models. Therefore, validation of these models against measurements is crucial [1].

Dealing with solar radiation and shadows is a common issue in Geographical Information Systems (GIS) that stems both from a wide range of interests that go from visualization and rendering to many potential applications such as predicting the radiation at crucial environmental like glaciers [2], the design of agroforestry exploitations [3] or, in the urban environment, project planning [4], the generation of solar radiation maps [5], which are essential for the deployment of concentrating solar power (CSP) plants [6], and for planning the optimal location and orientation of PV or solar thermal panels for distributed generation of electricity or heat for domestic applications at building roofs or facades [7,8,9,10,11,12,13,14,15,16]. The latter application demands reliable models for solar radiation and shadow casting at the local level that can deal with the complexities of urban architecture, which has a tremendous impact on blocking the solar radiation of nearby areas. Furthermore, this involves a high accuracy of the input so that it can guarantee to capture the mentioned architectural features. On top of all this, a contained computation cost is highly desirable.

The Sun radiation in the direct, diffuse and reflected forms, as well as the treatment of shadow casting is an active matter of research for which several methodologies have been proposed. Li and Liu [10] employed the Yallop algorithm to calculate the solar position and evaluate several methods for calculating diffuse radiation, but did not adress the issue of shadows cast by nearby obstacles. In [15], Toledo et al. report on the determination of the optimal PV panel tilting angle accounting for direct and diffuse radiation calculated by using the Erbs’s correlation and the effect of shadows cast by surrounding buildings parameterized with their height and distance. In [7,8], Dorman et al. proposed an open-source tool devised to modeling and shadow casting of a 2.5D vector-based environment. The shadow calculation is based on trigonometric relations between the sun’s rays, the ground, and any intervening obstacles. In [16] Tripathy et al. report on a system for computing the shadows as the intersection of parallelograms that arise from the transformation of a polygon-defined 3D city into a cloud of points with the terrain, facades and horizontal surfaces. Stendardo and coworkers [12] employed separate digital surface/elevation models with different accuracy levels for the terrain and for a 3D city in a system powered by ray-tracing algorithms implemented through parallelized GPU computations in order to estimate the solar potential in Geneva. Liang et al. [13] proposes a GIS-integrated open source tool that includes the r.sun package for the estimation of clear-sky, diffuse, and reflected solar radiation that is in turn coupled to a well-established OpenGL-based 3D-graphics engine to compute shadow casting. In [14] Toledo and coworkers employed the SOL algorithm coupled to the direct and diffuse solar radiation model from Perez [17]. The algorithm, implemented in MATLAB, employs the concepts of hyperpoints to calculate the height of shadows in facades in a rectangular grid that represents the 3D city. In [4] Kaynak et al. obtain the direct solar radiation the elevation angle constant (EAC) method, and use the local meteorological with the Anstrom-Prescott model to estimate the sunlight attenuation due to atmospheric conditions, while the finite element method, back-face detection and ray-tracing algorithms are used to make a precise calculation of cast shadows.

Regarding the accuracy level of the Digital Elevation Model (DEM), several works suggest that a finely detailed 3D representation of the geometric features is vital to accurately evaluate the shadows cast by nearby geometries. In [18] the authors used various GIS light projection models to analyse how the level of detail in a building’s 3D approximation affects the quality of its projected shadows. They found that finer details are relevant for local applications such as predicting solar panels’ potential. In [19] the solar yield was assessed based on the grid spacing of the point cloud, finding out that the error increases linearly with the grid spacing up to 4 m, that is generally larger for structured grids than for unstructured grids, and that considerably fine grids are necessary when considering roofs with complex shadowing artifacts.

In the context of GIS and the representation of architectural features in the urban environment, a common approach is the use of the so-called 2.5D polygons, which consist mainly of simplified representation of buildings as a composition of extruded shapes. The previous studies suggest that the polygon-based 2.5D urban mapping strategies might not be able to hold the needed detail when buildings that feature more complex shapes come into play. However, some works have extended the 2.5D geometry concept through the adaptation of more complex architectural shapes by using LiDAR data to reproduce the building geometry by means of a composition of arbitrarily and vertically (walls) oriented polygons within an adaptive mesh [20]. The present work uses a similar approach; we generate a high-definition DEM that accurately captures the terrain and the architectural geometry. This DEM is created from LiDAR data with a scan spacing of 0.30 m. Computing the shadows with this high-definition DEM is too computationally expensive, so instead, we use a triangular mesh adapted to it.

In this paper, we introduce a toolbox that accurately calculates solar radiation while considering the shadow casting of complex orography and architectural features. Section 2 describes this toolbox: the MAPSol model, the construction of the high-definition DEM, and the generation of the adapted triangular mesh. In Section 3, we explain the process for acquiring experimental data and how we validated our model. We use a case study where shadows originate from the intricate geometry of the tower of the Cathedral of Salamanca onto a nearby placed pyranometer. The study’s results confirm our model’s accuracy in accounting for the shadows that impact solar radiation.

2. Materials and Methods

This multidisciplinary work has required a broad methodology covering different areas of knowledge. It begins by recalling some basic concepts about solar irradiance that will allow a better understanding of the work. Then, the MAPSol 3D local-scale adaptive solar irradiance model for which the validation protocol described in this paper has been designed is presented. This model allows the detection of the shadows cast by complex orography, but this requires a high-resolution map of the study surface. Therefore, the methodology used to obtain high-resolution DEM is described. To detect the shadows, the model relies on a mesh of the surface defined by the high-resolution DEM, so the next step is to generate an adapted mesh in the most efficient way. Finally, the model is validated by means of actual solar irradiance measurements. For this reason, the corresponding procedure of experimental data acquisition is described in detail.

2.1. Solar irradiance

Irradiance is defined as the incident energy per unit of time and surface area, i.e., it is a power density [21]. The solar irradiance received on a point of interest depends on three key factors. The first factor is the Earth’s position concerning the Sun, determined by latitude, longitude, day, and time. This factor determines the angle at which the sunbeam strikes a given spot. The second factor is the orography and buildings, which impact shadows and the orientation between the sunbeam and the point of interest. The third factor is the atmosphere, which attenuates solar irradiance due to gases, solid and liquid particles, and clouds.

Depending on the factors considered, we can obtain two types of solar irradiance: real-sky irradiance (when all factors are considered) and clear-sky irradiance (when clouds are excluded) [22]. On the other hand, due to the scattering effects generated by the atmosphere and the effect of the angle of incidence, the global irradiance can be expressed as the sum of two components: the direct irradiance coming from the solar disk, and the diffuse irradiance, due to the scattering effects in the atmosphere.

The definition of the main irradiances is presented below [23].

• Diffuse Horizontal Irradiance, $DHI$: This is the solar irradiance collected on a horizontal surface from atmospheric scattering of light, excluding circumsolar radiation.
• Direct Normal Irradiance, $DNI$: It is the component of solar irradiance collected on a surface perpendicular to the Sun’s rays. The horizontal diffuse component, $DHI$, is neglected here. On clear days, this component is much larger than the diffuse component, while on days with high cloud cover, it is practically zero. As it is measured over the Earth’s surface, its values depend highly on atmospheric conditions and the time of the year.
• Global Horizontal Irradiance, $GHI$: This is the sum of all irradiance components collected over a horizontal surface. This includes the direct and diffuse components, as well as the reflected components, which are generally neglected because of their low value. The $GHI$ can be calculated from the following expression:

  $$GHI=DHI+DNI·sin\alpha$$

  (1)

  where $\alpha$ is the solar altitude angle, i.e, the complementary of the zenith angle of the Sun, $DHI$ is the horizontal diffuse component and $DNI$ is the normal component.
• Beam Horizontal Irradiance, $BHI$: It is the direct horizontal component of the irradiance, i.e., the direct irradiance on a plane perpendicular to the vertical of the site. It can be obtained as:

  $$BHI=GHI-DHI$$

  (2)

2.2. The MAPSol model

This section describes the solar irradiance model MAPSol developed at the University of Las Palmas de Gran Canaria [24,25,26]. MAPSol computes the clear-sky beam irradiance in every triangle of an adaptative mesh representing the orography surface. What sets MAPSol apart from other models is that it considers the orography shadows in the irradiance computation.

2.2.1. Clear-sky beam irradiance model

In this section, the main ingredients of the clear-sky beam irradiance model implemented in MAPSol are presented; a more detailed description can be found in [22,24].

The BHI equation for a horizontal surface is:

$$BHI=\left(I_0ϵ\right)\left(-0.8662T_{LK}m{\delta }_R\left(m\right)\right)\left(L_{f}sin\alpha \right).$$

(3)

In this equation, each of the terms can be related to the three factors impacting solar irradiance.

The first term $\left(I_0ϵ\right)$ refers to the Earth’s geometry. $I_0$ is the solar constant (the Sunbeam irradiance at the mean solar distance) [27], and $ϵ$ is a correction factor due to earth’s orbit eccentricity [22,24].

The second term $\left(-0.8662T_{LK}m{\delta }_R\left(m\right)\right)$ concerns to the atmospheric attenuation. Specifically, m and ${\delta }_R\left(m\right)$ are related to the attenuation by gas constituents, and $T_{LK}$ is related to the attenuation by solid and liquid particles [22]. The parameter m is the relative optical air mass [28]; $T_{LK}$ is the air mass 2 Linke atmospheric turbidity factor, and $\delta R\left(m\right)$ is the Rayleigh optical thickness at air mass m and both are calculated according to the improved formula by [29].

Finally, the third term $\left(L_{f}sin\alpha \right)$ is related to the shadows of both the orography and the buildings. Parameter $\alpha$ represents the solar altitude angle, while $L_f$ is the lighting factor that considers cast shadows. A Solar Position Algorithm computes the solar altitude angle; the present model uses the one developed at the Plataforma Solar de Almería: the PSA algorithm [30]. On the other hand, the light factor takes values between 0 and 1; $L_f=0$ means complete shadowing and $L_f=1$ total sun exposure. The following section describes the algorithm to compute this factor.

2.2.2. Shadow detection

The model computes the light factor using a triangular mesh of the terrain generated with the method shown in Section. Section 2.4. Therefore, $L_f$ relates directly to finding which triangles in the mesh are shadowed. This algorithm is described in detail in [24,26,31].

In short, this geometrical problem can be solved by finding those triangles that, looking at the mesh from the Sun, have at least another triangle that covers the former. This is done using a new reference system $x^{\prime }$, $y^{\prime }$ and $z^{\prime }$, where $z^{\prime }$ is in the direction of the beam irradiance (see Figure 1). Then, the mesh is projected on the plane $x^{\prime }{y}^{\prime }$. This is the view of the mesh from the Sun that allows finding out which triangles are overlapped, if any.

To transform the coordinates and align the direction of the $z^{\prime }$ axis with the beam irradiance, we need the solar azimuth $\theta$ and solar altitude $\alpha$ (shown in Figure 1). These angles are obtained from the PSA algorithm [30].

Examining intersected triangles at each time step and applying temporal refinement to each analyzed triangle, a shading level is allocated to shadowed triangles based on the count of warning points located within other triangles (refer to Figure 2). This allocation is the light factor $L_f$.

The drawback of this technique emerges when attempting its application to a sizable mesh, as the computational cost becomes high, making the process too slow for some purposes.

In response to this challenge, a technique was devised in [31]. This approach involves pre-filtering candidate triangles before examining the ultimate mesh overlap. The determination of whether a triangle is shadowed is contingent upon either self-shadows or shadows cast by other triangles.

Self-shadows occur when a surface casts shadows onto itself, for instance, when it is oriented away from the Sun. Determining this condition is straightforward by examining the incidence angle (${\delta }_{exp}$), defined as the angle between the solar vector and the normal to a triangle surface. The criterion for self-shadows is met when the incidence angle exceeds $\pi /2$. So, a Self Shadows Light Factor ($L_{fss}$) is defined with values:

$$L_{fss}=\left\{\begin{array}{cc}0\hfill & |{\delta }_{exp}|>(\pi /2)\hfill \\ 1\hfill & (\pi /2)\ge |{\delta }_{exp}|\ge 0\hfill \end{array}\right.$$

As the analysis of cast shadows demands substantial computational effort, an initial step involves circumventing the examination of all triangles oriented away from the Sun (self-shadowed). An effective strategy involves filtering triangles based on:

• In the absence of self-shadowed triangles (those facing away from the Sun), the entire mesh is illuminated, and no shadows are present.
• Only triangles oriented away from the Sun are capable of casting shadows. These are referred to as potential 1 triangles [31].

The shading analysis is finalised by examining the shadow projection on a specific set of points known as warning points (WP). These points are uniformly distributed within the triangle using the Rivara 4-T algorithm [32] and correspond to the geometric centres resulting from each subdivision of the triangle (see Figure 2).

A warning point is shaded when its projection onto the x’y’ plane falls within a triangle that can potentially cast a shadow on it and its z’ coordinate is smaller than that of its projection onto the said triangle (it is further from the sun than its projection on the triangle). So, the cast shadows light factor, $L_{fcs}$ is computed as:

$$L_{fcs}=\frac{n_{wp}-i}{n_{wp}}$$

(4)

where $i=0,1,\cdots ,n_{wp}$ is the number of shadowed warning points, and $n_{wp}$ the total number of warning points. The ultimate light factor, $L_f$ is:

$$L_f=L_{fcs}·L_{fss}$$

(5)

2.3. High-resolution DEM

There is a wide variety of technologies that allow a massive capture of terrestrial information [33]. Among the different geospatial data coming from different sensors that allow the massive capture of information, LiDAR data stand out [34]. LiDAR data, obtained from laser scanners on board different devices, are highly accurate three-dimensional data. Aerial LiDAR data are the most widely used data for massive data acquisition in urban environments [35]. Thanks to the INSPIRE (Infrastructure for Spatial Information in Europe) initiative [36] based on the Directive 2007/2 [37], many countries within the European Union provide LiDAR data to users. The National Geographic Institute (IGN) in Spain offers highly accurate and updated LiDAR data of the entire Spanish territory, captured periodically under the National Plan for Territory Observation (Plan Nacional de Observación del Territorio - PNOT) [38]. Previous studies have demonstrated the applicability of these LiDAR data in urban environments in Spain [39,40,41,42].

Although the IGN allows the user to download the LiDAR point cloud through the Download Centre [43], it also offers several three-dimensional products that are derived from LiDAR data: (i) Digital Terrain Model (DTM) (ii) Digital Elevation Model (DEM) (iii) Digital Slope Model (DSM). All of them have one feature in common: a spatial resolution between 2 and 25 metres. Although a spatial resolution of metres may be sufficient for some studies, in certain cases, a resolution of less than one meter is necessary to obtain accurate results. For the purpose of this research study, a high-resolution DEM was generated using the available LiDAR data from IGN over the study area, described in Section 3.2. The LiDAR data used belongs to the first coverage of the PNOT observation program and has the characteristics described in Table 1.

From this data, a high-resolution DEM is generated by applying the grid method of triangulation named as TIN (Triangular Irregular Network), which converts the network of individual points in a continuous surface of triangular facets. The vertices of the triangles are the points measured by the LiDAR sensor, forming the point cloud. The procedure behind the generation of the DEM is the Delaunay triangulation, where the length of the edges and the angles for the triangles can be limited in order to generate smoother surfaces with fewer peaks, avoiding sliver triangles. Thus, according to the Delaunay triangulation, triangles are generated for groups of three points, ensuring that there is no point inside the circum-sphere of any triangle. Once the triangulation is generated, the triangular faces are interpolated, using the values of elevation and slope in order to generate a smoother elevation layer.

2.4. Mesh generation

The MAPSol model described in Section 2.2 requires a discretisation of the orography provided by the high-resolution DEM described in Section 2.3. The uniform triangulation associated with the DEM-raster should be discarded, as it is not computationally feasible. As an alternative, the construction of an adapted mesh is proposed, which allows capturing the singularities of the study domain and computing the beam irradiance in a more efficient way.

The mesh generation is based on a local refinement driven by the error (vertical distance) between the approximation of the orography (adapted mesh) and the actual orography (DEM). The algorithm starts with a coarse mesh and locally refines [32] those triangles that do not satisfy the error criterion. It produces graded meshes that capture the features of the domain.

The process is fully automatic and only requires an initial discretisation, in this case the DEM, to capture the orography with a fixed accuracy. It has been implemented with the Neptuno library developed by L. Ferragut [44]. The entire mesh generation process, including the DEM construction described in Section 2.3, is summarized in Figure 3.

Although MAPSol does not require a high-quality mesh, the previous procedure could be combined with surface parameterization and/or smoothing techniques [45,46] to improve the mesh quality and reduce the number of triangles of the final mesh.

2.5. Experimental measurements of solar irradiance with pyranometers

The purpose of this section is to obtain a realistic picture of solar radiation in the area of the Trilingüe building of the University of Salamanca (Spain) through experimental data. The physical principles of the pyranometers used to collect solar irradiance data are presented below.

The main function of pyranometers is to convert a flux of incident photons into an electrical signal through a process of collection, selection and detection [47]. Depending on the magnitude to be measured, they require different adaptations. Here, the focus is on devices for measuring $GHI$ and $DHI$. In the case of GHI, the devices for measuring global radiation register very unequal fluxes, and they have also a strong dependence on the Sun position. In order for the electrical signal obtained to be representative of the incident flux, they have to mix the radiation inside, so that the resulting signal is isotropic. This can be achieved by using a diffuser sheet or an integrating sphere as receiving elements.

On the other hand, pyranometers dedicated to measuring DHI require a complement that eliminates the direct component of the radiation. A simple system widely used incorporates a shadow ring, which is nothing more than a ring of an opaque material that casts a direct shadow on the receiver. The ring is positioned so that its axis is parallel to the Earth’s axis in the equatorial coordinate system and passes through the receiver. Figure 4 schematizes the positioning of a shadow ring at a latitude of ${40}^{\circ }$ N, where the arrows on the right side represent the Sun’s rays. Once the ring has been correctly positioned, the receiver will be permanently shaded provided that the ratio of the receiver’s distance from the center of the ring to its radius is close to the tangent of the Sun’s declination. Due to the changing declination of the Sun over the course of the year, the shadow ring requires small adjustments in its position every two or three days, aiming to keep the receiver permanently shaded.

To compensate for the defect in the diffuse irradiance collection generated by the shadow ring, a correction factor is introduced. It depends on the latitude of the observation point as well as the time of the year. In reference [48], the ring manufacturer provides the full derivation of the correction factor, C:

$$C={\displaystyle \frac{1}{1-S}}$$

(6)

where S is obtained as:

$$S={\displaystyle \frac{2·0.185rad}{\Pi }}cos\delta (U_0·sin\lambda ·sin\delta +sin{U}_0·cos\lambda ·cos\delta )$$

(7)

where $\lambda$ is the angle of latitude, $\delta$ the angle of declination of the Sun, and $U_0$ the angle formed by the Sun’s position at sunrise and noon in the plane of the ring,

$$U_0=arccos(-tan\lambda tan\delta )$$

(8)

There is a wide variety of pyranometers that implement different technologies depending on the type of transducer they use. They are subdivided into two classes: indirect conversion detectors and direct conversion detectors [47].

• Indirect conversion detectors: they work by converting the incident photon flux into another type of flux (usually heat), but can also be a secondary photon flux. Heat flux detectors are widely used and their operation is relatively simple. To convert the photon flux into heat flux, a highly absorbing paint or varnish is applied to the detector, which causes its temperature to rise when the light beam is impinging on it. Knowing the temperature at two points, and assuming that the steady state is reached, the intensity of the flux is calculated, which will be proportional to the temperature difference. Figure 5 (a) shows a general scheme of the parts of an indirect heat flux conversion pyranometer. In the upper part there are two domes, the outer dome has the function of avoiding energy exchanges due to convective phenomena; as a whole, the domes act as an integrating sphere. As it can be seen, the detector is surrounded by an anti-radiation shield to prevent radiation penetrating from anywhere other than the dome. Figure 5 (b) shows the Pyranometer Kipp & Zonen SMP10, belonging to the Energy Optimization, Thermodynamics and Statistical Physics Group (GTFE), with which the Global Horizontal Irradiance measurements were performed.

  

  Figure 5. Heat flux sensing pyranometer: (a) Basic scheme [49] and (b) Pyranometer belonging to the Group of Energy Optimization, Thermodynamics and Statistical Physics (GTFE) of the University of Salamanca.

  Figure 5. Heat flux sensing pyranometer: (a) Basic scheme [49] and (b) Pyranometer belonging to the Group of Energy Optimization, Thermodynamics and Statistical Physics (GTFE) of the University of Salamanca.

  
• Direct conversion detectors: again there are two types. Photoemitter cells are based on the junction of an anode and a cathode, between which there is a large potential difference (in the range of kV) and an avalanche effect is produced. On the other hand, there are detectors based on PN junctions, the photodiodes, where the current generated is proportional to the incident flux. These types of detectors have better sensitivity than avalanche detectors and work with low voltage [50].

For this work, radiation data obtained in situ are used. These data are recorded with two Kipp & Zonen pyranometers, model SMP10. These are class A pyranometers (according to ISO 9060) that implement thermopiles as heat flux detectors. Thus, they work as indirect conversion detectors. With them, different magnitudes are collected: pyranometer 1 collects $GHI$ measurements, while pyranometer 2, equipped with a Kipp & Zonen CM 121 shadow ring, takes $DHI$ values. Table 2 shows the main characteristics of the model in question (detailed information about the pyranometers or the shadow ring can be found in the complete manufacturer’s manual [49]).

Figure 6 shows a comparison between the wavelength response of an SMP series pyranometer and the standard solar spectrum collected at sea level. As it can be seen, the pyranometer response is practically linear over the entire range of interest, from the infrared (around 2000 nm) to the ultraviolet limit of the visible (around 350 nm). In this region of the spectrum, most of the energy from the Sun is deposited.

3. Results

In this section, it is precisely described the case study with which the MAPSol model has been validated, including the description and reliability of the chosen irradiance data, the generation of the high-resolution DEM of the study area, as well as the adapted grid, and the simulation with the MAPSol model. The chosen case study provides the possibility to assess the effect of shadows caused by the intricate orography generated by buildings close to the location of pyranometers.

3.1. Experimental data acquisition

Experimental data were obtained from two pyranometers and a weather station, shown in Figure 7, installed on the rooftop of the Trilingüe building at the Faculty of Sciences (40.96062 N, 5.67075 W) of the University of Salamanca (Spain): A pyranometer Kipp & Zonen SMP10 (Figure 7 (a)), and a pyranometer SMP10 with shadow ring model Kipp & Zonen CM 121 (Figure 7 (b)).

During the day, the devices perform the corresponding measurements, storing $GHI$, $DHI$ values every five minutes. At the end of the day, the collected data are sent to the computer, generating a csv file containing the irradiance information. The next step is to classify and filter the files. To do this, a program developed in Mathematica^® identifies the files that belong to the same day and extracts the desired information from each of them. Finally, the program generates a new output file where all the information for the same day is merged.

The results presented in this work have been obtained with data between 07/30/2020 and 02/22/2023. In Figure 8, illustrative monthly averages have been performed on the irradiance, so that for each month $GHI$ and $BHI$ curves are obtained as a function of time (which is represented in UTC format). The monthly average is obtained by considering the data of all the days belonging to that month in the whole study period.

As previously mentioned, pyranometers take measurements of $GHI$ and $DHI$. To obtain the $DHI$ component, a correction factor which depends on the day of the year is required. For each of the measurements the correction factor is calculated as a function of the solar declination angle, hereafter $DHI$ is considered to include the correction factor. The $BHI$ and $DHI$ measurements are obtained, every 5 minutes, from equations (2) and (1), respectively. During the night, $GHI$ and $DHI$ values should be identically zero, but it is found that negative values appear. This is usual when performing this type of measurements. However, it does not make physical sense to work with negative irradiances. To solve this problem, the following filters are introduced into the program:

• If $DH{I}_i<0\Rightarrow DH{I}_i=0$
• If $GH{I}_i<0\Rightarrow GH{I}_i=0$
• If $(GH{I}_i-DH{I}_i)<0\Rightarrow (GH{I}_i-DH{I}_i)=0$
• If $sin{\alpha }_i<0\Rightarrow DN{I}_i=0$

The irradiance curves introduced below are intended to be as close as possible to the actual measurements, so no statistical treatment beyond averaging is done. Figure 8 shows the irradiance curves for four representative months of the year: March, June, September and December. In all these plots, a certain level of noise is observed within the irradiances, which evinces that the pyranometers are very sensitive to small changes. For example, on a summer day with isolated clouds, if a cloud comes between the Sun and the pyranometer for a few minutes, the direct components will become very small. When the cloud has passed, the direct component will again reach high values. When observing the set of figures, it can be noticed that the irradiance curves increase in width and height in the summer months. On the contrary, during winter, when the Sun is less time and at a lower height above the horizon, the curves become narrower and present lower peak values.

In order to assess the reliability of the measurements obtained, these values are compared with two independent external sources: AEMET [51] and Solargis [52] for solar irradiance. The Atlas of Solar Radiation in Spain, using data from the EUMETSAT climate FAS developed by AEMET [51], uses satellite data with a spatial resolution of $3\times 3$ km. Data from the period (1983 – 2005) were used to elaborate the atlas. In order to compare the results in a more convenient way, Table 3 shows the results of this work compared with those of AEMET. The relative differences obtained for 9 of the 12 months of the year, do not overcome 10 % in the daily accumulated energy for the global component ($GHI$). This number is reduced to 8 for the direct horizontal component ($BHI$).

Taking the experimental data of $GHI$, $BHI$ from Table 3, it is possible to obtain the annual cumulative values and also the average value per day of the energy received. The annual cumulative values are easily obtained by performing the operation:

$$E_x=\sum _{i=1}^{12}{x}_i{N}_i.$$

(9)

where x represents the irradiance component, the subscript i is the number of the month and $N_i$ is the number of days in month i.

The Atlas of Solar Radiation in Spain [51] using data from the EUMETSAT climate FAS provides daily average $GHI$ and $BHI$ values in the form of an irradiance distribution map. The values extracted from these maps can be seen in Table 4, together with the values recorded in situ.

The value for $GHI$ obtained in this work is within the AEMET interval, touching the upper limit. The experimental $BHI$ component is 0.01 kWh ${\mathrm{m}}^{-2}·{\mathrm{day}}^{-1}$ above the maximum value of the interval. AEMET does not offer the annual average accumulated values, but they can be calculated by multiplying the daily average by the number of days in the year (the result is also shown in the same table). Again, the experimental value of $GHI$ is within the range, while the value of $BHI$ is slightly above it.

As it is known, Sun chart is employed to present, at a specific location, the apparent position of the Sun, i.e., the height of the Sun at any hour of the day. From a Sun chart elaborated by the University of Oregon Program [53], together with a panoramic photo from the Trilingüe building at the University of Salamanca, Figure 9 is elaborated. This figure shows the Sun path chart from the Trilingüe building at the University of Salamanca (40.96062 N, 5.670759 O) between December 21^th and June 21^th, overlapped with a panoramic photo, taken from pyranometer that register $GHI$, to identify shadowing sources.

One of these shadowing sources on the location of the pyranometers is the 93-metre-high cathedral tower, located 296 metres to the east. As can be seen in Figure 9, the cathedral tower casts its shadow in this area during the early hours of some days of the year, specifically between 22 August and 11 September, and in its symmetrical months, February and March. In these months close to spring, cloudiness is higher, so it is more difficult to obtain irradiance measurements that allow the effect of the cathedral’s shadow to be clearly seen. However, in summer the sky is clearer and this phenomenon can be seen more clearly in the measurements obtained.

3.2. Area study, high-resolution DEM and adapted mesh

An area of $430m\times 176m$ has been selected which includes the Trilingüe building at the Faculty of Sciences of the University of Salamanca, where the pyranometers are located, and the cathedral tower located at $296m$ in a straight line to the east, which casts its shadow on the pyranometers affecting the reading of the irradiance data.

Provided the point cloud available in the area under study, the DEM generated with the process described in Section 2.3 presents a resolution of $33cm$ (Figure 10). This methodology performs the generation of the Delaunay triangulation using negligible computing time (3 seconds for $35.228$ points and $1.385.800$ triangles triangles), with a computer with Intel Core $i7-6700$ processor at $3.41GHz$, $64bits$, $32Gb.$ RAM and $931Gb$.

The size of the raster file is too large for the calculations involved in the MAPSol model, so the mesh is adapted using the procedure described in Section 2.4, achieving a reduction of $62\%$ in the number of cells. The adapted mesh for a tolerance of 1 m has $529.503$ triangles, and includes all the singularities of the study area, in particular the cathedral tower. This element holds significant importance in the DEM. It is not present in some of the raster files, such as those provided by the IGN (Digital Surface Model - Building and Digital Surface Model), due to their low resolution and dimensions.

The number of triangles can be further reduced without affecting the accuracy of the shadow and irradiance calculations, using a tolerance of 5 m. This means that in the mesh adaptation process, a triangle is refined if the distance from the triangle to the high-resolution DEM is greater than 5 m, or 1 m in the previous case. The number of triangles of the adapted mesh with a tolerance of 5 m is $44.313$, which represents a reduction of more than $95\%$ with respect to the original mesh, without, as we will show further, a significant loss of precision in the calculation of shadows and irradiance, but with a considerably lower computational cost.

Figure 11c displays the fine adapted mesh (1 metre) for the entire study area, and several zooms of both meshes, the fine adapted mesh (Figure 11b) and the uniform original mesh (Figure 11a), for the Cathedral building area. A 3D reconstruction of the study area is depicted in Figure 12, including the location of the pyranometers with a red point, which allows the complexity of the study area to be observed. Furthermore, this procedure demonstrates that with few free resources, it is possible to reconstruct a 3D image of a very complex area such as this one, in the old part of the city of Salamanca. For Figure 12, the most appropriate orientation has been chosen so that the orthophoto of the area projected on the fine adapted mesh provides the best visual result.

3.3. Simulation with MAPSol

Using the coarse adapted mesh (5 metres), the shadows and global irradiance in each triangle have been calculated for each day of a full year, with a time step of 5 minutes, with the MAPSol model. Computations have been performed in a computer equipped with two AMD EPYC 7313 CPUs, 128GB of RAM memory and Debian Linux version 11 operating system. MAPSol is written in Python and C++, using the latter for the core of the computation and the former mostly to deal with input/output. In this particular work, Python 3.9.16 version and GNU g++ 10.2.1 C++ compiler have been used. Simulations have been performed in parallel, using one core per month. The mean wall clock execution time was 3 hours and 56 minutes, including writing results to files for later analysis.

For each time step, a VTK file representing irradiance in the whole domain is written, as well as a CSV file with values of GHI in the location of the pyranometers. Monthly average values of global irradiance in the study area have been calculated, and in Figure 13 four of the most significant months of the different seasons of the year have been represented. In addition, the annual mean GHI has been depicted in a 3D image in Figure 14, notice that a view from the south has been chosen so that the facades most exposed to solar radiation can be appreciated.

Optionally, additional VTK files representing shadows cast in the domain for each time step can be obtained. Figure 15 represents the computed shadows cast in the study area on 4 September 2022 at 7.00 a.m., when the shadow of the cathedral tower, located in the blue dotted square, affects the reading of the pyranometers, located in the red dot.

3.4. Comparison of simulation results with experimental data

Taking into account the availability of data, due to cloud cover, possible technical difficulties and the reduced number of days on which the effect of the cathedral tower on the location of the pyranometers can be seen, the following dates have been chosen to compare the measured and calculated $GHI$: September 4^th and 11^th 2022, when the measurements will show the effect of the shadow of the cathedral, and March 15^th, 2021 and August 4^th, 2022, when the shadow phenomenon will not occur.

Both measured and calculated $GHI$ values have been taken every 5 minutes throughout the 24 hours of the day, with a total of 288 data. In the two lower graphics of Figure 16, it can be clearly appreciated that the model captures the effect of the shadow of the cathedral in the early morning. The cathedral tower shadows the pyranometers just at sunrise, causing an interruption in the increase of irradiance, either with a sharp and short decrease (Figure 16 bottom left) or with a sharp and delayed sunrise (Figure 16 bottom right). Both phenomenologies are well captured in the simulation, with the pyranometer measurement (purple line) matching the calculated $GHI$ (green line) very accurately. The upper graphics correspond to two dates where the shadow of the cathedral does not affect the irradiance data reading. The upper left graph shows the effect of cloud cover around 9 a.m. on pyranometer measurements (purple line).

To compare the accuracy of the results and the computational cost, the selected days were also simulated using the high resolution (1-metre) mesh. The results at the pyranometer’s locations are the same (with 6-digit accuracy) as those obtained using the 5-metre mesh, but the mean execution time is 68 times higher, so the 5-metre mesh was chosen for all irradiance simulations. However, the use of the high resolution mesh might be appropriate in case a very precise study of shadows over the whole domain is intended.

In order to evaluate the precision of the model for clear-sky, four representative statistical error indicators have been calculated for the selected dates.

• $MAE$: Mean Absolute Error

  $$MAE=\frac{{\sum }_{i=1}^n|{\widehat{GHI}}_i-GH{I}_i|}{n}$$

  (10)
• $NMAE$: Normalized Mean Absolute Error

  $$NMAE=\frac{{\sum }_{i=1}^n|{\widehat{GHI}}_i-GH{I}_i|}{n*GH{I}_{max}}·100$$

  (11)
• $RMSE$: Root Mean Square Error

  $$RMSE=\sqrt{\frac{{\sum }_{i=1}^n{({\widehat{GHI}}_i-GH{I}_i)}^2}{n}}$$

  (12)
• $NRMSE$: Normalized Root Mean Square Error

  $$NRMSE=\sqrt{\frac{{\sum }_{i=1}^n{({\widehat{GHI}}_i-GH{I}_i)}^2}{n*GH{I}_{max}^2}}$$

  (13)
• $R^2$: Coefficient of determination

  $$R^2=1-\frac{{\sum }_{i=1}^n{({\widehat{GHI}}_i-GH{I}_i)}^2}{{\sum }_{i=1}^n{({\widehat{GHI}}_i-\overline{GHI})}^2}$$

  (14)

where $GH{I}_i$ and ${\widehat{GHI}}_i$ represent respectively the measured and calculated value of the global irradiance at time $t_i$, corresponding to the $n=288$ time instants for which data are available. $GH{I}_{max}$ refers to the maximum measured $GHI$, and $\overline{GHI}$ is the mean of the measured GHI. Table 5 summarizes the statistical indicators described above for the four selected days, showing very small values, especially the normalised errors, resulting in a higher coefficient of determination ($R^2\ge 0.99$).

4. Discussion and conclusions

The combination of different technologies from diverse areas of knowledge has allowed the validation of a solar radiation model that is particularly useful in complex orographies such as those with buildings in an urban area. This is especially useful for the optimised design of solar installations that are increasingly being installed in cities.

The solar radiation model MAPSol allows to estimate the clear-sky beam irradiance and generate a very detailed shadow map, which is computed at a reasonable computational cost by pre-filtering the triangles of the surface mesh and using an adapted mesh.

For areas with complicated orography, a high-resolution DEM is required, which is calculated from LiDAR data with a resolution of $33cm$. Adapting the mesh associated with this DEM allows the number of elements to be reduced by more than $90\%$ without reducing the accuracy of the calculations, but at a much lower computational cost. The comparison of results and computational cost using two adapted meshes of different resolutions, one finer (1 metre tolerance) and the other coarser (5 metres tolerance) allows us to conclude that using the coarse mesh considerably reduces the computational cost without losing accuracy.

The accuracy of the MAPSol model has been validated by comparing the calculated $GHI$ data with data measured by two pyranometers. The calculated $GHI$ accuracy is close to $99\%$ on days with low cloud cover. It is worth noting that on some of these days, the sensors detect the shadow of the dome of the cathedral tower in the early hours of the day.

Therefore, the proposed procedure and model allow the calculation of irradiance in areas of very complex orography with a reduced computational cost and freely available data, achieving a great balance between efficiency and effectiveness.

In addition, the adapted mesh generated from the high-resolution DEM has proven to be useful for recreating 3D images of complex areas with limited resources and freely available data.

In the future, the described procedure and model will be used to elaborate a solar map of a larger area of the city of Salamanca. In addition, the comparison of the measured irradiance values, subject to variable cloud cover, and those calculated for clear skies, can be used to produce cloud cover indices, and progress in the use of the MAPSol model under real sky conditions.