A GIS-Based Method for Identification of Wide Area Rooftop Suitability for Minimum Size PV Systems Using LiDAR Data and Photogrammetry

: Knowledge of roof geometry and physical features is essential for evaluation of the impact of multiple rooftop solar photovoltaic (PV) system installations on local electricity networks. The paper starts by listing current methods used and stating their strengths and weaknesses. No current method is capable of delivering accurate results with publicly available input data. Hence a different approach is developed, based on slope and aspect using aircraft-based Light Detection and Ranging (LiDAR) data, building footprint data, GIS (Geographical Information Systems) tools, and aerial photographs. It assesses each roof’s suitability for PV deployment. That is, the characteristics of each roof are examined for fitting of at least a minimum size solar power system. In this way the minimum potential solar yield for region or city may be obtained. Accuracy is determined by ground-truthing against a database of 886 household systems. This is the largest validation of a rooftop assessment method to date. The method is flexible with few prior assumptions. It can generate data for various PV scenarios and future analyses.


Introduction
Solar energy has a vital role to play in the sustainable energy mix for achieving the 2020 deadline of the Paris Agreement [1]. This requires all cities to increase the supply of green energy and change consumption patterns. Knowledge of quantity and output timing of solar energy becomes increasingly important. It can help with managing self-consumption, pay-back times, profitability, carbon trading, and financial attractiveness of domestic PV systems [2]. There is some evidence that PV has less impact on the environment than other renewable technologies e.g., solar thermal [3], although use of flat plate collectors has reduced CO2 emissions in Greece [4]. Efficiency of PV systems and their simulation are improving [5].
Exact calculation of the solar exposure on pitched roofs is essential for modelling photovoltaic (PV) deployment in domestic scenarios. Two popular tools for calculating in-plane radiation and PV energy production are PVGIS (Photovoltaic Geographical Information System) [6] and PVWatts Calculator [7]. These are free, web-based, do not require installation, and are based on minimal inputs. In order to produce results, these tools either require the user to enter system tilts and azimuths, or employ default or optimal values. Commercial photovoltaic software and calculators [8][9][10] produce quick reliable information of PV potential, but also need the user to supply inclination and orientation of system. This research fills this gap. The UK lacks a national database capable of providing the necessary information. This is also a common omission in other countries. The EU Buildings Database [11] contains information on gross floor area and roof insulation type but has no details of roof inclination, orientation or pitched dimensions. The US is slightly better provided for with a report which assesses the rooftop solar PV potential of 23% of buildings nationwide [12]. In this paper, previous works [13][14][15] are revisited and expanded upon. That is, automated extraction of building roof plane features over extensive areas, shading techniques, and the influence of module orientation on yield are studied. An approach which calculates solar exposure on tilted roofs, based on grid squares or pixels, is proposed. Open data building footprint maps are combined with LiDAR (aircraft-based Light Detection and Ranging Data). These are analysed statistically within a Geographical Information Systems (GIS) environment. The new approach is ground-truthed against a database of about 2000 rooftop domestic solar systems deployed across the city of Nottingham, UK [16]. Up to the present, no other rooftop PV capacity studies have been as rigorously validated.
Urban PV is growing rapidly. Mono-and polycrystalline modules currently hold up to 90% of the market [17]. Recent developments in concentrating photovoltaic systems (CPV) to reduce weight and size have enabled this high efficiency technology for rooftops in areas receiving large amounts of direct sunlight [18,19]. To date in the UK, rooftop PV has been deployed faster than expected. 20% of total installed solar PV capacity (2.6 GW) comes from small scale 0 to 4 kW installations [20]. Net-metering to encourage self-consumption is taking place in many countries worldwide, including the UAE, Lebanon, Chile, and parts of India [21]. In the Netherlands household installations comprise 90% of PV capacity [22]. In China, rooftop distributed installations reached 15 GW in 2017, a rise of 400-500% compared to 2016 [23].
The modern urban environment delivers a variety of problems and advantages for PV deployment. Many neighbourhoods display complex amalgams of roof pitch (tilt) and building orientation (azimuth). Roof features such as aerials, chimneys, stack vents, velux, cross-gables, and dormers may further reduce possible system size.
Roof pitch is associated with building age and roofing material. These in turn are related to geographic regional variation in construction methods. Building direction (azimuth) is largely governed by road layout which is a reflection of local topography.
The tilt and azimuth of a PV installation has two major effects on energy output. Firstly, there is a larger or smaller amount of total annual yield: dependent on how well the roof slope and azimuth match the average sun position over the year [24]. Secondly, there is an impact on the seasonal or daily timing of peak energy generation [25]; although in future energy storage is expected to enhance local generation [26]. PV and battery components will need to be selected so as to minimize environmental impacts [3].
Module performance is also affected by variable operating conditions e.g., temperature and wind. Increased operating temperatures decrease efficiency. Convective heat loss is dependent on wind direction (i.e., PV system orientation), amongst other factors [27]. Solar PV modules may capture the maximum solar radiation by inclination at an optimal angle dependent on sun path (site latitude) and typical weather (diffuse fraction of solar radiation), which is about 38 degrees for most of the UK. Roofs which have higher or lower pitches than this optimal value receive less irradiation.
Current building stock does not always permit the use of this optimum and trade-offs in placement are necessary. Similar mismatches occur in other countries. Whereas the traditional UK roof pitch (40-50°) [28] tends to be steeper than the optimum angle for PV, in Mediterranean regions traditional tile roofs (20-25°) [29] are frequently shallower than the optimum 30-33°. Azimuth may also be less than ideal. Although total annual yield is lowered by non-optimal building orientation there may be positive side effects such as higher morning, evening and winter generation.
Analysis of roof characteristics and their impact on PV output and timing for an individual house is relatively straightforward. This research provides an efficient method for analysing areas too large to investigate manually.

Research Methodology
There has been extensive prior research into computerised recognition of 3D structural characteristics, including several excellent literature reviews [30][31][32]. The authors of a previous paper [30] divide rooftop area estimation methods into three. First, the constant value methods approach the problem collectively by scale-up [33,34]. Although quick and easy, they employ broad assumptions and produce generalised results. Second, manual selection [35] is time-consuming. Third, GIS-based methods [13] deliver detail and may be automated but require substantial computing power. No technique is widely accepted as definitive. This paper concentrates on GIS-based methods.
This research topic presents difficulties with regard to both noisy data and the vastness of LiDAR datasets. Recognition of three-dimensional structural features is nontrivial. This paper commences by testing current methods which take aerial photogrammetry and/or LiDAR data as inputs. The positive and negative results are discussed. Present methods include model-driven, peak searching, iterative voting, edge detection with LiDAR, edge detection of images, image recognition, and hill shading with ambient occlusion.
A LiDAR dataset comprises a grid of ground elevation values produced by an aircraft at constant altitude firing a pulsed laser to Earth and measuring the returns with a sensor. Data resolution is a resultant of the number of returns per square metre. LiDAR can give precise heights of objects (e.g., houses and trees), as well as surface undulation. In the UK, LiDAR data for limited areas and scales of resolution may be obtained from the Environment Agency [36]. One meter resolution was selected to accommodate both accuracy and availability. It covers approximately 70% of England. Aerial photography was obtained from Google Earth [37].
Since none of the existing methods was found to be adequate using the available input data, a new approach is elaborated. This is expedient for medium resolution LiDAR. Previous methods may result in imprecise values unless high resolution data is available. Rather than attempting to obtain exact roof areas, each roof is assessed for PV suitability. Suitability here is defined as a roof to which at least a minimum size photovoltaic system may be fitted (8 m 2 of roof area corresponding to a 1 kW system, assuming a module size of 1.6 m 2 [30]), with an azimuth east through south to west and tilt of between 15° and 60°. (Most UK homeowners install a 1 kW to 4 kW solar system on their roof [38].) Thus, the minimum potential solar yield for region or city may be obtained. The need to generate accurate roof areas and PV system sizes from inexact data is avoided because the goal is the minimum requirements for domestic PV only. Nor is this method intended to separately identify multiples of 1 kW. This is not possible with the available input data and domestic systems are almost entirely comprised of smaller multiples, e.g., 3.76 kW.
LiDAR heights on roof tops only are selected by clipping them out using building outlines (from OS MasterMap Topography Layer [39]) as patterns or "cookie cutters". The azimuth and tilt of each roof gridpoint is calculated by weighted least squares fit of a plane to a 3 × 3 pixel neighbourhood, centred on each LiDAR point [13]. Due to LiDAR inaccuracies, flight paths, and chance, this will result in a unique value for each pixel. Roofs slope evenly, so a statistical technique smoothing azimuth pixel values into groups is identified. Next the roof is divided into separate planes according to the grouping. The size, tilt, and azimuth of each plane can be deduced from the aforementioned calculation. A check is performed to ascertain whether a system of 8 m 2 can be mounted on any of the roof faces.
The methods are tested on two case studies. The Wollaton Park area of Nottingham, UK was chosen because its houses display many different architectural styles. The second case study is a database of approximately 2000 housing association residential rooftop installations in Nottingham. System sizes, placement, and installer records of the tilt and azimuth of individual systems have been gathered from a monitoring website.

Simple Roof Slope and Azimuth Extraction
Trigonometry may be applied to the LiDAR grid of roof point heights to calculate tilt and azimuth of every grid square or pixel (2 m, 1 m, 50 cm, 25 cm, or 15 cm, depending on the resolution available in the area of interest). The traditional approach is to group pixel values obtained from a 3 × 3 neighbourhood by compass point bands to produce realistic roof planes [40]. However, Figure 1 illustrates the problems which may occur. The building in the aerial photograph has a simple roof layout, comprising one north and one south-facing roof section. (Note: the overhead perspective images in this paper are of poor quality. This is a reflection of the data which is publicly available and is part of the problem which this paper seeks to address.) Whilst the azimuth diagram generated from LiDAR reproduces the two sections, there are many spurious small roof planes pointing in various directions. These result from the presence of chimneys, TV aerials, etc., as well as overhanging trees and surrounding structures such as garages. In the case of more complex roof structures, these pseudoplanes can be difficult to distinguish from genuine dormers and porches. Extraction of roof geometry from LiDAR has been the subject of extensive research over the last ten years. Existing solutions are categorised, reviewed and tested in the following paragraphs.

Model Driven Methods
This approach comprises the matching of the irregular roof segment shapes obtained from LiDAR to the best-fitting model in a library of basic building shapes. Jacques et al. (2014) [41] utilised it to classify small buildings in the city of Leeds, UK, using a restricted catalogue of common roof profiles (gabled, hipped, flat, complex, or unclassified).
Model-based methods do not work well for multifaceted roof shapes and intricate building construction. Looking at the topic from a country-wide perspective, there are numerous possible building types. Internationally, roof type is just as varied [42]. Some authors list as many as 50 categories with multiple subcategories. The Geograph Britain and Ireland Project [43], which collects representative photographs for every square kilometre of the nation, has captured examples of over 25 different roof profiles. Some have very different forms (e.g., flat, round, or hipped dormer). Due to the multiplicity of possible model shapes, this line of research was not pursued.

Histogram Discrimination/Peak Detection
This approach is perhaps the oldest and simplest. Peaks are searched for in elevation (above ground level), tilt, or azimuth histograms and used to segment the data. Spatial planes are fitted for each segment. Theoretically, simple gabled roofs should display a rectangular height histogram and that of hip roofs should resemble a trapezium [44]. In fact, these ideals cannot be achieved with real data, as Figure 2 explains. The Wollaton Park (Nottingham) house in the Figure 2 example is a complex but not unusual structure, comprising a hipped roof with a porch and dormers. Its elevation histogram should slope gently straight down (in the shape of a trapezium [44]) but in reality is concave (see dotted and solid lines in Figure 2 top right). The building in question is known to have a roof tilt of 38° (the same for each of the two major front and rear planes). Actually, it is barely possible to distinguish the 35-40 degree bin as the most frequently occurring in the tilt histogram. (This is obvious on a simpler roof form.) The azimuth histogram is a little clearer in as far as the major front (south-facing, 180°) and rear (north-facing, 360°) planes may be discerned. Then again, the east (90°), and west (270°) planes cannot be extricated from random noise in the data. Probably because of these kinds of problems, the peak-fitting method seems to have been largely replaced by other techniques. Furthermore, it analyses each height point in isolation. The spatial relationship between points is not considered, although all the points on a plane will have related values until an edge is reached. Newer methods refine peak detection with iterative voting (e.g., region-growing [45], random sample consensus algorithm (RANSAC) [46], and Hough Transform) [47]. These are all region dependent, i.e., they account for the spatial location of each height point with reference to its neighbour. The Hough Transform is the computationally fastest of these techniques.
Hough plane detection has two stages: edge detection, followed by grouping of the points inside the edges to generate the planes.

Edge Detection
Initially, a Canny edge detector [48] was applied to the 1 m LiDAR data for the Nottingham house in Figure 2 (GRASS software (Version 6.4.3); i.edge [49]). The Canny edge detector is well known and often used to process both LiDAR data and images. It works by marking local maxima in the LiDAR as edges. However, in the case of the Nottingham building in Figure 2, the Canny algorithm completely failed to discriminate any edges (roof ridges), due to noise and the relatively coarse resolution of the data. When tested on several of the smaller housing association properties, edges were detected but not all correctly ( Figure 3). On some homes the roof ridge is identified but on others an edge perpendicular to the expected position is located. With no clear or consistent pattern to these errors, further algorithms were trialled with the aim of improving reliability. These included the simple (moving window) filter of System for Automated Geoscientific Analyses (SAGA) GIS [50] and ArcGIS [51,52] low pass (3 × 3 cell area mean), majority (3 × 3 cell area mode), and high pass (3 × 3 cell area weighted) filters. There was no improvement in results. Roof ridges appeared too wide or were not detected.
The problem appears to be the resolution of the input LiDAR data. Roof features are too small to be easily perceptible in 1 m data. Figure 4 illustrates the LiDAR data for the example Nottingham house in the form of a simple graph. The larger the circle, the higher the roof elevation of the 1 m grid cell it represents. As may be seen, even with manual intervention, not all roof features are visible in 1 m LiDAR. Higher resolution LiDAR is only publicly available for small areas of the UK and not for the Nottingham test area. For this reason, tests were carried out with aerial photography instead of LiDAR.

Edge Detection Using Google Earth
Images captured from Google Earth were utilised because they are readily available and cover all areas. Several filters available in GIMP software were investigated [53,54], including the low pass, Sobel (horizontal and vertical moving windows), and Laplace (high pass). The basis of all of them is gradient calculation, with edges being defined when a threshold value is exceeded, similar to the Canny edge detector. The best results were achieved with the Laplace filter preceded by a 10 pixel blur to prevent false edges ( Figure 5).
When a wider area was investigated (the Wollaton Park suburb surrounding the example house), it became obvious that only two planes of four-plane roofs were being identified ( Figure 6). This image shows the south and west plane as one, and the north and east plane as one, for four-plane roofs. This is possibly due to the aerial photograph being taken in the afternoon and the filter merely distinguishing the sunny/less sunny sides of the roof. The next step was to investigate image recognition techniques.

Image Recognition
Initially, an unsupervised technique (i.e., image classification without the analyst's intervention) [55] was tried on the example Nottingham house. The ArcGIS software automatically groups image pixels with similar values into statistically distinct classes using iterative clustering around the mean (iso cluster algorithm). In this instance, the outcome was unusable. Almost every pixel in the image was treated as a separate roof plane, the exception being areas of shade which were well distinguished because they are much darker than the rest of the roof. Several supervised methods were then tested. That is, training areas representative of separate roof planes were created by manually digitising polygons. Next, these training samples were used to categorise all other pixels in the image via a classification algorithm. Training examples were digitised for all directions of the example Nottingham roof which may be identified manually (north, east, south, and west). Areas of shade on the roof and chimneys were also digitised for recognition as separate features. Two classification methods delivered reasonable results (Figure 7):


Class probability which employs Bayesian statistics to segment the image.  Maximum likelihood classification, which also uses Bayes theorem but weights classes if they are more likely to occur.
It may be seen that again there are problems with spurious features being identified. Having said that, the main difficulty is that the north plane is hard to distinguish from the east, and the south plane cannot be separated from the west. This is similar to findings from edge detection of Google Earth images. It appears that what is needed for accurate roof plane segmentation is aerial photographs taken at different times of the day. These images would show different roof directions in slightly different colours, as the sun lights each one in sequence on its daily path. A composite from several images would then deliver accurate results. However, multiple daily photographs are not available from Google Earth or any other freely available image source. Therefore the decision was taken to re-examine LiDAR as a data source.
(N.B. Google Earth's Voyager 3D Cities layer [56] is generated from multiple sources, including Sketch-up models and stereoscopic imagery. There is no 3D geometry currently accessible for download, which eliminates this resource at present.)

Hill Shading with Ambient Occlusion
The previous sections discovered a need for images captured at successive times during the day. This was achieved by applying the hill shading with ambient occlusion module from SAGA software [14,[57][58][59] to 1 m LiDAR data for the residential area in Nottingham. Hill shading models beam radiation from a single direction. Ambient occlusion adds the diffuse component of sunlight. It samples a hemisphere around each LiDAR height point and ascertains what proportion of that hemisphere is blocked by higher surrounding points. The pixel is shaded to suit. The combined technique was used to generate shading patterns on roofs at different times during the day. Preliminary results appear encouraging (Figure 8). North-and west-facing roof planes are shaded (darker) in the morning simulation (sun in southeast), north only at 2 pm, and north-and east-facing segments are in shadow slightly later in the afternoon (sun in southwest). Nonetheless, this technique has some shortcomings. It does not allow for beam reflection and transmission, e.g., through thin cloud, and therefore is not completely realistic. In addition, it is slow [58].

Review of Progress
All the techniques covered so far are based on grouping the unique values allocated to each LiDAR grid height point or Google Earth image pixel colouration to produce realistic roof planes. That is, roof segmentation traditionally precedes estimation of solar potential on building roofs. This may be the standard approach but, as illustrated above, there are many difficulties, summarised in Table 1. None of the above methods works well with the data resolution available in the UK (1 m for the most part). Existing methods cannot deal with the variety of roof styles, data noise, and temporal and spatial resolution of available data. Identification of suitable rooftops for PV is critical for research ranging from grid impacts to carbon emissions. No agreed method of obtaining this value has yet been published [30]. The following sections propose an alternative to conventional rooftop PV techniques. An accurate, sound rooftop PV assessment tool, applicable to many cities, is provided. Generates shading patterns at different times of day. Most promising of these methods but not completely realistic.

New Method
This method simply asks "Can PV be fitted to this roof?" There are no predefined templates, merely three basic technical requirements: (1) sloping roof with a pitch of between 15° and 60° (flat roofs are treated separately); (2) orientation in the southern half of the compass; and (3) roof area large enough to accommodate a 1 kW system, i.e., 8 m 2 roof plane. Figure 9 is a schematic representation of the stages involved in this new method. ArcGIS commands are provided in the Appendix A.
1. Select building elevation pixels only from the LiDAR grid, using OS Mastermap Topography layer [39] as a "pattern cut-out". Heights of other features e.g., vegetation, cars, salt bins, are eliminated ( Figure 9, Step 1). 2. ArcGIS Slope and Aspect tool are used to calculate the tilt and azimuth of each roof grid square ( Figure 9, Step 2). 3. The mean tilt of the entire roof is calculated. All main roof patch areas are assigned the same tilt.
Accuracy is maximised by taking the entire roof. Roofs with tilts of less than 15° and greater than 60° are excluded. Roofs with a pitch smaller than 15° are not steep enough for accurate azimuth estimations and are categorised as "flat". 4. Determine the mode of azimuth for all grid squares within each building footprint. 5. The initial assumption is that all houses have two major aspects. The area is manually checked using Google Earth and extra steps are carried out if more complex building forms are detected: a. Duo pitch roof: Considering due North to be zero degrees, if the mode is between 90° and 270°, add 180° to obtain the southern aspect for PV. In theory an "up and over roof" has two aspect "modes" but due to LiDAR flight paths, reflections, etc. one will predominate. Every pitched roof must have at least two diametric planes. (In contrast to standard peak finding, only one azimuth peak must be found, rather than two.) b. Square four hip roof: If the mode is between 90° and 180°, then add 90°. If the quantity of west-facing grid squares is greater than the amount of east-facing pixels, select the west-facing ones. These receive more insolation and cost will probably prohibit installation on both aspects. c. Three plane houses: Similar to four, but one direction will be missing. If no pixel values are detected in a 10° band when the mode is replaced with the 90° opposite value, this orientation is missing and process is abandoned. d. More than four main roof patches-complex roof formats are not suitable for PV and therefore not covered by this research.
6. Select roofs east through south to west (southern cardinal compass directions only). 7. Chose grid squares not outside half a standard deviation of the mode (Figure 9, Step 3). 8. Carry out a Rook's Case connectivity check to delete roof patches which are linked by the corners (touch diagonally). PV modules cannot be connected cornerwise (Figure 9, Step 5). 9. Apply a minimum 10 pixel filter (10 of 1 × 1 m grid squares) to the selected grid squares to eliminate small areas (Figure 9, Step 6). 10. Remove dangling pixels with a boundary clean. 11. Roof patch area may be calculated (see Appendix A) if desired, although this is not necessary. It goes beyond the main objective of this approach, which is to select roof planes which meet the three basic technical requirements for PV.
Note, for speed or in very large areas, the default of two plane roofs may be accepted. This is the most common roof type for houses of all ages (see photographs by Barrow [60]). Gabled (two plane) roofs are also found on terraced houses which comprise large areas of industrial cities. Minimum solar power system area was obtained using the roof direction rather than the pitch. This proved less variable than the tilt, making it easier to group grid square values around a statistical measure. By investigating ten houses where azimuth had been measured, the mode was found to be the most effective measure for categorizing pixels (rather than mean, maximum, etc.). There is less distortion from errors in the original LiDAR data.
In order to group roof pixels into areas which may be checked for minimum PV size requirements, the following mathematical formulae were tested. These all select azimuth pixels around the mode:


Equal interval plus or minus 45°.

 Jenks Natural Breaks [61] 
One quarter standard deviation of mode. This gathers almost one sixth of roof data (68% std/4 = 17%).  One-third standard deviation of mode. This gathers nearly a quarter of roof data (68% std/3 = 23%.  Half standard deviation of mode. This obtains about one-third of roof data (68% std/2).
These five classification methods were tested on a database of solar photovoltaic installation data maintained by a housing association (see Section 5). Eight-hundred-and-eighty-six of the homes have LiDAR coverage. System sizes are catalogued, so the area of individual systems may be estimated (1 kW = 8 m 2 ). The cosine rule (see Appendix A) was used to transform horizontal roof plane area to tilt area before checking for the minimum of 8 m 2 for PV. Results from the half standard deviation method most closely matched database values. Only 2.5% of housing association homes fitted with PV were missed. The other four classification methods missed about twice as many houses. Visual checking of the structurally more diverse houses in the Wollaton Park case study with Google Earth also found the half standard deviation method to be most accurate. Figure 9. Method to discover whether roofs are suitable for minimum size PV installation.
No specific software is required for this method, other than GIS. Processing speed is acceptable with a standard desktop PC and satisfactory results are possible with medium resolution LiDAR. The process relies on data processing. Automation is possible, but not essential. Sample results for the Wollaton Park area of Nottingham are illustrated in Figure 10. Forty roofs are identified as suitable for PV. Some complex roofs are wrongly identified in the top right of the image. These are inappropriate for PV installation because of dormers and cross gables. However, compared to the methods detailed in Section 3, this new approach performs very well. The next section describes its validation in detail. The new approach is, of course, limited by the availability of LiDAR data. It does need some manual checking in areas with complex architectures. Economic, social and market factors of PV uptake, competition with other renewable energies and rate of adoption are not investigated. On the other hand, it is low-cost, depending only on publicly available data. It is suitable for small buildings as well as medium and large, unlike some methods which tend to perform best on big buildings. The half-standard-deviation-of-mode method acts as a reference point in that it estimates minimum regional domestic PV deployment.

Results and Validation
Validation of the new method took place with a housing association database of PV installation data (Nottingham, UK). It was possible to obtain address (and therefore latitude/longitude), system size, and installers' values of tilts and azimuths for these systems. Eight-hundred-and-eighty-six of the homes have LiDAR coverage, so it is possible to compare modelled results to actual on-the-ground measurements.
The results are summarised in Table 2 and illustrated in Figures 11 and 12. Figure 11 graphs the percentage of actual systems within 5 degree bins of the modelled value of tilt/azimuth. Figure 12 charts the under/overestimation of roof plane size.

Tilt
Sixty-one per cent of tilts obtained with the half standard deviation of mode method were in 5° of the installation database (Table 2). Eighty-seven per cent of estimated tilts were in 10° of the database (Figure 11). The Mean Bias Error (MBE) is 4°, the Root Mean Square Error (RMSE) is 7°, and the most frequently occurring error is 4°. Houses with the same configuration vary by 3° [13], so these results are acceptable. In addition, the installers' figures are thought to be "rule of thumb" and not measured e.g., by inclinometer. A 10° tilt variation between the traditional UK roof pitches of 40-50° will only make a 1% difference to average annual plane-of-array irradiation received [15].

Azimuth
Thirty-three per cent of azimuths calculated with the new method were in 5° of the installation database (Table 2). Sixty-six per cent of estimated azimuths were in 15° of the database. All the estimated azimuth figures were in 45° of the recorded value ( Figure 11). Forty-five degrees azimuth variance impacts plane-of-array irradiation by 15%. The MBE is 14%, RMSE is 18%, and the most frequently occurring error 5°. Again, these figures are considered to be satisfactory. Tolerable results have been achieved despite the fact that difficulties were noted with the housing association dataset. Visual checks using Google Earth discovered cases where the LiDAR derived figure is correct and the installers' value is not. Tilts and azimuths appear to have been transposed in the database in some instances. Additionally, the installers appear to have estimated azimuth by the position of the sun without allowing for its annual path.

Roof Patch Area
The vast majority (97.5%) of established systems used in the validation process were correctly identified as being suitable for at least a minimum potential 1 kW system ( Figure 12). Eighty per cent had an area at least the size of the actual installed system. Comparing LiDAR derived values and values calculated from the system sizes, the MBE is 6% and RMSE 9%.   Table 2 and Figures 11 and 12 show that the majority of the modelled tilt and azimuth values fall within 15 to 20 degrees (6%) of the installers' values. However, the spread of size value differences is much wider. Only 12% of results are within 10% of the ground-measured value. Therefore this research concentrates on identifying suitability for a minimum size PV system only, with 98% accuracy.

Context of New Method
Comparable techniques have been developed by Boz et al. [32] and the National Renewable Energy Laboratory (NREL) [12,62,63]. This section discusses the similarities and differences between previous work and the current method.
Boz et al. [32] bin individual pixel values from LiDAR into seven slope (tilt) and five aspect (azimuth) classes. The results are smoothed by replacing pixel values based on the majority of a 3 × 3 neighbourhood to the present pixel (Majority Filter). Prior classification means that the fine detail of the original input LiDAR is lost during the roof segmentation. This work is validated visually by matching 150 rooftops in Philadelphia to aerial imagery.
Similar to Boz et al., NREL categorize each 1 × 1 m roof square into nine azimuth classes. For each distinct roof plane creating by azimuth classification, the mean tilt was determined (Zonal Mean). A PV installer's data set containing the location, tilt, and azimuth of 205 assembled PV arrays was used to validate the results of this analysis.
Boz et al. report that their method gives the most precise results when applied to simple roof structures. NREL's technique has virtually the same accuracy as the half standard deviation of mode method. Eighty-nine per cent of NREL modelled results were within 10° of the actual slope compared to 87% obtained by the current authors. Ninety-six per cent of NREL's modelled results have the same azimuth as the actual azimuth set against the 100% accuracy obtained by the current authors (allowing for categorization into compass bands to be compatible with NREL).
The techniques of Boz et al. (prior classification with Majority Filter) and NREL (single tilt value for each unique roof plane obtained from the azimuth) were tested with UK data. It was found that the Majority Filter did not add any accuracy. Straightforward categorization into compass band as previously carried out by the current authors [13] gives greater accuracy in the UK. This is due to the low number of LiDAR pixels bounded by the footprint of a typical UK home. Likewise, calculating the tilt value for every roof plane generated flawed results because segmenting small buildings gives imprecise results due to scarcity of data points.
The methods of Boz et al. and NREL are reported as working well for the larger homes of the US. The method presented here (entire building average for tilt and half standard deviation of mode for azimuth planes) is suitable for the smaller houses of the UK and other European and Asian countries. It is validated against more actual buildings' data than any previous method.

Summary and Discussion
The new method presented here is a competent investigative tool for the installation of residential solar systems. It is capable of estimating how much energy could be generated if a minimum size PV system was installed on all existing suitable roofs throughout an entire city, region, or country. This information is valuable to network operators and can be used for transmission planning.
This method has few prior assumptions. It just determines suitability for PV. The results represent a minimum boundary for potential PV deployment. There is no attempt to predict actual installation because this is influenced by a range of factors, e.g., affordability, attractiveness, and tenure.
The new method chooses grid squares within half a standard deviation of the azimuth mode. It categorises noisy data values to produce intelligible results. A map of roofs suitable for minimum size solar system installation is produced: size 8 m 2 or more and tilt 15-60°, southern compass direction.
A comprehensive validation has been carried out using two techniques. Firstly, by a check for wrongfully selecting inappropriate roofs as suitable for PV. This was carried out by manually matching 50 rooftops in Wollaton Park to Google Earth imagery ( Figure 10). Two roofs were incorrectly selected as opposed to 50 correctly categorised. Secondly, by a check for missing suitable roofs by comparison against the biggest installation database used by any analogous research to date.
The tool presented here can be used to calculate the ability of aggregations of households at various scales to offset their electricity consumption with PV. The majority of people live in cities, even in less developed countries. This proportion is set to increase to 68% by 2050. Overall population growth will add even more people to the urban environment. The greatest increase in city dwellers is expected in Africa and Asia [64]. Successful development will require sustainable urbanisation. Rooftop PV and accurate estimations of its yield have a significant role in urban expansion. The new method is suitable for both industrial and residential buildings, since it is capable of analysing all building sizes, and flat as well as tilted roofs. The detailed information it may provide allows for a more flexible grid and related infrastructure to accommodate renewable technologies.

Conclusions
This method is useful, effective and functions correctly with the data publicly available in the UK (predominantly 1 m LiDAR, Google Earth images and MasterMap building outlines, Section 2). In this respect, it provides a valuable contribution to the scientific field because the methods tested and reviewed in Section 3 require higher resolution input data than can be provided to produce usable results.
The unique attribute of the method presented here is that it is twice validated, by extensive ground-truthing against a database of 886 installations, and against aerial photography. This makes it the most thoroughly validated method to date.
The aim of any PV roof-area estimation method is to provide data for further analysis. This method is flexible. It allows individual houses as well as large numbers of properties to be examined, depending upon later requirements. Unlike some methods, it makes no assumptions when larger numbers are involved and does not rely on compass band classification. An example of use of the data generated would be a study of the relationship between azimuth and self-consumption. A range of PV-related research is enabled.
This method may be applied to estimate minimum rooftop area suitable for PV deployment in any country which has northern European-type architecture and/or equivalent building size, provided LiDAR data is available. In countries with larger house sizes, e.g., Australia, USA, and Canada, it may be better to use the NREL method [12]. For cities with predominantly flat roofs, e.g., Cairo, building footprint analysis will suffice. Following implementation of this method, photovoltaic yield may then be calculated using a free or commercial tool [6][7][8][9][10] or programming a thermal and an electrical model [65].
The knowledge of rooftop PV potential provided by this publicly available data and methods may assist planners, decision-makers, and network operators in analysing PV's facility to meet electricity demand. Future electricity generation portfolios in specific countries may also be modelled.
Currently the model provides information on minimum potential PV deployment. It is planned to identify or develop a technique to estimate maximum deployment, e.g. Boz et al. [32] or the NREL method [12,62,63]. These two sets of information will supply the range of possible solar installation scenarios. Future developments of the new method provided here will also include incorporation of meteorological data and performance models to deliver energy yield prediction, as commenced in [14,66]. It then becomes possible to investigate times and quantities of potential output and feed-in behaviour to the transmission system. The mass effect of many small-scale systems has, to date, mostly been studied using generalised data.

Conflicts of Interest:
The authors declare no conflicts of interest. The project funders were not involved directly in writing this article.

Appendix A
ArcGIS Commands for method to select suitable roofs for minimum size PV installation: 1. Preparation: Subtract Environment Agency digital terrain (DTM) LiDAR from the digital surface model (DSM) to obtain building height above ground level. Delete pixels with a lower height than 2 m. This removes rogue values whilst allowing for low eaves. 2. Cut out buildings only. Prepared LiDAR DSM-DTM grid and buildings from OS Mastermap Topography layer as inputs: ArcToolbox > Spatial Analyst > Extraction > Extract by Mask. 3. Spatial Analyst > Surface > Slope/Aspect. Create integer rasters to enable mode to be computed in the next step. 4. Spatial Analyst Tools > Zonal > Zonal Statistics > Mean 5. Spatial Analyst Tools > Zonal > Zonal Statistics > Mode 6. Mode may be in the north, so carry out some swaps: a. Swap mode raster values by 180: Raster calculator: i. Values between 90 and 270 are alright: Con(("aspectmode" ≤ 270) & ("aspectmode" ≥ 90), "aspectmode", 0) ii.