Recent vegetation phenology variability and wild reindeer migration in Hardangervidda plateau ( Norway )

More than others, arctic ecosystems are affected by consequences of global climate changes. The herbivorous plays numerous roles both in Scandinavian natural and cultural landscapes (Forbes et al., 2007). Wild reindeer (Rangifer tarandus L.) herds in Hardangervidda plateau (Norway) constitute one of the isolated populations along Fennoscandia mountain range. The study aims to understand temporal and spatial variability of intraand interannual home ranges extent and geophysical properties. We then characterize phenological variability with Corine Land Cover ecological habitat assessment and bi-monthly NDVI index (MODIS 13Q1, 250 m.). Thirdly, we test relationships between reindeer’s estimated densities and geophysical factors. All along the study, a Python toolbox (“GRiD”) has been mounted and refined to fit with biogeographical expectancies. The toolbox let user’s choice of inputs and facilitate then the gathering of raster datasets with given spatial extent of clipping and resolution. The grid generation and cells extraction gives one tabular output, allowing then to easily compute complex geostatistical analysis with regular spreadsheets. Results are based on reindeer’s home ranges, associated extent (MODIS tile) and spatial resolution (250m). Spatial mismatch of 0.6 % has been found between ecological habitat when comparing raw (100m2) and new dataset (250m2). Inter-annual home ranges analysis describes differences between inter-seasonal migrations (early spring, end of the summer) and calving or capitalizing times. For intra-annual home ranges, significant correlations have been found between reindeer’s estimated densities and both altitudes and phenology. GRiD performance and biogeographical results suggests 1) to enhance geometric accuracy 2) better examine links between estimated densities and NDVI.


Introduction
In the information age, « big data » constitutes a large part of the future of ecology and biogeography (Hampton et al., 2013).Transdisciplinary expertise for socioecological issues are nowadays relying on multi-factorial and multiscale analysis for environmental surveys.Massive datasets are easily computed and processed, thanks to tremendous computer hardware advances and trends for horizontal interactions between experts, citizens and institutions (Kambatla et al., 2014).That said, the era of data-intensive science needs some specific tools to collect and centralize the numerous and heterogeneous datasets related to a studied phenomenon, such as wildlife and herding issues, interacting with cultural and natural systems.Multi-temporal scalar information is another issue, long-term strong trends (climate change, biodiversity loss) interacting with interannual variability.Within the scope of environmental impact assessment, and particularly in the field of animal ecology, an equal spatial resolution of biotic and abiotic factors is required to describe, explain and forecast how a population is supposed to move in a limited range of space, according to future bioclimate and biocenoces interactions.Arctic and subarctic ecosystems are particularly affected by the consequences of global climate change (Uboni et al., 2016).Rangifer tarandus tarandus L. constitutes a key-component species, and its role in both natural and cultural sphere is increasingly important.With a wide circumpolar distribution, the capital-breeder and longdistance migratory ungulate moves through taiga and tundra biomes.Its particular geographical distribution and its ecological habitats selection are relevant to ecosystem structure and functioning.The ability of Rangifer to cope with its food needs and spatial requirements may represent a good opportunity to increase arctic and subarctic socio-ecosystems resilience to rapid effects of climate change (Post and Pedersen, 2008).At a global scale, populations of wild reindeer are classified as vulnerable by the International Union for Conservation of Nature (UICN), and currently red-listed by the institution.Warmer winter temperatures, rapid changes in the water cycle, increasing land-scape fragmentation and human disturbances contribute to deplete populations, biocenoces communities and individuals (Weladji et al., 2006).On the Hardangervidda plateau (Norway), the subpopulation of Rangifer tarandus tarandus L. is one of the last remaining wild reindeer population in Scandinavia.The population stock has been reduced sevenfold between early 1970s and today (Uboni et al., 2016).With Rondane and Snøhetta sub-populations, those three biogeographical isles could represent a genetic resource contributing to Scandinavian Rangifer resiliency to climate change and related cultural and natural landscapes.As such, it is important to improve our understanding of the capital-breeder ungulate migrations during the vegetation growth period and the critical times in its life cycle: response to spring snow melt, calving period and "fat accumulation" period in late summer (Klein, 1990;Skarin, 2008).Our hypothesis are the following: Rangifer tarandus tarandus L. has a seasonal behavior, choosing different habitats according to the key moment of its biological cycle during the early spring and the growing season (Skarin, 2008); strong estimated densities of Rangifer tarandus tarandus L. depend on the bioclimatic conditions during its migration, and on geographical features such as topography and habitats (Klein, 1990).Management of merged datasets, here by spatializing relationships between biocenoces and biotopes, may improve GIS toolbox as well as our understanding of Rangifer tarandus interactions.The study aims to present a GRiD-toolbox, combining some pre-existent tools in current versions of ArcGIS and QGIS software, applied to a multi-scalar reindeer home ranges study.

Material
Several free-access datasets, fully described in Table 1, are compiled, merged and analyzed in the study.Downloading datasets is the first step (url access in Table 1).The first basis is the reindeer GPS-tracks database used by Cagnacci et al. (2015) to assess and characterize most European's largest herbivorous seasonal residence times, associated home ranges and identifying migrational structures.7 wild female reindeers have been equipped and surveyed using GPS-radio collar.The second dataset is ASTER-GDEM, providing altitudes and other topographical in-formations (slopes, expositions).The third dataset involves average annual temperatures, interpolated from the Global Historical Climatological Network at high resolution (Hijmans et al., 2005).Corine Land Cover map is the fourth dataset.The Corine program assess ecological habitats and land uses at European scale, with remote sensing and photointerpretation observations.The downloaded datasets are available for 1990, 2000, 2006 and 2012.In our case, the 2006 edition is used in order to suit temporally with the locations of wild reindeer.The last dataset is represented by MODIS 13Q1 time series.Bi-monthly temporal acquisition has been chosen, from the beginning of 2000 to the end of 2015.Dates of imagery acquisition is expressed in Julian day.The acquired tiles extent to entire Scandinavia, with a spatial resolution of 250 meters.Since the native dataset is in*.HDR format, Normalized Difference Vegetation Index (NDVI) band has been chosen, and converted into GeoTiff format.Different softwares are used in the study, depending on processing steps.First of all, we used GIS software, such as ArcMap 10.2 and Qgis 2.18 for both preprocessing and processing.Then, R-studio and Ade-Habitat-HR (Home Ranges, Callenge, 2006) package allow us to calculate reindeer home-ranges, using Kernel Utilization Density technique, on a bi-monthly basis.Then home ranges have been exported in separate GIS files (*.SHP).These home ranges and associated quantification of reindeer's densities constitute the key-input for the generation of new datasets.Regular spreadsheets, such as Open Calc 4.1.3or Office Excel 2016 versions are used all along the study.Both R (package Rcmdr) and Xlstat can be used for complex statistical calculation.For this specific study, Xlstat has been privileged.The general approach requires to define the spatial and temporal unit to compare and to analyze, in this case bimonthly home ranges during the growing season (April to August).Grid extent will be generated at the regional scale and related cells used for the extraction will depend on bi-monthly home ranges distribution.The spatial and temporal resolution is forced here on the MODIS-13Q1 accuracy, but this can be specified in GRiD-toolbox parameters.

Monitoring wild reindeer (Rangifer tarandus L.) home ranges in Hardangervidda Plateau (Norway)
At the scale of the entire population equipped with GPScollars (seven individuals), the GPS-tracks has been sorted out by bimonthly period and exported in *.CSV format.Data have been sorted in two packages.The first sorting aggregates bi-monthly GPS-tracks during the whole available period (2007-2010) and represent the seasonal variability of the growing period (two bimonthly home ranges by month, from April to August, totaling ten files).The second sorting focuses on the inter-annual variability of the ten bi-monthly home ranges during three years (2007, 2008 and 2009), giving 30 files.
The GPS-tracks are then ready to estimate the reindeers' utilization distribution.Using the R-package ADE HabitatHR (Callenge, 2006), we compute the Kernel Utilization Distribution (KUD) for each fortnight.The bivariate normal kernel technique is specified, as well as the smoothing parameter calculated with « ad-hoc » method (Worton, 1989;Skarin et al., 2008).Grid size (~100 m²) and extent (h-value for estimated densities; 1km²) are also specified to fit with reindeer's large scale migration properties.The resulting convex hull has been vectorized, with 95% of the kernel home-range extent as boundary.Then datasets are exported, respectively into GeoTiff raster grid for Kernel Utilization Distribution densities, and into shapefiles for the 95% convex hull.
Finally, a grid is generated, intersected following reindeer KUD' shapes and geophysical descriptors systematically sampled.

GRiD -toolbox: ArcMap modelbuilding and Python programming (Qgis)
The figure 1 summarizes the processing chain used for GRiD (Grid Raster Information Dataset) tool-box creation.Two prototypes of GRiD are available, even if improvements are still on progress.Here the ArcMap Model-builder is preferred, but Qgis GRiD toolbox is ready to use for the third and last step shown in figure 1.
According to the goals of the study, the first three steps are flexible in the process chain.
A key step is to define a common projection in order to minimize geometrical error between datasets.The WGS84-UTM 32N projection has been chosen, initially used for locations of wild reindeers.After projecting every dataset, we generated a regional grid, based on MODIS 13Q1 tile extent and spatial resolution (250m).
Each cell is theoretically supposed to suit with MODIS 13Q1 coarse grid, and thus values of NDVI.The step 3 converts features (cells) to points located at the core of each cell.The resulting mesh of points still has a regional extent, such as the basic grid obtained from step 2. The generated bi-monthly home ranges of wild reindeers are put into the chain, in order to intersect the regional mesh of points.As output, each intersect processing produces one shapefile and one spatial unit for home range analysis.

2.2.3.1
Checking correspondences between raw and sampled datasets We quantify the possibility to lose spatial information, by upscaling geophysical datasets at MODIS -13Q1 spatial resolution (~250m).In the study, the comparison will be made between raw Corine Land Cover dataset (e.g.100m of spatial resolution) and new Corine dataset extracted using GRiD toolbox (250m).The analysis involves the calculation of particular and total error (%) between raw and new dataset.A threshold of 5% of total error has been retained.

2.2.3.2
Removing geophysical factors which are spatially autocorrelated The goal is here about deleting geophysical factors which could be redundant in the analysis between reindeer's home ranges and reindeer's kernel densities within bimonthly home ranges.At the scale of one of the most extent home range (2nd half of April), continuous values (annual temperatures, altitudes, exposition, slopes, NDVI) are crossed using Pearson's correlation.A threshold of 0.8 is retained for the significance of a relationship, allowing then to remove redundant factors.

2.2.4.1
Characterizing the seasonal variability of the chosen ecological habitats For the characterization of the inter-seasonal variability of the choice of ecological habitats, we summed every pixel belonging to each habitat class within the Kernel Utilization Distribution area.Adding up pixels allow us to compare habitat's areas (sum of pixel's absolute values, or percentage calculation) between seasonal home ranges, aggregated during the whole period.The statistical summaries and bar charts for every home range are then summed up by selecting the first and the second most frequent ecological habitat by bi-monthly reindeer's home range.A Khi²-test is then computed to verify the independence between the fortnight and the category of land cover (number of pixels).

2.2.4.2
Inter-annual variability of bioclimatic conditions: choice of particularly different years for spring migration, calving area and "fat accumulation" area When extracting MODIS-NDVI time series for annual home ranges, care was taken to sample the three Julian date acquisition relevant to bi-monthly home range.It gives three NDVI time series (columns) by bi-monthly home range.For example, when extracting geophysical raster based on 1st half of April, the 97th Julian day of 2007, 2008 and 2009 have been added to the 1st half of April home range.We then check the independency of these three inter-annual statistical distribution, using Kruskal-Wallis independency tests.Based upon the median of rank generated distributions, two different years in regards of NDVI inter-annual distributions have been selected.

Testing relationships between biotic and abiotic factors
The selected contextual bi-monthly home ranges correspond to three major moment in reindeer's annual biology (descent to summer pastures = 1st half of April, calving time = 2nd half of May, fat accumulation time = 2nd half of July).For those particular home ranges, we tested the relationships between contextual estimated densities and other extracted factors (altitudes, exposition, slopes or proxies of primary chlorophyll production such as NDVI).Spearman correlation test is preferred in this case.According to Spearman table of significance, with a confidence threshold of 95% and a number of observations >100, the given threshold for significance is above 0.197.

Differences in habitat choices
We compare the NDVI values distribution and the reindeer's densities probabilities distributions belonging to the same home range.Our aim is to verify the similarities or dissimilarities between these two kinds of distributions.Kruskal-Wallis independence tests has been conducted using habitat categories as sub-samples.Ecological habitat categories divide both densities estimation and NDVI values as many sub-samples as ecological classes the observations belong to.2) according to the spatial resolution, in one of the most extent home range (2nd half of April, e.g.17236 cells; ~4309 km²).Similar percentages can be found within coarse dataset (3rd line) compared to the percentages for GRiD-extracted dataset.

Results and suggestions
The last line subtracts the difference between GRiDextracted percentages to coarse Corine Land Cover.Each accounted error has been added up to the total, giving a total error about 0.59 %.Close to 99.4% of representativeness, the GRiD-toolbox extraction shows good results, in particular for Corine Land Cover ecological habitats.That said, some improvements could reinforce spatial accuracy of GRiD-toolbox.Firstly, the quantification of error between specific spatial resolution of coarse datasets and inputted spatial resolution using GRiDtoolbox has to be more thoroughly assessed.It could be interesting to repeat the same test for other datasets, such as raw ASTER-GDEM (spatial resolution e.g.30m) and extracted values (250m).At the scale of the analysis (here the 2nd most extent reindeer's home range), it is important to notice that the number of observations (e.g.pixels) can vary according to the considered spatial unit.
Since the number of geostatistical observations is very different between home ranges, checking spatial unit representativeness for each dataset would allow us to have a good overview of total error for a given spatial resolution as GRiD toolbox input parameter.
The second improvement is related to geometrical bias, probably generated by the chosen map projection (UTM-32N).Indeed, Mercator projection is more suited with equatorial region, and tends to elongate surfaces when dealing with polar or subpolar regions.Working with an official map projection clipped by zone (EUREF89-NTM in our case) will likely improve overall results for GRiDtoolbox and wild reindeer distribution analysis.Another parameter (choosing cartographic projection) could be added to the toolbox and automatically process the projection homogenization between the numerous datasets we are willing to extract.

Auto-correlated and redundant data
The auto-correlated factors are annual mean temperatures values  and altitudes.This can be explained by the altitudinal gradient, in particular between the high altitudes and low annual temperatures in the north of the studied area, low altitudes and higher temperatures in the south (Hardangervidda plateau).At the scale of the second most extent home range, the correlation is strong and significant (r = +0.83;p<0.0001).Then, it is statistically possible to express altitudes distribution as a good proxy for annual temperature distribution.That said, it could be interesting to reprocess the calculation with a more recent dataset (1980-2010 for example) and more accurate temperatures variables (minimal and maximal temperatures).

Ecological habitats and intra-annual (e.g. interseasonal) variability of home ranges
In absolute values of pixels, intermediary seasonal home ranges (early spring, late summer) have the greatest overall area, compared to calving period and full-season home ranges.For example, the most extent home range (1st half of April) counts 19736 geostatistical observations (~4480 km²) when the smallest one (2nd half of June) is about 2262 pixels (~565 km²).Indirectly, this fact is shown in the figure 2. Percentages of the two most frequent ecological habitats by fortnight and home ranges during the whole period are here displayed.Colors represent the type of ecological habitat, with a total number of 6.For both April and August fortnights, a relatively balanced distribution is noticed compared to June.Moors and heathland areas are widely represented, with no less than 35% as minima for every home range total areas.Water bodies are represented for the 1st July home range (~10%), since any internal differentiation has been made within bi-monthly home ranges.The Khi²-test (performed on the number of pixels of the 9 more represented habitats for the ten fortnights) shows a strong relationship between these variables (observed value: 9873.29,critical value: 92.80, p<0.0001).Some habitats are over-represented during early spring, particularly sparsely vegetated area.During spring and summer, moors and heathland and mixed forests are overrepresented.During late summer, coniferous forests and transitional woodlands are over-represented.Thus, due to the wider extent of reindeer's home ranges during intermediary seasons, habitat diversity appears to be higher when females start to go towards the calving areas.The same can be said during the summer, with mixed forests and later coniferous habitats emerging after the phenological peak of broad-leaved forests followed by the one of moors.Explanations can be found when looking at opportunistic reindeer feeding strategies, in particular in spring when annual herbaceous just start growing (Skarin et al., 2008).No clear correlation for "fat accumulation" home ranges (2nd fortnight of July) has been found.Here, the relative abundance of bio-mass, shortly after the phenological peak, can explain the lack of relationship.

Conclusions
Mounting GRiD toolbox has been incrementally encouraged by the needs of automation processing for the survey of wild reindeer migration and densities in Hardangervidda Plateau (Norway).The wide range spatial behavior of reindeer is an interesting case study to apply this method, more adapted to large territories than simple grid multi-data analysis previously realized by other authors in smaller areas (Jolivet 2014, Bortolamiol et al., 2016).Indeed, by choosing the core point of each grid cell as the common reference to cross all the datasets, the process is realized within a reasonable computation time while avoiding important bias according to the preliminary test realized in our study.Improvements still have to be made to ensure GRiD geometrical reliability, cartographic projection homogenization and minimizing error when extracting large and various raster datasets.Applying GRiD toolbox for wild reindeer survey during the three growth periods simplified pre-treatments by automatizing grid and creation of mesh of points; and finally datasets fusion into one attribute table.Based on geo-statistical observations (pixels) it allowed us to compute several complex statistical analyses, not usually present in GIS software.Based on Open Access datasets, the study firstly shows interesting links between reindeer intra-annual distribution and ecological habitat diversity; secondly relationships between inter-annual variability of estimated densities and abiotic factors, validating thus our second and third hypothesis.

Fig. 1 .
Fig. 1.General process line for GRiD toolbox, using ArcGIS Model-builder (credits: R. Courault, 2017) The final step is crucial: raster values are sampled and resulting attribute tables exported as a *.CSV file.Each generated CSV matrix blend spatial unit of analysis (here bi-monthly home ranges of reindeers) with different geophysical datasets we are willing to statistically sum up or establish relationships (columns).Each observation (line) is a pixel contained into a specific pre-calculated home range.These geostatistical observations are now characterized with one specific value both for biotic and abiotic factors, and are ready to analyze.

3. 1
GRiD toolbox: sufficient accuracy % of total area for each present class of habitat for raw Corine Land Cover area for each present class of habitat, extracted Corine Land Cover dataset,

Fig. 2 .
Fig. 2. Comparison between reindeer's fortnight home ranges along vegetational season (April to August) and the first two major frequent CLC habitat classes by home range (sources: Corine Land Cover ecological habitat, European Environmental Agency; AdeHabitat package, Callenge, 2006; GRiD-toolbox; credits: R. Courault, 2017) 3.2.3NDVI dissimilarities and selecting specific home ranges The Kruskal-Wallis test is used to verify if statistical distributions are significantly different.This is the case for the three annual home ranges tested (spring migration: 1st half of April; calving area: 2nd half of May; fat accumulation area: 2nd half of July), totaling 9 distributions.The NDVI values (2007-2008 and 2009) are statistically different (p<0.0001).Between the three years, it is now possible to determine what are the « cold

3. 3 . 2
Strong inter-annual variability in estimated densities and phenology between the inter-annual home rangesIn order to explore more thoroughly the relationship between reindeer estimated densities and NDVI values, Kruskal-Wallis tests have been computed.Statistical distributions have been clustered by sub-samples, e.g.ecological habitat classes.Every test calculated for the 9 specific home ranges has rejected the dependence hypothesis be-tween the geostatistical distributions of estimated reindeer densities and NDVI values (9 tests for the 9 specific home ranges).Such variability is displayed in figure3.The set of maps shows likeness in pattern distribution of estimated densities shapes (ethological/individual behavior, calving area), but also dissimilarities (environmental factors depending on contextual year, spring migration).

Table 2
compares the percentages of Corine Land Cover ecological habitats (1st line, table