Vulnerability Evolution and Interactive Mechanism Attribution of Vegetation Productivity in the Jinsha River Dry-Hot Valley Integrating Physics- Guided Neural Networks

Instructions for Data Processing and Code Usage
This archive contains all the data and code for this study. After extraction, it includes two folders: “data” and “code”.
1. Data Processing Workflow
The “data” folder contains all the raw data used in this study, including DEM, LAI, LST, LUCC, NIRv, Precipitation, RVI, SoilMoisture, Temperature, and VPD. Below is a detailed description of the processing methods and workflows for each dataset:
(1) DEM
The DEM data were obtained from the Landsat 8 OLI/TIRS satellite digital product (Geographic Data Cloud Platform, gscloud.cn). The data were resampled to a spatial resolution of 250 m.
(2) LAI (Leaf Area Index)
The LAI data were processed on the Google Earth Engine (GEE) platform: First, the MODIS 8-day composite 500-meter resolution product (MODIS/061/MOD15A2H) was called. The time range from January 1 to December 31 of the specified year was filtered, and the Lai_500m band was extracted. Next, the annual mean for this period was calculated and multiplied by a scale factor of 0.1 to restore the true physical values. Then, spatial clipping was performed using the study area boundary (roi), and the band was renamed to LAI. Finally, it was band-merged with the soil moisture data, resampled to a 250-meter resolution using bilinear interpolation, and reprojected to the EPSG:4326 coordinate system.
(3) LST (Land Surface Temperature)
On GEE, the processing workflow for LST is as follows: First, the ECMWF ERA5-Land monthly mean reanalysis dataset (ECMWF/ERA5_LAND/MONTHLY_AGGR) was called. The time range from January 1 to December 31 of the specified year was filtered, the image collection mean for this period was calculated, and spatial clipping was performed using the study area boundary (roi). Next, the skin temperature band (skin_temperature) was extracted. Since the unit of the original data is Kelvin (K), it was converted to Celsius (°C) by subtracting 273.15, and the band was renamed to LST. Finally, it was band-merged with the VPD data, resampled to a 250-meter resolution using bilinear interpolation, and reprojected to the EPSG:4326 coordinate system.
(4) LUCC (Land Use and Land Cover)
Land Use and Land Cover (LUCC) data were obtained from the 30-m China Land Cover Dataset (Yang and Huang, 2021), which provides consistent annual classifications spanning 1990 to 2025. The annual 30-m raster datasets corresponding to the target study years were clipped to the study area boundary. To ensure spatial alignment with the multi-source environmental variables, the data were then resampled to a uniform 250-m resolution.
(5) NIRv (Near-Infrared Reflectance of Vegetation)
On GEE, the detailed processing workflow for NIRv is as follows: First, a continuous array of years from 2000 to 2025 was automatically constructed based on the study area boundary (roi), and the MODIS 250-meter resolution, 16-day composite vegetation index product (MODIS/061/MOD13Q1) was called. Next, during per-scene image processing, strict quality control and cloud masking were performed using the SummaryQA band to select high-quality pixels with lte(1). Simultaneously, the NDVI and near-infrared bands (sur_refl_b02) were extracted and multiplied by a scale factor of 0.0001 to restore their true physical reflectance values. Then, the Badgley correction formula was applied to calculate the near-infrared reflectance of vegetation (\text{NIRv}=\left(\text{NDVI}-0.08\right)\times \text{NIR}), and physical constraints were implemented by masking out non-vegetated areas with NDVI>0.08, while forcing any negative NIRv values to zero (.where(nirv.lt(0), 0)). Finally, median compositing (.median()) was applied to the filtered image collection to extract annual steady-state features. Combined with spatial clipping for the study area, the spatial resolution was uniformly set to 250 meters and reprojected to the EPSG:4326 coordinate system.
(6) Precipitation
On GEE, the detailed processing workflow for Precipitation is as follows: First, an annual looping sequence from 2000 to 2025 was constructed based on the study area boundary (roi), and the UCSB-CHG/CHIRPS/PENTAD dataset (a 5-day pentad total precipitation dataset integrating satellite observations and ground stations) was called. Next, filterDate() was used to filter all precipitation rasters within the specified year (from January 1 to December 31 of that year), and the precipitation band was extracted. Then, by performing a sum operation (.sum()), the total annual precipitation for that single year was calculated (in single-year meteorological studies, annual average precipitation typically refers to the total cumulative precipitation of that year), and .float() was used to force-convert the data into 32-bit floating-point format. Finally, the precipitation data was band-merged with the temperature data (named Precip), clipped to the study area with missing values unmasked to zero (.unmask(0)), and uniformly set to a 250-meter resolution and the EPSG:32647 projection coordinate system.
(7) RVI (Ratio Vegetation Index)
First, a target array containing 6 key years (2000, 2005, 2010, 2015, 2020, and 2025) was defined based on the study area boundary (roi), and a vegetation threshold parameter of C=0.08, consistent with the NIRv calculation, was established. Next, the MODIS/061/MOD13Q1 product (MODIS 250-meter resolution, 16-day composite vegetation index product) was called, and image collections within the study area for the specified year (from January 1 to December 31 of that year) were filtered using filterBounds() and filterDate().
During per-scene image processing, strict quality control and cloud masking were performed using the SummaryQA band to select high-quality pixels with lte(1). Simultaneously, NDVI, the near-infrared band (sur_refl_b02), and the red band (sur_refl_b01) were extracted and uniformly multiplied by a scale factor of 0.0001 to restore their true physical reflectance values. Then, the ratio vegetation index was calculated using the formula RVI=NIR/Red, and a physical constraint was applied by setting a mask for NDVI>0.08 (C) to exclude non-vegetated areas, ensuring consistency with the input range of the PINN model. Finally, median compositing (.median()) was applied to the filtered image collection to obtain annual steady-state spatial features. Combined with spatial clipping for the study area, the spatial resolution was uniformly set to 250 meters and reprojected to the EPSG:4326 coordinate system.
(8) SoilMoisture
On GEE, the detailed processing workflow for Soil Moisture is as follows: First, the ECMWF ERA5-Land monthly mean reanalysis dataset (ECMWF/ERA5_LAND/MONTHLY_AGGR) was called. Using the specified target year, the time range from January 1 to December 31 of that year was filtered, and the first-layer volumetric soil water band (volumetric_soil_water_layer_1, with a depth of 0–7 cm), which is most sensitive to vegetation productivity responses, was extracted. Next, the annual mean (.mean()) of the image collection for this period was calculated, spatially clipped using the study area boundary (roi), and the band was renamed to SoilMoisture. Finally, it was band-merged with the LAI data, resampled to a 250-meter resolution using bilinear interpolation, and reprojected to the EPSG:4326 coordinate system.
(9) Temperature
On GEE, the detailed processing workflow for Temperature is as follows: First, an annual looping sequence from 2000 to 2025 was constructed based on the study area boundary (roi), and the ECMWF/ERA5_LAND/MONTHLY_AGGR dataset (ERA5-Land monthly aggregate meteorology and hydrology reanalysis dataset) was called. Next, filterDate() was used to filter the monthly mean image collection within the specified year (from January 1 to December 31 of that year), and the 2-meter temperature band (temperature_2m) was extracted. Then, the annual mean condition for that year was obtained by calculating the temporal mean (.mean()). Since the original data unit is in Kelvin (K), it was converted to Celsius (°C) by subtracting 273.15, and .float() was applied to force-convert the data into 32-bit floating-point format. Finally, the temperature band was renamed to Temp, band-merged with the precipitation data, clipped to the study area with missing values unmasked to zero (.unmask(0)), and uniformly set to a 250-meter resolution and the EPSG:32647 projection coordinate system.
(10) VPD(Vapor Pressure Deficit)
On Google Earth Engine (GEE), the detailed processing workflow for VPD is as follows: First, the ECMWF ERA5-Land monthly mean reanalysis dataset (ECMWF/ERA5_LAND/MONTHLY_AGGR) was called. The time range from January 1 to December 31 of the specified year was filtered, the image collection mean for this period was calculated, and spatial clipping was performed using the study area boundary (roi). Next, the 2-meter temperature (temperature_2m) and 2-meter dewpoint temperature (dewpoint_temperature_2m) bands were extracted and converted from Kelvin (K) to Celsius (°C) by subtracting 273.15, respectively. Then, using Tetens' formula via ee.Image.expression(), the vapor pressure deficit was calculated, and the resulting band was renamed to VPD. Finally, it was band-merged with the LST data, resampled to a 250-meter resolution using bilinear interpolation, and reprojected to the EPSG:4326 coordinate system.
2. Instructions for Code Usage
The code folder contains two subfolders: GEE Scripts and Python Scripts. The GEE Scripts folder stores the JavaScript codes used on the GEE platform to acquire raw data, where each script's name corresponds to the dataset it retrieves. The Python Scripts folder contains codes for analyzing the raw data, consisting of the following:
(1) Annual Mean NIRv.py
This script is used to batch-read multi-temporal 250-meter resolution NIRv (Near-Infrared Reflectance of Vegetation) raster images from a specified directory, calculate the annual regional mean, and perform time-series statistics, corresponding to Fig. 1 in this study.
(2) spatiotemporal_evolution_analysis.py, Trend Analysis.py, coefficient of variation (CV).py
First, running spatiotemporal_evolution_analysis.py yields three raster files: NIRv_Trend_Slope.tif, NIRv_Trend_Pvalue.tif, and NIRv_Stability_CV.tif. Building upon these, running Trend Analysis.py generates the Spatial distribution of vegetation productivity (NIRv) evolutionary trends (Fig. 5), while running coefficient of variation (CV).py produces Fig. 6.
(3) PGNN.py
PGNN stands for Physics-Guided Neural Networks. Running this script produces the model comparison results and physical parameter values presented in Table 1 of the paper.
(4) residual analysis.py
Running residual analysis.py yields the Comparison of residual kurtosis and skewness between PGNN and baseline models, corresponding to Fig. 7 in this study.
(5) SHAP-Summary Plot.py
Running SHAP-Summary Plot.py generates Fig. 8 in this study.
(6) SHAP-Dependence Plot.py
Running SHAP-Dependence Plot.py generates Fig. 9 in this study.
(7) SHAP-Interaction Plot.py
Running SHAP-Interaction Plot.py generates Fig. 10 in this study.