Mapping global eco-environment vulnerability due to human and nature disturbances

Graphical abstract

If applicable, include full bibliographic details of the main reference(s) describing the original method from which the new method was derived. The current method [1] can be considered as an extended work of our previous research outcomes [2,3] where we performed eco-environmental vulnerability assessment (EVA) at a provincial scale of Vietnam. The innovation in our current research [1] is to provide the first global-scale map of eco-environmental vulnerability. The six global environmental vulnerability (GEV) categories and their map of global "hotspots" are novel and will be proven highly useful to researchers and decision makers around the world working on the issues of sustainability, conservation, development, and vulnerability. Subsequently, three eco-environmental zones are introduced with functions and advices for activities and planning for the regions of concern. Considering both natural and manmade disturbances, GEV analysis can be of significant value in (i) Enriching the guidance of global and regional planning and construction, and protection of the ecological environment; (ii) Harmonizing information from the reports that employ different approaches or definitions; (iii) Providing a feasible framework template of EVA, which benefits environmental education; and (iv) Conveying information to the public for enhancing the role of communities in solving environmental problems. If applicable, include links to resources necessary to reproduce the method (e.g. data, software, hardware, reagent)

Method details
We propose a framework for evaluation of eco-environmental vulnerability at global scale by using GIS techniques and freely accessible datasets. However, our method and framework can be applied not only to any environmental issues due to any specific natural stresses such as forest fire, typhoon, flooding, drought or anthropogenic disturbances such as pollution, industrialization, urbanization, but also to any region of concern. There are many influential factors that can affect eco-environment. Choosing the right set of indicator and mapping procedure are crucial for achieving meaningful and reliable quantitative results.
In this study we extend our previous works conducted over a regional scale at Thua Thien-Hue Province of Vietnam [2,3] by upscaling the assessment framework. The upscaling was performed by Table 1 Indicators used to evaluate GEV including their sources, data description, and preparation, and a brief explanation of their roles.

Major disturbance determinants Indicators
Role in environment profile Sources This indicator shows average income of each country from high to low income (highly-developed countries to developing countries). In general, in the developing countries, the eco-environment is likely to be disturbed more than developed countries since they are on the fast growing processes of urbanization and industrialization. Income also reflects the education level as well as public awareness of eco-environmental protection.
https://data. worldbank.org [7] Distance from urbanized areas (C 7 ) This indicator determines the influence from the urban by spatial distance. Exposure from urban affected the eco-environment by the stress from the city like pollution from vehicles and air-condition, and trash from households, and wastewater. It is likely that the farther from the urban the better the eco-environment. considering influential factors retrieved from global datasets across five domains (natural hazards, society-economics, hydrometeorology, and topography and land resource) that aim to visualize the nature and human disturbances on eco-environmental vulnerability. Nguyen et al. (2016) proposed an assessment framework to evaluate the eco-environmental vulnerability in association with 16 variables with six of them constructed from Landsat 8 satellite image products. 16 variables were taken into account and organized into four domains: (1) hydro-meteorology, (2) socialeconomics, (3) land resource, and (4) topography. This framework is relied on remote sensing data, digital maps and in situ measurements with aid of AHP and GIS. In this framework, the weights are stressed on society-economics and topography (0.329) followed by hydrometeorology (0.2) and land resource (0.142). Details of weights and applied AHP process were given in Table 1  Globally speaking, it is almost infeasible to process all the data and variables at uniformed 100 m of resolution when, in particular, input variables cannot be derived from one kind of dataset or single satellite image as Liou et al. (2017) did. Thus, for a study over a globe scale, we further improved the assessment framework. Totally, 16 indicators were derived from ten global datasets (Table 1) and one global dataset of PM 2.5 was used for validation. In the global improved framework, we considered one more domain, i.e., natural hazards. In this domain, we took into consideration of four kinds of dominated hazards, including flood, drought, landslide and cyclones. The influence of the four kinds hazards were considered equally with weight of 0.25. Thus, the global weight was revised in comparison with the previous framework with criteria that can capture natural and human disturbance, as presented in Table 3 with summary given in the diagram in Fig. 2. All indicators are computed or converted and resampled to be of same GEOtiff format and resolution of a pixel size 0.0833 Â 0.0833 . The size of the globe for this resolution contains approximately 4320 columns and 1673 rows. All GEOtiff raster files were projected to spatial reference of WGS-84. Processing steps were mainly performed in ArcGIS 10.3 environment.
The selection of influential factors and spatial data processing here are optimized to identify vulnerability levels of eco-environment at global scale due to nature and human disturbances. The Slope constraint is a factor influencing land-use decision and the item "Land utilization possibilities". The influence of terrain on erosion is great important. Steeper slopes are also associated with shallower soils in general and with a higher risk for soil degradation and landslides [13][14][15].
http://www.fao.org/ geonetwork [12] Slope aspect (C 16 ) Slope aspect and topographic position contribute to define annual mean temperature, potential energy incoming and evapotranspiration. Resulting in vegetation structure, ground moisture, snow retention, plant communities and surface temperature are all characteristics influenced by aspect [16]. influential factors can be replaced and adjusted based on the condition of the region of interest. For instance, for the region with dominated disturbance of forest fire or earthquake, the natural hazards can be replaced by the named driving forces. The global weights are derived in detail steps as given in Tables 3-6 with summary of weight results as given in Table 7. Class weights are given in Table A1 in Appendix Supplementary material. Similarly, the local weights were derived with consideration of five domains. Results of local and global weights are all summarized in diagram (Fig. 2). Data are downloaded and processed in steps as depicted in Fig. 1 and detail processing steps are described in the following: * Step 1: Data preparation Details of data preparation is described one by one as seen below: (C 1 ) Soil moisture: Annual average soil moisture is retrieved from L-band microwave missions SMOS (Soil Moisture and Ocean Salinity) as output in NETCDF format, further converted to GEOTIFF format, and then mosaicked and classified into eight classes over a global scale of 0.25-degree grid. Soil moisture is provided by INRA (Institut National de la Recherche Agronomique) and CESBIO (Centre d'Etudes Spatiales de la BIOsphère). We interpolated the areas with insufficient data or data being removed due to unsuitable values of shallow effects.
(C 2 ) Precipitation: Global annual average precipitation (mm) with a resolution of 0.00833 from a GEOTIFF file is classified into eight classes. The precipitation data is interpolated using thin-plate splines with covariates including elevation, distance to the coast. Weather station data are used between 9000 and 60,000 stations with temporal range of 1970-2000 [5].
(C 3 ) Temperature: Global average temperature with a resolution of 0.00833 from a GEOTIFF file is classified into eight classes. The data are processed by using the same weather station and method as precipitation data [5].
(C 4 ) Distance from hydrological network: By using shape file formatted inland water surface data, we calculate the distances of interest by a Euclidean Distance tool in ArcGIS, which are then further classified into eight classes.
(C 5 ) Population: Population data are in excel file format, stored into a shape file, and further converted to GEOTIFF format. Finally, population is classified into eight classes.
(C 6 ) Income: Income data are in excel file format, stored into a shape file, and then further converted to GEOTIFF format. Raster of income is classified into five classes.
(C 7 ) Distance from urbanized areas: Urban areas are in shape file format. We calculate the distance by a Euclidean distance tool in ArcGIS, which is further classified into eight classes from near range to far range from the urbanized areas. (C 8 ) Land use/land cover (LULC): The MODIS Land Cover Type product in 2017 (Short Name: MCD12Q1) of 500 m SIN Grid in HDF format is download and further processed such as mosaicked and converted into GEOTIFF format and then classified into 16 classes following the legend and instruction of providers. References include USGS's website and following papers [18][19].
(C 9 ) Normalized Difference Vegetation Index (NDVI): Global MODIS vegetation indices, NDVI product Vegetation Indices 16-Day L3 Global 500 m are downloaded, mosaicked, and used to compute mean value NDVI of the year 2017, and then further classified into five classes.
(C 10 ) Drought; (C 11 ) Tropical cyclones; (C 12 ) Landslides; (C 13 ) Flood: These indicators are downloaded either in an excel file, shape file or raster GEOTIFF format, then further processed into same format of GEOTIFF at resolution of 0.0833 over a global scale, and aggregated by weighted sum function in ArcGIS (each type of natural hazard has the same weight). Final raster of natural hazards is classified into five classes.
(C 14 ) DEM: SRTM DEM in GEOTIFF file is downloaded with a resolution of 0.00833 and further processing is conducted such as filling holes, and classification into eight classes.
(C 15 ) Slope constraint: Slope constraint with 0.0833 resolution is downloaded and converted into GEOTIFF file and classified into 8 classes from low constraint to very frequent constraint.
(C 16 ) Slope aspect: Slope aspects are computed by using ArcGIS function with input SRTM DEM and then further classified into 10 classes. PM 2.5 -validation data: An annual global map of PM 2.5 is derived from MODIS data [20][21]. *Step 2: Derivation of global weights Analytical Hierarchy Process (AHP) AHP is an effective tool for dealing with a larger number of influential variables for making decisions in a structural way. Take for example risk, hazard, and eco-environmental impact in a multiindicator decision making process to derive vulnerability assessment [22][23], AHP offers the analysis of the spatial multi-indicator layers through the creation of a hierarchical structure by providing weighting and ranking with the views of expert opinions and users [23][24]. Here a weighted overlay technique is used to synthesize weighted and ranked spatial variables together. AHP has been widely used for mapping hazards, vulnerability, and risk of various natural disasters such as floods, landslides. Hence, it is considered suitable for global eco-environmental vulnerability assessment [25][26]. In this study, for weighting the selected indicators of eco-environmental vulnerability, the individual pairwise comparison matrix was established through the qualitative judgments of five experts, and a user, each of whom considered the influential factors and alternatives, which were integrated with every indicator. For prioritization of the indicators, experts were asked based on the pairwise comparison 9 point scales ( Table 2) developed by Saaty (2008) [22]. Five experts were selected at the international and national levels. The expert selection was based on their related depth knowledge and research experiences about the influential factors. The selected experts were from research institutions, and governmental and academic sectors.
The consistency ration (CR) was calculated to justify the consistency of comparisons given by experts and the user in the pairwise comparison matrix. The comparisons in the pairwise comparison Experience and judgement slightly favor one indicator over another 5 Strong importance Experience and judgement strongly favor one indicator over another 7 Very strong importance One decision indicator is favored strongly over another and its supremacy is established in practice 9 Extreme importance The evidence favoring one decision indicator over another is of the highest possible order of validity 2, 4,6,8 Intermediate values between the two adjacent judgements Compromise is needed  where Random Index (RI) represents the randomly generated average consistency index (CI), which is defined as follows: where l max refers to the largest eigenvalue of the matrix and n represents the order of the matrix (1).
In addition, the Natural Break Statistical method in ArcGIS environment was used during the spatial analysis to classify the indicator layers where it was required. All the classification of indicator layers and internal weights are explained in Table 7.  Table 7 Weightings of group indicators and indicators used for the calculation of global eco-environmental vulnerability (modified and adapted from (1,15). Consistency ratio of assessment is 0.007. Class weights and consistency ration of each indicator are provided in Table 6.

Derivation of global weights
For determining the weights of the global group variables, the individually pairwise comparison matrix was established by using the 9-point scale ( Table 2) of Saaty (2008) as shown in Table 3 that is further processed to obtain the pairwise comparison matrix of group variables as shown in Table 3. The table of normalized pairwise comparison matrix, weights, and consistency ratio (CR) is then given in Table 5 followed by Table 6 showing the consistency of judgments. Two types of consistency  methods, including geometric mean method and row average method, are used to check the consistency of our judgments, requiring the consistency ratio to be less than 0.1. From Table 6, it is seen that our judgments are acceptable. Under such circumstance, the derived global weights are adopted for further analysis. *Step 3: Classification of indicators and derivation of class weights Similarly, the class weights are derived as shown in Table 7. *Step 4: Mapping four major determinants and final global eco-environmental vulnerability maps  An analytical hierarchy process (AHP) and geographical information system (GIS) are implemented to combine multi-indicators in groups and then further aggregated to become one final indicator of GEV by using Eqs. (3) and (4):  where GEV denotes the global eco-environmental vulnerability (the higher the GEV value, the greater the vulnerability is likely to be), B i is the ith group determinant factor, W i is the weight of the ith group determinant factor, CJ is the jth indicator, wj is the weight of the jth variable, and nB i is the number of indicators in a group determinant factor B i introduced in Table 1 and Step 1 data preparation. Weights of 16 indicators and five groups are presented in the Table 7 and Figs. 3-7. To classify vulnerable intensity, the GEV is standardized and compared. In this study, we use histograms to reveal the statistical distribution corresponding to values of grid cells of eco-environmental vulnerability raster to classify GEV assessment into six categories, namely very low (< 1.66), low (1.66-2.13), medium (2.13-2.56), medium high (2.56-3.06), high (3.06-3.59), and very high (>3.59). *Step 5: Analysis of the result and validation After obtaining the maps of four major influential factors (Fig. 5) and global eco-environmental vulnerability map, the spatial analysis distribution is necessary to understand how it exhibits across to continents with influential factors. Spatial Analysis Tools including, Zonal Histogram, Tabulate Area and Zonal Statistic in ArcGIS are applied to investigate the statistical features of these maps. For example, Tabulate Area distribution of eco-environmental vulnerability over the land cover types is shown in Fig. 4. Individual indicators are presented in Figs. 6-10.  Validation is crucial to check the reliability of the results. A global map of PM 2.5 is derived from MODIS data [21] and considered as an independent variable for validation with global ecoenvironmental vulnerability map since it can be considered as an anthropogenic disturbance associated with nature and human-made influence. We choose 100 points randomly distributed over the globe by using the Random Points and Pixel Value to Points functions to get the attributing values of global eco-environmental raster and PM 2.5 raster. Analysis of correlations between EVA map results and PM 2.5 pollution distribution is conducted. It is found that the correlation coefficient reaches 0.82 as shown in Fig. 3.