A geographic information system-based large scale visibility assessment tool for multi-criteria photovoltaic planning on urban building roofs

Integration of photovoltaics (PV) into the urban environment will play a major role in the energy transition. However, installing PV systems on building roofs can be challenging, particularly for monumental buildings with strict architectural and social value restrictions. Assessing roof surface visibility is, therefore, key to finding as much permitted roof surface area as possible that may be used for PV installation. In this study, a GIS-based large-scale visibility assessment tool is developed that can assist in evaluating roof visibility, using LiDAR, road networks, and cadastral data as inputs. The tool delivers multi-level outputs, including maps of roof binary visibility, roof visual amplitude, roof PV system layout, roof PV system AC yield, and roof PV module visibility. After optimization, an average speed of 0.12 s/ m 2 is achieved. For each roof surface, an additional sensitivity analysis has been conducted. This step determines the optimal values for two visibility analysis parameters: assessment range and observer spacing, balancing the computational demand and result accuracy. Application of this workflow to the monumental buildings on the TU Delft campus revealed that approximately 2.68 GWh/year of electricity could be harvested from imperceptible PV modules, while an additional 0.42 GWh/year of energy is attributed to PV modules with medium visibility, and 0.37 GWh/year of energy is associated with PV modules with high visibility. This modeling workflow supports the multi-criteria decision-making process for urban roof PV planning.


Introduction
The integration of photovoltaic (PV) into urban environment is set to be the main contributor to the global transition to green energy.Building rooftops are particularly good hosts for PV system installations, considering their vicinity to energy consumption and the space scarcity in urban areas.In fact, studies have shown that rooftop PV systems have substantial potential for green electricity generation.It is stated in a report published in 2014 that the annual rooftop solar electricity production in the Netherlands can reach up to approximately 50 TWh, which comes down to almost 45% of the total national electricity consumption in 2021 [1,2].
Despite the promising energy potential, designing urban PV systems is often challenging due to sub-optimal irradiation caused by unfavorable roof tilt and orientation.Additionally, the complex urban morphology creates intricate shading patterns, which significantly reduce solar access and impede the deployment of PV systems.To overcome these challenges, several urban solar PV analysis tools have been developed.These tools use either aerial images or light detecting and ranging (LiDAR) data to detect the geometry of suitable roof surfaces for PV integration [3,4].Depending on the project objective and the availability of data sources, the result is delivered to one of the renewable energy potential hierarchies (physical, geographical, technical, and economic potential) [5][6][7].
However, the decision-making process in urban PV planning is not solely limited to technical or economic concerns.Social considerations may also play a crucial role in determining the viability of rooftop PV system installation [8].Monumental buildings, for instance, are often restricted from changing their appearance due to their architectural and historical values.Installing PV systems on such buildings requires additional visual impact analysis as the contrasting color of PV modules compared to traditional building materials can affect the building's identity in its context [9].This means that in some cases, identifying suitable roof surfaces for PV installation is not only dependent on the potential energy yield but also highly influenced by the visibility of these roof surfaces to the public domain and their social value.In the Dutch city of Delft alone, there are over 1500 monumental buildings.This makes large-scale visibility assessment of roof surface useful to find as much suitable roof surface area as possible for PV installation.There are three primary methods to perform the visibility assessment: (i) consulting the experts, (ii) surveying the general public, and (iii) using spatial modeling and analysis [10].These three methods are applied flexibly in the context of solar energy planning, with one or https://doi.org/10.1016/j.rser.2023 more methods being used depending on the specific requirement of the project and data source availability.
For the utility-scale solar power plant which is typically large and features centralized installations, it may pose a negative visual impact on the nearby residents by introducing contrasting artifacts into the natural landscape [11].Ana et al. have conducted a combined analysis on the aesthetic impact of solar power plants, where they proposed an expert model that derives four key indicators from photographs: visibility, color, fractality, and the local atmospheric condition.By calculating the weighted sum of these factors, the visual impact of the facility is determined and further validated with the public preference approach [12].The same model is adapted and utilized to study the aesthetic impact of four PV power plants, but the accuracy of obtaining the color indicator from photographs is constrained by the weather condition under which the photo is captured [13].Another method addressing both spatial and perceptual aspects was developed by [14].This approach delivers two visibility maps where the Boolean map indicates the binary visibility of the PV plant from a given location, and the visual perception map conveys the view angle of the PV plant from the affected location.Additionally, by taking into account the number of daylight hours, the potential observation hour (POH) can also be calculated which represents the aggregated value of the maximum number of hours in an average day that an object may be visible to each potential observer [15].In [16], a comprehensive methodology that incorporates all three visibility assessment methods is introduced, in which seven perceptual parameters are used to conduct the visual impact analysis in both quantitative and qualitative manners.
Rooftop solar PV systems, which are typically smaller in size and feature decentralized installations, pose a greater visual impact on the built landscape.The fact that these systems are usually spread across various locations on different buildings makes the visibility analysis more fragmented.A widely used approach to evaluate the visibility and integrity of PV placement is the LESO-QSV method [17].This method establishes two main categories: urban surface criticity including system visibility and context sensitivity, and the system integration quality, which is dependent on the factors such as system geometry, materiality, and pattern obtained from the photograph.It generates an acceptability grid that offers users an indication of the aesthetic impact of the PV system, thereby assisting in the decisionmaking processes.The system visibility is determined with a spatial model that is extensively explored in the research conducted by Florio et al. [10].Based on this model, a comprehensive method is proposed to estimate the energy generation potential of visually-accepted PVs, building energy consumption and economically viable microgrid operation [18].Lingfors et al. approached the visibility analysis from the perspective of vantage points on the ground rather than the building envelope itself [19].It is found that using vantage area instead of discrete vantage points leads to higher visibility for roofs.Drawing on alternative multi-criteria decision-making approaches, Thebault et al. utilized the ELECTRE TRI method to categorize the suitability of roofs for PV integration [20].In a related study, Sun et al. assessed the feasibility of building-integrated PV (BIPV) by examining both solar energy harvesting potential and visual impact [21].Additionally, they investigated the visual impact of hypothetical colored BIPV retrofits within the urban context using saliency mapping algorithm [22].
However, these methods necessitate the inputs such as photographs or vector-based building models, which are either weather condition dependent or difficult to acquire; thus, their applications are limited for scaling up and vast area assessment.The aim of this research is to develop a GIS-based large-scale visibility assessment tool that relies on minimal input datasets.In this work, a tool is developed that evaluates both the visibility of the surface and its visibility degree.This tool operates in a raster-based representation of the environment generated from freely available LiDAR data, along with cadastral data and road networks representing the public domain.The output contains the roof surface visual amplitude map which not only tells the binary visibility of the roof cell, but also shows how visible they are by categorizing the roof cells into low, medium and high visibility groups.The computational demand is reduced by using parallel computing and an observer filter based on the roof surface tilt and orientation.The visibility assessment target can be alternated between roof cells and observer points, depending on which has fewer points.A sensitivity analysis is also incorporated for every roof surface, and the assessment terminates when the optimal assessment range and observer spacing are found.Additionally, the in-house developed solar PV mapping tool is employed to perform the solar PV potential analysis [6].As a result, in addition to visibility maps for roof surfaces, this tool can also deliver results including the possible layout of the rooftop PV system, visibility categories for each module, and their annual AC yield.To demonstrate, this tool is applied to all the monumental buildings on the campus of Delft University of Technology (TUD).
The rest of the content is organized as follows: in Section 2, the methodology applied in this tool is elaborated in detail.In Section 3, the results from each step of the workflow are shown, where one building is selected as a demonstration.The hyperlink to a web scene is also provided, which visualizes the simulation outputs in this work.The limitations of this methodology are discussed in Section 4, together with the possible future improvements.Lastly, conclusions are drawn in Section 5.

Methodology
This section elaborates on the approach employed for visibility assessment in multi-criteria PV planning on building roofs in a large-scale urban environment.The main steps include (1) input data preparation; (2) roof plane detection and extraction; (3) roof surface visibility assessment and optimization; (4) PV module positioning, yield estimation, and visibility assessment.Each step is discussed in detail in its corresponding subsection.A complete flowchart is shown in Fig. 1, in which the steps mentioned previously are grouped into two clusters: input data preparation and visibility and yield assessment.

Input data preparation
The difficulties encountered when performing urban study mainly lie in the perplexing urban morphology.Therefore, the crucial initial step in an urban-centric analysis involves acquiring a precise and comprehensive representation of the urban environment.
LiDAR appears to be a promising solution to address this need, and it has been increasingly used to describe the urban context.It bases on a remote sensing technique that uses laser pulses to densely sample the surface of the earth [23].Generally, for a large-scale LiDAR collection, the sensors are mounted on airplanes which fly over the target area and collect the laser reflection signals as they go.The resulting dataset is presented as a geo-referenced point cloud, meaning that the 3D data is assigned geographic coordinates that allow it to be accurately located in space.In the Netherlands, several versions of LiDAR data are available, where the first batch of collection (AHN1) is dated back to 1997 and the most recent available version was released in 2022 (AHN4).The version employed in this project is the AHN3 dataset in which the data for the city of Delft was collected in 2014, having a resolution of three to five points per square meter [24,25].The original LiDAR data can be quite large in size and difficult to work with.Thus, to make it more manageable, it is often converted into a raster-based digital surface model (DSM) which is executed by the built-in geo-processing tool in ArcGIS Pro.This model serves as the input to the visibility and PV potential assessment tool, and it is essentially a grid that consists of equally sized cells, with each cell assigned a single height that represents the maximum height of that cell.Given the resolution of the AHN3 dataset, the cell size is chosen as 0.5 m, meaning that each cell covers an area of 0.25 m 2 .
Another important input is the public domain which establishes the observer points when performing the visibility assessment.While an infinite number of observer points could exist within the public domain, this is not feasible for practical applications.As such, a finite representation of the possible public domain must be determined.In this project, points that represent roads or streets are chosen, as these are the primary pathways people use for transit between locations; alternatively said, these are the regions where most of the traveling happens and possess a high probability of PV modules being visible.This road network is sourced from GeoFabrik, which provides data in the form of a set of lines that define roadways [26].Equally-spaced points are generated along the lines and subsequently mapped onto the digital terrain model (DTM).An additional 1.7 meters in vertical elevation is added to all the points to simulate the average human eye height.Fig. 2 examples two DSMs of part of the TUD campus with the public domain featured by red points.The primary distinction between these two DSMs is the presence or absence of vegetation: Fig. 2a displays the unfiltered environment and is utilized for conducting PV potential analysis, while Fig. 2b presents the vegetation-filtered environment and is employed for visibility assessment.Removing trees from the DSM is preferred in the visibility assessment because in winter the defoliated trees can be deemed as highly transparent.In fact, the visibility assessment is preferred to be done on a worst-case scenario.In other words, by employing the most open interpretation of the terrain, one can more accurately determine whether the PV system is visible or not.This approach ensures a more stringent analysis, reducing the possibility of overestimating the energy potential and underestimating the visibility of the PV system.
The last input required for the tool is the building footprints which can be obtained from the cadastral data of Basisregistratie Adressen en Gebouwen (BAG) [27].It consists of geo-referenced polygons, each corresponding to a registered household in the Netherlands.This dataset serves as a bounding box to pull out the building points for visibility and PV potential analysis.This procedure is country-independent but is obviously thought for and applicable to the Netherlands, as it is enabled by the digital data readily and publicly available for the country.In the case of other locations, more effort should be put by public and private stakeholders to also readily and publicly provide digital data to support a quicker implementation of photovoltaic technology in the urban environment.

Roof plane detection and extraction
The input data need to be processed to identify the roof surfaces as well as determine their tilt and azimuth angles.The algorithm to realize this step has been extensively discussed in the work carried out by J. Azanza [4].Here, some modifications are introduced to improve its accuracy and computational speed.Fig. 3 presents the modified steps to detect and extract the roof points.The starting point involves pulling out the building point cloud from the DSM using the building footprint as a clipping mask.Then, the roof points are segmented into distinct clusters based on their spatial distribution [28], where two neighboring points are grouped into the same cluster if their Euclidean distance is lower than a threshold or the angle between the LiDAR sensor and these two points is greater than a set value.This is exampled in Fig. 3a in which the red and blue points represent two different clusters.After that, the normal vector of each roof point is calculated and allocated to different bins considering the smallest polar angle difference.These bins essentially contain various combinations of tilt and azimuth angles, where tilt angles vary from 5 to 85 degrees in 10-degree increments, and azimuth angles are separated by 15-degree intervals, starting from 7.5 degrees and ending at 352.5 degrees.Based on the frequency count for each bin, they are sorted in descending order to identify the most frequently occurring orientation of the roof surface.The normal vector corresponding to the bin with the highest frequency count is selected as the primary reference vector for plane fitting [29].If the remaining number of points is greater than 1% of all the roof points within the cluster, the normal vector of the bin with the second highest frequency count is passed in as the reference vector.This process is iterated until there are insufficient roof points (less than 1%) remaining to fit a plane.Finally, the normal vector of the fitted plane (colored in green) is determined as shown in Fig. 3c, with which the tilt and azimuth angles are calculated and assigned to all the points found on the plane.By identifying the most frequently occurring roof orientations and using them as reference vectors for plane fitting, the modified algorithm manages to prioritize the reference vector selection, which enhances the computational speed and accuracy of the plane fitting process.

Roof surface visibility assessment and optimization
Prior to proceeding with the algorithm, a couple of assumptions need to be established given the limitation of data input.LiDAR data do not provide information on material properties or surface textures.Even though coupling the LiDAR data with surface context information is possible by means of satellite and street-view images [30], algorithmic processing of this information to obtain quantitative visibility results that are aligned with mathematical psychology is deemed infeasible within the scope of this study due to its complexity.Another limitation lies in the fact that roads are not perfectly represented by lines, as the latter do not account for the width of the roads, which can affect the final visibility results.Unfortunately, a suitable data source that provides such information on a large scale could not be found.Meanwhile, the social value of the observer points located at different zones of the urban area can also influence the visibility result, as human traffic is expected to be heavier around places with more social value, such as tourist attractions, compared to residential neighborhoods.Therefore, based on these limitations in data availability and implementation practicality, the following assumptions are made for the visibility assessment: • The material coherency, geometric coherency and modular pattern coherency for a rooftop PV system are not included in this algorithm.
• The observer points are generated based on the linerepresentation of road networks with various spacings.• Each observer is assumed to have equal social value in the public domain, meaning that there is no weight assigned to the observer points.
• The PV installations must adopt the tilt of the roof surface, regardless of flat or slanted roofs [31].

Binary visibility
By definition, if a photon reflected or originated from any object manages to enter a human's eye and evoke a visual perception in the brain, this object is deemed as visible.Therefore, to perceive an object, it is imperative to have an unobstructed line of sight between the eyes and the object.Based on this feature, it is possible to identify invisible roof surfaces by performing a line-of-sight (LOS) assessment between the roof surface cells and observers [32].Fig. 4 demonstrates the binary visibility results obtained from the proposed algorithm for the evaluation of visibility between observers and roof surface cells.The sub-figures depict two different observer points situated at different locations on the terrain, while the terrain and surface point remain consistent.The LOS assessment establishes a line between the surface point and the observer, subsequently calculating the angle formed between the observer and each of the points that are interpolated along the line.If there exists any angle larger than the one formed between the surface point and the observer, the binary visibility between these two points is marked as invisible.For the case in Fig. 4a, the observer is visible to the surface point, while in Fig. 4b, the observer is invisible to the surface point.In the mathematical definition, the binary visibility   between a pair of surface cell and observer point can be simply written as: By performing the LOS assessment between the roof cells and their respective observer points, the frequency with which the roof cell is visible from the public domain can be computed.Additionally, this analysis allows for the determination of the length of the public domain, which is essentially the product of the number of visible observers and the spacing between them.These insights are useful for identifying the low-hanging fruits for PV system installation, as they can reveal invisible areas that are free from regulatory restrictions.Furthermore, they provide an indication of the magnitude of the public domain that is visually impacted by PV modules.

Solid angle and visual amplitude
The detection of photons coming from an object is necessary but not sufficient to guarantee a successful perception.Failure to complete the visual task by evoking a visual perception in the human brain, even with the presence of detected photons, renders the object invisible.A way to incorporate this effect is through calculating the amount of field of view this object subtends to the observer, in other words, the solid angle as illustrated in Fig. 5 [14].It is defined as the surface's projection area   onto a sphere centered at the observer, divided by the square of the sphere's radius: The projection area   ranges from 0 to 2, where the latter corresponds to an observer standing directly on the planar surface, resulting in a hemispherical projection.In this work, an analytical solution is applied, which calculates the solid angle that a triangular surface subtends to an observer point [34].To compute the solid angle of a square roof cell, this algorithm divides the cell into two triangles and calculates their individual solid angle separately before adding them together.This process is visually described in Fig. 6 where  is expressed in units of steradians (sr).It should be noted that the inverse tangent function 2 is used in Eq. (3) to keep the calculated solid angle within the range of [0, ] [35].
To determine the vectors for the algorithm, the coordinates of the corner points of each visible roof cell are first generated.Fig. 7 shows a raster-based environment in which two points are depicted.The red point is the observer point, and the blue point represents the center of the visible roof cell.Given the raster resolution of , the vectors can be found as per: The  and  components of the vectors can be readily obtained, as both the raster resolution and the x,  coordinates of the roof cells and the observer points are known.The latter essentially represent the projection of the corner points at the x-y plane.To find the  ,,, component of the corner points, a vertical line is assumed to be extending from the projection and intersecting the roof surface.This simplifies the problem to determining the precise intersection between a line and a plane.A general illustration of this solution, including the interplay between the geometrical elements, is depicted in Fig. 8, and further detailed in Eqs. ( 9) and ( 10): where In the context of finding the z components of the corner points, point Q in Fig. 8 represents a point on the roof surface, with its corresponding normal vector denoted by ⃗ .Point  corresponds to the projection of the corner point on the x-y plane, while the direction vector ⃗  remains fixed as [0,0,1].Therefore, Eq. ( 9) can be rewritten as:  ,,, =  ,,, (11) The Euclidean length of vectors can be calculated by computing their magnitude.
To give a better insight into the amount of field of view a roof cell takes up for the average observer inside the visible public domain, the average solid angle is proposed.It is given in Eq. ( 12), assuming equal weight for all the observer points.The solid angles (  ) for a roof cell are calculated for each visible observer, added up, and the summation is subsequently divided by the number of visible observers   .
To incorporate the visual threshold used in psychophysics, visual amplitude (VA) is introduced [10,18].It is defined as a base-ten logarithm as shown in Eq. ( 13), where  is the average solid angle of the roof cell, and  0 is the minimum angle of resolution presented in steradians.

𝑉 𝐴
The latter is highly dependent on visual contrast, which refers to the difference in appearance between two objects due to their color, material, and texture variations [36].It determines the saliency of the target roof cell; in other words, it affects how easily the roof cell can be distinguished from the surroundings if this roof cell were to be covered by PV material.A significant variation in the material properties between PV and roof surface means a greater likelihood for the PV to be perceived given the same visual angle that subtends to the observers.This distinction creates a huge visual contrast, leading to a low minimum angle of resolution.Conversely, when the material properties of PV and the roof surface resemble, the visual contrast decreases and the minimum angle of resolution can be set larger.This adjustment brings down the VA given the same geometrical parameters.Consequently, the roof cell that was previously classified to the high-visibility group, may now be assessed as having medium or low visibility.However, accurately modeling the visual contrast of roof cells poses a significant challenge due to the lack of detailed material data, which can result in large uncertainties in the modeling process.To address this limitation, a fixed minimum angle of resolution of 1 square minute of arc, equivalent to 8.46 × 10 −8 [sr], is used in this project.This threshold corresponds to the scenario of a standard observer in conditions of infinite luminance contrast, which is aligned with the worst-case-scenario visibility assessment approach, ensuring that the visibility is not underestimated.This threshold can be adjusted if the material in the context is provided [36].For a planar surface, the VA ranges from minus infinite to 3.94 where the lower and upper limits can be found in the cases when the surface is infinitesimal, and the surface creates a hemispherical projection of 2 [sr], respectively.Fig. 9 lists some visual thresholds determined experimentally by assessing the visual response of different groups of volunteers [37].Three visibility groups are classified where a surface with a visual amplitude lower than zero is classified as having low visibility and can be barely detected by most healthy individuals; a surface with a visual amplitude between 0 and 1.4 is considered as medium visibility and can be perceived by the majority of people; surfaces with a visual amplitude falling within the higher end of the scale are classified as highly visible and are typically perceivable by everyone except those with visual impairments.These categories are adopted in this project to classify the visibility of the roof surface cell, assisting in better understanding and quantifying the visual impact.

Optimization
Performing many LOS assessments is time-consuming, and it is found to be the primary bottleneck of this algorithm regarding computational speed.Several approaches are adopted to make this algorithm more computationally efficient, including task allocation optimization, unnecessary LOS elimination, and adaptive observer selection.Additionally, the assessment range was increased and the observer spacing was reduced incrementally to further enhance the efficiency of the algorithm while maintaining its accuracy.
The most straightforward way to reduce computational burden is by using parallel computing [38].Parallel computing allows for the allocation of tasks among multiple cores and processes them simultaneously.Therefore, it optimizes efficiency by harnessing the collective power of multiple processing units, but its performance is ultimately limited by the physical constraints of the computer hardware.In this project, the simulations were performed on a laptop with an Intel(R) i7-8750H CPU (6-core) running at a base frequency of 2.2 GHz and 16 GB of DDR4 RAM.
The second approach used for speed improvement is by eliminating the observers that are lying behind the roof surface, considering the fact that the roof surface is only visible to the hemisphere that it is facing.This process is realized with Eq. ( 14), ( 15), ( 16) and illustrated in Fig. 10.
where  is roof tilt and  is the roof azimuth.If this condition is satisfied, the observer point is identified to be in front of the roof surface and a LOS assessment is performed.For a flat surface that is lying above the observer height, this observer filter eliminates all the observer points, drastically reducing the computational time without impacting the final results.Additional refinement to the filter is required if the PV module is mounted with a specific tilt on the flat roof.In such cases, PV modules close to the roof edge may still be visible from the public domain.
The third technique employed to reduce computation time is adaptive observer selection as illustrated in Fig. 11.LOS assessment evaluates the binary visibility between the point of interest and the public domain by constructing multiple lines connecting them.This process can be executed either based on the roof cell or the observer point.By adaptively selecting the assessment target between the roof cell and observer points, depending on which one has fewer points, the number of iterations is decreased.Consequently, this approach minimizes computational time by reducing the number of LOS assessments required.
To achieve an optimal balance between the accuracy of the visibility results and the computational time, it is important to select an appropriate assessment range.Both binary visibility and visual amplitude results are influenced by changes in the assessment range.Generally, increasing the assessment range results in a larger public domain, potentially leading to the identification of new visible roof surface cells.However, a convergence point will be reached, beyond which further increases in the assessment range yield no significant changes in visibility outcomes.Fig. 12 visualizes the method for determining the optimal assessment range.The algorithm operates on a roof-specific basis, meaning that the optimal assessment range is found for the entire roof surface rather than one single roof cell.This is achieved by calculating the total binary visibility   (the total number of visible roof cells) and the average VA of the roof surface within the assessment radius according to Eqs. ( 13) and (17).
where   is the average solid angle of each visible roof cell.While executing, the algorithm incrementally increases the assessment range (by 20 m) and calculates the change induced by the inclusion of new observer points.If the difference between the results for both   and    after th and ( − 4)th increments is lower than 5%, the algorithm terminates, and the optimal assessment range for the specific roof surface is determined.A minimum assessment radius of 100 m is set as default.Roof surfaces that are generally visible to distant observers will automatically be assessed for extended ranges.Conversely, roof surfaces that are visible only to nearby observers will stop further examination.With this approach, the assessment process for each roof surface is tailored, optimizing both accuracy and computational efficiency.Besides, another optimization process is conducted on the observer spacing, as it can also influence the visibility results.Given the assessment range determined from the previous step, the observer spacing is halved with each iteration, as illustrated in Fig. 13.The   and    of the roof surface are computed for the newly included observer points, and the results are cumulatively added to the previous iteration.If the difference between adjacent iterations is lower than 5%, the observer spacing is determined.In this project, the maximum observer spacing is set to 600 cm.

PV module positioning, yield estimation, and visibility assessment
In the multi-criteria decision-making process for roof PV planning, visibility constitutes only one aspect.As the primary purpose of installing PV modules is solar energy harvesting, evaluating the PV potential of roof-mounted PV systems serves as a vital piece of the puzzle to achieve a comprehensive understanding before concluding the roof PV planning.In this work, the PV potential analysis is conducted using the in-house developed workflow [3,6,39].This workflow fits the PV modules to the detected roof surface and calculates its annual AC yield considering the complex skyline where the PV module is located within [40].An economic threshold of 650 kWh/kWp is selected to exclude the PV modules whose installation is not economically viable, ensuring that only feasible options are considered in this analysis [6].Meanwhile, for each fitted PV module, the visibility is evaluated by generating a 2 by 2 grid on the PV module surface; thus, each point represents a quarter of the module surface area.These module surface points are used to first perform the PV potential calculation, followed by LOS analysis and VA calculation considering the previously established roof-specific assessment range and observer spacing.The final VA for each module is determined by averaging the sum of the four fractions, and further categorized based on Fig. 9.In general, a denser grid results in a more accurate visibility map for the PV module, but this increased accuracy comes at the expense of computational time.Eventually, the simulated AC yield for PV modules in each VA category is reported, providing useful insights for rooftop PV planning from both economic and environmental perspectives.

Results and discussion
In this section, the outputs from the visibility assessment tool for multi-criteria urban roof planning are presented.This includes the results for roof detection and extraction, roof binary visibility map, roof VA map, roof PV module layout, roof PV AC yield map, and roof PV visibility map.The computational speed for each step of the workflow is also presented.To better demonstrate the workflow in a consistent manner, the building of Electrical Engineering, Mathematics, and Computer Science (EEMCS) faculty on TUD campus is selected as a study case.An aerial image of this building can be seen in Fig. 14.

Roof detection and extraction
The outcomes of the roof detection and extraction algorithm applied to the building of EEMCS faculty are demonstrated in Fig. 15.The initial point cloud of the building, obtained from cropping the DSM using the building footprint, is illustrated in Fig. 15a.Subsequently, the extracted roof planes are depicted in Fig. 15b, with each roof plane assigned a unique color for differentiation.The minimum area for a roof plane to be registered is set to be 10 m 2 , which corresponds to a requirement of more than forty points for the identified roof plane.The majority of building points that are excluded due to this threshold are found to be rooftop artifacts or curved edges characterized by distinct normal vectors.Therefore, this exclusion has very little impact on the final results.

Roof binary visibility and VA maps
The binary visibility result is depicted in Fig. 16, which is essentially a binary map in which the yellow points correspond to visible roof cells and purple points denote invisible roof cells.As expected, flat roof surfaces remain obscured from the public domain.The analysis reveals two distinct visible sections: the first one is located at the saw-toothshaped roofs, with a corresponding Google Street View image displayed aside.The west part of this roof section sees the streets, leading to the majority of the western roof cells being visible.The other section includes the southern roof of the ESP lab, which also directly faces the streets, as evidenced by the corresponding Google Street View image.Although trees are present at the lower region of the roof, the visibility analysis was conducted using a DSM in which vegetation has been removed as introduced in Section 2.1, leading to an unobstructed LOS between the roof cells and the observer points.
Fig. 17 shows the normalized VA (nVA) result for the building of EEMCS faculty.It is calculated by dividing the VA by the maximum solid angle of a planar surface (3.94 sr).As the VA is solely calculated for the visible roof cells, the non-zero VA results are primarily distributed across the saw-tooth-shaped roofs and the southern roof of the ESP lab.The color scale indicates higher nVA values using lighter colors, illustrating that the west edge of the saw-tooth-shaped roofs exhibits a compatible visual impact on the public domain, predominately falling within the low to medium visibility range.Variations in this section are not significant due to the similarity in both their distance and the viewing angle relative to the public domain, meaning that the subtended solid angles are closely aligned.By contrast, the southern roof of the ESP lab demonstrates high visibility due to its shorter distance between the roof and the public domain.Additionally, a descending trend in visibility is observed from the bottom to the top edge of the roof.This can be mainly attributed to the distance and the viewing angle, which dictate the solid angle of the roof surface cell subtending to the observers.Considering the average eye height (1.7 m) of the observers, the bottom of the roof has a shorter distance as well as a smaller viewing angle relative to the observer points, resulting in a higher nVA.

Roof module placement and AC yield map
The Roof PV module placement result is shown in Fig. 18a, and the corresponding AC yield map for all the fitted PV modules is displayed in Fig. 18b.PV modules installed on slanted roofs adopt the tilt and orientation properties of the roof, while those placed on flat roofs are  configured to align with the longest edge and set to have a zero-degree tilt in compliance with local regulations as introduced in Section 2.3.Besides, this process uses the non-filtered DSM to account for potential shading effects caused by the surrounding urban environment, such as the trees in front of the roof of the ESP lab.For economic feasibility, modules with a simulated specific yield of less than 650 kWh/kWp are excluded.

Roof PV visibility map
To categorize the visibility of the fitted PV modules, the VA of each module is calculated and the result is presented in Fig. 19.Instead of using numerical categorization, the visibility of the PV modules is divided into low, medium and high visibility groups, using the standard classification introduced in Fig. 9.It is evident from the figure that all   low visibility, depending on their placement location.Generally, PV modules situated close to the west edge have higher visibility, which degrades gradually as they approach the east edge, eventually becoming completely invisible.Compared to the visibility results of roof surface cells in the same location, PV modules feature a larger area so that the solid angle subtending to the observers is more substantial, leading to an increased VA, namely a higher visibility.To gain more insights for multi-criteria decision making for urban roof PV planning, the annual AC yield of the PV modules from each visibility category is calculated, with the results shown in Table 1.The same workflow has been applied to the rest of the monumental buildings on TUD campus (refer to Fig. 20), and the results are listed in Table 2.It turned out that approximately 2.7 GWh annual AC yield can be expected from the installed PV modules that are imperceivable from the public domain, which accounts for 77% of the total potential on such monumental buildings.All the results are visualized and published as a webscene, with its screenshot shown in Fig. 21.
According to the local regulation in the city of Delft, a PV system can be installed on monumental buildings without permits as long as the installed PV systems are barely visible or invisible from the public domain [31].This simulation framework serves as a guideline for the decision-makers to identify these low-hanging fruits and helps streamline the initiation and execution processes of PV projects.Meanwhile, the regulation also emphasizes that installing PV systems on flat roof surfaces is preferred over placing them on tilted ones.Identifying the flat roof surfaces is one of the intermediate results from this simulation

Table 2
Annual AC yield results of the PV modules from different visibility categories for all the monumental buildings on TUD campus.The building number is in agreement with the numbering in Fig. 20.

Building [-]
Low visibility framework, thereby assisting in further expanding selection of preferred installation sites.Even for the roof surfaces with some degree of visibility, this workflow can still provide valuable information on their location and PV potential, setting the stage for any potential subsequent on-site visibility assessments.

Computation speed
The computation speed of each step of the workflow applied to all the monumental buildings on TUD campus is shown in Table 3.In most cases, the majority of time is spent in determining the optimal assessment range and observer spacing, primarily due to the timeconsuming nature of LOS assessment.For buildings B4, B7, B9, and B12, the visibility assessment steps take virtually no time, as these buildings only have flat roofs and do not require LOS analysis.On average, the workflow takes 0.12 s/m 2 , with the maximum and minimum values being 0.27 and 0.04 s/m 2 , respectively.This metric is obtained by dividing the time required for the entire workflow by the overall roof surface area processed.The difference in the computation time mainly lies in the level of complexity of the building layout, as those with a greater variety of roof types tend to include more observer points in the urban context, necessitating more LOS assessments.

Accuracy evaluation
The entire workflow has employed techniques that are broadly recognized for their accuracy, and the best practices have been followed when applying these methods.However, it is important to highlight that validating the specific implementation of these methods, especially on visibility outcomes, presents inherent challenges due to the subjective nature of visibility.The perception result can vary widely depending on the conditions of the observers; for instance, an individual who has been informed about the potential existence of PV modules on the monuments might perceive more PV modules than an uninformed observer.Therefore, the rigorous validation would have to implement the expert-based model or conducting questionnaires, where photos of the PV system from various observer points are taken and graded to establish a visibility scale.Such an approach, unfortunately, is not considered within the scope of this project.Instead, the simulation results from this work are serving as a guideline for decision-makers in identifying the low-hanging fruits.For roof surfaces with possible visibility concerns, this workflow still provides valuable location and PV potential information, which can be instrumental if a more detailed on-site visibility assessment is required.

Improvement in the public domain representation
In this study, the public domain is defined by a finite set of observer points generated from road networks, implying a one-dimensional representation.However, in reality, roads have a width dimension, and individuals standing at different positions on an identical road may lead to a different binary visibility result relative to the roof surface.Meanwhile, open spaces such as parks and squares cannot be accurately defined within the public domain because people can freely roam around this place.Therefore, it is necessary to explore the possibility of using a 2D representation (such as polygons) to define the public domain and perform visibility assessments based on the grid generated from this representation.A sensitivity analysis should be performed to determine the optimal grid density, striking a balance between computational time and result accuracy.

Improvement in the human traversal model
Observer points in the public domain are assumed to have equal weight, which may not accurately reflect real-world conditions.Streets with higher traffic, such as those near tourist attractions, are bound to have greater significance compared to the less-traveled ones.Therefore, to create a more realistic representation, it is necessary to incorporate an additional weight dimension for observer points by integrating human behavioral models and traffic data, thereby accounting for the varying importance of different locations in the public domain.

Improvement in the material information of urban elements
Materials can influence the visual contrast between the PV modules and their background.In the current analysis, an infinite visual contrast is assumed, which does not accurately reflect real-world outdoor urban environments.Integrating the material information into the large-scale LiDAR data allows for a roof-specific visual contrast identification, thereby enhancing the robustness of visibility assessment results.

Improvement in PV module mounting configuration
To comply with the local regulations in the city of Delft, PV modules are only permitted to be mounted horizontally with respect to flat roofs.However, as regulations may differ by the cities, investigating alternative PV configurations can bring more insights into the economic aspect of the multi-criteria decision-making process for urban roof PV planning.It is important to note that configuring the PV system with different tilts and orientations can lead to additional visible occasions, such as the visibility of tilted PV modules near the edge of the flat roofs.Adjustments to the workflow are needed to incorporate these scenarios to avoid underestimating the visual impact of the PV systems.

Fig. 1 .
Fig. 1.Workflow of the methodology divided into two parts: data preparation, and visibility and yield analysis.The color code for the blocks: blue for input data; yellow for intermediate output; green for operation; and red for final output.Acronyms: Digital Surface Model (DSM), Line-of-Sight (LOS), Visual Amplitude (VA).

Fig. 2 .
Fig. 2. DSMs of a section of the TUD campus, featuring the mapped public domain as red points: (a) trees included; (b) trees removed.

Fig. 3 .
Fig. 3. Steps for roof detection and extraction: (a) The roof points are extracted with the building footprint and then segmented into clusters.Red and blue points represent two different clusters; (b) For each cluster, the normal vector of each roof point is calculated and grouped into bins; (c) Finally, the roof plane is fitted using the normal vectors obtained in the previous step as reference vectors.The normal vector calculated from fitting the plane is assigned to all the points belonging to the plane and further used for tilt and azimuth determination.

Fig. 4 .
Fig. 4. Line-of-sight assessment for binary visibility assessment: (a) the line of sight is not blocked by the terrain or any artifacts and the observer is visible to the surface point; (b) the line of sight is obscured by the terrain and the observer point is invisible to the surface point [32].

Fig. 5 .
Fig.5.Depiction of solid angle subtended by a planar surface as viewed by a person centered at the gazing field with a radius of r, where ñ is the normal vector of the planar surface, and   is the projected area of the planar surface on the gazing sphere[33].

Fig. 6 .
Fig. 6.The surface cell is divided into two triangles whose solid angles  1 and  2 are calculated separately and summed up to determine the solid angle of the surface cell.

Fig. 7 .
Fig. 7.The calculation of the Euclidean distance between the observer point and the four corners of the surface cell in a raster-based representation of the environment, where the red point is the observer point, and the blue point is the center of the surface cell.

Fig. 8 .
Fig. 8. Illustration of finding the intersection between a line and a plane.Point Q lies on the plane defined with the normal vector ⃗ . is an arbitrary point outside the plane, defining line L with direction vector ⃗ .M is the intersection point between line L and the plane P.

Fig. 9 .
Fig. 9. Categorization of the VAs based on the thresholds determined experimentally by assessing the visual response of different groups of volunteers.Three groups are classified, indicating a low, medium, and high visual impact [37].

Fig. 10 .
Fig. 10.Visual depiction of hemisphere filter of a roof plane with a 135 • azimuth and 45 • tilt.The green points are in front of the roof plane which potentially can be visible; red points are behind the roof plane, and they are eliminated from LOS assessment.

Fig. 11 .
Fig. 11.Visual depiction of the adaptive observer selection method.The LOS can either be performed from (a) the observer point or (b) the roof cell, depending on which group has fewer points.

Fig. 12 .
Fig. 12. Visual depiction of the determination of the assessment range for each roof surface, with an incremental step of 20 m.The yellow point represents a roof surface cell.The green points are the observer points within the assessment range, while the red points are not falling within the assessment range.

Fig. 13 .
Fig. 13.Visual depiction of the determination of the optimal observer spacing for each roof surface.The yellow point represents a roof surface cell.The left figure shows the green observer points with spacing S, and the right figure halves the observer spacing by adding red observer points.

Fig. 14 .
Fig. 14.Orthogonal aerial image of the building of EEMCS faculty, which is highlighted with its footprint in yellow [41].

Fig. 15 .
Fig. 15.Roof detection and extraction from the point cloud of the building of EEMCS faculty, where (a) the unclustered original point cloud is shown; (b) the detected and extracted roof planes are displayed, and each roof is represented by a different color.

Fig. 16 .
Fig. 16.The binary visibility map for the building of EEMCS faculty, where the visible roof cells are in yellow and the invisible roof cells are in blue.

Fig. 17 .
Fig. 17.Normalized VA map for the building of EEMCS faculty.The lighter the color, the higher the visibility of the roof cell.

Fig. 18 .
Fig. 18.(a) PV module fitting results where the modules that failed in passing the economic threshold are excluded; (b) The AC yield map of all fitted PV modules.

Y
.Zhou et al.

Fig. 19 .
Fig. 19.Visibility categories of all fitted PV modules, where three categories are identified, namely low, medium and high visibility.

Table 1
Annual AC yield results of the PV modules from different visibility categories for the building of EEMCS faculty.
the PV modules on flat roofs have low visibility, signifying that most of the individuals moving in the public domain cannot perceive them.PV modules installed on the southern roof of the ESP lab exhibit high visibility, indicating that they are typically perceivable by everyone except those with visual impairments.As for the PV modules mounted on the saw-tooth-shaped roofs, their visibility ranges from high to Y.Zhou et al.