Remote Sensing Methods for the Biophysical Characterization of Protected Areas Globally: Challenges and Opportunities

: Protected areas (PAs) are a key strategy to reverse global biodiversity declines, but they are under increasing pressure from anthropogenic activities and concomitant effects. Thus, the heterogeneous landscapes within PAs, containing a number of different habitats and ecosystem types, are in various degrees of disturbance. Characterizing habitats and ecosystems within the global protected area network requires large-scale monitoring over long time scales. This study reviews methods for the biophysical characterization of terrestrial PAs at a global scale by means of remote sensing (RS) and provides further recommendations. To this end, we ﬁrst discuss the importance of taking into account the structural and functional attributes, as well as integrating a broad spectrum of variables, to account for the different ecosystem and habitat types within PAs, considering examples at local and regional scales. We then discuss potential variables, challenges and limitations of existing global environmental stratiﬁcations, as well as the biophysical characterization of PAs, and ﬁnally offer some recommendations. Computational and interoperability issues are also discussed, as well as the potential of cloud-based platforms linked to earth observations to support large-scale characterization of PAs. Using RS to characterize PAs globally is a crucial approach to help ensure sustainable development, but it requires further work before such studies are able to inform large-scale conservation actions. This study proposes 14 recommendations in order to improve existing initiatives to biophysically characterize PAs at a global scale.


Introduction
Protected areas (PAs) are one of the main conservation strategies to counter the current biodiversity crisis [1]. However, PAs are under ongoing social, economic and environmental threats, and so the conservation of biodiversity within PAs and the restoration of PAs constitute one of the main current sociopolitical challenges [2]. The long-term conservation benefits of PAs depend on timely management actions based on relevant data and models that can predict the responses of ecosystems to various stress factors [3,4].
Anthropogenic activities and the changes in land use they generate have an impact on how efficient PAs are at protecting biodiversity globally [5,6]. Moreover, climate change impacts severely affect PAs, including increased frequency of flooding, soil erosion and plant water stress [7]. It is increasingly recognized that even large, historically stable ecosystems (such as the Amazon) are threatened and could undergo regime shifts to alternative ecosystems within 50 years [8].
provide additional information on ecosystems contained within those ecoregions and have instead been used to prioritize the importance of conservation of larger regions [51].
This study seeks to provide recommendations for the biophysical characterization of terrestrial PAs at the global scale by means of RS. To this end, in Section 2, we discuss the importance of taking into account structural and functional attributes, as well as integrating a broad spectrum of variables, to account for the different ecosystem and habitat types within PAs, reviewing examples at the local and regional scales. In Section 3, we discuss potential candidate input variables at the global scale for the characterization of PAs, as well as the challenges and limitations of existing global environmental stratifications and the biophysical characterization of PAs, and offer recommendations. Computational and interoperability issues are also discussed, as well as the potential for cloud-based platforms linked to earth observations to support large-scale characterization of PAs. Finally, Section 4 provides a summary list of recommendations. Although focusing on terrestrial areas, we also mention a few examples of RS data used to characterize marine protected areas (MPAs).

Relevant Ecological Units and Descriptors
In order to comprehend the ecological complexity in PAs, biophysical characterization of PAs should take into account the different ecosystem and habitat types that are present within them and distinguish their ecological attributes as much as possible, including structural and functional attributes. To this end, a wide range of environmental descriptors should be included in the analysis, including drivers that ultimately shape ecosystems ( Figure 1). ISPRS Int. J. Geo-Inf. 2021, 10, x FOR PEER REVIEW 3 of 21 provide additional information on ecosystems contained within those ecoregions and have instead been used to prioritize the importance of conservation of larger regions [51]. This study seeks to provide recommendations for the biophysical characterization of terrestrial PAs at the global scale by means of RS. To this end, in Section 2, we discuss the importance of taking into account structural and functional attributes, as well as integrating a broad spectrum of variables, to account for the different ecosystem and habitat types within PAs, reviewing examples at the local and regional scales. In Section 3, we discuss potential candidate input variables at the global scale for the characterization of PAs, as well as the challenges and limitations of existing global environmental stratifications and the biophysical characterization of PAs, and offer recommendations. Computational and interoperability issues are also discussed, as well as the potential for cloud-based platforms linked to earth observations to support large-scale characterization of PAs. Finally, Section 4 provides a summary list of recommendations. Although focusing on terrestrial areas, we also mention a few examples of RS data used to characterize marine protected areas (MPAs).

Relevant Ecological Units and Descriptors
In order to comprehend the ecological complexity in PAs, biophysical characterization of PAs should take into account the different ecosystem and habitat types that are present within them and distinguish their ecological attributes as much as possible, including structural and functional attributes. To this end, a wide range of environmental descriptors should be included in the analysis, including drivers that ultimately shape ecosystems ( Figure 1). The assessment of structural attributes, such as vegetation height or heterogeneity by means of RS, helps distinguish characteristic ecosystems and habitats within PAs such as The assessment of structural attributes, such as vegetation height or heterogeneity by means of RS, helps distinguish characteristic ecosystems and habitats within PAs such as forests, wetlands, grasslands, shrublands, dunes and riparian habitats, among others. Furthermore, RS variables related to functional attributes, such as vegetation phenology or energy fluxes, have been proven to complement and improve habitat and ecosystem classifications based only on structural features by capturing the occurrence of natural disturbances and vegetation productivity, among others [52][53][54]. Several studies have reviewed the use of RS for assessing habitat and ecosystem structure, functions and condition in PAs at the local and regional scales [24,[55][56][57][58][59][60][61][62].
With regard to structural attributes, wetlands, riparian forests and dune habitats, for example, have been mapped by means of texture-and object-based RS data analysis and machine learning algorithms in order to characterize and monitor changes in PAs [63][64][65][66][67][68][69][70][71][72][73]. Grasslands have been accurately mapped using time series of RS data [74]. Forest and shrubland structure has been mapped by means of very high resolution imagery [75][76][77]. Tree species richness across the tropics has been mapped by means of full-waveform lidar data [78]. Vegetation structure has been mapped at the local and regional levels in PAs by means of manned and unmanned aerial vehicles carrying airborne LiDAR and multi-and hyperspectral sensors [79][80][81]. Chetan and Dornik [82] quantified changes in vegetation greenness and structures within Natura 2000 sites over 20 years. The vegetation heterogeneity and pattern was characterized by means of image texture measures (e.g., gray-level co-occurrence matrix) derived from RS data [83][84][85][86][87][88][89].
In relation to functional attributes, several studies have quantified vegetation productivity over time by means of remote sensing-derived indices and found correlation with biodiversity patterns [90][91][92][93]. Moreover, the effect of disturbances such as post-fire forest vegetation regrowth was studied by means of different RS vegetation indices [94,95]. For a recent review of methods, sensors and ecosystem structural and functional attributes assessed by means of RS in PAs, see [33].
Furthermore, given the inherent ecological complexity that can be found within PAs, their systematic characterization needs to extend specific habitat or ecosystem mapping and assessment methods so that all habitat and ecosystem components that are present within them are taken into account [96,97]. By stratifying the natural landscape into homogeneous regions defining ecological units, the complexity of PAs can be converted into something that is more manageable and understandable [98]. For example, if a protected landscape contains both a lake and mountains, separating both elements cartographically would help inform and support adaptive management. In this regard, methods to characterize PAs should rely on a comprehensive list of environmental quantitative descriptors based on RS data, which could be categorized into different topics: (1) vegetation, including structure, phenology and disturbances; (2) climate; (3) water budget; (4) energy exchanges; (5) terrain and (6) soil, among others.
As was previously mentioned, vegetation-related variables, such as the amount of woody and herbaceous biomass or different vegetation indices, can help us distinguish between broad ecosystem types (such as forests, grasslands or wetlands) by capturing their structure, phenology and productivity [99]. Climatic descriptors, such as precipitation and temperature, are also important variables to be included in biophysical assessments to represent seasonality, extremes and limiting climatic factors [100][101][102][103]. Topographic gradients drive many patterns and processes in hydrology and ecology and are key to understanding the variation of habitats and biodiversity [104,105]. Water-related variables are also a good proxy for plant water stress and the presence of aquatic ecosystems and can therefore supplement the information on climate and vegetation by distinguishing differing responses to available water [106][107][108]. Variables that describe the energy exchanges between the land surface and the atmosphere, as well as the partition of energy into ground and vegetation, are also essential for ecological assessment and modeling [109].
Soil data are often ignored when characterizing PAs, but more than 25% of the Earth's species live only in the soil [110]. Aside from that, soils form the foundation for many vegetation types and provide key supporting ecosystem services that are crucial for the maintenance of other types of services [111]. Given that soil biodiversity cannot be directly monitored by RS, soil descriptors that can be directly or indirectly monitored by RS and modeling can act as proxies [112,113]. In this regard, soil organic carbon appears to be one of the main drivers of soil microbial biodiversity at the global scale [114][115][116], particularly in extreme environments with low net primary productivity such as polar [117] and dryland regions [118]. Soil texture is also a relevant descriptor, since previous research has demonstrated that soil biota abundance and biodiversity, particularly soil microorganisms, increase with a decreasing soil particle size [119].

Global Input Variables and Data Sources
In the previous section, we reviewed the importance of taking into account the structural and functional attributes as well as integrating a broad spectrum of variables to account for the different ecosystem and habitat types within PAs. In this subsection, we give a list of potential candidate input variables mapped at global scale for global characterizations of PAs and discuss some limitations and recommendations.
Data sources presenting time series and regular updates at the global scale should be favored over single records in time to allow for the assessment of change over time and identify reference conditions. When correlated variables are used, principal component analysis can often be applied in order to compress them and use the resulting uncorrelated axes as input for the models to avoid redundant predictors [120].
Given that global RS data usually shows greater inaccuracies than local or regional datasets, the use of ensembles of different input data or models corresponding to the same variable might be advantageous, providing more accurate outputs as well as better conveying uncertainty [121][122][123][124][125]. Aside from that, many biophysical variables mapped at the local or regional scales are not available at a global scale, which might limit the relevance of global analyses for local scale management. Therefore, global characterization of PAs should be primarily aimed at informing larger scale conservation and management actions and plans, unless no better information is available at the local or regional scales. Table 1 lists a set of recommended variables that can be used at a global scale for the biophysical characterization of terrestrial PAs. The list is not exhaustive, but it provides a wide range of relevant variables, including potential data sources. A more comprehensive list of potential variables can be found at the Global Climate Observing System Programme (https://public.wmo.int/en/programmes/global-climate-observing-system /essential-climate-variables; accessed 1 June 2021) or the Copernicus Global Land Service (http://land.copernicus.vgt.vito.be/PDF/portal/Application.html#Home; accessed 1 June 2021). A table with additional information, including the URLs of data sources, can be found as supplementary material. For MPAs, previous studies have highlighted candidate variables measurable by RS which are relevant to characterizing marine habitats [126][127][128][129][130]. They include, among others, bathymetry, the concentration of chlorophylla, sea surface temperature and sea surface salinity. A comprehensive list of these marine variables-together with access to the RS measurements of these variables-can be also found at the Copernicus Marine Service (https://resources.marine.copernicus.eu /?option=com_csw&task=results; accessed 1 June 2021) and the Living Wales Geoportal (https://wales.livingearth.online/data/environmental-variables/marine/; accessed 1 June 2021).

Global Environmental Stratifications
There are several biophysical characterizations available at the global scale, partially or totally based on RS data and modeling. Metzger et al. [137][138][139] used a broad set of bioclimatic variables to stratify the world in 18 environmental zones in order to support global ecosystem research and monitoring. Ivits et al. [53] mapped global ecosystem functional types using vegetation phenology and productivity variables by means of principal components and cluster analysis. Sayre [140] developed a map of global ecological land units using bioclimate, landforms, lithology and land cover variables. Tuanmu and Jetz [141] developed 14 remote sensing-based metrics to characterize habitat heterogeneity at a 1 km resolution at a global scale based on textural information extracted from the enhanced vegetation index (EVI) [142] and found out that bird species richness was strongly associated with habitat heterogeneity. Jung et al. [17] developed a global map of terrestrial habitat types following the IUCN habitat classification scheme (https://www.iucnredlist.org/resources/habitat-classification-scheme; accessed 1 June 2021) based on land cover, climate and land use data. Sayre et al. [51] developed a global classification of world climate regions and world ecosystems based on environmental descriptors, such as landforms, moisture, temperature, vegetation type and land use. Finally, in [143], a global ecosystem typology, including indicative distribution maps based on a large set of different environmental descriptors, existing global occurrence maps of specific ecosystem types and previous global environmental characterizations, was developed. They used a hierarchical classification system that first characterized ecosystems by their ecological functions and then distinguished ecosystems with contrasting species assemblages.
These global stratification initiatives are not limited to PAs and are indeed useful for prioritizing the conservation importance of larger regions. However, RS and modeling efforts specifically aimed to systematically characterize PAs could provide more relevant information needed to inform several policy initiatives as well as support management applications in PAs at regional or global scales, such as the assessment of ecological representativeness, the prioritization of PAs, connectivity assessments and the mapping of new areas requiring protection.

Global Characterization of Protected Areas
In relation to global characterizations within PAs by means of RS and modeling, in [144], the EODHaM system for characterizing habitats in PAs and surrounding areas was developed using earth observation data and expert knowledge. They used a semiautomated statistical procedure based on data related to the terrain, vegetation, water balance and land use. As part of the Digital Observatory for Protected Areas (DOPA) [145], in [120], PAs were systematically stratified globally into different habitat functional types based on remote sensing data and modeling, allowing for the quantification of the similarity between a reference area (representing a habitat functional type) and the surrounding areas based on a set of ecological indicators [146][147][148]. The method also graphically compares the ecological features of each habitat functional type found in a PA to help identify their main characteristics and understand the main biophysical gradients that occur at the PA level ( Figure 2). The methodology uses a combination of several multivariate statistical analyses based on different global predictors that account for the climate, topography, vegetation and water exchanges. One of the advantages of this methodology is that the analysis is fully automated and can be performed at different spatial resolutions, which is especially important when dealing with smaller PAs. Furthermore, the similarity maps that are produced can also be used to identify new potential areas to be protected to strengthen ecological connectivity. When used in conjunction with forecasted bioclimatic data, the approach can further help identify new areas for conservation by considering current and climate change scenarios [147]. pography, vegetation and water exchanges. One of the advantages of this methodology is that the analysis is fully automated and can be performed at different spatial resolutions, which is especially important when dealing with smaller PAs. Furthermore, the similarity maps that are produced can also be used to identify new potential areas to be protected to strengthen ecological connectivity. When used in conjunction with forecasted bioclimatic data, the approach can further help identify new areas for conservation by considering current and climate change scenarios [147]. . NDVI stands for the normalized difference vegetation index, and NDWI stands for the normalized difference water index. A detailed description of the study variables and the methodology followed can be found in [120].
When prioritizing and ranking PAs, most studies have focused on species diversity to measure uniqueness [149,150]. However, biophysical characterizations have also been used along with biotic variables to perform gap and representation analyses in PAs [51,151]. Dubois et al. [147] proposed a methodology to assess the uniqueness of PAs based on biophysical variables. However, they lacked the means to decompose each analyzed area into areas with similar ecological features. The methodology proposed by [120] partially solved this issue by identifying habitat functional types and mapping similar areas at the ecoregion scale. This approach could be used to further create a composite indicator for each PA that reflects the biophysical richness of PAs and the uniqueness of their habitats. Coastal PAs in particular should be taken into account when developing . NDVI stands for the normalized difference vegetation index, and NDWI stands for the normalized difference water index. A detailed description of the study variables and the methodology followed can be found in [120].
When prioritizing and ranking PAs, most studies have focused on species diversity to measure uniqueness [149,150]. However, biophysical characterizations have also been used along with biotic variables to perform gap and representation analyses in PAs [51,151]. Dubois et al. [147] proposed a methodology to assess the uniqueness of PAs based on biophysical variables. However, they lacked the means to decompose each analyzed area into areas with similar ecological features. The methodology proposed by [120] partially solved this issue by identifying habitat functional types and mapping similar areas at the ecoregion scale. This approach could be used to further create a composite indicator for each PA that reflects the biophysical richness of PAs and the uniqueness of their habitats. Coastal PAs in particular should be taken into account when developing indices of this type, given their inherent complexity as ecotones and the higher pressures they are exposed to because of human developments that are often concentrated along the coasts [152][153][154][155][156].
Perhaps the main limitation of global biophysical assessments using RS is the lack of ground truthing and comparison maps in order to evaluate results [157]. In this regard, resulting habitat and ecosystem types based on RS methods could be classified according to existing global typologies in order to serve and support different initiatives of habitat and ecosystem monitoring globally. For example, a hierarchical classification framework could be applied to the ecological features resulting from the methodology developed in [120], in which some key variables guide the first broad set of typologies and other variables help distinguish more specific subclasses according to existing typologies. The previously mentioned recent global environmental stratification initiatives already provide potential comparison maps, such as the IUCN Global Ecosystem Typology [158] and the set of world climate regions and world ecosystems [51]. The proposed approach would allow for taking similar regional features into consideration as well as going deeper into a specific global ecosystem type (e.g., tropical moist forests and mangroves).
In relation to the marine realm, current efforts to globally characterize PAs by means of RS have focused on the use of bathymetry. As such, DOPA uses a model of global bathymetry that is partially based on RS data to compute a marine habitat diversity index for MPAs [145]. The fact that (1) most RS methods can only derive information from the upper layer of the ocean (with the exception of altimeters for coarse scale bathymetry), (2) the spatial resolution of available RS data may be too coarse to characterize MPAs and (3) RS-based management of MPAs requires large financial and human resources constitutes major impediments to the use of RS data for characterizing MPAs [130]. This may explain why global characterization of MPAs using RS is limited. However, initiatives to characterize PAs using a broader set of RS-measured variables are more numerous at the regional [130,159,160] and local scales [161,162]. Beyond the characterization of MPAs, RS data have been used to assess the connectivity of MPA networks [154,163] and to delineate bioregions that can be further used as a basis to inform the design of MPA networks [164][165][166][167].

Computing Infrastructures
Computational capacity is another important limitation when characterizing PAs at the global scale. Most models and processing workflows developed so far are limited by the fact that there is no direct integration with external data sources and models, with most of them being standalone desktop or server applications. In this regard, large computational advances have occurred in recent years based on cloud-based infrastructures that support remote sensing data acquisition and processing [168]. Several  Among other options, Google Earth Engine (GEE; [171]), ArcGIS Online (https://ww w.esri.com/en-us/arcgis/products/arcgis-online/overview; accessed 1 June 2021) and European Copernicus Data and Information Access Services (DIAS; https://www.copern icus.eu/en/access-data/dias; accessed 1 June 2021) offer data and services for cloud-based processing and remote sensing on large scales. Typical environmental applications include detecting deforestation, classifying land cover, estimating forest biomass and carbon or mapping the world's roadless areas [172]. The advantage of using those services lies in the easy data access (including time series), the possibility to create graphical user interfaces and their remarkable computation speed, as processing is outsourced to cloud servers. Moreover, OpenEO (https://openeo.org/; accessed 1 June 2021) allows interoperability with big earth observation cloud backends for several programming languages.

Concluding Remarks and Recommendations
While the methods for mapping and assessing habitats and ecosystems are equally useful within and outside PAs, integrated assessment methods that systematically characterize and measure the diversity of habitats and ecosystems within a region are especially relevant when applied within PAs at a global scale. The global characterization of PAs can provide multiple benefits and applications, such as (1) supporting short-, mediumand long-term management actions, especially at the regional and global scales, that can ensure the maintenance of biodiversity and maximize the provision of ecosystem services [173,174], (2) evaluating the effects of climate change in PAs [175] and (3) informing policy initiatives, such as the European Biodiversity Strategy or the post-2020 Global Biodiversity Framework, on how to develop monitoring tools and indicators to promote sustainable management of PAs [176]. These kinds of analyses not only need to be done at a global scale, but also repeated if possible (e.g., annually) to document the changes that occur [177]. In this regard, the use of variables representing longer-term periods is also useful for capturing the presence of potential habitats and ecosystems, which can then be used as a reference for monitoring and condition assessment purposes. Furthermore, although locally derived variables are better descriptors of the ecosystems, global data sources are needed in order to systematically compare PAs across the globe and inform larger scale conservation actions.
In the last decade, cloud-based infrastructures have greatly improved the access to time series of relevant earth observation variables, which are crucial to the proper monitoring and assessment of ecosystems [178], bringing new opportunities for the global characterization of PAs. However, it is also necessary to translate the results from the global characterization of PAs into information that can be used in the real world, such as by sharing all data and models generated using online interoperable tools [179][180][181][182]. As an example of this, DOPA provides access to various global datasets and indicators that can inform decision-making and PA management [148], such as climate and topographic statistics, information about pressures, the occurrence of extreme events, land cover, land degradation and fragmentation, ecosystem services and species. Moreover, the Protected Planet website allows for exploring the World Database on Protected Areas (WDPA), maintained by the UN Environment Programme World Conservation Monitoring Centre (UNEP-WCMC). The CBD-mandated WDPA is the key reference dataset for any global protected area analysis, and it includes both spatial (mapped boundary or point location) and nonspatial (e.g., name, type, size, age and status) information for over 230,000 protected areas worldwide [183]. Despite accelerated efforts to improve the global PA data, the quality of the WDPA data still varies greatly between countries and regions, and this should be acknowledged in any analysis using the WDPA. Only limited information related to the systematic global biophysical characterization of PAs can be found online thus far, such as the terrestrial habitat diversity index in DOPA [145].
Systematic information related to the uniqueness or the importance of PAs based on biophysical variables could, among other things, further support the ranking and prioritization of PAs based on the diversity of their habitats and ecosystems. Biophysical studies also allow us to study the role of habitats and ecosystems in maintaining biodiversity in the context of climate change, since species populations can adapt to changes by moving to new areas that meet their ecological requirements [146]. Several applications of habitat models have shown a high correlation between biodiversity and the diversity of habitat types, and they can help identify potential new areas that should be protected in order to maintain species protection into the future [120,184,185]. Table 2 gives an overview of applications of different environmental descriptors, including methods and data, which are relevant for the biophysical characterization of PAs, highlighting the importance of taking into account structural and functional attributes as well as integrating a broad spectrum of environmental descriptors in the global biophysical characterization of PAs. Table 2. Summary table with example applications of different environmental descriptors, including data and methods, that are relevant for the biophysical characterization of PAs. Acronyms used: object-based image analysis (OBIA); normalized difference vegetation index (NDVI); normalized difference water index (NDWI); machine learning (ML); principal components analysis (PCA); light detection and ranging (LiDAR); digital elevation model (DEM); normalized difference blue/red ratio (NDBR); wide dynamic range vegetation index (WDRVI); soil-adjusted vegetation index (SAVI); green-red vegetation index (GRVI); plant senescence reflectance index (PSRI); and water band index (WBI).

Application Environmental Descriptors RS and Ancillary Data Methods
Wetland and dune habitat mapping [63][64][65][66][67][68][69][70]73,178] • Finally, we give a summary of the recommendations proposed to improve the global biophysical characterization of PAs in relation to different aspects, which are listed below.
Environmental attributes and descriptors: • The structural and functional attributes of ecosystems and habitats within PAs should be addressed; • A broad set of variables representative of key biophysical quantitative descriptors should be used to produce integrated assessments, potentially including vegetation, energy, climate, water, terrain and soil.
Data sources and processing: • Global data sources presenting time series and regular updates should be preferred; • Dimensionality reduction techniques are often used to deal with correlated input variables; • The use of ensembles of different input data or models corresponding to the same variable is recommended to provide more accurate outputs and deal with uncertainty. Methods: • The use of interoperable RS cloud-based infrastructures is recommended for largescale processing; • Analyses should be regularly repeated to document changes; • The analysis should extend beyond specific habitat or ecosystem mapping and assessment methods so that a variety of habitats and ecosystem types can be identified; • The resulting habitat and ecosystem types within PAs should be, to the greatest extent possible, comparable with existing global typologies; • There is a clear need and potential to develop methodologies for assessing the biophysical uniqueness of PAs that could support prioritization analyses; • Methods should allow the prediction of climate change impacts on ecosystems by using forecasted bioclimatic data.
Application in policy and practice: • Translate the results into information that can be used by policy and decision-makers; • Ensure transparency and reproducibility by sharing all data and models generated using online interoperable tools; • Global characterization of PAs should be specifically aimed at informing larger scale conservation and management actions and plans, unless no better information is available at the local or regional scales.