Phenology-Based Biomass Estimation to Support Rangeland Management in Semi-Arid Environments

Livestock plays an important economic role in Niger, especially in the semi-arid regions, while being highly vulnerable as a result of the large inter-annual variability of precipitation and, hence, rangeland production. This study aims to support effective rangeland management by developing an approach for mapping rangeland biomass production. The observed spatiotemporal variability of biomass production is utilised to build a model based on ground and remote sensing data for the period 2001 to 2015. Once established, the model can also be used to estimate herbaceous biomass for the current year at the end of the season without the need for new ground data. The phenology-based seasonal cumulative Normalised Difference Vegetation Index (cNDVI), computed from 10-day image composites of the Moderate-resolution Imaging Spectroradiometer (MODIS) NDVI data, was used as proxy for biomass production. A linear regression model was fitted with multi-annual field measurements of herbaceous biomass at the end of the growing season. In addition to a general model utilising all available sites for calibration, different aggregation schemes (i.e., grouping of sites into calibration units) of the study area with a varying number of calibration units and different biophysical meaning were tested. The sampling sites belonging to a specific calibration unit of a selected scheme were aggregated to compute the regression. The different aggregation schemes were evaluated with respect to their predictive power. The results gathered at the different aggregation levels were subjected to cross-validation (cv), applying a jackknife technique (leaving out one year at a time). In general, the model performance increased with increasing model parameterization, indicating the importance of additional unobserved and spatially heterogeneous agro-ecological effects (which might relate to grazing, species composition, optical soil properties, etc.) in modifying the relationship between cNDVI and herbaceous biomass at the end of the season. The biophysical aggregation scheme, the calibration units for which were derived from an unsupervised ISODATA classification utilising 10-day NDVI images taken between January 2001 and December 2015, showed the best performance in respect to the predictive power (Rcv = 0.47) and the cross-validated root-mean-square error (398 kg·ha−1) values, although it was not the model with the highest number of calibration units. The proposed approach can be applied for the timely production of maps of estimated biomass at the end of the growing season before field measurements are made available. These maps can be used for the improved management of rangeland resources, for decisions on fire prevention and aid allocation, and for the planning of more in-depth field missions.


Introduction
The livestock sector is economically important in Niger (for a map see Figure 1), contributing on average 10% to the gross domestic product (GDP) of Niger during the period 2009-2013 [1]. The agriculture and livestock census of 2005/2007 [2] estimated the total number of livestock to be around 31 million, composed mainly of cattle, sheep, and goats. The geographical distribution of livestock is not homogeneous in the country, and the largest numbers of livestock are in the regions of Zinder, Tahoua, Maradi, and Tillabery. Three different livestock systems exist in Niger, which are adapted to the agro-ecological conditions in the different zones of the country. A sedentary livestock system (accounting for 66% of the livestock) is practised together with the cultivation of crops in farms in the agricultural zone of the south (see Figure S1). Low-distance transhumance (i.e., seasonal movement of livestock) to pastoral enclaves during the rainy season is performed to avoid crop damage. Nomadic herding (18% of the livestock) is the main or only activity of herders in the pastoral zone or in the transition between the agricultural zone and the pastoral zone to its north. Livestock movements depend on the availability of water and pasture. The transhumance system (involving 16% of total livestock) is practised by herdsmen belonging to the Fulani ethnic group and is characterised by the seasonal movement of large herds between two distinct areas through well-defined corridors. Generally, the movement from north to south occurs in late winter, and transboundary transhumance to Benin, Nigeria, or Burkina Faso is common [2].
The rangelands of Niger are mainly located in the Sahel zone, a semi-arid region between the Sahara desert in the north and the Sudanian savannah in the south and are highly vulnerable as a result of the high inter-annual precipitation variability [3,4] and, hence, rangeland production variability. In addition to the inter-annual weather variability, the Sahel region of Niger shows great spatial precipitation heterogeneity, making rangeland production highly variable in space in a single year. Taking this variability into account would allow better management of rangeland resources and animal movement. The availability of a biomass production map at the end of the growing season represents a key tool to enhance sustainable and efficient rangeland management and improve food security in this region as a large number of households depend on livestock to sustain their livelihoods. Together with the number of livestock estimated by administrative units, this map could be used to compile a forage balance, identifying areas potentially affected by forage surplus or deficit that create a risk for fires or livestock mortality, respectively.
There is a long tradition of remote sensing (RS)-based herbaceous biomass estimation in Niger. The first studies in the late 1980s [5][6][7][8] applied linear regression between maximum standing biomass and normalised difference vegetation index (NDVI)-based indices (time-integrated and maximum NDVI) derived from the National Oceanic and Atmospheric Administration (NOAA) Advanced Very High Resolution Radiometer (AVHRR) imagery. The Ministry of Livestock in Niger currently uses a similar method based on NDVI data from the Satellite Pour l'Observation de la Terre-Vegetation (SPOT-VGT) and, since 2014, data from the Meteorological Operational satellite programme (MetOp) AVHRR. The ministry fits a linear regression model between ground measurements of the current season and a NDVI-based metric (maximum NDVI, mean NDVI, or cumulative NDVI are tested, and the best-performing metric is selected) to create a map of rangeland biomass. In 2014, Nutini et al. [9] estimated the extent of rangeland biomass at three sites in Niger with a radiation use efficiency model utilising cumulative dry matter productivity derived from SPOT-VGT [10] over a fixed time period, corrected with the evaporative fraction and derived from Moderate-resolution Imaging Spectroradiometer (MODIS) albedo and thermal measurements.
In addition to Niger, aboveground biomass has been estimated in several other semi-arid regions using RS data, including optical, radar, and combined multi-sensor imagery and modelling approaches (see [11] for a review). The majority of studies reviewed by Eisfelder et al. [11] applied medium-and low-resolution optical or radar data to derive a RS indicator that was used in an empirical relationship with field biomass measurements. The review by Ali et al. [12] on satellite RS of grasslands in general (not restricted to semi-arid regions) includes methodologies to retrieve grassland biophysical Remote Sens. 2017, 9,463 3 of 16 parameters. The authors concluded that regression models based on vegetation indices were the predominant approach for biomass estimation, while machine learning algorithms and the fusion of multisource data in biophysical simulation models still require further research to better utilise their potential [12].
Recent studies in the Sahel zone have used linear regression to estimate aboveground biomass from satellite-derived phenological metrics such as maximum value, the start-of-season and end-of-season values, and the cumulative value of the RS index during the growing season. These studies used time series of different RS indices such as the NDVI [13], the fraction of absorbed photosynthetically active radiation (FAPAR) [14,15], and the vegetation optical depth (VOD) [16]. In addition to simple linear regression, recent studies have also applied multi-linear models of phenological metrics [14] and machine-learning models combining FAPAR seasonal metrics and agro-meteorological data [17] to estimate annual (herbaceous) biomass production in Sahelian ecosystems of Senegal.
High spatial heterogeneity due to topography, cover, or relatively small area like in the Mediterranean rangelands challenge RS-based vegetation mapping and monitoring [18]. This heterogeneity can be addressed with higher spatial resolution imagery (for a review see [11]) or spectral unmixing modelling as a substitute for higher spatial resolution [18]. Landsat and very high-resolution data (e.g., IKONOS, QuickBird, HyMap) have been used for biomass estimation in Botswana, the Sahel and Sudan, West Africa, South Africa, and Zimbabwe [11]. However, biomass estimation has been mostly limited to a single year or very few years and relatively small target areas. An example of the application of an unmixing method for the retrieval of herbaceous biomass is the study by Svoray and Shoshany [19] that merged synthetic aperture radar (SAR) images with unmixed Landsat TM images.
This study aims to support existing rangeland management activities in Niger by developing an RS-based approach for mapping rangeland biomass. In contrast to the current method used by the Niger Ministry of Livestock, which is calibrated with the field observation of the current year, we intended to develop a method that does not require field measurements of the current season to estimate the herbaceous biomass at the end of the season. The proposed approach uses the phenological timings of the start and end of season extracted from the RS time series to retrieve a proxy of rangeland biomass at the end of the growing season; namely, the cumulative value of NDVI during the growing season (cNDVI). A linear regression model is calibrated with multi-annual field measurements of herbaceous biomass. The calibration of the regression model using multi-annual and multi-site observations can be performed with different aggregation levels, ranging from the site level (i.e., site by site) to the global level (i.e., all sites pooled together). The site-level calibration can be used to gain insights into the robustness of the relationship by analysing the slope and intercept across sites [9]. However, some level of aggregation is needed to map the biomass. Pooling all observations to derive a single regression model may represent an optimal use of the data if the relationship between measured biomass and the RS variable is constant in space. However, location-specific factors such as herbaceous species type and background reflectance may generate differences in the relationship between the RS variable and biomass. This spatial heterogeneity in the relationship can be modelled by grouping the field observations into more homogeneous sets. Different spatial aggregation schemes are proposed and tested to select the design providing the highest prediction power in a similar way to that suggested by Meroni et al. [20]. The method that has been developed can also be used in a predictive fashion to estimate the biomass at the end of the current season before field survey results are made available. The model can also serve as a backup solution in the event that field surveys are not carried out in a specific year or a specific region.

Study Area
The West African state of Niger features a pronounced climatic gradient. According to the Köppen-Geiger climate classification [21], the northern part is characterised by a hot arid desert and the southern part by a hot arid steppe climate. Mean annual precipitation ( Figure 1) ranges from <100 mm in the north to 600-700 mm in the south of the country (calculated from TAMSAT rainfall estimates from the TARCAT v2.0 dataset for the period 1986-2015 [22,23]). Most of the precipitation falls during the rainy season from June to September and is associated with the West African monsoon season [3]. The mean annual surface temperature ranges from about 23 • C in the northernmost part of the country to 27-29 • C in the central and southern parts, with the highest mean monthly values occurring during the summer (e.g., range across the country in June: 31-35 • C) and the lowest during the winter (e.g., January: 13-25 • C) during the period 1981-2010 [24,25].
Remote Sens. 2017, 9,463 4 of 16 <100 mm in the north to 600-700 mm in the south of the country (calculated from TAMSAT rainfall estimates from the TARCAT v2.0 dataset for the period 1986-2015 [22,23]). Most of the precipitation falls during the rainy season from June to September and is associated with the West African monsoon season [3]. The mean annual surface temperature ranges from about 23 °C in the northernmost part of the country to 27-29 °C in the central and southern parts, with the highest mean monthly values occurring during the summer (e.g., range across the country in June: 31-35 °C) and the lowest during the winter (e.g., January: 13-25 °C) during the period 1981-2010 [24,25]. . Precipitation totals were computed from 10day TARCAT data [22,23]. Dark blue dotted lines refer to 100 mm isohyets.
Following the precipitation gradient, the vegetation also changes from north to south. The main ecoregions in Niger include the Sahara desert, the south Saharan steppe and woodlands, the Sahelian acacia savannah, and the west Sudanian savannah [26]. The focus of this study is the pastoral zone of Niger, which is mainly located in the area with a mean annual precipitation of 100-300 mm (see Figure 1) and generally belongs to the Sahelian acacia savannah. As a precise delineation of the pastoral area is not available, the following limits were adopted: the 380 mm mean annual precipitation isohyet marks the southern limit, while the extent of vegetated area (as defined in Section 2.3.1) marks the northern limit (the final outline is shown in Section 3). In this way, we excluded the more agricultural landscape in the south while including all field measurement sites.
The mean annual precipitation sum of the sampling sites ranges from about 107 mm to 372 mm (calculated from TARCAT data [22,23]), while the mean annual temperature ranges from about 25.6 to 29.3 °C (calculated from ECMWF ERA INTERIM model data [27]). The dominant species in the herbaceous layer at the sample sites are Cenchrus biflorus, Aristida spp., Schoenefeldia gracilis, Dactyloctenium aegyptium, and Alysicarpus ovalifolius. In general, the species composition does not vary significantly annually at any particular site unless there is a bush fire. In contrast, the species composition varies between sites and is mainly linked to the soil type (A. Boureima, unpublished data).   , location of sample sites used for the regression model, and administrative regions of Niger (dark grey). Precipitation totals were computed from 10-day TARCAT data [22,23]. Dark blue dotted lines refer to 100 mm isohyets.

Data
Following the precipitation gradient, the vegetation also changes from north to south. The main ecoregions in Niger include the Sahara desert, the south Saharan steppe and woodlands, the Sahelian acacia savannah, and the west Sudanian savannah [26]. The focus of this study is the pastoral zone of Niger, which is mainly located in the area with a mean annual precipitation of 100-300 mm (see Figure 1) and generally belongs to the Sahelian acacia savannah. As a precise delineation of the pastoral area is not available, the following limits were adopted: the 380 mm mean annual precipitation isohyet marks the southern limit, while the extent of vegetated area (as defined in Section 2.3.1) marks the northern limit (the final outline is shown in Section 3). In this way, we excluded the more agricultural landscape in the south while including all field measurement sites.
The mean annual precipitation sum of the sampling sites ranges from about 107 mm to 372 mm (calculated from TARCAT data [22,23]), while the mean annual temperature ranges from about 25.6 to 29.3 • C (calculated from ECMWF ERA INTERIM model data [27]). The dominant species in the herbaceous layer at the sample sites are Cenchrus biflorus, Aristida spp., Schoenefeldia gracilis, Dactyloctenium aegyptium, and Alysicarpus ovalifolius. In general, the species composition does not vary significantly annually at any particular site unless there is a bush fire. In contrast, the species composition varies between sites and is mainly linked to the soil type (A. Boureima, unpublished data).

Biomass Data
In this study we used a dataset of measured aboveground herbaceous biomass (B) at the end of the growing season for the period 2001-2015, provided by the Ministry of Livestock of Niger. The dataset included 90 sites within the pastoral zone, with varying levels of information available for the study period, and comprised 926 records of dry matter production (kg·ha −1 ) in total. The fresh and dry leaf matter are considered and are here assumed to represent the total aboveground biomass of the concluded season. The date of sampling varied but was usually between mid-September and the end of October. Each sampling site represents a (relatively) homogeneous area of 3 km × 3 km. The sampling of herbaceous biomass was performed along a transect, applying a double-sampling approach combining weighed biomass and visual estimates that finally yields estimates of dry matter production. The sampling method is described in more detail in Maidagi et al. [5] and Wylie et al. [8].
Biomass measurements should be made in the initial stages of the senescence phase, after the time of maximum vegetation development, typically between mid-August and mid-September. In this way, the measured biomass should not be affected by the natural decay of leaf matter and possible grazing that may take place during senescence. In addition to the tendency to underestimate the total biomass, late biomass measurements are less correlated to the RS-derived biomass proxies, which are dominated by the signal recorded during the period of maximum vegetation cover. To exclude measurements that are poorly representative of the seasonal biomass production we thus discarded biomass measurements from sampling carried out after 25 October (date proposed by local experts).
The quality of data from the database is known to be variable and field operator-dependent; therefore, we analysed the dataset and screened it for outliers. We inspected the site-level relationship between biomass and the RS indicator (see Section 2.3.1). As a relatively good site-specific relation (i.e., positive correlation) between the biomass and the RS indicator is expected, we opted to exclude from further analysis those sites with a coefficient of determination (R 2 ) of less than 0.2 for the linear regression between the two time series (number of sites excluded = 13). Sampling sites with fewer than four observations in the study period were discarded from the dataset (number of sites excluded = 21). Finally, 56 sites with a total of 616 records were retained for model parameterisation.
The biomass data at the end of the growing season of the rangelands in Niger showed high variability in space and time during the period 2001-2015. The multi-annual average biomass (dry matter production) of single sites ranged from 331 kg·ha −1 to 1778 kg·ha −1 , with a mean value of 760 kg·ha −1 . The annual spatial average of all sites varied over the study period from 246 kg·ha −1 (2004) to 1228 kg·ha −1 (2007), with an overall mean value of 732 kg·ha −1 . This high spatial and temporal variability of end-of-season biomass is mainly due to the high precipitation variability e.g., [28,29] and is in accordance with observations in other parts of the Sahel region e.g., [14,30,31].

Remote Sensing Data
This study used the eMODIS NDVI product, based on MODIS data acquired by the Terra satellite and provided by the United States Geological Survey (USGS; data portal: http://earlywarning.usgs. gov/fews/). The eMODIS product is a 10-day maximum value NDVI composite [32], temporally smoothed with the Smets algorithm [10]. Composite images are produced every five days; hence six temporally overlapping composites are generated per month. Here we only used the composite images for days 1-10, days 11-20, and days 21 to the last day of each month for the period January 2001 to December 2015. More information on the eMODIS NDVI product can be retrieved from the product documentation [33].

RS Proxy for Biomass Production
The seasonal cumulative NDVI (cNDVI) is defined as the integral of NDVI during the growing season subtracted from the area under the baseline (i.e., the NDVI minimum level before the start of the growing season), as proposed by Meroni et al. [34] for the seasonal cumulative FAPAR. The time interval for integration is defined by the start of season (SOS) and the end of season (EOS), estimated for each pixel and season. The model-fit approach of Meroni et al. [35], with further modification as in Vrieling et al. [36], was used to calculate the phenology parameters (SOS, EOS, maximum value of NDVI). SOS was estimated as the point at which the fitted NDVI model for the season exceeded 20% of the local growing amplitude (i.e., between minimum NDVI before green-up and maximum NDVI of that season), and EOS was estimated as the point at which the model falls below 80% of the decay amplitude (i.e., between maximum NDVI of the season and the following minimum NDVI after decay). Pixels are considered non-vegetated when the variability of the NDVI time series, as measured by the difference between the 95th and 5th percentiles, is below 0.05 NDVI units. Other methodological details can be found in Vrieling et al. [36]. The value of cNDVI is thus controlled by the integration limits, the baseline value, and the amplitude and shape of the NDVI seasonal trajectory.

Linear Regression Model and Spatial Aggregation Levels
Linear regression models are built by matching the field measurements at a given site and time with the corresponding cNDVI. The cNDVI values for the site were calculated as the mean of a kernel of 11 × 11 pixels centred on the site coordinates to roughly match the 3 km × 3 km area for which the field measurements are representative. A smaller kernel size (2750 m) was set to avoid border effects.
Relationships between cumulative vegetation indices and dry biomass are empirical and have only local values [37], meaning that a regression established in a given agro-ecological context may not be successfully extrapolated to a different region. One approach to tackle this problem would be to define smaller geographical units with similar physical characteristics and to treat them separately. However, improved performance is not guaranteed as the total number of observations available is constant, and the sample size (over which the single regression is calibrated) is reduced by subsampling. Thus, a trade-off between model specificity and data availability exists. We explored this trade-off empirically to find the modelling solution providing the best performance in prediction.
Different levels of spatial data aggregation (and thus model parameterisation) were considered for model calibration. In the simplest considered model, the whole set of measurements was pooled, and a single linear regression model was established (hereafter referred to as the 'global' model). Such a model assumes a constant relationship between biomass and cNDVI in both space and time. However, this relation might not hold in all circumstances, in particular with respect to spatial variation. The relation between the two variables may, in fact, vary spatially as a result of different species composition or background reflectance affecting NDVI estimates. Such spatial heterogeneities may be accounted for by limiting the spatial domain at which the model is calibrated. The stratification of the study area into smaller and more homogeneous units is thus expected to reduce the estimation error in the temporal domain. We therefore attempted a number of stratification strategies (Figure 2) with the aim of selecting the one providing the best performance in prediction. This selection was done empirically on the basis of a cross-validation (cv), applying a jackknife technique (leaving out one year at a time) to assess the predictability of biomass based on cNDVI. For each year left out, the linear regression coefficients were estimated based on the remaining dataset and used to predict the biomass values that were not used in the calibration. Then, all single jackknifed predictions were compared with the actual measurements to compute the cross-validated R 2 (R 2 cv ). Based on the predicted biomass values from the jackknifing, a cross-validated root-mean-square error (RMSE cv ) was calculated as follows: whereB i is the predicted biomass derived from the jackknifing, B i is the measured biomass, and n is the number of samples. The RMSE cv gives an indication of the magnitude of the errors in prediction.
Remote Sens. 2017, 9, 463 7 of 16 more biologically meaningful models. For a better parameterisation of the relationship between biomass and cNDVI, different stratifications with a more biophysical rationale were tested. First, we used the Global Agro-Environmental Stratification (GAES [39]). The GAES is derived from the segmentation of 13 input layers with 1 km spatial resolution, aiming to stratify agricultural production zones according to the region's agro-environmental characteristics, including climate, altitude, irrigation, production, phenology, growing cycles, crop type, and field size parameters. The GAES is produced at four hierarchical levels (level 4 being the most detailed) using eCognition software on a continental basis for the segmentation. For our study, we used GAES level 4 to stratify the sample sites (referred to as the "GAES" model), which are represented by six GAES strata. The GAES stratum represented by just one site (ID 177) was assigned to another GAES stratum based on spatial proximity and agro-environmental type description, resulting in a total of five GAES units of the model. Second, soil types from the African soil map [40] were used for stratification (referred to as the 'soil' model). The basic information for this harmonised continental scale soil map was derived from the Harmonized World Soil Database [41], and the naming of the soil types followed the World Reference Base for Soil (WRB) classification and correlation system [42]. The sample sites were characterised by nine different soil types. Soil types represented by just one sample site were assigned to other soil types based on spatial proximity and similar characteristics, i.e., eutric fluvisol was merged with haplic vertisols (both soil types are influenced by water) and vertic cambisols were merged with eutric cambisols (both are cambisols), resulting in a total of seven soil units.
Third, we created a stratification by intersecting GAES and soil type units (referred to as the "GAES + soil" model) to further refine the segmentation of the study area. This intersection resulted in a total of 13 combinations. Combinations represented by just one sample site were merged with another GAES + soil combination with the same GAES unit, i.e., "GAES unit 1 + soil unit 2" was assigned to "GAES unit 1 + soil unit 1" and "GAES unit 5 + soil unit 2" was merged with "GAES unit 5 + soil unit 3", resulting in a total of 11 GAES + soil units.
Fourth, we generated a stratification layer based on the unsupervised classification of NDVI images from eMODIS (see Section 2.2.2), applying an unsupervised clustering algorithm, namely the Iterative Self-Organizing Data Analysis Technique (ISODATA), similar to the approach of de Bie et al. [43]. The aim was to obtain a biophysical stratification purely relying on RS data. All 10-day NDVI images from January 2001 to December 2015 (n = 540; spatial subset defined by the mask of the study area) were used as input for the ISODATA clustering. We ran the ISODATA clustering with a fixed number of clusters, testing all options from five to 15. We then calibrated the regression model for each of the obtained 11 stratifications and selected the one with the highest R 2 cv as the 'biophysical model' An initial attempt to include local effects in the regression was the use of spatial proximity for stratification of available field measurements. For that we used administrative units (departments) (referred to as the 'department' model) taken from the Global Administrative Unit Layers (GAUL) level 2 [38], as administrative units are often used as stratification layers in official statistics. Note that the department with just one sample site (Goure) was assigned to the closest department, leading to a final total of 10 units. However, the departments are administrative units with arbitrary boundaries that do not guarantee any internal affinity in a biophysical sense. That said, the department model can still serve as a benchmark, a kind of 'arbitrary reference', to evaluate its performance relative to more biologically meaningful models. For a better parameterisation of the relationship between biomass and cNDVI, different stratifications with a more biophysical rationale were tested.
First, we used the Global Agro-Environmental Stratification (GAES [39]). The GAES is derived from the segmentation of 13 input layers with 1 km spatial resolution, aiming to stratify agricultural production zones according to the region's agro-environmental characteristics, including climate, altitude, irrigation, production, phenology, growing cycles, crop type, and field size parameters. The GAES is produced at four hierarchical levels (level 4 being the most detailed) using eCognition software on a continental basis for the segmentation. For our study, we used GAES level 4 to stratify the sample sites (referred to as the "GAES" model), which are represented by six GAES strata. The GAES stratum represented by just one site (ID 177) was assigned to another GAES stratum based on spatial proximity and agro-environmental type description, resulting in a total of five GAES units of the model.
Second, soil types from the African soil map [40] were used for stratification (referred to as the 'soil' model). The basic information for this harmonised continental scale soil map was derived from the Harmonized World Soil Database [41], and the naming of the soil types followed the World Reference Base for Soil (WRB) classification and correlation system [42]. The sample sites were characterised by nine different soil types. Soil types represented by just one sample site were assigned to other soil types based on spatial proximity and similar characteristics, i.e., eutric fluvisol was merged with Remote Sens. 2017, 9, 463 8 of 16 haplic vertisols (both soil types are influenced by water) and vertic cambisols were merged with eutric cambisols (both are cambisols), resulting in a total of seven soil units.
Third, we created a stratification by intersecting GAES and soil type units (referred to as the "GAES + soil" model) to further refine the segmentation of the study area. This intersection resulted in a total of 13 combinations. Combinations represented by just one sample site were merged with another GAES + soil combination with the same GAES unit, i.e., "GAES unit 1 + soil unit 2" was assigned to "GAES unit 1 + soil unit 1" and "GAES unit 5 + soil unit 2" was merged with "GAES unit 5 + soil unit 3", resulting in a total of 11 GAES + soil units.
Fourth, we generated a stratification layer based on the unsupervised classification of NDVI images from eMODIS (see Section 2.2.2), applying an unsupervised clustering algorithm, namely the Iterative Self-Organizing Data Analysis Technique (ISODATA), similar to the approach of de Bie et al. [43]. The aim was to obtain a biophysical stratification purely relying on RS data. All 10-day NDVI images from January 2001 to December 2015 (n = 540; spatial subset defined by the mask of the study area) were used as input for the ISODATA clustering. We ran the ISODATA clustering with a fixed number of clusters, testing all options from five to 15. We then calibrated the regression model for each of the obtained 11 stratifications and selected the one with the highest R 2 cv as the 'biophysical model' for comparison with the other aggregation schemes. Table S1 and Figure S2 show the model results for the different ISODATA clustering-based stratifications. Note that not all the clusters obtained were represented by sample sites. This does not affect the calibration of the model (there will just be no regression for the uncovered clusters), but it does need to be addressed in the creation of the biomass maps (see Section 2.3.3).
The biophysical model that performed best in prediction was obtained with 11 clusters (clusters 1 and 11 were not represented by samples; therefore, nine calibration units were used) and is shown in Figure 3. In the northern part of the study area, south of the Sahara desert, there occur some well-defined bands of clusters that are linked to the precipitation gradient (see isohyets in Figure 3). In the south, the spatial pattern of clusters is more complex, especially in the southwest of Niger. Here, a small-scale heterogeneous pattern occurs, attributed to higher variability of morphological features (e.g., rivers, depressions) leading to different plant growth conditions and different land uses (e.g., presence of agriculture). Figure S3 shows the ISODATA clusters compared with the elevation (derived from Shuttle Radar Topography Mission data [44]) of Niger, highlighting that some spatial patterns of the clusters are clearly linked to elevation differences. for comparison with the other aggregation schemes. Table S1 and Figure S2 show the model results for the different ISODATA clustering-based stratifications. Note that not all the clusters obtained were represented by sample sites. This does not affect the calibration of the model (there will just be no regression for the uncovered clusters), but it does need to be addressed in the creation of the biomass maps (see Section 2.3.3). The biophysical model that performed best in prediction was obtained with 11 clusters (clusters 1 and 11 were not represented by samples; therefore, nine calibration units were used) and is shown in Figure 3. In the northern part of the study area, south of the Sahara desert, there occur some welldefined bands of clusters that are linked to the precipitation gradient (see isohyets in Figure 3). In the south, the spatial pattern of clusters is more complex, especially in the southwest of Niger. Here, a small-scale heterogeneous pattern occurs, attributed to higher variability of morphological features (e.g., rivers, depressions) leading to different plant growth conditions and different land uses (e.g., presence of agriculture). Figure S3 shows the ISODATA clusters compared with the elevation (derived from Shuttle Radar Topography Mission data [44]) of Niger, highlighting that some spatial patterns of the clusters are clearly linked to elevation differences. Figure 3. Results of the best-performing ISODATA clustering (11 clusters, of which nine were represented by sample sites). Dark grey lines represent mean annual isohyets (in mm) computed from 10-daily TARCAT data from 1986-2015 [22,23]. Areas outside the study area (as defined in 2.1) are masked out in grey.
In summary, a set of different stratification schemes with varying spatial detail was tested for the model calibration, ranging from the more detailed models at the GAES + soil (number of calibration units = 11), department (n = 10), and biophysical (n = 9) levels to those at the soil (n = 7), GAES (n = 5), and, finally, global levels (n = 1). Given the linear model used, the total number of parameters to be estimated (with the same amount of data) is 2n (n being the number of stratification units). The final selection of the most suitable model is made in cross-validation to address the trade- Figure 3. Results of the best-performing ISODATA clustering (11 clusters, of which nine were represented by sample sites). Dark grey lines represent mean annual isohyets (in mm) computed from 10-daily TARCAT data from 1986-2015 [22,23]. Areas outside the study area (as defined in 2.1) are masked out in grey.
In summary, a set of different stratification schemes with varying spatial detail was tested for the model calibration, ranging from the more detailed models at the GAES + soil (number of calibration units = 11), department (n = 10), and biophysical (n = 9) levels to those at the soil (n = 7), GAES (n = 5), and, finally, global levels (n = 1). Given the linear model used, the total number of parameters to be estimated (with the same amount of data) is 2n (n being the number of stratification units). The final selection of the most suitable model is made in cross-validation to address the trade-off between the increased accuracy resulting from the increased model parameterisation and the decreased robustness and predictive power related to a reduced sample size on which each model is calibrated through an increased number of calibration units.

Map of Estimated Biomass
The estimated biomass (B e ) at the end of the season is then mapped using the developed regression model: where a and b are the model coefficients and cNDVI is the cumulative NDVI over the growing season of the year of interest. The maximum value between zero and the linear regression estimate is taken to discard negative biomass estimations that may be originated by models with negative intercept when the cNDVI is close to or equal to zero. The same coefficients are applied for all pixels when using the global regression model. In case of a stratified model, unit-specific coefficients are applied for each pixel within a certain unit.
In the case of the biophysical model, two clusters were not represented by sampling sites (i.e., cluster 1 and cluster 11). To obtain an estimation of biomass for those clusters, the coefficients of the next closest ISODATA clusters were used (cluster 2 parameters for cluster 1 and cluster 10 parameters for cluster 11).

Results and Discussion
Spatial patterns of the mean and coefficient of variation of the cNDVI over the period 2001-2015, show a rough north-south gradient (Figure 4), following the annual precipitation gradient (see Figure 1). The north of the study area (at the border to the Sahara) is characterised by a low mean cNDVI and high coefficient of variation values, the latter indicating high inter-annual variability. The south shows a higher mean cNDVI and lower inter-annual variability.
Remote Sens. 2017, 9,463 9 of 16 when the cNDVI is close to or equal to zero. The same coefficients are applied for all pixels when using the global regression model. In case of a stratified model, unit-specific coefficients are applied for each pixel within a certain unit.
In the case of the biophysical model, two clusters were not represented by sampling sites (i.e., cluster 1 and cluster 11). To obtain an estimation of biomass for those clusters, the coefficients of the next closest ISODATA clusters were used (cluster 2 parameters for cluster 1 and cluster 10 parameters for cluster 11).

Results and Discussion
Spatial patterns of the mean and coefficient of variation of the cNDVI over the period 2001-2015, show a rough north-south gradient (Figure 4), following the annual precipitation gradient (see Figure 1). The north of the study area (at the border to the Sahara) is characterised by a low mean cNDVI and high coefficient of variation values, the latter indicating high inter-annual variability. The south shows a higher mean cNDVI and lower inter-annual variability.

Regression Models
The R 2 value between measured biomass and cNDVI at each site can serve as a primary indication of the strength of the relationship between these two parameters. R 2 varies notably among the sites (ranging from 0.25 to 0.99), as indicated in Figure 4, without showing any clear spatial pattern. This variability might be due to data quality issues and/or other processes affecting the measured biomass (e.g., grazing) that are not considered in the simple linear regression with cNDVI. Mbow et al. [45] showed in their study in the Sahel of Senegal that species composition significantly

Regression Models
The R 2 value between measured biomass and cNDVI at each site can serve as a primary indication of the strength of the relationship between these two parameters. R 2 varies notably among the sites (ranging from 0.25 to 0.99), as indicated in Figure 4, without showing any clear spatial pattern. This variability might be due to data quality issues and/or other processes affecting the measured biomass (e.g., grazing) that are not considered in the simple linear regression with cNDVI. Mbow et al. [45] showed in their study in the Sahel of Senegal that species composition significantly affects the relationship between NDVI and biomass, concluding that temporal and spatial variation in species dominance could add noise to this relationship. In their study, the effect was less pronounced when applying the seasonal integral of NDVI instead of the peak NDVI value. Olsen et al. [13] demonstrated on the basis of a long-term grazing trial in north Senegal that grazing-induced variations in the composition of annual plants changed the relationship between the seasonal integral of NDVI and the end of the season biomass.
Considering the regression models with the different aggregation levels, all stratified models show a notably better performance than the global model (Table 1). There is a tendency for increasing R 2 in fitting values with increasing complexity of the model (increasing number of coefficients to be adjusted). While the global single-unit model has an R 2 of 0.33, the models with five to seven units have R 2 values in the range 0.42-0.44, and the models with 10 or 11 units show R 2 values in the range of 0.50-0.51. However, it is not strictly the case that an increase in calibration units always causes an increase in R 2 . The biophysical model is characterised by the highest R 2 without having the highest number of calibration units ( Figure 5). This indicates that the type of the stratification plays an important role in determining the performances of the model. number of calibration units ( Figure 5). This indicates that the type of the stratification plays an important role in determining the performances of the model. The performance of the different models in prediction is indicated by the R 2 cv and the RMSEcv. The global model has by far the lowest R 2 cv and the highest RMSEcv. The models building on existing thematic stratifications show increasing R 2 cv and decreasing RMSEcv, with increasing numbers of calibration units. The biophysical model deviates from these relationships and shows the highest R 2 cv (0.47) and the lowest RMSEcv (393 kg ha −1 ), while having nine calibration units. The second highest R 2 value for the department model appears to be due to overfitting as this model shows the largest drop when changing to the prediction mode (R 2 compared with R 2 cv).  The performance of the different models in prediction is indicated by the R 2 cv and the RMSE cv . The global model has by far the lowest R 2 cv and the highest RMSE cv . The models building on existing thematic stratifications show increasing R 2 cv and decreasing RMSE cv , with increasing numbers of calibration units. The biophysical model deviates from these relationships and shows the highest R 2 cv (0.47) and the lowest RMSE cv (393 kg·ha −1 ), while having nine calibration units. The second highest R 2 value for the department model appears to be due to overfitting as this model shows the largest drop when changing to the prediction mode (R 2 compared with R 2 cv ).
Note that the regressions over elementary units of the models vary in performance. For example, the biophysical model shows a range of R 2 values from 0.34 to 0.73, and, in cross-validation, the values range from 0.25 to 0.67. Elementary model units with a very low R 2 indicate that it is difficult to establish a good relationship between measured biomass and the RS-based biomass proxy in this area. This is not necessarily related to the small number of sites/samples in these units (e.g., the elementary model units of the biophysical model with the lowest R 2 are represented by eight and 10 sites, respectively). These low R 2 values could be the result of unsuitable stratification or the reasons mentioned for the single sites (data quality issues of certain samples, grazing, changing species composition).
In summary, the stratification of the study area in several calibration units increases the performance of the global regression model. An improved model performance for the estimation of biomass based on seasonal RS-derived indicators and in situ measurements through the subdivision of the study area into biophysically meaningful units was also shown for the Sahel in Senegal. Diouf et al. [14] stratified their study area by ecoregions and increased the R 2 from 0.68 to 0.77, applying a multi-linear model of three phenology-related parameters.
The biophysical model appears to be the most suitable among the tested stratifications for the estimation of end-of-the-season biomass in Niger as it shows the highest R 2 cv and the lowest RMSE cv . Furthermore, this model presents a biophysically meaningful stratification of the study area that can be easily retrieved from the same RS data that are used to build the model and does not require ancillary data.
Compared with field surveys, our approach provides continuous biomass estimates across the study region and not just point information like other RS-based methods. The advantage of the proposed method compared with the current RS-based approach of the Niger Ministry of Livestock is that it can be applied before the field measurements are made available. This translates to a time-gain of two to four weeks that the field trips normally last. Therefore, the produced maps of estimated biomass could be used for the planning of more in-depth field missions and for the better management of rangeland resources. Additionally, our model can also serve as a backup solution in the event that field surveys are not carried out in a specific year or a specific region. For example, the number of available measurement sites per year varied between 20 and 54 during the study period (see Figure S4) and the frequency of measurements at a single site between four and 14. However, we strongly recommend to continue the field measurements as they provide additional information like species composition and allow the recalibration of the model, yielding to an improved statistical reliability of the model.
Three ways to improve the model can be tested: first, incorporating more seasonal metrics in a multi-linear model as proposed by Diouf et al. [14]; second, utilising a seasonal metric or a combination of parameters that can be retrieved earlier in the season and therefore allow a biomass estimation before the end of the season; and, third, using a combination of NDVI-or FAPAR-derived seasonal metrics with meteorological data [17].

Estimated Biomass Map
Sample maps of estimated biomass at the end of the growing season derived from the biophysical model are presented in Figure 6 for the years 2004 and 2007. As shown by field measurements (see Section 2.2.2), 2004 was characterised by a very low spatial average of measured biomass, while 2007 showed very high measured biomass values. The biomass maps highlight the high inter-annual but also spatial variability of biomass production in the region, underlining the importance of flexible and production-adapted rangeland management. estimation before the end of the season; and, third, using a combination of NDVI-or FAPAR-derived seasonal metrics with meteorological data [17].

Estimated Biomass Map
Sample maps of estimated biomass at the end of the growing season derived from the biophysical model are presented in Figure 6 for the years 2004 and 2007. As shown by field measurements (see Section 2.2.2), 2004 was characterised by a very low spatial average of measured biomass, while 2007 showed very high measured biomass values. The biomass maps highlight the high inter-annual but also spatial variability of biomass production in the region, underlining the importance of flexible and production-adapted rangeland management.  Figure 7 shows the estimated biomass map for 2015 as the most recent example, indicating one possible use of such biomass maps. In addition to the biomass for 2015, the biomass anomaly for this year was also calculated as the difference between the estimated biomass of 2015 and the long-term average (in our example, the mean of the years 2001-2014). The anomaly map highlights areas with above or below average biomass production. With this information at hand, the relevant authorities could make certain rangeland management decisions. In the 2015 example, there is a large negative  Figure 7 shows the estimated biomass map for 2015 as the most recent example, indicating one possible use of such biomass maps. In addition to the biomass for 2015, the biomass anomaly for this year was also calculated as the difference between the estimated biomass of 2015 and the long-term average (in our example, the mean of the years 2001-2014). The anomaly map highlights areas with above or below average biomass production. With this information at hand, the relevant authorities could make certain rangeland management decisions. In the 2015 example, there is a large negative anomaly in the Eastern Tahoua and Southern Agadez regions, while there is a significant positive anomaly in the Eastern Zinder and Western Diffa regions (circled in Figure 7). Together with information about the availability of water for the animals and the location and size of herds, the authorities could then assess if it is necessary to reduce the livestock stocking rates or to provide additional feed (e.g., enriched alfalfa, wheat bran, cotton seed cake) in areas with large negative anomalies. In areas with pronounced positive anomalies authorities could evaluate if it would be appropriate to channel additional animals to these areas (simultaneously strengthening the livestock safety in view of the high concentration of animals) or if some fire protection measures are required (e.g., firewalls and awareness-raising campaigns). The report on the pastoral campaign from the Ministry of Livestock [46] from November 2015 reported variable forage production in the pastoral areas of Niger, ranging from good to mediocre production. Large pockets with low watering resulting in low biomass production were observed in South Ingall, Aderbissenat, North-West Abalak, Tchintabaraden, and Bermo, which are areas of large negative anomalies in Figure 7. Having considered the biomass production and the number of livestock (resident livestock and exceptional cases of livestock from Malian refugees), the report concluded with a note of concern for the 2015/2016 dry season as all regions were experiencing an overall forage deficit. areas of Niger, ranging from good to mediocre production. Large pockets with low watering resulting in low biomass production were observed in South Ingall, Aderbissenat, North-West Abalak, Tchintabaraden, and Bermo, which are areas of large negative anomalies in Figure 7. Having considered the biomass production and the number of livestock (resident livestock and exceptional cases of livestock from Malian refugees), the report concluded with a note of concern for the 2015/2016 dry season as all regions were experiencing an overall forage deficit.

Conclusions
In this study a phenology-based predictive model to estimate rangeland biomass in Niger was developed. The cNDVI during the growing season, derived from 10-day MODIS NDVI data applying a model-fit phenology retrieval method, was used as a remote sensing proxy for biomass production. The relationship between this variable and the measured herbaceous biomass at the end of the growing season was modelled by linear regression. Different spatial aggregation levels for the model calibration were tested to find the most suitable one for biomass prediction.
A general, but not strict, tendency for increased estimation performance with increased model complexity (in terms of number of parameters to be calibrated) was observed. This indicates the importance of additional and spatially heterogeneous agro-ecological unobserved effects (which

Conclusions
In this study a phenology-based predictive model to estimate rangeland biomass in Niger was developed. The cNDVI during the growing season, derived from 10-day MODIS NDVI data applying a model-fit phenology retrieval method, was used as a remote sensing proxy for biomass production. The relationship between this variable and the measured herbaceous biomass at the end of the growing season was modelled by linear regression. Different spatial aggregation levels for the model calibration were tested to find the most suitable one for biomass prediction.
A general, but not strict, tendency for increased estimation performance with increased model complexity (in terms of number of parameters to be calibrated) was observed. This indicates the importance of additional and spatially heterogeneous agro-ecological unobserved effects (which might relate to grazing, species composition, optical soil properties, etc.) in modifying the relationship between cNDVI and herbaceous biomass at the end of the season. The biophysical model showed the best performance with respect to the predictive power (R 2 cv = 0.47) and RMSE cv value (398 kg·ha −1 ), without having the highest number of parameters to be calibrated. All models that were tested using several calibration units outperformed the simple model, using a unique linear relationship for the entire study area.
The presented approach can be applied for the timely production of estimated biomass production maps at the end of the growing season and before the field measurements are made available. This would mean a time gain of two to four weeks, which is the length of time the field trips normally last. Therefore, such maps could be used for the planning of more in-depth field missions, for the better management of rangeland resources, and for timely decisions on aid allocation and fire prevention. Additionally, the approach can serve as a backup solution in the event that field surveys are not carried out in a specific year or a specific region.
In summary, the presented approach should be used as a supplement to and not as a replacement for field measurements. We highly recommend that field measurements are continued for the following reasons:

•
They provide additional information such as species composition that currently cannot be retrieved by RS data. • The presented model to estimate biomass was calibrated with field data obtained over a 15-year time period. However, the full variability in biomass production in the study area may not be captured by this 15-year period, as this variability may also be affected by longer-term climate variability. As longer satellite image time series become available in the coming years, future research should analyse this temporal variability. In general, whenever a new field survey becomes available, the model should be recalibrated as more data improve the statistical reliability of the model.

Supplementary Materials:
The following supplementary data are available online at www.mdpi.com/ 2072-4292/9/5/463/s1. Figure S1: Agro-ecological zones of Niger; Figure S2: Performance of different ISODATA-clustering-based stratifications with respect to R 2 cv and RMSE cv ; Figure S3: Results of the best-performing ISODATA clustering compared with elevation from SRTM data; Figure S4: Available number of biomass measurements during the study period; Table S1: Model results for ISODATA clustering with a fixed number of classes.