Spatio ‐ Temporal Response of Vegetation Indices to Rainfall and Temperature in A Semiarid Region

: In this research, vegetation indices (VIs) were analyzed as indicators of the spatio ‐ temporal variation of vegetation in a semi ‐ arid region. For a better understanding of this dynamic, interactions between vegetation and climate should be studied more widely. To this end, the following methodology was proposed: (1) acquire the NDVI, EVI, SAVI, MSAVI, and NDMI by classification of vegetation and land cover categories in a monthly period from 2014 to 2018; (2) perform a geostatistical analysis of rainfall and temperature; and (3) assess the application of ordinary and uncertainty least squares linear regression models to experimental data from the response of vegetation indices to climatic variables through the BiDASys (bivariate data analysis system) program. The proposed methodology was tested in a semi ‐ arid region of Zacatecas, Mexico. It was found that besides the high values in the indices that indicate good health, the climatic variables that have an impact on the study area should be considered given the close relationship with the vegetation. A better correlation of the NDMI and EVI with rainfall and temperature was found, and similarly, the relationship between VIs and climatic factors showed a general time lag effect. This methodology can be considered in management and conservation plans of natural ecosystems, in the context of climate change and sustainable development policies.


Introduction
To evaluate the effect of temperature and rainfall on vegetation, it is necessary to understand the temporal and spatial dynamics of the vegetation. Furthermore, land use patterns play a very important role in the evaluation of environmental processes for the sustainable management of natural resources and ecosystems conservation [1]. It has been stated that monitoring vegetation dynamics and land cover changes is fundamental for modeling climatic and hydrological scenarios [2,3].
In this sense, the development of remote sensing (RS) has simplified the generation of new methods based on the analysis of satellite images and the use of geographic information systems (GIS). Remote sensing, is information acquisition from a distance, where satellite data are used to infer topography, vegetation cover, system productivity, and other landscape attributes. Remote sensing instruments accumulate the electromagnetic energy reflected by the sun or reflected energy emitted by some device, such as a radar. This information is detected and collected in different ranges of wavelengths (bands), there is a spectral resolution (number of bands and range of band wavelength width) associated to each sensor, and at the same time, a spatial resolution that indicates the area of minimum ground it represents (cell or grid) [4,5]. Remote sensing and geographic information system area. UWRL regression model was used to assess the correlation between the VIs and the climatic variables.

Study Area (SA)
The study area corresponds to the state of Zacatecas, Mexico, with an area of 11141.56 km 2 . It is considered semiarid with an estimated average monthly temperature for the study area of 16.5 °C and average monthly rainfall of 47.45 mm (period 2014-2018) (Figure 1). . The study area is well described by Anuard et al. 2019 [30]. In the study area, the maximum rainfall occurs in the summer. The rainy season extends from May to October, with July, August, and September being the months with the highest rainfall. The dry season begins in November, the month in which the rain is very similar to that of April. The driest months are usually February and March. In December and January there is a slight increase in precipitation due to the cold waves and western circulation. The warm season is from May to September in most of the study area. The warmest months are May and June. October and April are transition months, as the cold season is from November to March, with the coldest being December and January. The current land cover is explained in Table 1.

Image Classification and Accuracy Assessment
The Landsat 8 image classification was created using a maximum likelihood classification (MLC), commonly used for thematic mapping with satellite multispectral imagery [34]. The samples were created by selecting about 80% of the sample class data for the satellite image. Seven thematic classes were selected to represent the land cover; the classes were forest, irrigation agriculture, rainfed agriculture, scrub, urban area, water bodies and grassland.
The supervised classification in ArcGIS 10.2.2 allowed to classify the raster by means of each pixel of a cell, for the seven classes, one works with training samples to obtain a final classification of the image. For the irrigation agriculture the classification was made by analyzing the irrigation calendars in order to choose the month that allowed us to differentiate the irrigation and rainfed agriculture areas. April was chosen and different Landsat 8 images of this month were analyzed over the years of study (2014-2018), resulting in the integration of a single classification, which is presented in this study. In addition, the location of the groundwater extraction wells was considered, as sampling points to irrigation agriculture. As well, the other classes were categorized by the type of vegetation and ecoregions, according to the terrestrial ecosystems in Mexico, from the outlook of National Commission for the Knowledge and Use of Biodiversity (CONABIO) [35].
An accuracy assessment was performed to evaluate the performance of the classification algorithm. It was conducted by creating confusion matrices between the land cover maps and the test data, calculating assessment parameters, such as overall classification accuracy (OCA), producer's accuracy (PA), user's accuracy (UA), and overall kappa coefficient (OKC). The OKC is a measure of overall statistical agreement of an error matrix, which takes non-diagonal elements into account [36]. The Kappa coefficient is calculated from Equation (1) [37]: where N is the total number of samples, ∑ is the total correct samples for all classified phenomena, and ∑ * is the total product determined by multiplying the total number of samples and the number of correct samples per phenomenon.

Satellite Imagery
The Landsat 8 images were processed. Landsat 8 and the operational land imager (OLI) provide seasonal coverage of the global landmass at a spatial resolution of 30 meters (visible, NIR, SWIR), 100 meters (thermal), and 15 meters (panchromatic) [38]. These images used the Worldwide Reference System (WRS), a global notation system for Landsat data. It enables a user to inquire about satellite imagery over any portion of the world by specifying a nominal scene center designated by path and row numbers.
The WRS has been proven valuable for the cataloging, referencing, and day-to-day use of imagery transmitted from the Landsat sensors. The Landsat path/row (29/44) covered the study area; one scene per month from 2014 to 2018 was acquired, and the images were processed to remove clouds and were re-projected into the Universal Transverse Mercator (UTM), Datum WGS 84; zone13N.
The Landsat 8 OLI image was obtained from USGS Earth Explorer [39] as a Level-1 Surface Reflectance product. The reflectance calculation was performed in the ArcGis 10.2.2 software. A top of atmosphere (TOA) planetary reflectance value with correction by solar angle (Equation (2) and (3)) [40,41] was applied.
where ρλ´= TOA planetary reflectance, without correction for solar angle, Mρ = band-specific multiplicative rescaling factor from the metadata, Aρ = band-specific additive rescaling factor from the metadata, and Qcal = quantized and calibrated standard product pixel values. The TOA reflectance with a correction for the sun angle is shown by Equation (3).

Vegetation Indices (VIs)
The most used VI is the NDVI, a quotient that represents the functional characteristics of the active plant that contrasts the reflectance of the near infrared and red bands [42] (Equation 4): where NIR corresponds to the near infrared band, and R corresponds to the red band. By definition, the NDVI values vary between +1 and −1 with higher values for dense vegetation and very low (or negative) values for snow, water and clouds [43,44].
The EVI is similar to the NDVI and can be used to quantify vegetation greenness. The EVI corrects for some atmospheric conditions and canopy background noise and is more sensitive in areas with dense vegetation. It incorporates an "L" value to adjust for canopy background, "C" values as coefficients for atmospheric resistance and values from the blue band. These enhancements allow for index calculation as a ratio between the R and NIR values while reducing the background noise, atmospheric noise, and saturation in most cases (Equation 5). * * (6) Here, is the atmospherically corrected or partially atmospherically corrected surface reflectance, L is the canopy background adjustment that addresses nonlinear and different NIR and R radiant transfer through the canopy, and C1 and C2 are the coefficients of the aerosol resistance term, which uses the blue band to correct for aerosol influences in the red band. The coefficients adopted in the EVI algorithm are L = 1, C1 = 6, C2 = 7.5 and G (gain factor) = 2.5 [45].
The SAVI is used to correct NDVI for the influence of soil brightness in areas where vegetative cover is low. Landsat surface reflectance-derived SAVI is calculated as a ratio between the R and NIR values with a soil brightness correction factor (L) defined as 0.5 to accommodate most land cover types (Equation 6). * 1 The MSAVI minimizes the effect of bare soil on the SAVI. MSAVI is calculated as a ratio between the R and NIR values with an inductive L function applied to maximize reduction of soil effects on the vegetation signal (Equation (7)) [46].
The NDMI is used to determine vegetation water content. It is calculated as a ratio between the NIR and SWIR values in the traditional fashion (Equation 8): where SWIR is the short-wave infrared [14].

Climatological Time Series
A geostatistical analysis was carried out with the climatological database of the National Water Commission (CONAGUA), provided by the State Hydrometeorology Department. The average monthly rainfall and temperature values for the period of 2014-2018 from 25 weather stations were analyzed. Inverse distance weighting (IDW) was used. This is a widely used interpolation method, both for small and large datasets. The unknown value is calculated based on all known points and is inversely proportional to their distances [47].

Bivariate Analysis
The BiDASys (Bivariate Data Analysis System) program was used for the application of ordinary and uncertainty least squares linear regression models to experimental data from the response of vegetation indices to climatic variables. It applies recursive discordancy tests to both models (studentized residuals) for the detection of discordant outliers. The uncertainty weighted linear regression (UWLR) were applied. The UWLR should is considered as a better alternative because the use of uncertainty has a probability connotation with a strict confidence level of 99%, which was used in the present study. The regressions were acquired correlating the average monthly values of temperature and rainfall (obtained by IDW) and the average values of each VI classification [48,49].

Image Classification
According to the classification, the current soil use corresponds to 792.82 km 2

Accuracy Assessment
The classification accuracy is presented in Table 2; the accuracy of the satellite image classifications was approximately 92.68%. The Kappa coefficient was 90.33%. To evaluate the results of land cover, the confusion matrix for the land cover map was created (Table 3). In the table, the column of ground truth represents the accuracy along the matrix and the total numbers of points for each classification. An accuracy of 83.33% was found for the classification of irrigation agriculture, 93.66% for scrub, 91.71% for rainfed agriculture, 96.30% for forest, 95.59% for grassland, 100% for water bodies, and 95.24% for the urban area.

Spatio-Temporal Patterns of VIs and Climatic Variables
The spatio-temporal changes in the VIs were highly seasonal; the evolution of the VIs in the monthly time series showd a parallel dependence on the behavior of vegetation and changes in rainfall and temperature. In 2014 and 2015, the highest VI values were observed in summer and autumn, the lowest in spring and winter. In 2016, the highest values were in autumn and winter, and the lowest were in summer and spring. In 2017, the highest values were in autumn and summer, and the lowest were in winter and spring, and in 2018, the highest were in autumn and spring, and the lowest were in winter and summer.
In some cases, where the rain exceeds 200 mm (for example, June 2015), values close to 0 were obtained, reflected low values in the VIs. The NDMI values indicated a low canopy cover and dry cover (Figures 3-12).          In the rainfed agriculture classification, the VIs showed a very low cover, with values decreasing from the NDVI, EVI, SAVI to MSAVI, in this order of values nearest to zero. The NDMI showed a mid-to-low cover with high water stress (Figures 9 and 10).
In the scrub classification, the VIs showed a very low cover with values decreasing from the NDVI, EVI, SAVI to MSAVI, in this order of values nearest to zero. The NDMI showed a low cover with low water stress (Figures 11 and 12).
There was a positive correlation between the indices and climatic variables which implies that, as the vegetation begins to turn green, the temperature and precipitation also increase and vice versa. The highest values of VIs, rain and temperature correspond to the wet season (summer) and the lowest to the dry season (winter) with spring and autumn as transition months.
The greener classes, such as irrigation agriculture (Figures 7 and 8) and forest (Figures 3 and 4), had the highest VI values, and the lowest VI values were observed in scrub (Figures 11 and 12), grassland (figure 5 and 6) and rainfed agriculture (Figures 9 and 10). However, the temporal variability of the indices, explained by climatic conditions, should be considered.
In the forest classification, the VIs showed mid-to-low canopy cover, with values decreasing from the NDVI, EVI, SAVI to MSAVI, in this order of values nearest to zero. The NDMI showed a mid-to-low canopy cover, with high water stress (Figures 3 and 4).
In the grassland classification, the VIs showed a very low cover, with values decreasing from the NDVI, EVI, SAVI to MSAVI, in this order of values nearest to zero. The NDMI showed a mid-to-low cover, with low water stress (Figures 5 and 6).
In the irrigation agriculture classification, the VIs showed a very low cover, with values decreasing from the NDVI, EVI, SAVI to MSAVI, in this order of values nearest to zero. The NDMI showed a mid-to-low cover, with low water stress (Figures 7 and 8).

Regression Coefficients
Uncertainty-based weighted least-squares linear regression (UWLR) were performed between VIs for each classification versus rainfall and temperature. The overall quality of the regression was evaluated from the linear correlation coefficient values.
In general, in all classifications the best correlations of the VIs of each classification was obtained for the NDMI and EVI, and the lowest correlations were for the SAVI, MSAVI, and NDVI, versus rainfall and temperature. Also, the UWRL regressions with a lag of one month is presented (Table 4).

Discussion
The classification accuracy was acceptable, similar to investigations carried out with comparable methodology [36], and the Kappa coefficient was 90.33%, which indicates a precise classification. The behavior patterns of the VIs in the forests respond to summer rains, where the highest values were observed and that correspond to the growing season, NDVI, SAVI, MSAVI, and EVI increase appreciably during the beginning and the end of the season of growth [50]. Forests have a longer response of lag to other types of vegetation, when they face extreme weather; their own responses to stress and their adaptability make them less influenced by climatic anomalies and have a lag response to prolonged drought due to their deep root system [51]. In the case of grasslands, there is a tendency to decrease in the values of the VIs, for the dry seasons, seasonality of precipitation, and more specifically, the proportion of total precipitation in summer, is related to the productivity of grasses [52], the lack of water can cause the grass to appear dry and straw, which affects the sensitivity of the VIs to the color.
For agriculture classifications a seasonal behavior was found, in rainfed agriculture, it depends totally on the climatic conditions to maintain its greenery, the highest values of the VIs correspond to the summer, where there are the most suitable rain and temperature conditions in the vegetation development; in irrigation agriculture usually the greenery conditions are maintained throughout the year. The scrubs have a resilience to water conditions in these semiarid regions and generally maintain constant VI values, also as a result of greater sensitivity to water stress and greater dependence on annual rainfall, it is notable low VI values in classifications of grassland and scrubland [53].
According to Liu [54]., vegetation tends to regulate the impact of climate change on soil-water dynamics, so that in times of growth, vegetation always grows under a naturally optimized environment. Furthermore, most of the rainfall usually occurs during the growing seasons, and therefore groundwater recharge occurs mainly during these seasons. The results show that the dynamics of vegetation are driven by climatic variability. In winter, low temperatures limit biomass, vegetation growth, resulting in low values [23]. In very particular cases and in conditions of high humidity the VI values decreased, most of the correlations were positive in rainfall and temperature [55].
The spatio-temporal correlations were higher for NDMI and EVI vegetation indices in this study area. The NDVI shows higher values than EVI, SAVI and MSAVI, however it does not have a good correlation with the temperature and the rainfall, according to the results obtained. The EVI was the second best index better correlated, one of the reasons for these results, is due the NDVI has been reported to be susceptible to the spectral influence of soil texture and moisture in gaps between vegetation for desert grass and desert shrubland [56]. The MSAVI, SAVI and NDVI had a weak correlation with climatic variables.
For this research, NDMI represented the best short-term response to climatic variables [57] in the study area and reflects the temporal patterns of water stress in which these ecosystems are found. The NDMI in UWLR regressions, showed the highest levels of correlation (R 2 ) with climatic variables. Similar to Lin et al. (2009), the results indicated that the vegetation dynamics are highly correlated with land surface moisture in an arid environment. This is most likely because the NDMI is based on water absorption, which implies that the land surface moisture content is greatly influenced by vegetation dynamics [58].
Meanwhile, the EVI presented better correlations with the climatic variables than the SAVI, MSAVI, and even NDVI. The EVI represents an improvement on the NDVI, showing a reduced saturation in high vegetation cover regions, a reduction in atmospheric influences and a de-coupling of the canopy background signal. These improvements are based on the introduction of the blue band to reduce the effects of the atmospheric aerosols in the red band, and on some correction coefficients to reduce the effect of soil reflectance. According to these differences, the NDVI is more sensitive to the chlorophyll content, whilst the EVI is more sensitive to the structural characteristics of the vegetation cover [59]. Some studies used the NDVI and EVI in semiarid climates and reported that the EVI can provide greater dynamic range than the NDVI [60]. In this study, the EVI is more useful to predict the spatio-temporal response to rainfall and temperature.
Likewise, most VIs are designed based on the concept of soil-line, which represents a linear correlation between bare soil reflectance at the red and near-infrared (NIR) bands. Unfortunately, the soil-line can only suppress brightness variation, not color differences of bare soil. Consequently, soil variation has a considerable impact on vegetation indices, such as the SAVI and MSAVI [61]. According to Ren and Feng [62] the SAVI and MSAVI may not be adequate to describe green vegetation information in arid and semi-arid ecosystems with low green vegetation cover.
The NVDI is a vegetation index that has demonstrated its usefulness in many studies. However, in some situations, other vegetation indexes might be more appropriate [63] as in this study, where NDVI is limited in its discrimination potential because of the sensitivity to the reflectivity of the soil, for example, in an area with low vegetation density, the reflectivity of a pixel in the infrared band and in the red band would be determined primarily by the ground, with a small variation due to the presence of vegetation. The result is that an index of vegetation in that area would give very similar results to those of bare soil and it would be impossible to detect the presence of vegetation. In fact, this problem is quite serious when the vegetation cover is less than 50%.
The VIs had the highest values in the wet season. In the dry season, they showed a tendency for very low values throughout the study area. This reflects the importance of a relatively cold climate in late winter and early spring, where vegetation most likely performs one of the most important strategies in arid and semi-arid environments and undergoes lethargy or dormancy before the favorable conditions for germination take place.
A better correlation was observed between the NDMI and EVI with the climatic variables than with the MSAVI, SAVI, and NDVI. Furthermore, the R and R 2 were greater within a lag of one month, implying that the vegetation has a lag response to rain and temperature in semi-arid ecosystems. There was abundant evidence that corroborated the spatio-temporal existence of time-lag effects of the antecedent temperature on vegetation growth. The vegetation had short time responses to asymmetric warming. The long lags in semi-arid and arid regions can be attributed to the mechanisms possessed by green plants that allow them to tolerate water shortages resulting from warming temperatures [64].

Conclusions
In this research, the spatio-temporal response of VIs to rainfall and temperature was studied. The relationship between VIs and climatic factors showed a general time lag effect and seasonal patterns in VIs were documented. Apart from the high values in the indices that indicate good health, the climatic variables that have an impact on the study area should be considered. For this investigation, a better correlation of NDMI to rainfall and temperature was found, and similarly, the EVI index showed a good correlation. The MSAVI, SAVI, and NDVI may not be feasible to study the spatio-temporal response to rainfall and temperature in semi-arid ecosystems, as the correlations in this study were not representative.
It is important to consider the unique foliar distribution of vegetation in semi-arid areas, seasonality, and hydraulic capacity of the soil when selecting vegetation indices. Future work should take into account the precipitation as a possibly the more relevant factor for the time lag, as abrupt variation in precipitation between dry and rainy seasons is much more pronounced compared to the gradual changes in temperature during the annual course.
Likewise, in this study, remote sensing and GIS techniques were appropriate tools to understand the VI response to the climatic conditions of the vegetation in a semi-arid area.