Influence of LiDAR Point Cloud Density in the Geometric Characterization of Rooftops for Solar Photovoltaic Studies in Cities

The use of LiDAR (Light Detection and Ranging) data for the definition of the 3D geometry of roofs has been widely exploited in recent years for its posterior application in the field of solar energy. Point density in LiDAR data is an essential characteristic to be taken into account for the accurate estimation of roof geometry: area, orientation and slope. This paper presents a comparative study between LiDAR data of different point densities: 0.5, 1, 2 and 14 points/m2 for the measurement of the area of roofs of residential and industrial buildings. The data used for the study are the LiDAR data freely available by the Spanish Institute of Geography (IGN), which is offered according to the INSPIRE Directive. The results obtained show different behaviors for roofs with an area below and over 200 m2. While the use of low-density point clouds (0.5 point/m2) presents significant errors in the estimation of the area, the use of point clouds with higher density (1 or 2 points/m2) implies a great improvement in the area results, with no significant difference among them. The use of high-density point clouds (14 points/m2) also implies an improvement of the results, although the accuracy does not increase in the same ratio as the increase in density regarding 1 or 2 points/m2. Thus, the conclusion reached is that the geometrical characterization of roofs requires data acquisition with point density of 1 or 2 points/m2, and that higher point densities do not improve the results with the same intensity as they increase computation time.


Introduction
There has been an exponential increase in the global population in recent years [1]. This population growth implies an increase in the needs for transport, electricity, and air conditioning [2]. Another characteristic of the demographic changes is the abandonment of rural areas in favor of cities [3], in such a way that cities become the main points of consumption. According to the United Nations [4], it is estimated that 66% of the world population will live in urban zones in 2050. Regarding cities, buildings are the origin of 40% of the total energy consumption and 36% of CO 2 emissions in Europe [5]. Regulations state that these values must be reduced through the increase in energy efficiency and the incorporation of sustainable energy resources, such as renewable energies [6].
Among renewable energies, solar energy has suffered the highest increase, 28% increase between 2017 and 2020 [7], due to its adequacy for generation at small-scale (self-consumption). This increase is a result of the effort made for the integration of solar energy in cities, which has been accompanied by a high number of scientific publications regarding photovoltaic solar energy in cities [8]. Existing studies

Materials and Methods
In order to evaluate the influence of the density of LiDAR point clouds in estimating the area, data from the IGN have been used. As case studies, different cities in Spain with different point cloud densities have been selected. However, the study is performed in such a way that the location of the case studies is not relevant to the results, so that these can be considered as a guide for those countries where LiDAR data are not available, in such a way that the geometric requirements are clear for the optimal investment in equipment when geometric data are to be provided.

Geospatial Data: Spanish National Geographic Institute
The IGN is a public agency from the Ministry of Transports, Mobility and Urban Agenda (Government of Spain), responsible for measuring and compiling data of the Earth's surface [34], which is made publicly available [35]. The compilation includes updated geographic data in all the territory of Spain, following the INSPIRE European Directive: topographic data in vector and raster format at different scales and thematic geographic data such as land cover/use [36].
The main product offered by IGN is aerial photography. This product is used as the basis for the generation of most of the cartographic products offered. The PNOA Project (which stands for National Plan for Aerial Orthophotography) [37] consists of the periodic acquisition of the territory of Spain with aerial photography. It has been active from the year 2004. PNOA implies the acquisition of images with a temporal resolution of 2-3 years. The main product derived from the aerial images is the aerial orthophotography, which has the main characteristic of spatial resolution: 25 or 50 cm (Table 1). Aerial orthophotographs from PNOA are georeferenced to the Spanish Reference Coordinate System [38]: ETRS89 in the Iberian Peninsula and REGCAN95 for the Canary Islands. The aerial orthophotographs are used in this work as a reference and ground-truth for the evaluation of geometry, due to its high spatial resolution. This work uses orthophotographs from the year 2017. Table 1. Main characteristics of data derived from the National Plan for Aerial Orthophotography (PNOA) Project. Data offered by the Spanish National Geographic Institute on their web page.

Product
Flight GSD 1 (cm) GSD Orthophoto (cm) Planimetric Accuracy in Orthophoto (m) PNOA  In 2008, the PNOA Project included the acquisition of geometric data with LiDAR sensors, with the aim of increasing the accuracy of Digital Terrain Models. This led to the creation of the PNOA-LiDAR Project [39]. LiDAR data are acquired across all territory every 2 or 3 years. Currently, two different LiDAR coverage sets coexist ( Figure 1). The first coverage started in 2008 and finished in 2015. Data acquisition had a minimum point cloud density of 0.5 point/m 2 and elevation accuracy (RMSEz o Root Mean Square Error in the elevation axis (the z-axis)) between 0.20 and 0.40 m ( Table 2). From 2015 up to now, data are being acquired for the second coverage, with a minimum point cloud density of 1 point/m 2 and elevation accuracy (RMSEz) between 0.15 and 0.20 m ( Table 2). PNOA-LiDAR data have the same coordinate reference system as the aerial orthophotographs.  In addition, the points in the point clouds are classified according to the classification from the American Society for Photogrammetry and Remote Sensing (ASPRS) [40], following a methodology presented in [41] that shows that the performance of the classification for points of the class "buildings" is over 90%. For this reason, the classification provided by the dataset is exploited in this study for the analysis of rooftop area.

Methodology
The methodology proposed consists of computing the area of the rooftops in point clouds with different densities, with the aim of evaluating the influence of the density in the accuracy of the area estimated. For this purpose, it is necessary to use aerial orthophotography as reference data. This kind of data is selected due to its planimetric precision. For both cases, the methodology consists of two main stages ( Figure 2). The first one to prepare the buildings ( Figure 3) and the second one to compute the area.   In addition, the points in the point clouds are classified according to the classification from the American Society for Photogrammetry and Remote Sensing (ASPRS) [40], following a methodology presented in [41] that shows that the performance of the classification for points of the class "buildings" is over 90%. For this reason, the classification provided by the dataset is exploited in this study for the analysis of rooftop area.

Methodology
The methodology proposed consists of computing the area of the rooftops in point clouds with different densities, with the aim of evaluating the influence of the density in the accuracy of the area estimated. For this purpose, it is necessary to use aerial orthophotography as reference data. This kind of data is selected due to its planimetric precision. For both cases, the methodology consists of two main stages ( Figure 2). The first one to prepare the buildings ( Figure 3) and the second one to compute the area.  In addition, the points in the point clouds are classified according to the classification from the American Society for Photogrammetry and Remote Sensing (ASPRS) [40], following a methodology presented in [41] that shows that the performance of the classification for points of the class "buildings" is over 90%. For this reason, the classification provided by the dataset is exploited in this study for the analysis of rooftop area.

Methodology
The methodology proposed consists of computing the area of the rooftops in point clouds with different densities, with the aim of evaluating the influence of the density in the accuracy of the area estimated. For this purpose, it is necessary to use aerial orthophotography as reference data. This kind of data is selected due to its planimetric precision. For both cases, the methodology consists of two main stages ( Figure 2). The first one to prepare the buildings ( Figure 3) and the second one to compute the area.  computation of the area. With the aim at estimating the area contained in a set of points, an intermediate function between all the previous methodologies has been used. Specifically, the "Boundary" function, which is implemented in MATLAB ® . Unlike the convex hull, the boundary can shrink towards the interior of the hull to envelop the points [51]. Thus, the user can specify a shrink factor "s". This is a scalar factor between 0 and 1. If s = 0, this function gives the convex hull. On the other hand, if s takes the value of 1, the result is a compact boundary. For this case, the value of the shrink factor is 0.5 (default value).

Case Studies
Industrial and residential buildings were selected for the computation of the area of their rooftops ( Figure 5) with the methodology proposed in Section 2 and for the evaluation of the influence of the point cloud density on the accuracy of the computation. The study analyzes industrial and residential buildings separately because the first present regular geometry, mostly rectangular, while the second usually present complex and irregular geometries, with more than four edges and angles ( Figure 6; Figure 7). To make the study as homogeneous as possible, a total of 120 industrial

Aerial Orthophotography
A reference value of the area is required in order to calculate the error of the estimation of the area from the point cloud. Although the most accurate procedure would imply the measurement of the areas with classical surveying, the high number of rooftops included in the study and their separate locations make the use of geospatial data necessary. In particular, in this study, the aerial orthophotography with 25 cm spatial resolution is established as reference, since it is the highest resolution data available from all the case studies.
The computation of the area of roofs from orthophotographs requires a digitalization process. In this case, manual digitalization of the roofs is performed using QGIS software [42], where the rooftops are delineated into a shapefile (Figure 4a). Thanks to the different tools of QGIS, it is possible to calculate automatically the area of the shapefile elements. It should be highlighted that the orthophotography is a planimetric product, and consequently the area calculated is the projected area (2 dimensions).
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 21 In computational geometry, there are different methods of boundary extraction that allow the automatic computation of the area contained in a set of points [47]. The most known technique is the convex hull [48]. There are many methods derived from this technique. Some of them are Alpha shape [49] and concave hull [50]. The main difference between these methodologies and the convex hull methodology is the closure form. While the convex hull is a technique that allows a convex enclosure on a set of points, alpha shape and concave hull allow a non-convex enclosure. For the first of the cases, this involves an over-estimation while the second involves an under-estimation in the computation of the area. With the aim at estimating the area contained in a set of points, an intermediate function between all the previous methodologies has been used. Specifically, the "Boundary" function, which is implemented in MATLAB ® . Unlike the convex hull, the boundary can shrink towards the interior of the hull to envelop the points [51]. Thus, the user can specify a shrink factor "s". This is a scalar factor between 0 and 1. If s = 0, this function gives the convex hull. On the other hand, if s takes the value of 1, the result is a compact boundary. For this case, the value of the shrink factor is 0.5 (default value).

Case Studies
Industrial and residential buildings were selected for the computation of the area of their rooftops ( Figure 5) with the methodology proposed in Section 2 and for the evaluation of the influence of the point cloud density on the accuracy of the computation. The study analyzes industrial and residential buildings separately because the first present regular geometry, mostly rectangular, while the second usually present complex and irregular geometries, with more than four edges and angles ( Figure 6; Figure 7). To make the study as homogeneous as possible, a total of 120 industrial

LiDAR Point Cloud
In order to calculate the area using the LiDAR data of each of the buildings selected for the study, it is necessary to know the points that belong to the roofs ( Figure 3). Although the LiDAR point cloud offered by IGN is classified according to the ASPRS, further processing is required. This processing has been automated using MATLAB ® , and consists of the following steps, that are further detailed in [25]:

•
Removal of points belonging to the ground and part of vegetation. The first step relies on the classification of the point cloud in ASPRS classes. Those points with a single LiDAR return and with a different classification than class 2 (ground) continue in the process. In this way, the process Remote Sens. 2020, 12, 3726 6 of 21 not only ensures the elimination of points belonging to the ground but also part of the points belonging to vegetation ( Figure 3b).

•
Removal of remaining points belonging to vegetation areas. A previous computation of the Normalized Vegetation Index (NDVI) [43] is required, using the red and near infrared information collected by the LiDAR and applying Equation (1).
A NDVI threshold of 0.3 [43,44] was established in such a way that points with higher NDVI value are considered vegetation and consequently removed ( Figure 3c).

•
Removal of points belonging to building façades. This segmentation is based on the Z-component of the normal vector (Nz) of each 3D point. Each normal vector is estimated through the adjustment of a 2D subspace that is tangent at the point of interest based on pairwise point relationships [45]. Specifically, a surface is adjusted to the 8 nearest neighbor points of each 3D analyzed point (search resolved by k-d tree) within a maximum radius of 4 m. In this calculation, the covariance matrix of the 8 points is analyzed and the normal vector of the surface adjusted is determined as the normal vector of the analyzed point. A threshold of 0.15 was established in such a way that any point with a Nz value below the threshold is considered a point corresponding to façades ( Figure 3d).

•
Removal of noise points. Possible residual points from the previous steps (such as ground points with a single LiDAR return that pass the first filter) are deleted by applying a statistical outlier removal (SOR) filter [46]. A threshold of 1 m deviation in 30 neighbor points was established. After this step, the resulting point cloud corresponds only to roof areas ( Figure 3e).
Once the point cloud only contains points belonging to the roofs of the buildings (Figure 3e), the area is calculated. Notice that LiDAR point data are a 3D product while the reference data (orthophotography) are a 2D product. Due to this difference, the comparison of the results between both products requires the dismissal of the third dimension from the point cloud (projected area) ( Figure 4).
In computational geometry, there are different methods of boundary extraction that allow the automatic computation of the area contained in a set of points [47]. The most known technique is the convex hull [48]. There are many methods derived from this technique. Some of them are Alpha shape [49] and concave hull [50]. The main difference between these methodologies and the convex hull methodology is the closure form. While the convex hull is a technique that allows a convex enclosure on a set of points, alpha shape and concave hull allow a non-convex enclosure. For the first of the cases, this involves an over-estimation while the second involves an under-estimation in the computation of the area. With the aim at estimating the area contained in a set of points, an intermediate function between all the previous methodologies has been used. Specifically, the "Boundary" function, which is implemented in MATLAB ® . Unlike the convex hull, the boundary can shrink towards the interior of the hull to envelop the points [51]. Thus, the user can specify a shrink factor "s". This is a scalar factor between 0 and 1. If s = 0, this function gives the convex hull. On the other hand, if s takes the value of 1, the result is a compact boundary. For this case, the value of the shrink factor is 0.5 (default value).

Case Studies
Industrial and residential buildings were selected for the computation of the area of their rooftops ( Figure 5) with the methodology proposed in Section 2 and for the evaluation of the influence of the point cloud density on the accuracy of the computation. The study analyzes industrial and residential buildings separately because the first present regular geometry, mostly rectangular, while the second usually present complex and irregular geometries, with more than four edges and angles ( Figure 6; Figure 7). To make the study as homogeneous as possible, a total of 120 industrial buildings and 95 residential buildings were selected (Table 3) with a wide range of areas, from 27 m 2 to 3000 m 2 ( Figure 5).  In the experiments, LiDAR datasets with density of 0.5, 1, 2 and 14 point/m 2 were used. The equipment used for the data collection is different and therefore impacts the RMSE. For the first coverage, the equipment used in two of the three cases is the same and the RMSEs (RMSE xy and RMSE z) are the same (Table 4). Although for the second coverage of LiDAR data the equipment is the same, for case 1 and case 2 the RMSEs are not the same (Table 5).  The orthophotographs used for the computation of the projected area used as reference present the same technical specifications for all the cases of study (Table 2). Although flight pixel size is 35 cm for case 1 and 22 cm for cases 2 and 3, all the orthophotographs generated present the same pixel size of 25 cm. All the orthophotographs were acquired under favorable meteorological conditions, as shown in Figure 6; Figure 7, which include an example of the buildings subjected to the study.

Building Boundary Extraction from Different Point Densities
The methodology for the computation of the projected area (of reference and from point cloud) is applied to all the cases of study (Figure 8). The different values of the projected area were computed from the different point cloud densities towards the analysis of the influence of the density on the

Building Boundary Extraction from Different Point Densities
The methodology for the computation of the projected area (of reference and from point cloud) is applied to all the cases of study (Figure 8). The different values of the projected area were computed from the different point cloud densities towards the analysis of the influence of the density on the accuracy of the projected area calculated.   43 36 In the experiments, LiDAR datasets with density of 0.5, 1, 2 and 14 point/m 2 were used. The equipment used for the data collection is different and therefore impacts the RMSE. For the first coverage, the equipment used in two of the three cases is the same and the RMSEs (RMSE xy and RMSE z) are the same (Table 4). Although for the second coverage of LiDAR data the equipment is the same, for case 1 and case 2 the RMSEs are not the same (Table 5). The orthophotographs used for the computation of the projected area used as reference present the same technical specifications for all the cases of study (Table 2). Although flight pixel size is 35 cm for case 1 and 22 cm for cases 2 and 3, all the orthophotographs generated present the same pixel size of 25 cm. All the orthophotographs were acquired under favorable meteorological conditions, as shown in Figure 6; Figure 7, which include an example of the buildings subjected to the study.

Building Boundary Extraction from Different Point Densities
The methodology for the computation of the projected area (of reference and from point cloud) is applied to all the cases of study (Figure 8). The different values of the projected area were computed from the different point cloud densities towards the analysis of the influence of the density on the accuracy of the projected area calculated.
The error of the area computed from each point cloud density is calculated, using as reference the area value from the orthophotographs: the error of the area is the result of the intersection between the reference area (green line in Figure 8) and the extracted contour (red line in Figure 8 The error of the area computed from each point cloud density is calculated, using as reference the area value from the orthophotographs: the error of the area is the result of the intersection between the reference area (green line in Figure 8) and the extracted contour (red line in Figure 8 (Figure 9). This result appears for both residential and industrial buildings, and for all point cloud densities considered in the study. When the average RMSE is calculated for the buildings under study (Table 6; Table 7), the results are opposite to the relative error. This is because the relative error shows the relative difference expressed in % while RMSE shows the area difference in m 2 .  The error of the area computed from each point cloud density is calculated, using as reference the area value from the orthophotographs: the error of the area is the result of the intersection between the reference area (green line in Figure 8) and the extracted contour (red line in Figure 8 (Figure 9). This result appears for both residential and industrial buildings, and for all point cloud densities considered in the study. When the average RMSE is calculated for the buildings under study (Table 6; Table 7), the results are opposite to the relative error. This is because the relative error shows the relative difference expressed in % while RMSE shows the area difference in m 2 . When the average RMSE is calculated for the buildings under study (Table 6; Table 7), the results are opposite to the relative error. This is because the relative error shows the relative difference expressed in % while RMSE shows the area difference in m 2 .  The results for point cloud density of 14 points/m 2 present lower error than for lower point cloud densities. In this case, the relative error value does not present a difference with the area of the building: the relative error is 6% for residential buildings and 3% for industrial buildings.
For an average building of each type (residential and industrial), the mean computation time for each point cloud density using a notebook with Intel Core i7-6700 processor at 3.41 GHz, 64 bits, 32 Gb. RAM and 931 Gb. disk is shown in Tables 8 and 9. In addition to the relative error in the roof surface computed with LiDAR data, there are other metrics to be taken into account for the evaluation of the validity of LiDAR for roof characterization. These metrics are precision, recall, and F-score [52] defined as follows: where TP, FP and FN are the number of true positives, false positives, and false negatives, respectively. The mean values for Recall and F-score are shown in Figures 10 and 11 according to the different point cloud densities and the use of the building. In the conducted experiments, the precision metrics always reach 100%, since the methodology avoids false positives in boundary extraction.

Influence of the Error in Area on the Computation of PV Energy Production
The design of a solar PV installation on a rooftop requires the knowledge of geometric parameters of the rooftop, among which the area is included [25]. For this reason, the error in the computation of the area of the rooftop will have an influence on the computation of the total solar potential of the rooftop.
Provided that the buildings selected as cases of study are located in Spain, the influence of the area on the solar potential is computed for this location. In all cases, the solar radiation is computed for panels with South orientation, and optimal tilt angle. In this way, the only cause of difference is the error in the area computed. In order to correlate incoming solar radiation with solar potential, PV panels of 330 W are considered. It should be highlighted that this and all further assumptions are equally applied to all the cases of study, in a way to validate the relative comparison.
According to the Spanish Technical Construction Code [53], Spain has five solar climate zones ( Figure 12); their average daily solar radiation values are shown in Table 10.

Influence of the Error in Area on the Computation of PV Energy Production
The design of a solar PV installation on a rooftop requires the knowledge of geometric parameters of the rooftop, among which the area is included [25]. For this reason, the error in the computation of the area of the rooftop will have an influence on the computation of the total solar potential of the rooftop.
Provided that the buildings selected as cases of study are located in Spain, the influence of the area on the solar potential is computed for this location. In all cases, the solar radiation is computed for panels with South orientation, and optimal tilt angle. In this way, the only cause of difference is the error in the area computed. In order to correlate incoming solar radiation with solar potential, PV panels of 330 W are considered. It should be highlighted that this and all further assumptions are equally applied to all the cases of study, in a way to validate the relative comparison.
According to the Spanish Technical Construction Code [53], Spain has five solar climate zones ( Figure 12); their average daily solar radiation values are shown in Table 10.

Influence of the Error in Area on the Computation of PV Energy Production
The design of a solar PV installation on a rooftop requires the knowledge of geometric parameters of the rooftop, among which the area is included [25]. For this reason, the error in the computation of the area of the rooftop will have an influence on the computation of the total solar potential of the rooftop.
Provided that the buildings selected as cases of study are located in Spain, the influence of the area on the solar potential is computed for this location. In all cases, the solar radiation is computed for panels with South orientation, and optimal tilt angle. In this way, the only cause of difference is the error in the area computed. In order to correlate incoming solar radiation with solar potential, PV panels of 330 W are considered. It should be highlighted that this and all further assumptions are equally applied to all the cases of study, in a way to validate the relative comparison.
According to the Spanish Technical Construction Code [53], Spain has five solar climate zones ( Figure 12); their average daily solar radiation values are shown in Table 10.  Five cities distributed by the different climatic zones were selected (Table 11) for the performance of the analysis of the effect of the error in the computation of the usable area of the rooftop, in such way that the effect is analyzed for a wide range of solar intensities. The solar radiation value was obtained through the "Solar radiation tool" offered by PVGIS (PhotoVoltaic Geographic Information System) [54] considering the average altitude of each of the cities (Table 11), a South orientation and an optimal tilt angle. The error analysis was performed for both residential and industrial buildings with a projected rooftop area of 100 m 2 , representing buildings under 200 m 2 , and a projected rooftop area of 400 m 2 as the example of buildings with an area over 200 m 2 . In the case of industrial buildings, which present mostly gable roofs, the computation was performed, taking into account the installation of the solar PV panels only in the slopes with South orientation. Thus, 50% of the total area of the rooftop is considered for the computation of the number of panels and of the solar potential. In the case of  Five cities distributed by the different climatic zones were selected (Table 11) for the performance of the analysis of the effect of the error in the computation of the usable area of the rooftop, in such way that the effect is analyzed for a wide range of solar intensities. The solar radiation value was obtained through the "Solar radiation tool" offered by PVGIS (PhotoVoltaic Geographic Information System) [54] considering the average altitude of each of the cities (Table 11), a South orientation and an optimal tilt angle. The error analysis was performed for both residential and industrial buildings with a projected rooftop area of 100 m 2 , representing buildings under 200 m 2 , and a projected rooftop area of 400 m 2 as the example of buildings with an area over 200 m 2 . In the case of industrial buildings, which present mostly gable roofs, the computation was performed, taking into account the installation of the solar PV panels only in the slopes with South orientation. Thus, 50% of the total area of the rooftop is considered for the computation of the number of panels and of the solar potential. In the case of flat roofs, which are also common among industrial buildings, the reduction in the usable area is applied to consider the tilt angle of the panels and the space needed between panels in order to avoid the projection of shadows.
Tables 12-15 present the solar power losses produced by the errors in the estimation of the area of the rooftop in industrial buildings with gable or flat roof.  Zone 1  14  8  10  3  Zone 2  15  8  11  3  Zone 3  16  9  12  3  Zone 4  17  9  12  3  Zone 5  19  10  14  3 Gable or hipped roofs are the most common in residential buildings. For this reason, only 50% of the projected area has been considered as usable area for gable roofs, considering that one slope is South-oriented and that is the slope where the PV panels are installed. In the case of hipped roofs there is no standard regarding area of each slope and the total area of the rooftop, in such a way that the study was performed installing the PV panels on the 33% of the total area. Tables 16-19 show the losses in the solar potential due to the error in the computation of the area of the rooftop from LiDAR data with different point cloud densities for residential buildings, which present gable or hipped roofs. Table 16. Average daily PV energy production loss (kWh) due to the error in the estimation of the rooftop area for different point cloud densities: residential building with a gable roof of 100 m 2 .

Discussion
The point cloud density that produces the highest error in the computation of the projected rooftop area is 0.5 point/m 2 ; in this case, the error for buildings smaller than 200 m 2 is 30% of the total Remote Sens. 2020, 12, 3726 15 of 21 area. The improvement on the computation of the projected area with double point cloud density (1 point/m 2 ) is an 8% reduction in the error, although the next doubling, from 1 to 2 points/m 2 , does not imply any great improvement: the error increases of 2%. The use of a point cloud density seven times higher, with 14 points/m 2 , implies a reduction in the error of 12% for buildings smaller than 200 m 2 and 6% for bigger buildings. However, the improvement is less intense than the improvement in the difference between point cloud densities. This difference in intensity of improvement is shown in Figure 13: the error reduction for the doubling of point cloud density from 0.5 to 1 point/m 2 has a mean value of 58% in a graphical representation (30 • slope) for all building cases, while the doubling from 1 to 2 points/m 2 presents a flat or an uphill difference in the value of the relative error; thus, no reduction in error is associated with this improvement of point cloud density. The density increases from 2 to 14 points/m 2 generate a mean reduction in the relative error of 1.12%: the mean slope value in the graphic representation of the error is 4% (2 • ).
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 21 higher, with 14 points/m 2 , implies a reduction in the error of 12% for buildings smaller than 200 m 2 and 6% for bigger buildings. However, the improvement is less intense than the improvement in the difference between point cloud densities. This difference in intensity of improvement is shown in Figure 13: the error reduction for the doubling of point cloud density from 0.5 to 1 point/m 2 has a mean value of 58% in a graphical representation (30° slope) for all building cases, while the doubling from 1 to 2 points/m 2 presents a flat or an uphill difference in the value of the relative error; thus, no reduction in error is associated with this improvement of point cloud density. The density increases from 2 to 14 points/m 2 generate a mean reduction in the relative error of 1.12%: the mean slope value in the graphic representation of the error is 4% (2°). A trend can be extracted regarding residential and industrial buildings: the first result in a higher error than the second, probably caused by the higher angle of the roof slope. The exception is found for buildings smaller than 200 m 2 and point cloud density of 0.5 point/m 2 . In this case, the error in the computation of the projected area is smaller for residential buildings than for industrial buildings, with 5% difference. However, not only the computation but the relative error in the estimation of the area shows that the lowest point cloud density (0.5 point/m 2 ) presents a higher error than the rest of the densities. Buildings with an area less than 200 m 2 have recall and f-score lower than buildings with an area greater than 200 m 2 . Factors that, for the case of point cloud of 1 and 2 points/m 2 , barely show variations. Point clouds with the higher density (14 points/m 2 ) present a recall and f-score value very close to 100%. As the methodology does not classify points as false positives, the precision in all the cases is at 100%.
With respect to PV energy production per square meter loss, due to the computation of the area with point cloud data, Figures 14-17 show the RMSE daily PV energy production per square meter loss.
Industrial buildings with flat roofs suffer a PV energy production per square meter loss between 613 and 870 Wh/m 2 . The PV energy production per square meter loss for industrial buildings with gable roof is between 397 and 846 Wh/m 2 . A trend can be extracted regarding residential and industrial buildings: the first result in a higher error than the second, probably caused by the higher angle of the roof slope. The exception is found for buildings smaller than 200 m 2 and point cloud density of 0.5 point/m 2 . In this case, the error in the computation of the projected area is smaller for residential buildings than for industrial buildings, with 5% difference.
However, not only the computation but the relative error in the estimation of the area shows that the lowest point cloud density (0.5 point/m 2 ) presents a higher error than the rest of the densities. Buildings with an area less than 200 m 2 have recall and f-score lower than buildings with an area greater than 200 m 2 . Factors that, for the case of point cloud of 1 and 2 points/m 2 , barely show variations. Point clouds with the higher density (14 points/m 2 ) present a recall and f-score value very close to 100%. As the methodology does not classify points as false positives, the precision in all the cases is at 100%.
With respect to PV energy production per square meter loss, due to the computation of the area with point cloud data, Figures 14-17 show the RMSE daily PV energy production per square meter loss.
Industrial buildings with flat roofs suffer a PV energy production per square meter loss between 613 and 870 Wh/m 2 . The PV energy production per square meter loss for industrial buildings with gable roof is between 397 and 846 Wh/m 2 .
Residential buildings with gable roofs suffer a PV energy production per square meter loss between 453 and 914 Wh/m 2 . The PV energy production per square meter loss for residential buildings with hipped roofs is between 582 and 924 Wh/m 2 . Remote Sens. 2020, 12, x FOR PEER REVIEW 16  Residential buildings with gable roofs suffer a PV energy production per square meter loss between 453 and 914 Wh/m 2 . The PV energy production per square meter loss for residential buildings with hipped roofs is between 582 and 924 Wh/m 2 .  Residential buildings with gable roofs suffer a PV energy production per square meter loss between 453 and 914 Wh/m 2 . The PV energy production per square meter loss for residential buildings with hipped roofs is between 582 and 924 Wh/m 2 .  Residential buildings with gable roofs suffer a PV energy production per square meter loss between 453 and 914 Wh/m 2 . The PV energy production per square meter loss for residential buildings with hipped roofs is between 582 and 924 Wh/m 2 .  If the focus is set on the different solar climate zones, the results are grouped: values of solar power loss are almost equal for solar climate zones 1 or 2 (mean annual solar radiation between 3.8 kWh/m 2 and 4.2 kWh/m 2 ), and for solar climate zones 3 and 4 (ranging between 4.2 kWh/m 2 and 5 kWh/m 2 ). In general, there is a difference of 100 Wh/m 2 loss between each of the groups: 600 Wh/m 2 (zone 1 and zone 2), 700 Wh/m 2 (zone 3 and zone 4) and 800 Wh/m 2 (zone 5).

Results Regarding Other Density Analysis
Drešček has also performed a quality analysis of building outline extraction for different point cloud densities [32]. For this case, the analyzed point cloud comes from a photogrammetric flight carried out with UAV (Unmanned Aerial Vehicles) with 1, 2.5, 5, 10 and 20 cm. point spacing (that is, point cloud densities of 10,000, 1600, 400, 100 and 25 points/m 2 , respectively). Table 20 shows the results obtained from a direct approach. The variation of the three quality parameters evaluated is practically visible. Comparing the results of Table 20 with those obtained for the point clouds evaluated in this article ( Figure 10; Figure 11), it can be observed that the point cloud with 14 points/m 2 has a very similar recall and f-score values to the point cloud with a 20 cm point spacing (25 points/m 2 ). From this comparison, it can be deduced that a significant improvement of the results does not require a significant increase in the point cloud density.

Conclusions
The exhaustive analysis of the point cloud density allows the determination of the accuracy of the geometric characterization of the rooftops performed from point cloud data, and of the subsequent simulation of the energy production in the roof. The results of the analysis have shown that the highest accuracy range is obtained from point cloud densities of 1 or 2 points/m 2 , especially in buildings with an area bigger than 200 m 2 , both residential and industrial. Lower densities present If the focus is set on the different solar climate zones, the results are grouped: values of solar power loss are almost equal for solar climate zones 1 or 2 (mean annual solar radiation between 3.8 kWh/m 2 and 4.2 kWh/m 2 ), and for solar climate zones 3 and 4 (ranging between 4.2 kWh/m 2 and 5 kWh/m 2 ). In general, there is a difference of 100 Wh/m 2 loss between each of the groups: 600 Wh/m 2 (zone 1 and zone 2), 700 Wh/m 2 (zone 3 and zone 4) and 800 Wh/m 2 (zone 5).

Results Regarding Other Density Analysis
Drešček has also performed a quality analysis of building outline extraction for different point cloud densities [32]. For this case, the analyzed point cloud comes from a photogrammetric flight carried out with UAV (Unmanned Aerial Vehicles) with 1, 2.5, 5, 10 and 20 cm. point spacing (that is, point cloud densities of 10,000, 1600, 400, 100 and 25 points/m 2 , respectively). Table 20 shows the results obtained from a direct approach. The variation of the three quality parameters evaluated is practically visible. Comparing the results of Table 20 with those obtained for the point clouds evaluated in this article ( Figure 10; Figure 11), it can be observed that the point cloud with 14 points/m 2 has a very similar recall and f-score values to the point cloud with a 20 cm point spacing (25 points/m 2 ). From this comparison, it can be deduced that a significant improvement of the results does not require a significant increase in the point cloud density.

Conclusions
The exhaustive analysis of the point cloud density allows the determination of the accuracy of the geometric characterization of the rooftops performed from point cloud data, and of the subsequent simulation of the energy production in the roof. The results of the analysis have shown that the highest accuracy range is obtained from point cloud densities of 1 or 2 points/m 2 , especially in buildings with an area bigger than 200 m 2 , both residential and industrial. Lower densities present an important decrease in accuracy (50% when using point cloud density of 0.5 point/m 2 ). Higher densities increase the accuracy, but the improvement is not as intense as the increase in computation time (1.12% increase in accuracy vs. mean increase in computation time from 13 s to 90 s (14.45% time increase) working with point cloud density of 14 points/m 2 instead of 2 points/m 2 ). A density greater than 14 points/m 2 does not improve the quality of the final results to the same extent as the increase in the number of points of the point cloud. In addition, it can be seen that the shape of the roof (gable or hipped roofs) does not influence in the accuracy of the determination of the area.
Although the computation of the rooftop area contains the whole building, solar PV panels are usually not installed in all the usable area. Panels are mostly installed in the slopes of the rooftop with the highest incoming solar radiation. That is, South orientation in the North hemisphere, and North orientation in the South hemisphere. Taking this fact into account, the percentage in error on the computation of the area and on the simulation of energy production per square meter is equal for the complete roof than for part of it. In this case, the essential factor is the incoming solar radiation. An analysis of the five solar climate zones in Spain has shown that the higher the incoming solar radiation on the rooftop, the higher the loss on solar power provoked by the error in the computation of the area. This is an important factor when translating the results of the study to countries with lower solar radiation than Spain: the same point cloud densities for the computation of the rooftop area will provoke lower error in the simulation of the solar power per square meter.
Future studies will deal with the analysis of point cloud densities different from those analyzed in this study. Effort will be put on point cloud densities between 2 and 14 points/m 2 , with the aim of reducing the gap of knowledge between these values, and to determine the highest point cloud density that allows the geometric characterization of the roof optimizing accuracy and computation time. This study will be possible when data from the last flights performed by IGN are made available, with a point cloud density of 4 points/m 2 .