# Methods

The geospatial and ecological datasets presented in this report were collated from a wide range of sources. This section summarises the key datasets used and the data analysis methods. More detail for a technical audience may be found at the references provided. 

## Data

### Urban extent 
Analyses of urban ecosystems are highly sensitive to the definition of the borders of each urban area. As this report is part of a systematic analysis that can be applied anywhere in the world, we used an existing definition of urban boundaries from 2020. This was the GHS-SMOD R2023A dataset (Schiavina et al 2023a). For all maps, these urban areas are identified with a dark grey dashed line. A buffer of 1000 m surrounding the urban area boundary is shown for reference.

### Land cover and vegetation 
Urban land and water cover types were mapped in 2020 at a spatial resolution of approximately 10 m by 10 m pixel size. These data were were extracted from the ESA WorldCover dataset 2020 (Zanaga et al. 2021). Tree height was mapped in 2020 at a spatial resolution of approximately 10 m by 10 m pixel size. We extracted vegetation heights only for regions mapped as tree cover, to ensure internal consistency. These data were extracted from the global canopy height database (Lang et al. 2023). The dominant plant functional types in each region in 2020 were extracted from the Global Plant Functional Types dataset v2.0.8 (Harper et al. 2023). 


### Topography and soils 
Elevation data were extracted for the urban area from the Amazon Web Service Terrain Tiles (https://registry.opendata.aws/terrain-tiles/) using the elevatr R package (Hollister et al. 2023). The elevation map was masked to only show terrain that is not normally covered by open water. Soil hydrologic type data were extracted at a 250 m resolution from the Global Hydrologic Soil Groups dataset (Ross et al. 2018). Soil erodibility (K) value was extracted at a 25 km resolution from Version 1.1. of the JRC/UNiBasel RUSLE-based Global Soil Erosion Modelling Platform (Borrelli et al. 2017).


### Climate
Long-term average annual temperature and precipitation data were extracted at a 30 arc second resolution from the WorldClim dataset (Fick et al. 2017). To indicate the magnitude of extreme precipitation events, the 100-year return values of daily precipitation were extracted at a 1 degree resolution from a series of ensemble predictions (Ruff and Pfahl, 2023). Rainfall-runoff erosivity values (R) were extracted at a 25 km resolution from Version 1.1. of the JRC/UNiBasel RUSLE-based Global Soil Erosion Modelling Platform (Borrelli et al. 2017). 

The azimuth and zenith of incoming solar radiation on five daylight occasions during a summer month was calculated using the insol R package (Corripio 2003; Corripio 2021). The indicative summer month was taken as June in the Northern Hemisphere, and December in the Southern, and the five occasions were sampled randomly within daylight hours of the month. The clear-sky ultraviolet index in 2020 was extracted from the Tropospheric Emission Monitoring Internet Service UV Index data product v2.2 (Van Geffen et al. 2017). The 95th percentile of the values recorded in 2020 for each pixel were taken, to represent a high UVI value for the location. Monthly Universal Thermal Comfort Index (UTCI) values from 2020 were extracted at a resolution of 1km from the GloUTCI-M dataset (Yang et al. 2024). The month with the highest and lowest mean monthly UTCI were extracted to represent the hottest and coldest thermal conditions respectively. 

Future climate projections were extracted from CMIP6 multi-model ensemble averages (Government of Canada 2024). Projected average annual precipitation and temperature for the period 2060-2069 were extracted for each of five scenarios (SSP1-1.9, SSP1-2.6, SSP2-4.5, SSP3-7.0, SSP4-6.0, SSP5-8.5) at a 1 degree resolution. 

### Sociological
The total population and population density within the urban area at a resolution of 100 m was extracted from the GHS-SMOD R2023A dataset (Schiavina et al 2023b). Average building height above the ground (ANBH)was extracted at a 100 m resolution from the Global Human Settlement building height dataset GHS-BUILT-H R2023A (Pesaresi et al. 2021; Pesaresi and Politis 2023). Place names including the name of the urban area and suburbs within were retrieved from the GeoNames database (Rowlingson 2019; GeoNames 2024). The difficulty in traversing space as a pedestrian was extracted from a global dataset of walking travel friction (Weiss et al. 2018; Weiss et al. 2020). 

### Biodiversity records 

Geolocated biodiversity records were extracted from the Global Biodiversity Information Facility (GBIF) using the rgbif package for R (Chamberlain et al. 2017; Chamberlain et al. 2024). The downloaded data for this urban area can be found using the download-specific reference (GBIF 2024). 



## Data analysis

All data processing and analyses were conducted in R v4.1. (R Core Team 2024), with geospatial processing conducted using the terra package (Hijmans 2024). All maps were created using tmap (Tennekes 2018). The methods collated here are an extension of previous systematic global comparisons of urban ecosystem services (Richards et al. 2022).


### Biodiversity
Biodiversity records from GBIF are largely derived from citizen science reports of animal and plant occurrences. These records are susceptible to biases in the taxa that are sampled and reported, particularly in urban areas. Therefore quantitative approaches were not used to analyse these data, they were discussed qualitatively as an indicator of likely species presence. Birds are typically better represented in citizen science datasets than other taxa, so for the bird records an analysis of species richness was performed at the scale of the entire urban extent community. The occurrence of bird species across sampling day/ observer combinations was analysed using the Chao2 estimator, which is based on species accumulation curve theory (Chao et al. 2016; Chao and Chiu 2016). Modelled species richness was reported if the estimated coverage estimate for the dataset is greater than 0.95.  

Ecological connectivity was modelled for two model species that represent (1) a mobile species with relatively high tolerance of urban and disturbed econsystem conditions and (2) a less mobile species specialised to forest and woodland ecosystems with low tolerance of urban and open ecosystems. These two model species are hypothetical and are used to demonstrate how the landscape configuration around the urban extent influences the suitability for different types of organisms. The relative friction of each ecosystem for the two model species was parameterised according to Table 1, with this score modified by the area of habitat within patches of suitable ecosystem cover (Richards et al. 2024). This modification represents the value of larger landscape patches for urban biodiversity (Beninde et al. 2015). 


### Regulating ecosystem services
Runoff retention by ecosystems was modelled using a curve number approach (USDA 1986), with curve numbers for combined soil, land cover, and plant functional types assigned following the GCN250 dataset (Jaafar et al. 2019). Runoff retention was modelled under an extreme (100-year return period) daily rainfall (Borrelli et al. 2017). Coastal risk values and the protection provided by coastal ecosystems was extracted from an existing published dataset (Chaplin-Kramer et al. 2019). The contribution of vegetated ecosystems in reducing soil erosion was modelled using a Revised Universal Soil Loss Equation (RUSLE) approach applied to quantify vegetation contributions (Guerra et al. 2014). Parameters for soil erodibility, rainfall erosivity, and the modification of erosion by vegetation were applied following a similar global modelling analysis (Borrelli et al. 2017). Rainfall was modelled under average antecedent conditions using the hbrc package (Richards 2023a; Richards et al. 2023b). 

Thermal discomfort under hot weather conditions was extracted as the maximum UTCI value for the urban area in 2020 (Yang et al. 2024). The shade provided by urban vegetation was modelled using a mechanistic 3D modelling approach (Richards 2023; Richards et al. 2023). Solar radiation was modelled as it travelled through a 3D environment using the raytracer R package (Morgan-Wall 2024), under the summertime sun locations described above. Vegetation and elevation data were extracted as described above. The role of elevation and buildings in providing shade was modelled and thus tree shade benefits were accounted as additional shade provided (Richards 2023; Richards et al. 2023). Protection from ultraviolet radiation at the neighbourhood scale was modelled using a general approach proposed by Na et al. (2014), implemented in the hbrc package (Richards 2023; Richards et al. 2023). Tree cover, height, and UV exposure were parameterised using the datasets described above.

Carbon stocks and potential sequestration (net primary productivity, NPP) were modelled based on the area of tree cover present in the urban area. The area of tree cover was multiplied by a modelled estimate of per ha carbon stock and potential sequstration given the climate conditions. Potential NPP was modelled as a function of rainfall and temperature based on the NCEAS equation (Del Grosso et al 2008). A modification value for urban trees and conversion factor between urban tree NPP and carbon stock was taken from Jenerette et al (2024).


### Socio-cultural ecosystem services
Accessibility to urban ecosystems was indicated by the travel time between each person's home and the closest urban ecosystem or water ecosystem, using a travel friction modelling approach (Weiss et al. 2018; Weiss et al. 2020; Richards et al. in review). Urban ecosystems were defined as patches of vegetation or bare ground at least 1 ha in contiguous patch size, following World Health Organisation guidelines for suitable recreational spaces (Richards et al. 2022). Water was defined as all open water and wetlands present, regardless of patch size (Richards et al. in review). Population distribution and pedestrian friction were parameterised following the datasets outlined above. Other cultural ecosystem services were assessed qualitatively through a literature review.


### Provisioning ecosystem services
Provisioning ecosystem services were assessed qualitatively through referring to the land cover map, which indicates some land covers that may be relevant for food production. This was supplemented with a literature review.



# References

Richards DR, Belcher RN, Carrasco LR, Edwards PJ, Fatichi S, Hamel P, Masoudi M, McDonnel MJ, Peleg N, Stanley MC (2022). Global variation in contributions to human well-being from ecosystem services provided by urban vegetation. One Earth 5: 522-533.

Morgan-Wall T (2024). rayshader: Create Maps and Visualize Data in 2D and 3D. R package version 0.37.3, <https://CRAN.R-project.org/package=rayshader>.

Richards DR, Belcher R, Schindler M. (In review). Quantifying travel time to access urban ecosystems globally. 

Richards DR, Herzig A, Abbott A, Ausseil A-G, Guo J, Sood A, Lavorel S (2023b). Diverse contributions of nature to climate change adaptation in an upland landscape. Ecosystems and People 19, 2225647.

Richards, D., Etherington, T. R., Herzig, A., & Lavorel, S. (2024). The Importance of Spatial Configuration When Restoring Intensive Production Landscapes for Biodiversity and Ecosystem Service Multifunctionality. Land, 13(4), 460.

Richards DR (2023a). hbrc: Functions for estimating ecosystem services indicators using high-resolution spatial datasets. <https://github.com/manaakiwhenua/hbrc>.

Chaplin-Kramer, R., Sharp, R. P., Weil, C., Bennett, E. M., Pascual, U., Arkema, K. K., ... & Daily, G. C. (2019). Global modeling of nature’s contributions to people. Science, 366(6462), 255-258.

Jaafar, H. H., Ahmad, F. A., & El Beyrouthy, N. (2019). GCN250, new global gridded curve numbers for hydrologic modeling and design. Scientific data, 6(1), 145.

USDA (1986). Urban Hydrology for Small Watersheds. In United States Department of Agriculture Technical Release 55 (TR-55)).

Chao, A., & Chiu, C. H. (2016). Nonparametric estimation and comparison of species richness. eLS, 1-11.

Chao A, Ma KH, Hsieh TC, Chiu C (2016). SpadeR: Species-Richness Prediction and Diversity Estimation with R. R package version 0.1.1, <https://CRAN.R-project.org/package=SpadeR>.

R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. <https://www.R-project.org/>.

Hijmans R (2024). terra: Spatial Data Analysis. R package version 1.7-78, <https://CRAN.R-project.org/package=terra>.

Tennekes M (2018). “tmap: Thematic Maps in R.” Journal of Statistical Software, 84 (6), 1-39. doi:10.18637/jss.v084.i06 <https://doi.org/10.18637/jss.v084.i06>.

Rowlingson B (2019). geonames: Interface to the "Geonames" Spatial Query Web Service. R package version 0.999, https://CRAN.R-project.org/package=geonames.

GeoNames. (2024). GeoNames Application Programming Interface. Retrieved August 2024, from http://www.geonames.org

Schiavina M., Melchiorri M., Pesaresi M. (2023a): GHS-SMOD R2023A - GHS settlement layers, application of the Degree of Urbanisation methodology (stage I) to GHS-POP R2023A and GHS-BUILT-S R2023A, multitemporal (1975-2030)European Commission, Joint Research Centre (JRC) PID: http://data.europa.eu/89h/a0df7a6f-49de-46ea-9bde-563437a6e2ba, doi:10.2905/A0DF7A6F-49DE-46EA-9BDE-563437A6E2BA

Schiavina M., Freire S., Carioli A., MacManus K. (2023b): GHS-POP R2023A - GHS population grid multitemporal (1975-2030).European Commission, Joint Research Centre (JRC) PID: http://data.europa.eu/89h/2ff68a52-5b5b-4a22-8f40-c41da8332cfe, doi:10.2905/2FF68A52-5B5B-4A22-8F40-C41DA8332CFE

Zanaga, D., Van De Kerchove, R., De Keersmaecker, W., Souverijns, N., Brockmann, C., Quast, R., Wevers, J., Grosu, A., Paccini, A., Vergnaud, S., Cartus, O., Santoro, M., Fritz, S., Georgieva, I., Lesiv, M., Carter, S., Herold, M., Li, Linlin, Tsendbazar, N.E., Ramoino, F., Arino, O., 2021. ESA WorldCover 10 m 2020 v100. https://zenodo.org/records/5571936

Lang, N., Jetz, W., Schindler, K., & Wegner, J. D. (2023). A high-resolution canopy height model of the Earth. Nature Ecology & Evolution, 7(11), 1778-1789.

Yang, Z., Peng, J., Liu, Y., Jiang, S., Cheng, X., Liu, X., ... & Yu, X. (2024). GloUTCI-M: a global monthly 1 km Universal Thermal Climate Index dataset from 2000 to 2022. Earth System Science Data, 16(5), 2407-2424.

Harper, K. L., Lamarche, C., Hartley, A., Peylin, P., Ottlé, C., Bastrikov, V., ... & Defourny, P. (2023). A 29-year time series of annual 300 m resolution plant-functional-type maps for climate models. Earth System Science Data, 15(3), 1465-1499.

Hollister J, Shah T, Nowosad J, Robitaille A, Beck M, Johnson M (2023). elevatr: Access Elevation Data from Various APIs. doi:10.5281/zenodo.8335450, R package version 0.99.0, https://github.com/jhollist/elevatr/.

Ross, C.W., L. Prihodko, J.Y. Anchang, S.S. Kumar, W. Ji, and N.P. Hanan. 2018. Global Hydrologic Soil Groups (HYSOGs250m) for Curve Number-Based Runoff Modeling. ORNL DAAC, Oak Ridge, Tennessee, USA. https://doi.org/10.3334/ORNLDAAC/1566

Fick, S.E. and R.J. Hijmans, 2017. WorldClim 2: new 1km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37 (12): 4302-4315.

Van Geffen, J., Van Weele, M., Allaart, M. and Van der A, R. TEMIS UV index and UV dose operational data products, version 2. Royal Netherlands Meteorological Institute (KNMI), Dataset, 2017. doi:10.21944/temis-uv-oper-v2.

Borrelli, P., Robinson, D. A., Fleischer, L. R., Lugato, E., Ballabio, C., Alewell, C., Meusburger, K., Modugno, S., Schütt, B., Ferro, V., Bagarello, V., Van Oost, K., Montanarella, L & Panagos, P. (2017). An assessment of the global impact of 21st century land use change on soil erosion. Nature communications, 8(1), 2013.

Ruff, F., & Pfahl, S. (2023). Global estimates of 100-year return values of daily precipitation from ensemble weather prediction data. EGUsphere, 2023, 1-21.

Corripio, J. G.: 2003, Vectorial algebra algorithms for calculating terrain parameters from DEMs and the position of the sun for solar radiation modelling in mountainous terrain, International Journal of Geographical Information Science 17(1), 1-23.

Corripio JG (2021). insol: Solar Radiation. R package version 1.2.2, commit 4fe6f26351af22121a5319b4a8b4950b66ca6226, https://github.com/cran/insol/.

Weiss, D. J., Nelson, A., Vargas-Ruiz, C. A., Gligorić, K., Bavadekar, S., Gabrilovich, E., ... & Gething, P. W. (2020). Global maps of travel time to healthcare facilities. Nature medicine, 26(12), 1835-1838.

Weiss, D. J., Nelson, A., Gibson, H. S., Temperley, W., Peedell, S., Lieber, A., ... & Gething, P. W. (2018). A global map of travel time to cities to assess inequalities in accessibility in 2015. Nature, 553(7688), 333-336.

Chamberlain S, Barve V, Mcglinn D, Oldoni D, Desmet P, Geffert L, Ram K (2024). rgbif: Interface to the Global Biodiversity Information Facility API. R package version 3.8.0, https://CRAN.R-project.org/package=rgbif.

Chamberlain S, Boettiger C (2017). “R Python, and Ruby clients for GBIF species occurrence data.” PeerJ PrePrints. https://doi.org/10.7287/peerj.preprints.3304v1.

Government of Canada. (2024). CMPI6 climate scenarios. Retrieved August 2024, from https://climate-scenarios.canada.ca/?page=cmip6-scenarios

Guerra, C. A., Pinto-Correia, T., & Metzger, M. J. (2014). Mapping Soil Erosion Prevention Using an Ecosystem Service Modeling Framework for Integrated Land Management and Policy. Ecosystems, 17(5), 878–889. 

Na, H. R., Heisler, G. M., Nowak, D. J., & Grant, R. H. (2014). Modeling of urban trees’ effects on reducing human exposure to UV radiation in Seoul, Korea. Urban forestry & urban greening, 13(4), 785-792.


Del Grosso, S., Parton, W., Stohlgren, T., Zheng, D., Bachelet, D., Prince, S., Hibbard, K., Olson, R., 2008. Global potential net primary production predicted from vegetation class, precipitation, and temperature. Ecology 89, 2117–2126. https://doi.org/10.1890/07-0850.1

Jenerette, G.D., Herrmann, D.L., 2024. A global synthesis of reported urban tree carbon production rates and approaches. Front Ecol Evol 11. https://doi.org/10.3389/fevo.2023.1244418

Beninde, J., Veith, M., & Hochkirch, A. (2015). Biodiversity in cities needs space: a meta‐analysis of factors determining intra‐urban biodiversity variation. Ecology letters, 18(6), 581-592.

Pesaresi, M.; Politis, P. (2023). GHS-BUILT-H R2023A - GHS building height, derived from AW3D30, SRTM30, and Sentinel2 composite (2018).European Commission, Joint Research Centre (JRC). doi:10.2905/85005901-3A49-48DD-9D19-6261354F56FE

Pesaresi M., Corbane C., Ren C., and Edward N. (2021). Generalized Vertical Components of built-up areas from global Digital Elevation Models by multi-scale linear regression modelling. PLOS ONE 16, e0244478. doi.org/10.1371/journal.pone.0244478
