Early Identification of Land Degradation Hotspots in Complex BioGeographic Regions

The development of low-cost and relatively simple tools to identify emerging land degradation across complex regions is fundamental to plan monitoring and intervention strategies. We propose a procedure that integrates multi-spectral satellite observations and air temperature data to detect areas where the current status of local vegetation and climate shows evident departures from the mean conditions of the investigated region. Our procedure was tested in Basilicata (Italy), which is a typical bio-geographic example of vulnerable Mediterranean landscape. We grouped Landsat TM/ETM+ NDVI and air temperature (T) data by vegetation cover type to estimate the statistical distributions of the departures of NDVI and T from the respective land cover class means. The pixels characterized by contextual left tail NDVI values and right tail T values that persisted in time (2002–2006) were classified as critical to land degradation. According to our results, most of the critical areas (88.6%) corresponded to forests affected by erosion and to riparian buffers that are shaped by fragmentation, as confirmed by aerial and in-situ surveys. Our procedure enables cost-effective screenings of complex areas able to identify raising hotspots that require urgent and deeper investigations. OPEN ACCESS Remote Sens. 2015, 7 8155


Introduction
Land degradation is a complex problem that refers to the persistent reduction or loss of biological and economic productivity of lands caused by human activities, sometimes exacerbated by natural phenomena [1].The detrimental impacts of these processes on land resources, which are widely recognized as environmental problems of global concern [2][3][4][5], become really dramatic in arid, semi-arid and dry sub-humid regions (drylands), where land degradation is more properly defined as desertification.
In fragile areas where drylands are widespread and the anthropogenic impact is currently increasing, as occurs in the Mediterranean basin, continuous and accurate monitoring strategies are required to assess vulnerability and to single out emerging degradation [6][7][8][9].Because of the inherent bio-geographic complexity of these environments, which is coupled with a large degree of diversity in socio-economic conditions, scientists and decision-makers usually adopt land degradation indices and procedures to synthesize the most important information about the systems involved (e.g., [10][11][12]).As an example, a widely adopted approach in operational contexts (e.g., by EEA-European Environment Agency) for the definition of such composite indices is the Environmentally Sensitive Areas (ESA) methodology [13].Such an approach is very useful for assessing structural vulnerability to land degradation and for identifying the prevailing predisposing causes that can induce degradation.
However, such estimates of vulnerability are not conclusive for the assessment of actual degradation.As an example, within the ESA approach, the Vegetation Quality Index (VQI) takes into account risk factors (e.g., fire risk, erosion protection) that are only based on the land cover type.However, even within regions characterised by similar land covers and VQI, we can find areas where vegetation shows different levels of degradation as well as areas in healthy condition [14].Therefore, the development of rapid methodologies for the early identification of degraded areas and the evaluation of land degradation levels is still a topic of particular interest for both scientists and decision-makers [15][16][17][18][19].
In such a context, remote sensing techniques represent a valuable tool for assessing the actual status of land cover [20].Moderate resolution multispectral time series (AVHRR, MODIS, VGT, MERIS) are generally exploited to assess vegetation response to land degradation at broad scale [16,[21][22][23][24][25][26].Normalized Difference Vegetation Index (NDVI) obtained from these data has also proven to be reliable for identifying vegetation stress in complex Mediterranean landscapes [14,[27][28][29][30][31][32].The high temporal frequency of these series is suitable for characterizing vegetation dynamics and monitoring trends, but not appropriate for detailed scale applications.
On the other hand, high-resolution information on vegetation and soil features can be retrieved from airborne hyperspectral sensors, such as HyMap, AVIRIS, MIVIS, CASI.Plant stress and/or degradation connected to prolonged water deficit, soil erosion and compaction, high salt concentrations, and other forcing can be detected with a high level of accuracy by adopting specific spectral indices (e.g., pigment indices, leaf water content, leaf chemical indices) or Spectral Mixture Techniques (SMA) [33][34][35][36][37]. Similarly, soil characteristics can be retrieved for establishing erosion and salinity levels, carbon, clay or carbonate content [34,38,39].The most recent airborne laser scanners (ASL-Lidar), providing a 3D target reconstruction, can be used jointly with multispectral data and provide detailed information on vegetation structure (e.g., tree height and diameter) and fragmentation [40][41][42].They can also be used for characterizing terrain morphology (e.g., DTM, erosion patterns, gullies) through more consolidated applications [43][44][45].The high accuracy provided by airborne hyperspectral and Lidar sensors is fundamental for a detailed characterization of the current degradation processes at local scale, but inadequate for a continuous monitoring of large areas mainly because of the high costs and the large amounts of data to be processed.
Currently, medium resolution multispectral satellites, particularly the Landsat series, still represent the main data source used for land degradation assessments over large regions and long term periods, especially where no information is available on the specific causes of degradation (e.g., [28,46,47]).Many studies using Landsat TM and ETM+ data have focused on drylands and on soil processes such as erosion and crusting [46,[48][49][50][51]. Very complex procedures, such as SMA, have been adopted to map degraded forests [52][53][54] or to estimate trends of fractional vegetation cover [55].Unfortunately, such approaches are quite difficult to be adopted in operational contexts for supporting local administrations because of the complexity of the procedures that require expert analyses.
For this reason we propose a quite different, straightforward approach with the aim of developing a methodology which should be: (1) useful for monitoring large areas; (2) applicable in operational contexts to identify priorities in complex vegetated regions; (3) based on low cost and accessible data; (4) based on simple mathematical-statistical tools providing easily readable results.
More importantly, it should also account for bio-geographic variability by discriminating intrinsic spatial fluctuations from anomalies due to stress factors.
Our rationale is founded on the contextual effects of processes driving degradation both on vegetation cover and local climate.Irrespective of the trigger causes, the primary sign of degradation is reduced biological productivity (e.g., [18,21,23,56]).Consequently, degradation in vegetated lands should enhance alterations in vegetation properties, such as health, density, etc.By looking at the spatial distribution of NDVI values, we should be able to identify anomalies in the structural and functional characteristics of the vegetation cover and to detect the areas where a given vegetation cover is less productive [48,57].In addition, in order to better understand if a given NDVI anomaly can be interpreted as an actual emerging contribution to degradation in the wider context, we can look at the local climatic scenario where the anomaly is located.Thus, within the target region, the detection of atypical local climates that are not able to support given land covers or are altered by land degradation should help us to select the right hots posts in the region.
Air temperature (T) was selected as a proxy for climatic conditions.Rainfall, which is the main factor for vegetation life, was included as ancillary information because of the limited availability of good quality data for the study area at local scale.It is well known that both repeated drought and land degradation usually result in local warming phenomena [58][59][60] that can exert negative pressure on the existing vegetation thereby triggering a positive feedback toward desertification.Data from ground stations were used in our analyses.We discarded the possibility of using the Landsat thermal band since it is strictly related to the surface conditions at the acquisition time and is not representative of the local climate that is of interest in our context.Moreover, thermal data are highly correlated with NDVI and do not provide substantial additional information to NDVI when used alone.
In our procedure, the bio-geographic complexity of the investigated area and the consequent correlation between NDVI and T are minimized by splitting the investigated region on the basis of the land cover class and by removing effects of geographic factors from temperature spatial variability.The analysis of the distributions of the departures of NDVI and T from their respective mean values estimated per cover class enables us to detect contextual negative NDVI and positive T extremes, which identify areas exhibiting both signs of vegetation stress and atypical climatic conditions that could be associated with emerging land degradation.
Our procedure was specifically designed for early warning.In fact, it focuses on green areas that are still productive.Degradation, if any, is in the emerging phase.In addition, we do not look at trends suited for detecting long lasting processes.We pick up the worst cases for productivity and climate in the area at a given year and in a single NDVI image.Few additional observations are used just to verify persistence in the most negative part of the regional statistics to avoid the detection of casual fluctuations.Within a continuous monitoring strategy, we can repeat the analysis in shifted time windows so as to pick up the time when an area moves toward a critical status.
In this work, we analysed three Landsat TM/ETM+ images acquired on 2002 (June), 2004 (May), and 2006 (July).We searched for common anomalies in images from different years and slightly shifted phenol logical periods in order to separate actual sustained signs of degradation from fluctuations.The low-cost and easy identification of critical areas across extended regions can be used to detect areas that require further investigation (e.g., by airborne hyperspectral and Lidar surveys, field measurements) devoted to the study of the specific processes underlying incoming degradation.

Study Area
The Basilicata region (Figure 1) is one of the areas in Mediterranean Italy that is most vulnerable to land degradation [61,62].In particular, 30.4% of the territory is characterized by a general structural vulnerability to desertification [62] that is common to the whole Mediterranean region due to the region's climate, terrain, soil properties, and ecologically unbalanced human pressure [63][64][65].
Land cover in this region is highly heterogeneous with a distribution that roughly reflects orography (Figure 2).In general, natural forested areas are present at high elevation, whereas anthropogenic land covers largely characterize level areas.The Pollino Massif (up to 2248 m a.s.l.) and the Lucanian Dolomites, in the southern-western part of the region, are covered by dense forests and woods, whose natural values have led to the creation of several protected areas, including the Pollino National Park.The north-east part, bordering on the Apulia region, is characterized by cultivated uplands, where cereals are the main cultivation; in the volcanogenic Vulture Basin vineyards dominate.The central-eastern part of the region predominantly comprises bare clay hills largely affected by erosion [49,64,66].The Metaponto plain, linking the central areas to the Ionian Sea, derives from the expansion of the main rivers (Basento, Bradano, Agri, Sinni, Cavone) over coastal deposits left by sea regression, where salinization phenomena are frequent [67].This plain represents the main agricultural area of the region, with intensive farming that contributes to a gradual impoverishment of the soil.Along the Apennine chain, the thermoregulatory effect of the sea is negligible; thus western and south-western areas of Basilicata experience mountainous climatic conditions, characterized by cold snowy winters and mild summers.In the hilly central areas, the daily temperature range is very high both in summer (18-22 °C) and in winter (10-15 °C), and drought is severe.Along the coast, the climate is typically Mediterranean, with cool wet winters and hot dry summers.The marked bi-seasonality and the large daily temperature range expose land to erosion and are predisposing factors for vegetation stress.Historical data  on rainfall support the dry-sub humid character of the eastern part of the region (e.g., [68]) where most of the vulnerable areas are located, but to date reliable rainfall databases that could be used for high resolution studies are not available.

Landsat TM/ETM+ Data
Three images of the Landsat TM/ETM+ multispectral sensors acquired on 14 June 2002, 26 May 2004, and 19 July 2006 were used for the analysis performed in this paper.
These images were selected in the period of maximum vegetation cover activity, falling late in spring and early in summer, after verifying a scarce presence of clouds.
The imagery was pre-processed by transforming the digital numbers into radiance units and then into reflectance according to the NASA-GSCF calibration coefficients [69,70].To minimize for changes in viewing and illumination geometry and atmospheric conditions the No-Change regression (NC) normalization technique was applied to make the reflectance of the three images comparable.It proved to be the best performing for the heterogeneity of the study area [71].The master image was selected by using the Dense Dark Vegetation (DDV) approach [72]: the image identified as the moststable was the scene of 2004.For the normalized images, the NC provided a correlation of 0.99 and 0.97 for RED and NIR channels, respectively.Then, the corrected channels were used to compute the NDVI as the ratio (NIR-RED)/(NIR+RED).

Meteorological Data
We used meteorological data collected by regional networks distributed in Basilicata.Monthly averages were obtained from daily temperature time series (2000)(2001)(2002)(2003)(2004)(2005)(2006) recorded at 54 ground stations (Figure 3) of the regional networks Lucanian Agency of Development and Innovation in Agriculture (ALSIA) and Regional Agency for Environmental Protection of Basilicata (ARPAB).

Elevation and Land Cover Maps
A Digital Elevation Model (DEM) with a spatial resolution of 20 m × 20 m was provided by the Basin Authority of the Basilicata region (Figure 2a).
Land cover information was extracted from the 2006 CORINE Land Cover (Coordination of Information on the Environment Land Cover, CLC).The map was downloaded from the European Environment Agency (EEA) (http://www.eea.europa.eu/data-and-maps/data/corine-land-cover-2006raster-2).For the purpose of this study, the original CORINE levels 1-3 were recoded to group artificial surfaces and water bodies and to preserve detailed information on the vegetation cover groups.The revised map had 10 classes(Figure 2b In order to overlay both the elevation and land cover maps with the satellite data, we reprojected and resampled all the data to the same projection (UTM 33-WGS84) and spatial resolution of the satellite images (30 m).

Auxiliary Data
Field surveys and aerial photos were used to verify the implemented methodology.Many field surveys were available, which were performed in the framework of the MILDMAP-MEDIA (Methodology Integration of EO Techniques as Operative Tool for Land Degradation Management and Planning in Mediterranean Areas) project from 2006 to 2008 in areas affected by strong erosion.Additional inspections were performed in forested and riparian areas mainly during 2009.Photographic reports on the landscape status were carried out and linked by GPS position and thenimported into a GIS environment.Aerial photographs (1:10,000) were used to interpret and validate the results of the proposed procedure.In particular, a vector map of areas affected by erosion and vegetation fragmentation was digitized from othophotos and then rasterized and overlain on the satellite-based map of critical areas.

Spatialization of Air Temperature Records
Temperature surface sat the Landsat resolution (30 m × 30 m) were derived from point data of air temperature collected at ground meteorological stations.A combined deterministic/stochastic descriptive model, separating the structural component of temperature linked to geography from stochastic fluctuations, was developed.In our model, the contribution of the structural part is expressed by two deterministic components, which model annual mean Ta and seasonal cycle Ts with elevation and distance from the coastline as independent variables, whereas a stochastic component ɛ accounts for local scale fluctuations according to the equation: where t is time (1÷12 months), z is elevation, d is distance to sea.Such a model was previously devised at the Southern Italy scale [73] and subsequently adapted to a smaller scale.Mean annual temperature Ta was described by means of linear regression on elevation (Ta= −0.0053 z) thereby explaining more than 80% of the temperature spatial variance.The monthly temperature seasonal trend Tswas approximated by a Fourier series, estimated from seven years of data for each station; in particular, we retained only the annual component of the Fourier series, because it provided a good approximation of the seasonal periodicity.Amplitude A and phase φ of this component were related to distance to sea 2π ( , ) with A(d)= −0.013d+0.83 and φ(d)= −0.012d − 0.77; other parameters were found to not improve the model.In particular, variations in latitude do not introduce significant variability on scale ranges less than 1 degree [74], which is about the latitudinal extension of the whole Basilicata region.We removed such deterministic components from the monthly temperature values recorded on the acquisition months of the satellite images.Once these components were removed, the variance of the observational data was reduced by 90%.This reduction provides an estimation of the significant contribution of geography to the spatial variability of the data.The residuals obtained (hereafter T) were finally interpolated by ordinary kriging [75].Note that the parameters in Equation ( 1) have to be estimated just once for any given region, after which they can be applied directly in the context of continuous monitoring plans.

Identification of Critical Areas
Our procedure is based on the statistical analysis of NDVI and T per cover class.Once the influence of elevation and distance to sea on temperature was minimized (see Section 4.1) and the residuals were spatialized at the same resolution of the satellite data, the departures of NDVI and T from the means of the relevant land cover class were estimated separately.Since we focused on low values of NDVI, we had to minimize all those spurious effects that could lower NDVI values.A preliminary analysis was performed to remove possible residual effects due to cloud and bare soil contamination (Section 4.2.1).Finally, we proceeded with the identification of critical pixels in three steps: estimation of the distributions of the variables normalized; definition of distribution tails to establish ranges of anomalous values (low NDVI and high T); and definition of a score variable ψ (range [0, 1]) for land degradation (Section 4.2.2).

Removal of Residual Non-Vegetated Effects
We paid attention to negative outliers within the sample distributions of NDVI.Mixed pixels contaminated by residual clouds, water, or bare soil could be erroneously confused for anomalies and therefore they were removed from the statistical distributions.The occurrence of outliers on the extreme left tail (less than the 0.5% threshold) of the NDVI distribution was checked for each land cover and, where they were present, the location of the corresponding pixels was evaluated in the imagery.As an example, some of these outliers corresponded to dam borders and some pixels were very close to areas flagged as cloudy.Thus, we detected water bodies and clouds by applying a threshold (NDVI ≤ 0) to the three images separately.Then, a mask was obtained by merging these three non-vegetated maps.

Statistical Procedure
The core of the procedure can be summarized in the following steps: 1. Estimation of normalized values of NDVI and T for each pixel according to the mean μ and standard deviation σ of the relevant land cover class i: The transformation of NDVI mostly filters out the average greenness of each land cover and the amplitude of its characteristic variability; the transformation of T removes the corresponding mean and variability range from the temperature values observed in all of the areas where the considered land cover is present within the region; 2. Definition of cut-off thresholds for fixing the anomaly ranges.
In order to homogenize the statistical analysis, which should focus on the left tails for NDVI (low values) and on the right tails for T (high values), we changed the sign of NDVI*: So as to define anomalies by focusing on the right tails of both the distributions of −NDVI* and T*.In order to make the method independent of the data distribution, we computed the cut-off by dividing the cumulative sample distributions in percentiles; we selected the cut-off at the 85% probability level, which corresponds to values greater than σ in the case of a Gaussian distribution; according to a partition in percentiles, the lower limit for * and * anomalies will be the 85th percentile (p85 (T) and p85 (NDVI), respectively); 3. Estimation of the product map P between the anomalies of * and which enhances those areas where both NDVI and T deviate significantly from the average behavior, and reduction of the variability range of P to the interval [0,1]: This transformation returns an easily readable score variable, which associates values close to ψ=1 to the highest anomalies of both NDVI and T in the region.Pixels characterized by high values of ψ that persist in time are classified as critical and are the best candidate for the identification of land degradation hotspots.
The scheme of the procedure is illustrated in Figure 4.This assumes that a unimodal distribution is associated to each land cover.Of course, bi-modal or multi-modal distributions could also be found.In this case, some other factors are likely to influence NDVI and/or T thus determining an additional separation within a given land cover class.In this case, we can examine them separately according to the procedure described above.If instead the additional distributions account for few/isolated pixels, their successive processing depends on the target of the requested information; for our purpose, which looks at collective behaviors at the regional scale, they can be labeled as "not-classified" and eliminated from the statistics.Finally, we would stress that a local misclassification (e.g., due to land use or land cover change) in one or two maps should do not affect the final result because anomalies are classified as critical only if they are present in all the investigated images.

Critical Areas
According to the analysis approach described in Section 4, we filtered out leftover spurious (non-vegetated) pixels and focused on the distributions of the residuals of NDVI and T by cover class.The actual need of performing the analyses within the masks of the different land covers to remove bio-geographic effects is well highlighted by Figure5, which reports the maps of the deviations of NDVI and T from the respective means over the whole region.It is possible to point out two macro-areas thatroughly reflect orography and land cover distribution over the region (see elevation and land cover maps in Figure 2).The western part is mainly mountainous and covered by natural vegetation whereas the eastern areas are predominantly plains covered by anthropogenic vegetation (permanent crops and arable lands).Thus, strength and sign of departures are largely due to differences among the average greenness levels of the vegetation classes, to the different phenol logical stages on the acquisition date of the satellite images, and to the different average climatic conditions that characterize these two macro-areas.The W-E separation also reflects the rain availability; by following the UNEP classification, historical rainfall data  show the western part of the region as humid/hyper-humid, with cumulative annual rainfalls generally higher than 1000 mm, and the eastern part as sub-humid/dry, with cumulative annual rainfalls generally lower than 700 mm.As a consequence, we estimated the normalized departures ( * and * ) by cover class after having removed the geographical structure from the temperature field (Section 4.1).Figure 6a,b shows the map of * over the whole region and a detail of the spatial texture of temperature that accounts for both the smooth transitions between slightly different climates and the patchiness due to different land covers.A comparison of Figure 6a with Figure 6c, which illustrates the FAO-UNEP aridity index estimated for the period 1999-2008, shows that most of the anomalous areas are located in drylands.The northernmost area is located in the high-Bradano basin, which is the driest in the region.The central area is next to badlands and the southernmost one is located in the central-eastern part of the Sinni basin, which is dry in the easternmost part.More in general, the good agreement between the aridity map and the map of * confirms the validity of using temperature as a proxy of climate in studies on land degradation in our region.In order to verify the minimization of the bio-geographic influence, we evaluated the correlations between the two measured variables and between * and * .As we expected, it was found that the measured variables were inversely-correlated (r= −0.71) whereas the correlation strength was reduced to r= −0.37 by our filtering.This correlation reduction implies that the filtering mostly removes information that are common to NDVI and temperature and that the joint use of both the variables actually represents an added value to the mere use of NDVI.Figures 7a,b show the sample distributions of * and * under the Forest class for 2002.Both distributions exhibit long tails to the right; values greater than the 85% cut-off correspond to the values we classified as anomalous.Figure 7c shows the anomalous pixels for all land covers identified by considering only the thresholds on * (in blue) and by taking also those on * (in red) into account; more than 63% of pixels anomalous for NDVI were excluded by the inclusion of T. Under the overlain mask (red pixels), the correlation between * and * was further reduced (r= −0.27).
For such an area, we computed the product map (P) and finally the score ψ-map (Figure 8).The analysis pointed out that most of the Basilicata territory does not show signs of degradation.Nevertheless, there is a very large aggregation of anomalous pixels in the south-eastern area (red and dark red pixels in Figure 8), which includes the badland formations of the Matera Hills.Many of these pixels represent arable lands that may have been included within the anomalies by chance (harvesting, meteorology, etc.).In order to detect actual degraded areas, we examined the intersection of the ψ-maps obtained from the three images (2002, 2004, and 2006) that represent slightly shifted phenological periods corresponding to the maximum of photosynthetic activity for the main land covers.Then we identified pixels that were selected as anomalous in all the cases (Figure 9a).The majority of these areas (about 88.6%, Table 1) is covered by forests and a small percentage (6.8%)by natural sparsely vegetated areas.These critical pixels were compared with a map obtained by separating the vegetated classes from the original CORINE shape file and digitizing fragmented areas (used as reference polygons) from orthophotos (Figure 9b,c).The digitizing process was conducted manually by one expert in photo-interpretation, with a deep knowledge of the study area.The aerial photo we used was an orthophoto with a resolution less than 1 m.Possible errors are negligible because the photo resolution is much higher than that of the Landsat image (30 m).The presence of extended erosion is an evident sign of degradation that affects the equilibrium of surrounding vegetation; also fragmentation of structured riparian vegetation can be an important sign of land degradation because of its fundamental role as natural filter to preserve the river ecosystems from the potential polluted water coming from the surrounding agricultural areas [76,77].
In order to better understand the origin of these anomalies, we analyzed orthophotos in detail and examined reports from field surveys performed inside and outside the critical areas (Figure 10a).Degradation in riparian vegetation was found both in the northern (mainly along the Bradano river) and south-eastern parts (along the Sinni, Agri, and Cavone rivers) of the region; the fragmentation of vegetation buffers is mainly present close to agricultural areas (see e.g., Figure 10b).The large patches found in the south-eastern part mainly involve areas covered by non-native vegetation.In the past, such areas were the focus of many reforestation interventions performed to increase the slope stability of badland formations [78].Thinly-planted coniferous-trees (mainly Pinus halepensis, Pinus radiatae, Pinus pinea) were introduced to restore the original natural ecosystem conditions through the formation of organic matter in order to reinvigorate broad-leaved forests.Nevertheless, such interventions were not effective and the conversion to native forest species has been scheduled in the last two plans [79,80].
These mixed coniferous forests are sparse and fragmented by the presence of thin ridges and deep valleys (Figure 11a) caused by strong erosion [64,81].Similar fragmentation processes are also detected within more natural vegetation covers.Here, also broad-leaved forests are broken by the presence of badland formations (Figure 11b).Critical conditions in these areas seem to be quite strictly confined to outcropping badland structures (see, e.g., Figure 12).The fragility of such fragmented forests in the south-eastern areas is exacerbated by rainfall regimes.Local data showed a significant decrease in annual rainfall for the period 1977-2006 compared to the historical mean from 1921 to 2000 (from 600 to 493 mm in Pisticci station; from 578 to 370 mm in Ferrandina station; from 635 to 385 mm in Pomarico station).The reduced water availability combined with the presence of badland formations (high runoff) contributes to an increasingly local warming phenomena.Land cover fragmentation in areas that are not subjected to these erosion processes is generally not sufficient for producing high ψ scores.In fact, many anomalies in NDVI, e.g., simply due to different vegetation density, are not critical because values of temperature are not anomalous, as previously shown in Figure 7c.

Accuracy, Consistency, and Exportability
A correct estimation of the accuracy of our analysis is very difficult because of the lack of spatially detailed information on actual land degradation (i.e., extensive in situ observations) at the scales investigated.The reference polygons obtained from orthophoto (see Sections 3.4 and 5.1) circumscribe areas that show general signs of fragmentation.
The majority of critical pixels (96%) actually falls into these fragmented areas along rivers or in bad land areas.The 4% error comes from mixed pixels mainly located in transition areas between forests and agricultural lands.Although largely acceptable, this commission error could be further reduced by introducing specific heterogeneous land cover classes to estimate the statistics of transition zones.
The scale difference between orthophoto and Landsat image (900 orthophoto pixels approximately match 1 Landsat pixel), determined an overestimation of the omission error.About 19% of non-critical pixels were located within fragmented areas.In all the cases, vegetation largely dominated within the Landsat pixels thereby testifying that the presence of bare soil detected at the micro scale was not sufficient to make biomass production anomalous at the Landsat scale.In spite of the non-critical character of these pixels when evaluated alone, we would note that just this complex alternation of critical/non-critical pixels observed within the reference polygons should be considered an alarming result because this configuration is coherent with the typical spatial patterns of desertification, which generally emerges as a patchy process.
Even neglecting that the estimated error was likely inflated, the overall accuracy of the analysis was rather good (about 87%).
In order to verify the consistency/exportability of the approach, we repeated the procedure by considering the two province administrations (equivalent to the NUTS-3 level, see Figure 1) separately and found that the same core areas were identified (~80% overlap); the main difference consisted in the magnitude of the ψ values.In addition, the procedure was also repeated in random rectangular subsamples of different sizes throughout the region.Our results showed that within sub-samples located in the western part (healthy vegetation), no particular criticalities were identified since very few pixels satisfied the coincidence between NDVI and T anomalies.For boxes located in the southeastern part, the dimension of more critical areas was a function of the size of the sub-sample.Therefore, if we are merely interested in identifying the most critical areas within a region, the procedure is able to locate the relative priority; whereas, if we are interested in the identification of critical areas in absolute terms, weshould focus on an area which is sufficiently large to pick up a complete representation of land cover properties of the region, healthy vegetation included.

Conclusions
We used NDVI (Normalized Difference Vegetation Index) maps from Landsat sensors and interpolated air temperature patterns to search for degraded areas in a vulnerable region of Mediterranean Italy (Basilicata).Differently from usual methodologies developed to analyse large areas, which focus on general vulnerability, we looked for the presence of stressed/sparse vegetation cover in atypically warm areas that could be a consequence/cause of complex interactions between stressed land cover and local climate.Our procedure, applied to a heterogeneous biogeographic area, was able to single out critical pixels (hotspots) showing clear signs of incoming degradation.
A comparison with reference polygons for fragmented areas digitalized from orthophoto showed that the majority of critical pixels (96%) actually falls into fragmented vegetation along rivers or in bad land areas.The 4% error comes from mixed pixels mainly located in transition areas between forests and agricultural lands.Although largely acceptable, this commission error could be further reduced by introducing specific heterogeneous land cover classes to estimate the statistics of transition zones.An apparently larger omission error (19%) was found because of the occurrence of fragmentation that is recognisable at the microscale but is not sufficiently significant to emerge at the Landsat scale as a factor of anomalous biological productivity.Altogether the overall accuracy obtained was 87%.
The general picture described by our results accounts for areas affected by a complex alternation of critical/non-critical pixels that is coherent with the typical patterns of desertification, which generally emerges as a patchy process.The support of rainfall data and aridity maps confirmed the effectiveness of temperature anomalies as a proxy for climate in studies on land degradation in the area.According to these supplementary data, critical pixels are mostly located in drylands; in particular, major clusters were detected close to the bad lands of the Matera Hills, which are still able to sustain productivity in cultivated areas but seem to be not able to guarantee ecological stability in some structured vegetation types (e.g., in forests).Differently, in areas where temperature was anomalous but climate was substantially humid, vegetation did not show significant signs of stress, and very few critical pixels were detected.
Because of its ability to identify anomalous patterns on the basis of accessible data and basic statistical and image processing tools, our procedure can be particularly useful in operational contexts for the systematic monitoring of large vegetated regions and can help to establish priorities in areas characterized by general vulnerability.Here, local detailed investigations based on more expensive tools, such as aerial hyperspectral imagery or field surveys, can be planned for understanding the underlying local phenomena.In order to have a more complete picture of the actual degradation levels, it is recommended to have a reference for vegetation in healthy condition, such as protected areas, within the study region.
The application of this approach in the context of a continuous monitoring strategy could be a useful tool for assessing the evolution of critical areas in land degradation studies.Such a tool could provide valuable and cost-effective support for the evaluation of the efficacy of intervention strategies as well as to assess the impacts of specific legislation promoted to regulate forestry and agricultural activities (e.g., rationale planning of cutting, protection of the buffer areas along rivers, etc.) whose efficacy can result in the reduction of vegetation fragmentation [82].
Our method can be used with any satellite and/or climate data, which is a very appealing feature, especially in the perspective of the next availability of the new Sentinel-2 products that will provide continuity of Landsat-type data [83] at higher resolution.

Figure 1 .
Figure 1.(a) The Basilicata region in Southern Italy; (b) RGB of the region with the separation of two administrative provinces (yellow line) and principal cited toponyms.

Figure 3 .
Figure 3. Localization of the 54 ground stations within Basilicata region (red-coloured dots).Drylands are located in the eastern part of the region.

Figure 4 .
Figure 4. Flowchart of the proposed procedure for the identification of critical areas.

Figure 5 .
Figure 5. (a) Departures of NDVI from the mean value over the whole regional territory (ΔNDVI); (b) departures of T from the mean value over the whole regional territory (ΔT).The figures refer to data acquired in 2002.Similar results were found for 2004 and 2006.Data are rescaled to 8 bit.

Figure 6 .
Figure 6.(a) Map of * for the Basilicata region for the year 2002; white areas within the region's boundaries indicate water bodies or clouds masked according to the procedure illustrated in Section 4.1; critical pixels are reported in purple; (b) zoom in on the black box in (a) to better highlight the spatial texture of the temperature deviations from the mean values of the relevant land covers (coordinates of the centre of (b): 587584E; 4501629N); (c) FAO-UNEP aridity index (1999-2008).

Figure 7 .
Figure 7. Example of 2002 data of the sample distribution of * (a) and * (b) under the forest class (the red line marks the 85% confidence level) and (c) map of anomalous pixels for all the classes: blue pixels correspond to the threshold only on the * , whereas red pixels correspond to the thresholds on both * and * .

Figure 8 .
Figure 8.(a) ψ-map of the anomalous pixels obtained from data acquired on June 2002.A detail on riparian vegetation (coordinates of the centre of the image: 591090E; 4515293N) (b) and on forests (coordinates of the centre of the image: 623501E; 4480242N) (c) highlighting the variability of the function.Red and dark red pixels correspond to the highest anomalies in the region of both NDVI and T variables.

Figure 9 .
Figure 9. (a) Map of critical pixels obtained from the intersection of the ψ-maps (2002, 2004, and 2006); (b) an example of digitalization for the validation procedure: green polygons represent digitalized forests fragmented by badland formations; grey colors represent other Corine classes (principally cultivations); ortophoto is below (coordinates of the centre of the image: 635893E; 4478008N); (c) comparison between the critical pixels identified applying the proposed methodology (in red) and the reference polygons (in green)used as truth map: very few pixels are located outside the reference polygons, they mainly correspond to borders between forests and cultivations (see yellow arrows).

Figure 10 .
Figure 10.Identified critical pixels overlapped on Landsat-TM RGB (a) with the indication of the sites of the performed field surveys (blue stars); (b) example of critical area related to fragmentation of riparian vegetation along the Sinni River (circled area n.2 in 10a) on orthophotos (coordinates of the centre of the image: 640153E; 4447721N).

Figure 11 .
Figure 11.Examples of forested areas fragmented by the presence of outcropping badlands for (a) forest having a low-naturality level (reforestation patterns are clearly visible; coordinates of the centre of the image: 626322E; 4481670N) and (b) mixed natural forests (coordinates of the centre of the image: 630165E; 4457011N), corresponding to circled areas n.3 and n.4 in Figure 10a, respectively.

Figure 12 .
Figure 12.Examples of badland formations affecting the Basilicata territory (photos extracted from the reports of the field surveys).

Table 1 .
Percentage of critical pixels for each land cover type.