Research on remote sensing ecological livability index based on Google Earth Engine: a case study from Urumqi-Changji-Shihezi urban cluster

Greenness (NDVI).

NDVI has become the most widely applied vegetation index in remote sensing, used for monitoring surface green vegetation cover (Rouse et al., 1974). In this paper, this index is utilized to represent the greenness within the RSELI. The formula is as follows: (1)${\displaystyle \text{NDVI}=\left({\rho }_{\mathrm{NIR}}-{\rho }_{\mathrm{red}}\right)/\left({\rho }_{\mathrm{NIR}}+{\rho }_{\mathrm{red}}\right).}$

In the formula: ρ_NIR represents the near-infrared band; ρ_red represents the red band.

Wetness (WET).

WET can be calculated through the Tasseled Cap Transformation (K-T Transformation) (Crist et al., 1985; Baig et al., 2014), extracted from MODIS satellite imagery data. In this paper, this index is utilized to represent wetness within the RSELI. The formula is as follows: (2)${\displaystyle \mathrm{WET}={\rho }_{\text{blue}}\times 0.0315+{\rho }_{\text{green}}\times 0.2021+{\rho }_{\mathrm{red}}\times 0.3102+{\rho }_{\mathrm{NIR}}\times 0.1594-{\rho }_{\text{SWIR}1}\times 0.6806-{\rho }_{\text{SWIR}2}\times 0.6109.}$

In the formula: ρ_blue represents the blue band; ρ_green represents the green band; ρ_SWIR1 represents short-wave infrared 1; ρ_SWIR2 represents short-wave infrared 2.

Dryness (NDBSI).

NDBSI is effective for monitoring environmental dryness, calculated as the average of the Built-up Index (IBI) and the Bare Soil Index (SI) (Xu, 2008; Essa et al., 2012). In this paper, this index is used to represent dryness within the RSELI. The formula is as follows: (3)${\displaystyle \mathrm{IBI}=\frac{\frac{2{\rho }_{\text{SWIR}1}}{{\rho }_{\text{SWIR}1}+{\rho }_{\mathrm{NIR}}}-\frac{{\rho }_{\mathrm{NIR}}}{{\rho }_{\mathrm{NIR}}+{\rho }_{\mathrm{red}}}-\frac{{\rho }_{\text{green}}}{{\rho }_{\text{green}}+{\rho }_{\text{SWIR}1}}}{\frac{2{\rho }_{\text{SWIR}1}}{{\rho }_{\text{SWIR}1}+{\rho }_{\mathrm{NIR}}}+\frac{{\rho }_{\mathrm{NIR}}}{{\rho }_{\mathrm{NIR}}+{\rho }_{\mathrm{red}}}+\frac{{\rho }_{\text{green}}}{{\rho }_{\text{green}}+{\rho }_{\text{SWIR}1}}}}$ (4)${\displaystyle \mathrm{SI}=\frac{\left({\rho }_{\text{SWIR}1}+{\rho }_{\mathrm{red}}\right)-\left({\rho }_{\mathrm{NIR}}+{\rho }_{\text{blue}}\right)}{\left({\rho }_{\text{SWIR}1}+{\rho }_{\mathrm{red}}\right)+\left({\rho }_{\mathrm{NIR}}+{\rho }_{\text{blue}}\right)}}$ (5)${\displaystyle \text{NDBSI}=\frac{\mathrm{IBI}+\mathrm{SI}}{2}}$

Heat (LST).

Heat (LST) involves the inversion of land surface temperature from satellite imagery using an atmospheric correction method. The emissivity (REf) of the surface is calculated through the Fractional Vegetation Cover (FVC). The original blackbody radiance values B(Ts), after radiometric correction, provide the radiance at the satellite, which is then atmospherically corrected to remove the effects of water vapor. Finally, the blackbody radiance is converted into land surface temperature using the Planck function (Jimenez-Munoz et al., 2009; Paolini, 2004; Weng, 2009; Nichol, 2005). The formula is as follows: (6)${\displaystyle \mathrm{FV\; C}={\text{NDVI}}_{\left(<0.05\right)}\times 0+{\text{NDVI}}_{\left(>0.7\right)}\times 1+{\text{NDVI}}_{\left(0.05\sim 0.7\right)}\times \frac{{\text{NDVI}}_{\left(0.05\sim 0.7\right)}-0.05}{0.7-0.05}}$ (7)${\displaystyle E_{\left(\text{water}\right)}=0.995E_{\left(\text{building}\right)}=0.9589+0.086\times \mathrm{FV\; C}-0.0671\times {\mathrm{FV\; C}}^2{E}_{\left(\text{natural}\right)}=0.9625+0.0614\times \mathrm{FV\; C}-0.0461\times {\mathrm{FV\; C}}^2}$ (8)${\displaystyle \mathrm{REf}={\text{NDVI}}_{\left(\le 0\right)}\times E_{\left(\text{water}\right)}+{\text{NDVI}}_{\left(0\sim 0.7\right)}\times E_{\left(\text{building}\right)}+{\text{NDVI}}_{\left(\ge 0.7\right)}\times E_{\left(\text{natural}\right)}}$ (9)${\displaystyle \mathrm{B}\left(\mathrm{Ts}\right)=\frac{{\rho }_{\text{TIRS}1}-L_{\mathrm{up}}-t\times \left(1-\mathrm{REf}\right)\times L_{\text{down}}}{t\times \mathrm{REf}}}$ (10)${\displaystyle \mathrm{LST}=\frac{{\mathrm{K}}_2}{ln\left(\frac{{\mathrm{K}}_1}{\mathrm{B}\left(\mathrm{Ts}\right)}+1\right)}.}$

In the formula: E_(water) represents the emissivity of water body pixels; E_(building) represents the emissivity of urban pixels; E_(natural) represents the emissivity of natural surface pixels; L_up represents the atmospheric upward radiance, L_down represents the atmospheric downward radiance, t represents the atmospheric transmittance in the thermal infrared band, and K₁ and K₂ are preset constants for satellite emission. All the above data are automatically acquired from the GEE platform.

Elevation (DEM) and slope (SLOPE).

DEM and SLOPE data used in GEE are from the USGS/SRTMGL1_003 dataset, jointly measured by the National Aeronautics and Space Administration (NASA) and the National Imagery and Mapping Agency (NIMA) of China. The data were collected using the Shuttle Radar Topography Mission (SRTM) system aboard the Endeavour space shuttle, leading to the creation of digital terrain models, which are the current SRTM terrain product data (Aizizi et al., 2023).

Aerosol optical depth (AOD).

AOD is selected to represent air quality conditions, reflecting the total amount of particulate matter suspended in the atmosphere. The MCD19A2 AOD data within the China region are chosen, which are the version 6 data products of the MODIS Terra and Aqua combined, using the Multi-Angle Implementation of Atmospheric Correction (MAIAC) land aerosol optical thickness (AOD) at grid level 2. These can be directly accessed and processed for annual average data compilation using JavaScript code on the GEE platform (Duan et al., 2021).

Population density (PD).

The WorldPop project employs machine learning methods to analyze the relationship between population density and a series of geospatial covariate layers, disaggregating the most recent census-based population counts matched to corresponding administrative units into approximately 100x100m grid cells. The mapping method is based on a random forest algorithm for dasymetric redistribution (Lloyd, Sorichetta & Tatem, 2017), allowing for direct extraction and utilization through the GEE platform.

RSELI model.

Due to the uneven dimensions of the eight factors, directly applying them in PCA would result in unbalanced indicator weights. Therefore, it is necessary to normalize these eight factors before performing PCA, converting each indicator value into a dimensionless value within the range of 0 to 1. This normalization ensures that all factors are on a comparable scale, allowing for a more accurate and balanced PCA calculation (Amani et al., 2020). The formula for this normalization process is as follows: (11)${\displaystyle {\text{XI}}_{\text{i}}=\frac{{\mathrm{I}}_{\mathrm{i}}-{\text{I}}_{\min }}{{\text{I}}_{\max }-{\text{I}}_{\min }}.}$

In the formula: XI_i represents the value after normalization, I_i represents the value before normalization, and I_max and I_min represent the maximum and minimum values before normalization, respectively.

After normalization of the eight factors, the first principal component (PC1) is calculated using the band synthesis and principal component analysis modules within the GEE platform. The formula for calculating PC1 is as follows: (12)${\displaystyle \text{RSELI}=\mathrm{PC}1\left[\text{NDVI},\mathrm{WET},\mathrm{DEM},\text{SLOPE},\mathrm{AOD},\text{NDBSI},\mathrm{LST},\mathrm{PD}\right].}$

The model described above can be compiled and processed using JavaScript code in the Google Earth Engine (GEE). GEE allows for the direct provision of the final results in TIFF format, which can be downloaded for analysis and use. Leveraging the GEE platform significantly reduces processing time due to its cloud computing capabilities, enabling the handling of multi-year remote sensing imagery in one go. This efficiency is particularly beneficial for extensive datasets covering long time periods, making GEE a powerful tool for environmental and ecological research.