High-Resolution Bioclimatic Surfaces for Southern Peru: An Approach to Climate Reality for Biological Conservation

1. Introduction

The use of bioclimatic models to study the ecological niche of species has been consolidated as a research field, since its application can contribute to knowing the potential distribution of living organisms, bringing us closer to macroecology [1,2]. The use of climate models related to the distribution of species can be applied to learn about habitat changes from various points of view (conservation of species, dispersion of disease vectors, colonization of organisms foreign to new ecosystems, effects of climate change, etc.) [3,4,5], which contributes to biodiversity management, or making appropriate decisions for the sustainability of ecosystems. On the other hand, although bioclimatic models are used for the present time, they can also be projected to the past (for instance, [3,6,7]), or to the future (for instance, [8,9]), so that they can be related to other aspects, providing information to understand certain phylogenetic, biogeographic or adaptability aspects of the treated species.

The information that determines the distribution of an organism is provided by its environmental tolerances, intrinsic relationships, and dispersion potential [10]; however, the last two are somehow hard to be known and, therefore, the use of only environmental variables has become extensive [11]. Information such as precipitation, maximum temperature, minimum temperature, humidity, solar radiation, evaporation, among other variables, help to determine the distribution of an organism and to understand biological and evolutionary aspects [5]. However, much of the information collected with specialized instruments is limited to only some variables (which is common in meteorological stations of a certain age), so that some of the environmental variations of certain localities are unknown [12]. Nevertheless, the production of climate models, in relation to some topographic covariates, can help to fill information gaps, which is possible using interpolations through the specific knowledge of variables close to the area and even the use of some satellite indices to correlate climate events [12,13,14,15,16,17,18].

At present, there are worldwide (or global) bioclimatic models that have been calculated with data from ground-based meteorological stations [15], satellite information [19,20] or a combination of both [14]. In some cases, comparisons were made between models generated with ground and satellite information, showing that the latter has a better performance, although it should be applied to sectors with low density of ground stations (such as tropical rainforests); in addition, for sectors with many topographic irregularities, data validation should be performed [21,22,23]. Regional bioclimatic models were produced in order to improve the local climate representativeness, taking into consideration a greater amount of climate information, so as to achieve a more “realistic” representation of the local climatic variability and, therefore, the ecological niche of the species [12,24,25]. Usually, the maximum resolution of bioclimatic models is found at approximately 1 km² of the surface, but can also be found at lower resolutions [3,14,19].

In the case of Peru, and specifically the southern area (Arequipa, Moquegua, and Tacna), no specific climate surfaces have been developed. However, the National Meteorology and Hydrology Service (SENAMHI) has disclosed climate information contained in raster layers with a resolution of 5 km² for Peru. This information, obtained from the combined analysis of data from meteorological and satellite stations [26], though no topographic covariates were used, represents a breakthrough for obtaining these data in Peru. On the other hand, the so-called “global” climate models, although representing the south of Peru, can be considerably biased given the scarce data used for their interpolation, and even the topographic complexity of the area makes the data presented by the “global” climate models even more variable [4,13,14,19,27,28].

Finally, given the current knowledge in the application of the bioclimatic models and the lacking thereof in our environment, we need to understand more “accurately” the climate variability in the south of Peru. Preparing species niche models with the right information will guarantee that their predictions are optimally made and reflect the environmental suitability of the species. Producing ecological niche models that reflect the local reality and that are based on reliable climate information will guarantee better decision making in terms of conservation or even health and environmental safety. Therefore, this study seeks to produce surfaces of maximum and minimum temperature, precipitation, and bioclimatic layers of the southern area of Peru for the conservation of biodiversity.

2. Materials and Methods

2.1. Area of Study

It includes the departments of Arequipa, Moquegua, and Tacna, from approximately 14° 39′ to 18° 18′ south latitude, and from 0 to 6427 m.a.s.l. (Figure 1). The topography is variable because the coast holds the old coastal mountain range, the peaks of which do not exceed 1500 m.a.s.l.; immediately after, to the east, there is the desert coastal plateau, characterized by extensive hyperarid plains, which sometimes reach 2500 m.a.s.l., and which are delimited by the foothills of the volcanic mountain range, the peaks of which easily surpass 5500 m.a.s.l. Finally, crossing the volcanic mountain range, there is the so-called Andean plateau or altiplano, which is crossed by the deepest canyons on the planet and at the bottom of which run turbulent rivers that flow into the Pacific Ocean.

The climate is typically dry, with only two distinguishable seasons: Wet and dry. The areas corresponding to the foothills of the volcanic mountain range and the altiplano have rainfall during the austral summer, but are dry during the rest of the year. On the other hand, in the coastal area, the wet season happens during the winter and the humidity coming from the Pacific Ocean enters the continent as stratocumulus clouds trapped in the spurs of the coastal mountain range or, in some cases, in the last foothills of the Andes (including the valleys that end in the sea). Lastly, the area of the desert plateau receives almost no precipitation and has high solar radiation throughout the year.

2.1. Collecting climate information

The climate information was collected from the National Meteorology and Hydrology Service (SENAMHI) (https://www.senamhi.gob.pe/?&p=descarga-datos-hidrometeorologicos), the National Water Authority (ANA) (https://snirh.ana.gob.pe/ObservatorioSNIRH/), Explorador Climático (https://explorador.cr2.cl/) and privately owned meteorological stations. The data collected were: Maximum temperature, minimum temperature, and precipitation, covering from 1964 to 2018. On the other hand, since not all the stations have the previously indicated time interval, the information thereof was taken into account that they provide at least one full year of records [24]. Given that the climate information in the area is poor, it is necessary to consider that having more records will help to better represent the climate variability of the area, improving the recreation of the produced models [12,24]. Finally, additional information was collected from stations located close to the limits established for the area of study so as to avoid the edge effect [24].

For the coastal area (lomas communities), 20 “virtual” meteorological stations were used and randomly located within the lomas ecosystems [29]. The monthly values of the “virtual” stations were assigned taking into consideration the ratio between the monthly precipitation of actual stations and the NDVI (Normalized Difference Vegetation Index) values of the lomas ecosystems. The NDVI values were obtained from the MOD13Q1 V6 product corresponding to the MODIS satellite, which has a spatial resolution of 250 m every 16 days. This product was selected because of its sensitivity in the desert of western South America [30]. In fact, it is the most suitable, since its products mask clouds, dense aerosols and cloud shadows [30,31]. Images were obtained from https://earthexplorer.usgs.gov, covering all the months of the year, from 2000 to 2018. The images were processed in QGIS v.3.26.2 [32] and monthly NDVI averages were obtained for an 18-year range. Once obtained, the averaged surfaces for each month were cut with a shapefile corresponding to the lomas ecosystem [29]. The localization of the weather stations located in the lomas (Atiquipa and Sama Grande stations) was plotted on the monthly NDVI surfaces, and the corresponding pixel values were extracted. The monthly precipitation and NDVI values extracted from the stations were respectively averaged. The localization of the “virtual” stations was plotted and the respective monthly NDVI was extracted. Finally, the precipitation values of the “virtual” stations were found by proportionally and directly relating the known precipitation according to the corresponding NDVI, all this with respect to the NDVI extracted for the localization of the “virtual” station. The following mathematical relationship was used:

$$\mathit{P}{\mathit{p}}_{\mathit{d}\mathit{m}}=\frac{\mathit{N}\mathit{D}\mathit{V}{\mathit{I}}_{\mathit{e}\mathit{v}}\mathit{x}\mathit{P}{\mathit{p}}_{\mathit{e}\mathit{r}\mathit{m}}}{\mathit{N}\mathit{D}\mathit{V}{\mathit{I}}_{\mathit{e}\mathit{r}}}$$

where:

$\mathrm{Ppdm}$ = Unknown monthly precipitation (mm).

$\mathrm{Pperm}$ = Monthly precipitation of the real station (mm).

$\mathrm{NDVIev}$ = NDVI extracted for the virtual station.

$\mathrm{NDVIer}$ = NDVI extracted for the real station.

It should be noted that the calculation is independent for each monthly NDVI value obtained.

2.2. Georeferencing weather stations

The georeferencing of stations was initially performed with information contained in the SENAMHI and ANA websites. However, some errors (wrong localization) were undetected and, considering that some stations were rather old, it was decided to verify their position taking some considerations into account. First, for the stations that are currently operating, localization was correct in most cases and only a small adjustment was made. Second, localization for old stations was set towards the central square of the village indicated in their information sheets. All georeferencing was performed with high-resolution satellite images (0.3 m) contained in Google Earth Pro v3.4.

2.3. Treatment of climate information

The climate information from every station was stored in spreadsheets. Six daily series of information were obtained and preliminarily treated with filters to eliminate extreme or atypical values. Then, using pivot tables, the data sequence, information gaps, and monthly averages for each year on record (temperatures in °C) or the total sum of monthly values for each year (precipitation in mm) were verified. After verifying the data, homogenization and autocompletion (analysis to detect and correct atypical data, as well as to complete missing data through correlation taking into account latitude, longitude and altitude) were performed using the Climatol v3.1.2 package (https://github.com/cran/climatol) [33] installed in RStudio [34].

2.4. Covariates, modeling of climatic surfaces and validation

We used an ASTER digital elevation model (DEM) from the Ministry of the Environment of Peru (https://geoservidor.minam.gob.pe/), with a resolution of 0.0083°, equivalent to approximately 1 km² (30 seconds), with a WGS84 reference system of longitude, latitude and datum. Also, as recommended some authors [13,14,24,35], given the topography of the area, the covariates used were slope, orientation, cloud cover, additionally including the Topographic Wetness Index (TWI) and the distance to the coast. Slope and orientation were obtained using the DEM, through the modeling implemented in QGIS v3.26.2 [32]. The TWI was obtained using the methodology of [36] and [37]. Cloud cover data were obtained from http://www.earthenv.org/cloud [35]. Finally, the distance to the coast was obtained from https://oceancolor.gsfc.nasa.gov/docs/distfromcoast.

The climatic surface were modeled with ANUSPLIN 4.4 [38]. ANUSPLIN uses data from the different covariates through the thin plate smoothing spline algorithm, producing surfaces with smooth values and without drastic changes. As for the evaluation of the reliability of the produced surfaces, and given the number of meteorological stations, it was decided to use a 10-fold (k-fold) cross-validation. To that effect, first, all the stations were randomly arranged and, then, a percentage thereof was removed so as to finally obtain 10 data lists where a 20% of information was missing (it should be noted that this was done independently, both for maximum and minimum temperature, and for precipitation). Then, the models were produced with each information list obtained; 10 iterations were applied, choosing among them that which presented the lowest value of residues and showed a linear behavior thereof. After choosing the corresponding surfaces, these were evaluated with the data of the stations that were left out of the model, using the root mean square deviation of cross-validation (RMSEcv). Finally, the surfaces produced for each month (corresponding to the 10 folds made) were averaged. A similar process was performed for the RMSEcv. Additionally, the mean absolute deviation of error (MAD) was calculated.

2.5. Comparison with other climatic surfaces for the area of study

The produced surfaces, i.e., accumulated annual precipitation and mean temperatures, were compared to those of WorldClim 2 and CHELSA [14,19] (they were not compared to [20] because their layers had a lower resolution and, instead of using precipitation, they used specific humidity). In this case, in order to better understand the area of study regarding precipitation, it was subdivided into two areas. The first included the coastal area, i.e., from 0 to 1500 m.a.s.l., and the second included the Andean area, above 1500 m.a.s.l. In order to obtain rainfall data, 2000 random points were superimposed in each area and the corresponding information was extracted at them. As for mean temperature, 4000 random points in total were used in the area of study. The Kruskal–Wallis statistical test was used with a 95% level of reliability, so as to detect potential significant differences.

2.6. Production of bioclimatic layers

The 19 bioclimatic surfaces were created from the produced climatic surfaces (12 of maximum temperature, 12 of minimum temperature, and 12 of precipitation) in [34] and using the Dismo package [39]. All of them were developed according to WorldClim and ANUNCLIM [14,15,40]. The mathematical details to create each of the surfaces can be found in [19,40].

3. Results

3.1. Climate information

Information from 118 stations was used. Eighty provided information on precipitation and temperature, 30 on precipitation only and eight on temperature only (Figure 2).

Mean temperature behavior shows that the warm months are January, February, and March (corresponding to the austral summer), while the cold months are June, July, and August (corresponding to winter). The highest temperature is recorded in February (24.8 °C) and the lowest in July (-0.6 °C) (Figure 3A). On the other hand, an individual analysis of the maximum and minimum temperatures shows that the average maximum temperatures in the areas under 1500 m.a.s.l. (coast, desert plateau, and arid valleys) take place in the summer months, February being the warmest with an average of 27.4 °C, while the colder months are June, July, and August, July recording the lowest values with an average of 20.8 °C. As for the areas above 1500 m.a.s.l. (western slopes and altiplano), the warmer months are September, October, and November, the warmest being November with an average of 20.1 °C, while the colder months were June and July, the latter being the coldest with an average of 18 °C. Regarding the minimum temperature, a single pattern is observed, the warmer months corresponding to summer and February being the warmest with an average of 9.7 °C, while the lower values are recorded during the winter, July being the coldest with an average of 3.8 °C.

Precipitation in the coastal area predominantly occurs during the winter months, i.e., in the months of July, August, and September, the latter presenting the highest values (an average of 1.8 mm). It is worth noting that the stations close to hill ecosystems, such as Aquitipa (Arequipa) and Sama Grande (Tacna), have the highest precipitation values (13 and 7.1 mm, respectively). March, April, and May are the driest months in the coastal area, April being the lowest with 0.02 mm (Figure 3B).

In the highlands (above 1500 m.a.s.l.) the rainy season corresponds to the summer months of January, February, March, including December. In contrast, the rest of the year, from April to November, are months with little precipitation. In the Andean area, the month with the most precipitation is February (with a mean of 89 mm). On the other hand, July is the “driest” month, with a mean precipitation of 1 mm (Figure 3C).

3.2. Modeling of climatic surfaces

Thirty-six climate surfaces per month were obtained, corresponding to maximum temperature (12 surfaces), minimum temperature (12 surfaces) and precipitation (12 surfaces) (Figure 4). As for their average maximum temperature values, it can be noted that the warmer months are November and December, November being the one with the highest value. For the same variable, June, July, and August were the colder months, July being the coldest (Table 1).

Thirty-six thousand surfaces were produced, which, after the evaluations and procedures performed, were reduced to 36 corresponding to maximum temperature (12 surfaces), minimum temperature (12 surfaces) and precipitation (12 surfaces) (Figure 4). As for their average maximum temperature values (Figure 5), it can be noted that the warmer months are November and December, November being the one with the highest value. For the same variable, June, July, and August were the colder months, July being the coldest (Table 1).

Regarding the minimum temperature, the warmer months were January, February, and March, February being the warmest. On the other hand, June, July, and August were the months with lower temperatures, July recording the lowest value (Table 2).

As for precipitation, January, February, and March had the higher values, January and February having the highest with a mean of 65.0 mm and 67.3 mm, respectively. On the other hand, the months with the lower values were May, June, and July, the driest being June, with a mean of 2.2 mm, and July, with a mean of 2.0 mm.

On the other hand, the patterns of the produced surfaces in terms of altitude, ocean–inland distance and orientation are recognizable (Figure 6). The temperature surfaces are clearly closely related to altitude (Figure 6A,B) in that it lowers with altitude, and there is no marked relation to the ocean–inland distance and orientation. Precipitation (Figure 6C) shows a certain relation to altitude, although it is only marked for the western part and the vicinity of the great mountains (precipitation increases with altitude). However, the greater the ocean–inland distance (which goes hand in hand with orientation), precipitation also increases. It should be noted that maximum temperature and precipitation in the coastal area have particularities. It can be seen that the maximum temperature notoriously decreases in the area of formation of lomas and then increases (which does not happen with the minimum temperature). Precipitation in the coastal area increases (almost matching rain values around 40 km inland) and then decreases to very low values.

3.3. Validation of temperature and precipitation surfaces.

Values were obtained for each test (Table 4). The RMSEcv range for maximum temperature is between 1.4 and 1.7, with an average of 1.5. In the case of MAD, values are between 1.0 and 1.3, with a mean of 1.1.

The RMSEcv range for minimum temperature was between 0.58 and 1.24, the average being 1.2. The MAD interval was 0.6 to 1.2, with a mean of 0.9.

The RMSEcv for precipitation had a 1.3 to 12.2 interval, the mean being 5.8. The MAD interval was 1.0 to 8.6, with a mean of 4.0.

3.4. Comparison with other surfaces and production of bioclimatic layers

In general, the obtained surfaces indicate more precipitation and higher temperature values compared to those of [6] and [12]. However, in the case of coastal precipitation, on average, an “intermediate” value was obtained compared to [6] and [12]. This was not the case for the Andean Area, where there was a higher value (Table 5). After applying the corresponding tests, it was found that there are statistical differences for the evaluated surfaces. In some cases, there were highly significant differences (p < 0.01). Globally, it can be indicated that the surfaces produced herein are statistically closer to those of [6] (Figure 5), since they presented significant differences (p < 0.05) in two cases.

In general, the obtained surfaces indicate more precipitation and higher temperature values compared to those of [14] and [19]. However, in the case of coastal precipitation, on average, an “intermediate” value was obtained compared to [14] and [19]. This was not the case for the Andean Area, where there was a higher value (Table 5). After applying the corresponding tests, it was found that there are statistical differences for the evaluated surfaces. In some cases, there were highly significant differences (p < 0.01). Globally, it can be indicated that the surfaces produced herein are statistically closer to those of [14] (Figure 7), since temperature did not show a statistical difference and was merely significant regarding Andean precipitation.

Finally, the 19 bioclimatic layers were produced from the 36 climatic surfaces, and are available for free in Zenodo, at: https://zenodo.org/record/6975996#.YvJpcXbMLtR (the 36 produced surfaces are also available for consultation at the same link). The denomination of each is mentioned below: bio1 = mean annual temperature (°C), bio2 = mean daytime range (°C), bio3 = isothermality, bio4 = temperature seasonality (standard deviation), bio5 = maximum temperature of the warmest month (°C), bio6 = minimum temperature of the coldest month (°C), bio7 = annual temperature range (°C), bio8 = mean temperature of the wettest quarter (°C), bio9 = mean temperature of the driest quarter (°C), bio10 = mean temperature of the warmest quarter (°C), bio11 = mean temperature of the coldest quarter (°C), bio12 = annual precipitation (mm/year), bio13 = precipitation of the rainiest month (mm/month), bio14 = precipitation of the driest month (mm/month), bio15 = precipitation seasonality (coefficient of variation), bio16 = precipitation of the wettest quarter (mm/quarter), bio17 = precipitation of the driest quarter (mm/quarter), bio18 = precipitation of the warmest quarter (mm/quarter), bio19 = precipitation of the coldest quarter (mm/quarter).

4. Discussion

4.1. Patterns of the climatic information obtained

As already known, temperature shows some patterns, like its decrease related to the a higher topographic altitude [41] or higher values during the summer season or a drop at higher geographical latitudes [42,43]. The analysis of the temperature data included herein shows the common patterns of temperature variation based on altitude; however, an analysis per season shows that the maximum temperature, unlike the minimum (with high values during the summer and low values during the winter), in the Andean areas (> 1500 m.a.s.l.) takes place in September, October, and November. The aforementioned pattern for maximum temperature was also observed by [44], also indicating that the month with the highest values is November. On the other hand, the reason why September, October, and November are the warmer months in the Andean area would be related to seasonal rainfall. The summer months (including part of December) in the Andean area correspond to the rainy season; therefore, the presence of cloudy masses is much more frequent, thus limiting solar radiation towards the ground surface. This “shadow” created by the clouds, together with the winds created by low pressure, would be responsible for lower maximum temperatures during the summer months. It is worth noting that the variation of the maximum temperature in the Andean areas is not very sharp throughout the year. Like [42], this research shows that this variation is around 2 °C. As for the temperature and longitude relation, as [42] and [45] indicate, it decreases from north to south, but the maximum temperature is sometimes contradictory. In these cases, the location of the station is worth noting; e.g., if it is located in plains or inside valleys. Topographic features model some temperature and precipitation characteristics [13], which, in turn, would potentially interact with wind speed and direction. In the case of valleys, a higher temperature could be caused by the limited movement of air masses.

Precipitation in our environment determines seasons throughout the year, i.e., there are only two seasons: Wet and dry [27,30,31]. As for patterns of precipitation, the latter seems to be related to altitude; however, it is possible to find a certain idiosyncrasy in it [12]. The general precipitation patterns for the western Andes of Peru indicate that precipitation increases with altitude [30] and decreases from north to south [27]; here, however, we will analyze the coastal and Andean areas, specifically. The coastal area of Peru is influenced by three factors defining its wet season: The Humboldt current, the trade winds and the thermal inversion layer in the lower troposphere [30,32], which, during the austral winter, produce very fine rain known as drizzle or misty rain that comes from the stratocumulus clouds and does not go higher than 1000 m.a.s.l. [30,33]. Based on what was obtained, precipitation in the stations located close to the coast is scarce and it increases in those located further inland, where there are hill ecosystems (precipitation occurs even inside the valleys). Probably, the increased precipitation is caused by the clouds “compressing” and producing rain when bumping into the western areas of the coastal mountain range. It is worth noting that the regular wet periods occur during winter. Nevertheless, in December, January, and February there is also a fair amount of precipitation that is due to ENSO events, when the Humboldt current debilitates and, therefore, the ocean becomes warmer, causing heavy rains [32]. As for the Andes, the rainy season is quite noticeable, being limited to summer. Rainfall in the Andean area is produced by geographic and physical factors, the most outstanding of which are the Intertropical Convergence Zone (ITCZ) and the Amazon basin [27,30], which regulate precipitation, decreasing from north to south and from east to west (the latter being related to altitude). All the aforementioned patterns are visible for the area of study, where the higher the elevation, the greater the amount of precipitation (the reference being the constant geographic longitude), while the northern area of Arequipa receives the most precipitation, unlike Tacna (this, considering the same altitude).

Precipitation in our environment determines seasons throughout the year, i.e., there are only two seasons: Wet and dry [42,45,46]. As for patterns of precipitation, the latter seems to be related to altitude; however, it is possible to find a certain idiosyncrasy in it [19]. The general precipitation patterns for the western Andes of Peru indicate that precipitation increases with altitude [45] and decreases from north to south [42]; here, however, we will analyze the coastal and Andean areas, specifically. The coastal area of Peru is influenced by three factors defining its wet season: The Humboldt current, the trade winds and the thermal inversion layer in the lower troposphere [45,47], which, during the austral winter, produce very fine rain known as drizzle or misty rain that comes from the stratocumulus clouds and does not go higher than 1000 m.a.s.l. [45,48]. Based on what was obtained, precipitation in the stations located close to the coast is scarce and it increases in those located further inland, where there are lomas ecosystems (precipitation occurs even inside the valleys). Probably, the increased precipitation is caused by the clouds “compressing” and producing rain when bumping into the western areas of the coastal mountain range. It is worth noting that the regular wet periods occur during winter. Nevertheless, in December, January, and February there is also a fair amount of precipitation that is due to ENSO events, when the Humboldt current debilitates and, therefore, the ocean becomes warmer, causing heavy rains [47]. As for the Andes, the rainy season is quite noticeable, being limited to summer. Rainfall in the Andean area is produced by geographic and physical factors, the most outstanding of which are the Intertropical Convergence Zone (ITCZ) and the Amazon basin [42,45], which regulate precipitation, decreasing from north to south and from east to west (the latter being related to altitude). All the aforementioned patterns are visible for the area of study, where the higher the elevation, the greater the amount of precipitation (the reference being the constant geographic longitude), while the northern area of Arequipa receives the most precipitation, unlike Tacna (this, considering the same altitude).

4.2. Modeling of climatic surfaces and validation

An important aspect to be taken into consideration in modeling is the representativeness of the area, i.e., the quality, quantity and, in this case, the correct geographical position of the input data to perform the interpolations. Some of the produced models use only information from weather stations, satellites or a combination of both [12,14,15,19,20,24]. In tropical countries, the problem of the lack of climate information, either due to the lack of meteorological stations, to the difficult access to information or to incomplete records is highlighted [15,27,28,48]. In this sense, it is possible that the models produced by other authors for the area of study of this research will not properly reflect the climate reality, because only a few ground stations were used or because data were mainly collected by satellite [14,15,19,20], with biases related to topographic complexity and to the density of meteorological stations [27,48]. For this case, the use of as much ground-based climate data available, considering the recommendations by Cuervo-Robayo [24], allowed to obtain “truer” surfaces, better reflecting the temperature and precipitation patterns. It should be noted that the coastal area, unlike the Andean area, has fewer meteorological stations (27% of the total) and that most of them are located near the coast. It is known that the coastal area of Peru has peculiar weather characteristics, with drizzle during the winter, which intensifies inland, towards the windward side of the lomas, creating the lomas ecosystems [42,49]. As previously indicated, the meteorological stations for this area are close to the coast and only two are close to the lomas ecosystems (Atiquipa station, in Arequipa, and Sama Grande, in Tacna).

Given the peculiarities of the lomas ecosystems, and in order to better represent them in the bioclimatic models, it was decided to produce “virtual” stations in the sense that they would have correlated information between the real weather stations (Atiquipa and Sama Grande) and the NDVI index, somehow covering the information “gaps”. The NDVI index, as such, represents health and vegetation cover; therefore, in natural situations, this index is mainly influenced by the amount of moisture in the soil [50,51] and, therefore, of precipitation. In fact, the NDVI-precipitation relation has already been written about, making it clear that they are closely correlated [52,53,54]; although temperature is also related [55,56,57], it seems that, in places where thermal fluctuation is not very wide, this variable does not play a preponderant role [57]. On the other hand, it is also mentioned that the NDVI response, timewise, depends on the nature of rainfall [57]. Thus, in areas where the wet season is sudden, there will be a lack of synchronization between precipitation and NDVI. In turn, in places with a progressive increase of humidity, NDVI peaks and precipitation will be quite close. It should also be noted that, in areas with a marked seasonality, the NDVI-precipitation relation is also close [57]. Then, taking into account the characteristics of precipitation in the formations known as Lomas, it evidently is the marked seasonality and a progressive appearance of vegetation [30,49,58]. In this research, apparently the “virtual” stations created based on the NDVI and precipitation values from real stations was satisfactory, since it shows very significant differences compared to other models. Also, the areas predicted by the model almost perfectly match the lomas vegetation maps [29,30]. In fact, the results show how precipitation occurs in certain sectors of the coastal mountain range (at a certain altitude, slope and orientation), recalling that these ecosystems are specific oases of life surrounded by one of the driest deserts in the world [59]. Anyhow, this is just an applicable solution for these particular ecosystems with closely related climate and a vegetation [60]. Although the performed mathematical operations seem simple, they are a way towards a solution regarding lack of information in the lomas ecosystems. Though in the coastal area of southern Peru there are weather stations which are very close to the coast, these do not reflect the precipitation occurring a few kilometers inland, in the lomas ecosystems. Even during the ENSO (El Niño-Southern Oscillation) events, in stations located close to the coast, precipitation is considerably lower compared to that in the lomas [61]. The veracity of the values found for the “virtual” stations could not be verified here, since the only two stations present in these ecosystems were used for data generation. However, the intention here is that the precipitation conditions in these unique ecosystems are adequately reflected, since other global models only show “flat” values or precipitation mostly restricted to the coastal valleys. The lomas ecosystems are unique in the world, the greatest number of endemic species being found towards the south of Peru [62]. Thus, for a good representation of the niche and distribution models of the species, it is necessary to show them as close as possible to what happens in reality. Further studies are needed to understand the relationship between precipitation and NDVI in lomas ecosystems.

The use of geographic (latitude and longitude) and topographic covariates is common in the modeling of climatic surfaces [14,15,24]. Nevertheless, some others also consider the use of other covariates, such as the cloud cover or the distance from ocean coastline to inland [13,14,63]. Generally, the most used covariate in modeling is altitude [14,15,24] and it is closely related to temperature variations [64]. On the other hand, precipitation, in very broad terms, can show a relation with altitude [45]. However, locally, it is idiosyncratic [19] and will depend on other factors, such as wind currents, topography and the diurnal cycle of solar radiation, these being the most complex in mountain areas [13,28]. In this case, as recommended by some authors [13,24] and considering the topographic complexity, slope, land orientation, TWI, distance to the coast, and cloud cover, were used for modeling, besides altitude, producing good results. In fact, towards the coastal area, where other surfaces show a poor representativeness of precipitation in lomas ecosystems, here they are clearly shown to occur towards the western slopes of the hills (windward).

The use of parameters to measure errors in models is common and necessary, because these parameters indicate de degree of bias or the difference between data collected on site and data interpolated by some model [65,66]. On the other hand, the use of validation values is common to verify the interpolated results in climatic layers with real values collected on site [12,14,24,67]. In this case, a k-fold cross-validation was used, since the reduced amount of available data made it difficult to take a certain percentage for a single training and test [68], which is why it was decided to perform 10-folds with random data and produce models with each of the “sets” of data obtained, finally evaluating the surfaces obtained through the RMSEcv. According to RMSEcv and MAD, the precision of the produced models is adequate and presents relatively low values. Nevertheless, there are some considerations regarding the interpolated variables. Regarding temperature, the RMSEcv and MAD values are higher for the maximum temperature, which indicates more variability, and this is reasonable since the maximum temperature could be influenced by different factors (as explained above), while the minimum temperature is more homogeneous and responds to seasonality [42]. Precipitation, unlike temperature (maximum and minimum) has slightly higher values for RMSEcv and MAD, a pattern also observed in [14,24]. As mentioned before, precipitation is idiosyncratic, with different patterns according to the combination of different factors [28]. This is why the variability of the precipitation models is somewhat high, even more so in places with a complex topography and a poor number of meteorological stations [14]. In this sense, probably for this study, integrating the use of covariates related to topography, humidity and the use of the largest amount of data from meteorological stations favored the low RMSEcv and MAD values. Finally, it is worth noting that the higher RMSEcv and MAD values related to precipitation were obtained in the months of the wet season in the coast (winter) and the Andes (summer), showing once more that it is during these periods when the most variations take place according to the precipitation patterns of each place.

4.3. Comparison with other models produced for the area of study

Statistically, the models produced throughout this research are different from those prepared for the same area, this is, WorldClim and CHELSA [14,19]. The data from the produced layers present highly-significant differences with regard to the CHELSA layers. It is likely that the data obtained for CHELSA [19], satellite-obtained, do not accurately show the temperature and precipitation patterns for the area of study. Usually, specific algorithms related to cloud cover and altitude are used to calculate climate patterns from satellite images, taking into account the analysis of the infrared spectrum in terms of temperature. In this sense, this being a “calculated” process, it needs to be validated with instruments located on the ground, even more so in areas with a complex topography [23]. When compared to WorldClim [14], significant differences were found for Andean precipitation and temperature, and highly significant differences with regard to coastal precipitation. The fact that [14] have obtained more approximate values than those presented here is probably related to the origin of the data, since the authors refer that they collected them from different sources. As for precipitation in the coastal area, the WorldClim and CHELSA models presented highly significant differences compared to the one presented here. This is probably due to the fact that the models produced in this work integrated real and “virtual” values for precipitation in the lomas sectors, which is not very well reflected in [14,19]; in fact, in [14], which is basing its data on [15], it can be observed that no meteorological station is present in the lomas sectors of southern Peru.

4.4. Expectations about the use of bioclimatic layers

The purpose of the produced 19 bioclimatic layers is to provide support for the conservation of flora and fauna species distributed in the departments of Arequipa, Moquegua, and Tacna. Having data that are more representative of the climatic reality of this area will guarantee the production of ecological niche models closer to reality, better outlining the environmental requirements of the potentially studied species, and this not only for endemic species, but also for “problem” species (for example, introduced or parasitic species), or even to detect evolutionary relationships, through the concept of niche conservatism. On the other hand, the surfaces produced here could also be used, for instance, to evaluate the historical biogeography of the species or how climate change could affect them (this, through models projected to the past or to the future, considering the processes or methodologies applied to obtain them). The modeling of the ecological niche (ENM) of the species with layers that are more similar to the climatic reality will allow us to better know the ideal environmental conditions for the species, and, with this, after its projection in the geographic space, we will be able to obtain, in a more approximate way, its geographic and potential distribution (SDM). Finally, it is expected that his work is a source of inspiration to better recreate the environmental and bioclimatic space of Peru, since, as we know, our territory, having a complex topography, needs on-site data for a better climatic representation of its surface and, thus, better understand the ecoclimatic requirements of its species.

5. Conclusions

The analyzed meteorological data present a 54-year amplitude, covering from 1964 to 2018. The patterns obtained regarding the variation of temperature and precipitation are as expected, i.e., that temperature presents an inverse ratio to altitude, a direct ratio to geographic latitude, and values are higher during the summer and lower during the winter. On the other hand, precipitation is idiosyncratic. However, it follows general patterns, which indicate a decrease from north to south and an increase from west to east, which is also related to the increase in altitude. Regarding some specific characteristics, the maximum, minimum and mean temperatures in the coastal and Andean areas are partially as expected, since, in the Andean area, the maximum temperature reaches its highest values in September, October, and November. As for precipitation, in the coastal area, precipitation appears during the winter, while in the Andean area it occurs during the summer. With the climate data, 36 high-resolution climatic surfaces were modeled (12 for maximum temperature, 12 for minimum temperature, and 12 for precipitation, each corresponding to a month of the year), which, according to the RMSEcv and MAD data, are optimal compared to the original data. The comparison between the models produced here and those produced by the other authors, who worked in the same area of study, shows significant and highly significant differences, a greater difference with the precipitation in the coastal area being specified. Finally, based on the climatic surfaces, 19 high-resolution (30 sec ~1 km2) bioclimatic layers were produced for the departments of Arequipa, Moquegua, and Tacna. The 19 produced surfaces will be used to model the ENM and SDM of the species of interest, thus allowing a “true” approach to their environmental and geographical reality, leading to a better understanding of their biology and thus allowing to devise conservation plans or to conduct new researches, based on the nature of the objective or the hypothesis of study.