Overlap and Segregation among Multiple 3D Home Ranges: A Non-Pairwise Metric with Demonstrative Application to the Lesser Kestrel Falco naumanni

Simple Summary Modeling animal space use in 3D is more realistic than confining research to 2D methods and can greatly increase our understanding of interspecific and intraspecific competition, predator–prey relationships, habitat selection and use. In particular, home range overlap/segregation is a fundamental property of animal interactions with deep implications for biodiversity conservation and management. In order to solve the issue of measuring the degree of overlap/segregation among an arbitrarily large number of 3D volumetric home ranges, we introduced the novel non-pairwise index MVOI (Multiple Volumetric Overlap Index) and its complement to 100 MVSI (Multiple Volumetric Segregation Index). Results show that traditional 2D spatial analyses can significantly overestimate the overlap between the individuals, population and species that occupy a habitat with a strong vertical component. Both the MVOI and MVSI can also be used to quantify how the 3D home ranges change over time (i.e., 4D home ranges) and the robustness of 3D home range assessment through the degree of overlap among different 3D estimators. We applied the MVOI and MVSI to birds, but they can be readily applied to any animal species, in particular those with a significant vertical component to their space use. Abstract In this study we solved the issue of measuring the degree of overlap/segregation among an arbitrarily large number (n ≥ 2) of 3D volumetric home ranges (i.e., x, y, and hg; where hg is height above ground level) for the first time. For this purpose, we introduced the novel non-pairwise index MVOI (Multiple Volumetric Overlap Index) and its complement to 100 MVSI (Multiple Volumetric Segregation Index). Regardless of the number of 3D volumetric home ranges, the MVOI and MVSI generate a single score of overlap/segregation between 0 and 100, making ecological interpretation much easier and more meaningful when compared to n × n pairwise overlap indices. As a case study, we applied the MVOI and MVSI to 12,081 GPS points of five lesser kestrels (Falco naumanni) during the nesting period at Santeramo in Colle (Apulia region; Italy) in an area with the most elevated density of lesser kestrels in urban colonies worldwide. The 3D volumetric home ranges ranged between 1.79 km3 and 8.19 km3. We found that the tracked birds had different vertical profiles, possibly to limit intraspecific competition, resulting in a 3D home range overlap that was only 61.1% of the 2D overlap and 52.8% of the probabilistic one.

One main measure of space use is the home range (HR hereafter), i.e., the space used by an individual to obtain resources and meet its requirements for survival and Biology 2023, 12, 77 2 of 12 reproduction [7]. Although the HR has been described for many taxa almost exclusively in two dimensions (x, y), for species that have a strong vertical component (z) to their movement (i.e., height above ground level for flying species, depth for fish species and elevation for species moving in highly variable terrain), such two-dimensional (2D) representation neglects some key components of their ecology, such as the actual size of their HRs and overlap between the HRs. Several authors [8,9] showed that a 2D representation can underestimate the surface area of HRs relative to three-dimensional (3D) estimates for animals that occupy habitat with a strong vertical component. The interspecific and intraspecific vertical stratification of the HR has been well documented in birds [9]; however, these studies quantified space use by comparing heights directly and ignored the other two dimensions, whereas others created 2D utilization distributions and then separately analysed vertical data. Regrettably, by assuming that individuals utilize the same vertical space, traditional 2D spatial analyses computed through the horizontal plane can likely overestimate the overlap between individuals that stratify vertically [10].
To date, very few studies examined the multidimensional space use of animals, arguably because of the lack of tools adapted to 3D data. Recently, however, two methods have been proposed to estimate 3D HRs. A study by [8] first stressed the inadequacy of existing modelling techniques to draw advantages from three-dimensional datasets and proposed the use of 3D kernel density estimators [11] to calculate 3D probabilistic HRs (i.e., x, y, z and p; where p measures the probability density that an animal is found at a given point within a certain 3D space). A study by [12] first proposed a 3D volumetric (i.e., x, y and h g ; where h g is the height above ground level) HR estimator based on initial 2D HR estimation, followed by square (or hexagon) tessellation and 3D extrusion that assembled the 3D volumetric HR in the form of n adjacent parallelepipeds with different heights above the topographic surface.
If the issue of calculating 3D HRs has been adequately approached to date, the same is not true for their overlap. As yet, only studies concerning the pairwise overlap of 3D utilization distributions (i.e., 3D probabilistic HRs) have been proposed. A study by [13] first suggested a methodology to estimate the pairwise (i.e., n = 2) overlap between 3D probabilistic HRs by the adaptation of 2D overlap indices. A study by [14] solved the issue of computing the non-pairwise overlap of an arbitrary number (i.e., n ≥ 2) of probabilistic HRs for the first time, but the solution was limited to 2D utilization distributions (i.e., x, y and p; where p measures the probability density that an animal is found at a given point within a certain 2D space). Although 2D probabilistic HRs can be considered as a particular case of 3D volumetric HRs (HR 3D hereafter), where z = p, they pose two strong simplifications for overlap estimation when compared with volumetric ones. Firstly, in 2D probabilistic HRs the sum of p values (that corresponds to the volume under the probabilistic surface) is equal to 1 (or 100%) by definition [15]. By contrast, in the HR 3D the volume is measured in cubic kilometres and can assume any possible value ≥ 0 [12]. Secondly, in the 2D probabilistic HRs, the lower reference surface is simply a two-dimensional plane (x, y) where animal locations are placed, which can be operatively considered as a flat surface at z = 0 m a.s.l. Instead, HR 3D estimation raise the challenging issue of taking into account the topographic surface over which the HR 3D lays. In fact, the lower reference surface here is a (sometimes very complex) topography where each pair of coordinates is associated to a particular topographic height a.s.l. that can assume any possible value, even ≤0. These differences make the index developed for multiple overlap of 2D probabilistic HRs [14] useless for the HR 3D .
Accordingly, in this study we solved the issue of measuring the degree of overlap/segregation among an arbitrarily large number (i.e., n ≥ 2) of HR 3D for the first time.
To this aim, we upgraded the general overlap index [16] that was recently introduced for the computation of non-pairwise 2D HR overlap among n individuals, populations or species. Using spatial location data from the lesser kestrels Falco naumanni in southern Italy, we demonstrated the properties of our new index and discussed its potential applications to behavioural and movement ecology.

Study Species and Study Area
The study species is a small raptor present among the Annex I species of EU Wild Birds Directive 2009/147/EEC. This species relies on natural cliffs or rural and urban buildings as nesting sites and breeds in steppe-like grasslands and non-irrigated arable crops where it forages predominantly on invertebrates (in particular, grasshoppers [17]). The study area corresponded to the lesser kestrel colony of Santeramo in Colle (Apulia region, southern Italy; Figure 1), i.e., an agricultural landscape with an elevation gradient from 348 to 509 m a.s.l., located within the SPA (Special Protection Area) "Murgia Alta" IT9120007 and included within the Important Bird Area "Murge". This area has the most elevated density of lesser kestrels in urban colonies worldwide [18]. The lesser kestrel population of Santeramo in Colle is characterized by an elevated within-colony overlap for both 2D [16] and probabilistic [14] HRs and elevated segregation with the adjacent lesser kestrel colonies [16]. Italy, we demonstrated the properties of our new index and discussed its potential applications to behavioural and movement ecology.

Study Species and Study Area
The study species is a small raptor present among the Annex I species of EU Wild Birds Directive 2009/147/EEC. This species relies on natural cliffs or rural and urban buildings as nesting sites and breeds in steppe-like grasslands and non-irrigated arable crops where it forages predominantly on invertebrates (in particular, grasshoppers [17]). The study area corresponded to the lesser kestrel colony of Santeramo in Colle (Apulia region, southern Italy; Figure 1), i.e., an agricultural landscape with an elevation gradient from 348 to 509 m a.s.l., located within the SPA (Special Protection Area) "Murgia Alta" IT9120007 and included within the Important Bird Area "Murge". This area has the most elevated density of lesser kestrels in urban colonies worldwide [18]. The lesser kestrel population of Santeramo in Colle is characterized by an elevated within-colony overlap for both 2D [16] and probabilistic [14] HRs and elevated segregation with the adjacent lesser kestrel colonies [16].

Tagging of Birds
We tracked five lesser kestrels between 13 and 29 June 2017 during the nesting period (Table 1). We fitted the lesser kestrels with data loggers at their nest boxes. We employed TechnoSmart GiPSy-4 and GiPSy-5 data loggers (23 mm × 15 mm × 6 mm, <5 g weight) in order to collect spatio-temporal information (date, time, latitude, longitude, height a.s.l. and speed). The weight of the loggers in relation to that of the tracked individuals was <4%. All devices were tied dorsally using a 2 mm large Teflon tape knotted with a triple knot, and two tapes were crossed without a knot at the height of the sternum. We

Tagging of Birds
We tracked five lesser kestrels between 13 and 29 June 2017 during the nesting period (Table 1). We fitted the lesser kestrels with data loggers at their nest boxes. We employed TechnoSmart GiPSy-4 and GiPSy-5 data loggers (23 mm × 15 mm × 6 mm, <5 g weight) in order to collect spatio-temporal information (date, time, latitude, longitude, height a.s.l. and speed). The weight of the loggers in relation to that of the tracked individuals was <4%. All devices were tied dorsally using a 2 mm large Teflon tape knotted with a triple knot, and two tapes were crossed without a knot at the height of the sternum. We downloaded the data from the data loggers after the birds were recaptured at their nest boxes. Data acquisition occurred every three minutes following deployment.

Data Preparation
About 93% of the GPS points (i.e., 12,081 out of 12,984 points) were linked to 6-8 satellites. In these favourable conditions, the error associated to the measurement of the height a.s.l. (h s ) was in the order of ±4-6 m [19]. We discarded the remaining 903 GPS points, linked to 5 or fewer satellites, because the measurement error associated to h s was in the order of ±20-25 m [19]. For each individual we also excluded the GPS points with h s greater than the 95th percentile (critical threshold) of the frequency distribution, so as to remove vertical outliers (possibly due to anomalous fine-scale weather conditions or GPS issues, e.g., atmospheric interference, low battery, multi-path effects) that could have determined an overestimation of the HR 3D . GPS data were transferred to the ArcView GIS [20] and superimposed on the topographic surface of the study area ( Figure 1) digitized at a 1:2000 scale by the authors from the available topographic maps of Apulia Region. For each GPS point, the height above ground level (in meters) was calculated as where h t is topographic height a.s.l. In order to make pairwise comparisons between the statistical distributions of bird heights above ground level, we used a two-sample z-test. Given the 5 individuals, we performed 5 × (5 − 1)/2 = 10 pairwise tests that were considered to be significant for (two sided) p < 0.05.

3D Volumetric Home Ranges
Within GIS, we calculated the HR 3D by using the methodology developed by [12]. Firstly, we computed the 2D polygonal HR of each individual. We employed the minimum convex polygon algorithm [21] with fixed arithmetic mean of all x (longitude) and y (latitude) coordinates and retained 95% of points closest to the arithmetic mean point. The 95% isopleth is widely used in the literature to exclude possible horizontal (i.e., in the X-Y plane) outliers due to fine-scale weather conditions or errors affecting GPS accuracy [22]. Secondly, we computed the smallest possible convex polygon (SCP) around the detected 2D HRs so as to delimit the study area over which the HR 3D and overlaps had to be calculated. Thirdly, we tessellated the SCP by using a regular square grid where the square size was determined through a trial-and-error computation. We used a set of possible candidates ranging from 1 hectare (i.e., 0.01 km 2 ) to SCP incremented by 0.1 hectares, and the square size corresponded to the minimum size so that at least 70% of the squares contained at least 1 GPS point inside (i.e., k ≥ 1).
Fourthly, for each individual we assigned the maximum value of h g (i.e., h g max ) to each square. In case of squares with no GPS points inside, we set h g max = 0 m a.g.l. In mathematical terms, Fifthly, for each individual we extruded each square upon the topographic surface by an elevation equal to h g max . This step assembled a polyhedron in the form of n adjacent parallelepipeds above the topographic surface, therefore the HR 3D of each lesser kestrel was Lastly, for each lesser kestrel, the volume of the HR 3D was calculated as in case ∆x ∼ = 0 and ∆y ∼ = 0. Although the 2D HR of each lesser kestrel was a sub-region of the SCP, we used the SCP as the domain of Equations (3) and (4) because h g max was equal to 0 outside the 2D HR, therefore V 3D assumed the same value over these two domains.

3D Home Range Overlap
We upgraded the recently introduced general overlap index (GOI; [16]) that allows for the computation of non-pairwise overlap among an arbitrarily large number (n ≥ 2) of 2D HRs. The GOI follows a simple idea: given n 2D HRs, it is always possible to calculate the extent of two spatial configurations, perfect segregation and perfect overlap. The GOI simply measures the distance of the observed overlaps from these two extremes. In the case of perfectly disjointed (i.e., perfect segregation) 2D HRs, the total area (T A ) covered by the 2D HRs is simply the sum of their extents ∑ A i . In the case of perfectly nested (i.e., perfect overlap) 2D HRs, T A is the extent of the largest HR (i.e., max(A i )). In the intermediate case (i.e., partially overlapping HRs), T A corresponds to the spatial union of the HR polygons, ∪A i . The larger the overlap, the smaller the ∪A i . The GOI is simply the ratio between the observed (Dist OBS ) and maximum (Dist MAX ) distances from the perfectly segregated (i.e., non-overlapping) situation, calculated as If Dist OBS = 0 (i.e., perfect non-overlap), then the GOI = 0; if Dist OBS = Dist MAX (i.e., perfect overlap), then the GOI = 100. In the intermediate cases, then 0 < GOI < 100. A general segregation index (GSI) can be computed as the complement to 100 of the GOI. A geometrical elucidation of these two indices can be found in [16].
Both the GOI and GSI only consider the 2D spatial domain of the individual HRs and ignore the third dimension (i.e., h s , h g and h t ) associated to each GPS points. In terms of volumetric HRs, in case of perfect segregation ∑ A i becomes the sum of the HR 3D .
In case of perfect overlap, max(A i ) simply becomes the largest HR 3D In the intermediate case (partially overlapping HR3 D ), ∪A i corresponds to the spatial union of the HR 3D , i.e., the volumetric HR where each grid square assumes the maximum value among all the h g max of the HR 3D , calculated as By inserting Equations (6)-(8) into Equation (5), the non-pairwise overlap among an arbitrarily large number (n ≥ 2) of HR 3D reads as In case ∆x ∼ = 0 and ∆y ∼ = 0, the MVOI must be calculated by using double integrals as The MVOI can also be calculated by directly using h s and h t as follows Equations (10) and (11) return the same result, but Equation (10) is probably more straightforward. Finally, a volumetric general segregation index (MVSI) can be computed as

Results
The five 2D HRs ranged from 22.76 to 111.48 km 2 , with lower values for the female lesser kestrels. The SCP was 142.4 km 2 ( Figure 2). After several trial-and-error attempts, the square size was set to 0.09 km 2 (i.e., ∆x = ∆y = 300 m); therefore, the SCP was covered with 1584 squares.   Table 1. The tracked lesser kestrels showed dissimilar vertical profiles ( Figure 3). All the pairwise differences between the statistical distributions of bird heights above ground level were significant, except for the two pairs of individuals F18-M18 and F24-M18 (Table 2).   Table 1. 9.55 ** M24 * p < 0.05, ** p < 0.01.
The median hg of the tracked individuals ranged from 16 to 29 m a.g.l., while the 95th percentiles of hg were between 283 and 448 m a.g.l. (Table 3). The female F18 presented  Table 1. The median h g of the tracked individuals ranged from 16 to 29 m a.g.l., while the 95th percentiles of h g were between 283 and 448 m a.g.l. (Table 3). The female F18 presented many GPS points with h g > 150 m and up to 250 m a.g.l., though the most frequent h g were those around the median (20 m a.g.l.). The female F24 preferred lower h g (the 25th percentile of h g was close to 0 m a.g.l.), and the probability density of h g decreased constantly as h g increased and was almost null for h g > 100 m a.g.l. The male M4 showed two peaks of the probability density of h g , the first peak close to 0 m a.g.l. and the second peak close to the median (29 m a.g.l.). In addition, it showed an elevated amount of GPS points up to 200 m a.g.l. The male lesser kestrel M18 preferred lower h g (the 25th percentile of h g was close to 0 m a.g.l.), but it showed a fair amount of GPS points up to 200 m a.g.l. The male lesser kestrel M24 exhibited the highest values of h g for both the 75th (71 m a.g.l.) and 95th percentile (448 m a.g.l.).
The five HR 3D s ranged between 1.79 km 3 and 8.19 km 3 (Table 3). Max(V 3D ) and ΣV 3D were 8.19 km 3 and 21.89 km 3 , respectively. The volume of ∪V 3D was 15.08 km 3 (Figure 4)  200 m a.g.l. The male lesser kestrel M18 preferred lower hg (the 25th percentile of hg was close to 0 m a.g.l.), but it showed a fair amount of GPS points up to 200 m a.g.l. The male lesser kestrel M24 exhibited the highest values of hg for both the 75th (71 m a.g.l.) and 95th percentile (448 m a.g.l.). The five HR3Ds ranged between 1.79 km 3 and 8.19 km 3 (Table 3). Max(V3D) and ΣV3D were 8.19 km 3 and 21.89 km 3 , respectively. The volume of 3D V  was 15.08 km 3 ( Figure   4), thus the MVOI was 100 × (21.89 − 15.08)/(21.89 − 8.19) = 49.71% and the MVSI was 100% -49.71% = 50.29%. Table 3. For each lesser kestrel, the 2D and 3D properties of space use are shown. a.g.l. stands for 'above ground level'. IDs are the same as those in Table 1.

Discussion
The issue of estimating the degree of overlap among multiple 3D volumetric HRs has remained unsolved to date. Accordingly, in this study we first introduced a non-pairwise metric of overlap/segregation among multiple HR 3D , whose ecological interpretation is much easier and more meaningful if compared to n × n pairwise overlap matrices computed through pairwise indices. In addition, one overlap index is more effective if estimates of the overlap are to be meaningfully compared across several studies. This fulfils the demand by many authors [23,24] who argued that any overlap index should be intuitive and easy to interpret.
Incorporating the vertical component into the representations of animal space use can provide novel ecological insights with benefits for conservation management [8]. Although 2D analyses are informative about the locations of the tracked individuals, volumetric analyses provide the benefit of combining 3D information into metrics that represent the actual space use of animals [10]. In fact, when animal behaviour includes movements with a substantial vertical component, like in birds, simplifying the assumptions of 2D HRs can affect ecological inferences by overestimating interactions between individuals, population and/or species [25]. By contrast, three-dimensional spatial analyses enable the accurate description of the patterns of spatial overlap/segregation [10]. Results from our analyses support this argument based on the increased amount of detailed information gained from incorporating the vertical dimension with animal location information. In fact, in our previous studies we found that the non-pairwise overlap index (GOI) among the 2D HRs of the lesser kestrels used in this study was 81.38% [16], and the overlap was 94.016% by using a non-pairwise probabilistic HR index (PGOI; probabilistic general overlap index; [14]). In both cases, the degree of overlap among the individuals of this colony was very elevated. In this study, we found that their 3D HR overlap was only 61.1% of the 2D overlap and 52.8% of the probabilistic one. This difference in the estimation of overlap occurred because the 3D analysis allowed separation between individuals that occurred in the same location but showed different vertical profiles above ground level. The lesser kestrel colonies present in the study area (e.g., Altamura, Cassano Murge, Gravina, Santeramo in Colle) are characterized by elevated segregation with the adjacent colonies [16,26]. Because the colony-specific home ranges result in mutually exclusive areas, each colony has only limited space available, which raises intra-colony competition due to increased 2D overlap among individuals. Accordingly, we hypothesize that the detected segregation along the vertical dimension could be an adaptive behaviour employed by the lesser kestrel population of this colony to decrease the elevated intra-colony competition. Alternatively put, we suggest that, as a general rule, the larger the intra-colony 2D overlap is, the larger the 3D segregation will be.
We argue that these three non-pairwise metrics of HR overlap/segregation (i.e., GOI, PGOI, MVOI) quantify the different properties of the animal space use, thus none is exhaustive if considered separately, and instead they should be combined into a vector of three synthetic metrics capable of taking into account all the fundamental properties of animal HR overlap/segregation.

Applications to Tracking Studies
The methodological approach proposed here enabled the description of animal distribution in the same number of dimensions as the environment in which they live, thus providing a realistic description of their space use. The most basic and intuitive application of our HR 3D overlap estimator is to gain ecological insights into the volumetric requirements of birds in order to endure in their distribution areas. For example, the smallest HR 3D found in this study (i.e., female lesser kestrel ID F18; HR 3D = 1.79 km 3 ) involved a 3D space above ground level equal to: (a) 716,000 Olympic-sized swimming pools (i.e., 50 m × 25 m × 2 m), (b) 1356 Colosseums or (c) 737 Great Pyramids of Giza. Instead, the 3D HR overlap indicates where most birds can find the most suitable conditions in the 3D environment, which includes (a) low disturbance due to anthropogenic structures with considerable vertical dimension (e.g., wind farms and power lines [27,28]); (b) favourable vertical thermal gradients (bands of warm and cool air) and air currents (headwinds, tailwinds, crosswinds); and (c) non-negative interactions with other avian species. The application of 3D analyses to the study of biotic interactions could provide the ability to better understand these interactions over current 2D approaches. For example, in the study area, the lesser kestrel is predated by the magpie (Pica pica), the lanner falcon (Falco biarmicus) and the peregrine falcon (Falco peregrinus) that occupy the same 2D space as that of lesser kestrels [29]. The application of our overlap/segregation metrics to the HR 3D of these species could possibly reveal that these predators have different degrees of 3D overlap with the lesser kestrel, and thus the different degree of biotic interactions. Overall, our results suggest that airspace is a habitat with as much potential to influence birds' movements and space use as other environments, and it should not be considered only as something birds move through between land or water habitats.
Although the basic utilization of the MVOI and MVSI deals with 3D space use and overlap at individual, population or species level, we suggest two further possible applica-tions. Time has rarely been considered when estimating HR size and overlap [30], and it has never been incorporated as a fourth dimension (i.e., 4D volumetric home ranges; x, y, h g and t). Space use can vary in time and such variations may have implications for habitat selection and use, trophic interactions and predator-prey relationships [9]. Accordingly, the MVOI and MVSI can be used to quantify how the individuals' HR 3D s overlap changes over time in a certain population (e.g., between life history stages or before and after experimental manipulations) or to evaluate the degree of overlap between the HR 3D of the same individual at different phenological stages (e.g., nesting and post-nesting periods). Another important application of our new metrics is the evaluation of the robustness of the HR assessment through the degree of overlap of several HR 3D estimators: if they are in good agreement, then the MVOI will be close to 100.
By combining overlaps into a single metric, one could lose useful information regarding the extent of overlap between animals with different attributes (e.g., sex, size, diet preference, nest location, etc.), from which one could obtain statistics. Accordingly, although we conceived the MVOI and MVSI as non-pairwise metrics among multiple 3D home ranges, nothing prevents researchers from applying them in a pairwise manner by splitting the dataset of the tracked animals based on the attribute of interest (e.g., individual i versus individual j, females versus males, youngs versus adults etc.).

Mathematical Properties of the Proposed Metrics
The proposed metrics have the same mathematical properties as the GOI and the PGOI: (1) whatever the number of HR 3D s under study, the MVOI and MVSI return a single overlap measure; (2) in case of perfectly segregated HR 3D , the MVOI is equal to 0 and the MVSI to 100; (3) in case of perfectly overlapping HR 3D , the MVOI is equal to 100 and the MVSI to 0; and (4) in any other case, the MVOI and MVSI return a value between 0 and 100. Although apparently complex, the MVOI corresponds to the linear equation Y = 100 × (b − X)/(b − a), where a is the volume of the largest HR 3D , b is the sum of the HR 3D and X is the volume of the union of the HR 3D , which varies depending upon the degree of overlap. The first derivative of the MVOI with respect to X is equal to −100/(b − a), thus every unitary (e.g., 1 km 3 ) increase/decrease of ∪V 3D (due, for instance, to changes of the HR 3D over time) determines a decrease/increase in the MVOI that is constant and independent of the initial value assumed by ∪V 3D . This assures that small/big changes to the HR 3D overlaps proportionally determine small/big changes to our overlap indices, regardless of the initial degree of overlap.

Conclusions
While it is evident that animals live in a 3D environment, to date, biotelemetry data have mostly been analysed and modelled in two dimensions. However, advances in biotelemetry are increasingly providing the opportunity to gather additional data beyond 2D location coordinates.
Accordingly, in this study, we proposed a new approach that integrates mathematical concepts with telemetry data to provide the opportunity to define 3D space overlap among n individuals, populations or species. The usefulness of our metrics is not limited to avian studies and can be applied to any dataset that includes 3D coordinates. A further strong point of our new metrics is that they can be calculated using standard GIS operations and can be obtained by using any free GIS software.

Institutional Review Board Statement:
The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Institutional Review Board of ISPRA (protocol code 38403/T-A 31, dated 28 June 2016).

Informed Consent Statement: Not applicable.
Data Availability Statement: The dataset used in this study is available from the first author on reasonable request.