Iskander Benhadj
email: benhadj@cesbio.cnes.fr

Benoît Duchemin

Philippe Maisongrande

Vincent Simonneaux

B. DUCHEMIN ‡, P. MAISONGRANDE ‡, V. SIMONNEAUX ‡ † S Khabba

Ghani Chehbouni

A Chehbouni

Automatic unmixing of MODIS multitemporal data for inter-annual monitoring of land use at a regional

Keywords: NDVI, land use, semi-arid, linear unmixing, endmembers, MODIS

Introduction

Changes in Land Use and Land Cover (LULC) is a major issue in Environmental Science, interconnected with many question concerning climate change, carbon cycle and biodiversity [START_REF] Aspinall | A Land Use and Land Cover Change science strategy (Summary of a Workshop Organized by the US Climate Change Science Program (CCSP) Land Use Interagency Working Group (LUIWG)[END_REF][START_REF] Lepers | A synthesis of information on rapid land-cover change for the period 1981-2000[END_REF]. The monitoring of LULC is also vital for managers and policy makers to make informed decisions regarding the sustainability of agriculture and provision of safe drinking water, especially in semi-arid areas. Remote sensing is very well-suited to achieve this monitoring since it allows observations regularly distributed in space and time [START_REF] Rogan | Remote sensing technology for mapping and monitoring land-cover and land-use change[END_REF]Chen 2004, Prenzel 2004).

Multi-temporal images are widely investigated for mapping and monitoring land-cover and land-use changes. At the present time, time series of images can be obtained at a high spatial resolution by programming a series of SPOT or FORMOSAT-2 acquisitions. These images with both high spatial resolution (~10 m) and high temporal repetitivity (a few days) offer strong opportunities to monitor land surfaces over small areas: 25x25 km² for FORMOSAT-2, 60x60 km² for SPOT. However, constraints related to acquisition, cost and processing often prevent the use of high spatial resolution data. Multi-temporal data acquired by low or moderate spatial resolution sensors such as NOAA-AVHRR, SPOT-VEGETATION or TERRA-MODIS are thus preferred for regional and continental studies (e.g. Lambin andEhrlich 1997, Hansen et al. 2000;[START_REF] Lunetta | Land-cover change detection using multi-temporal MODIS NDVI data[END_REF][START_REF] Matsuoka | Land cover in East Asia classified using Terra MODIS and DMSP OLS products[END_REF][START_REF] Stibig | A land-cover map for South and Southeast Asia derived from SPOT-VEGETATION data[END_REF]). Indeed, they offer a costless global coverage of the Earth on a daily basis. However, the spatial resolution of large field of views sensors -from 250m for MODIS to 1 km for VEGETATION and AVHRR -is generally much higher than the size of homogeneous areas (units) at the Earth surfaces. These sensors generally provide images with pixels that include a mixture of different units (mixed pixels). Consequently, the use of low spatial resolution data for a directly monitoring of LULC is not straightforward. Furthermore, conventional classification approaches based on signature clustering (like maximum likelihood, [START_REF] Richards | Remote Sensing Digital Image Analysis[END_REF] are not suitable since they aim to identify an unique class for each pixel.

For these reasons, the linear unmixing model has been developed [START_REF] Adams | Spectral mixture modelling: A new analysis of rock and soil types at the Viking Lander Site[END_REF][START_REF] Smith | Vegetation deserts: I. Regional measure of abundance from multispectral images[END_REF][START_REF] Elmore | Quantifying vegetation change in semiarid environments: Precision and accuracy of spectral mixture analysis and the Normalized Difference Vegetation Index[END_REF] based on the following assumption: the signature of a mixed pixel results from a linear combination of the distinctive signatures (endmembers) that are representative of the various land surfaces included in the study area. These typical signatures must describe as well as possible a pure component having meaningful features for an observer [START_REF] Strahler | On the nature models in remote sensing[END_REF]. Knowing these signatures is a prerequisite for applying the linear unmixing model [START_REF] Cross | Subpixel Measurement of Tropical Forest Cover Using Avhrr Data[END_REF][START_REF] Quarmby | Linear mixture modelling applied to AVHRR data for crop area estimation[END_REF][START_REF] Foody | Sub-pixel land cover composition estimation using a linear mixture model and fuzzy membership functions[END_REF][START_REF] Milton | The identification of reference endmembers using high spatial resolution multi-spectral images[END_REF]. Unmixing approaches can be divided into two categories depending on how the endmembers are estimated:

• Supervised approaches use the spectral signatures of endmembers as a priori information. These typical spectra can be collected at field or laboratory to define predefined library endmembers [START_REF] Adams | Classification of Multispectral Images Based on Fractions of Endmembers -Application to Land-Cover Change in the Brazilian Amazon[END_REF][START_REF] Roberts | Mapping chaparral in the Santa Monica Mountains using multiple endmember spectral mixture models[END_REF][START_REF] Smith | Vegetation deserts: I. Regional measure of abundance from multispectral images[END_REF]. They can also be derived from high spatial images using a training data set (small region where the land use is known). The use of predefined libraries may be not appropriate since differences in the acquisition conditions (e.g. sun-target-sensor geometry, atmospheric effects) may occur between endmembers and the data to be unmixed [START_REF] Song | Monitoring forest succession with multitemporal Landsat images: Factors of uncertainty[END_REF].

• In unsupervised approaches (see [START_REF] Plaza | A quantitative and comparative analysis of endmember extraction algorithms from hyperspectral data[END_REF] for a review), the identification of endmembers is automated. The common point in unsupervised algorithms is that they search endmembers directly from images [START_REF] Atkinson | Mapping sub-pixel proportional land cover with AVHRR imagery[END_REF][START_REF] Elmore | Quantifying vegetation change in semiarid environments: Precision and accuracy of spectral mixture analysis and the Normalized Difference Vegetation Index[END_REF][START_REF] Ridd | Exploring a V-I-S (Vegetation-Impervious Surface-Soil) Model for Urban Ecosystem Analysis through Remote-Sensing -Comparative Anatomy for Cities[END_REF][START_REF] Wessman | Detecting fire and grazing patterns in tallgrass prairie using spectral mixture analysis[END_REF]. In this case, the endmembers are retrieved at the same scale and conditions than the data to be unmixed.

The temporal variability of the observations is generally not considered in the abovementioned studies, though it is also an important source of information. In particular, the time courses of vegetation indices such as the Normalized Difference Vegetation Index NDVI allow to monitor the phenology of vegetation [START_REF] Gutman | Global Land Monitoring from Avhrr -Potential and Limitations[END_REF][START_REF] Justice | The Moderate Resolution Imaging Spectroradiometer (MODIS): land remote sensing for global change research[END_REF][START_REF] Duchemin | Potential and limits of NOAA-AVHRR temporal composite data for phenology and water stress monitoring of temperate forest ecosystems[END_REF]). This may be very useful for discriminating land classes. Differences in phenology depicted by vegetation indices can be used to map land surfaces using low spatial resolution data (e.g. [START_REF] Kerdiles | Noaa-Avhrr Ndvi Decomposition and Subpixel Classification Using Linear Mixing in the Argentinean Pampa[END_REF][START_REF] Cardot | Functional approaches for predicting land use with the temporal evolution of coarse resolution remote sensing data[END_REF][START_REF] Ballantine | Mapping North African landforms using continental scale unmixing of MODIS imagery[END_REF], Knight et al. 2006). These studies showed that: 1) land use maps are more accurate when vegetation indices are used instead of reflectances; 2) the use of NDVI with a linear approximation for its combination results in minor inaccuracies; 3) linear unmixing provides satisfactory results when the number of endmembers is limited. These considerations, which are of prime importance in unmixing procedure, are accounted for in this study.

In this context, the primary objective of this study is to evaluate the potential of MODIS data for monitoring the land use on the semi-arid Tensift/Marrakech plain. A secondary objective is to analyse the space-time variability of land classes in relation with water availability. The methodology is based on the unmixing of MODIS multi-temporal NDVI images. Land use maps are evaluated using ground data and high spatial resolution images, and their space-time variability is analysed together with information on irrigation water.

Research Design

The methodology is an unsupervised unmixing approach based on a statistical analysis for identifying endmembers directly from MODIS multi-temporal images at 250 spatial resolution (MOD13Q1 product, i.e. 16-day NDVI composite images by CV-MVC algorithm, [START_REF] Huete | Overview of the radiometric and biophysical performance of the MODIS vegetation indices[END_REF]. The algorithm first extracts typical NDVI profiles, then selects the endmembers amongst these profiles based on their ability to reproduce the space time variability of MODIS NDVI time series. The approach requires the two following assumptions: (1) pure pixels can be identified at the 250m resolution and

(2) endmembers are stationary over the Tensift-Marrakech plain.

The approach is set up to retrieve the fractions (surface covered by homogeneous units within each pixel) of three classes: orchard, non-cultivated areas and annual crop. These classes are predominant in the study area, they display distinct phenological features and they encompass the range of crop water needs: non-cultivated areas (no needs), annual crops (water needs ~ 400 mm/y) and orchards (water needs ~ 1000 mm/y). MODIS estimates are quantitatively evaluated against ground truth collected on a 9 km² area and a reference land use map derived from a time series of high spatial resolution images (SPOT and Landsat). These data were collected during the 2002-2003 agricultural season. The evaluation is based on classical statistical variables (correlation R², efficiency EFF, RMSE and bias) computed between land use fractions estimated with MODIS and derived from the validation data sets at 1 km resolution. In order to test the robustness of the algorithm, the performance of the algorithm is also discussed from the results obtained with the whole MODIS data set (2000-2006 period). Here we analyse the inter-annual variability of both endmembers and land use maps using rainfall and irrigation data as an indicator of water availability and vegetation growth.

Materials and Methods

In this section, we present the study area, the ground and satellite data, and the linear unmixing algorithm.

Study area and ground data

The study area is the eastern part of the semi-arid Tensift plain, a 3000 km² region located in center of Morocco (figure 1). The climate of this region is arid, with annual rainfall around 250 mm/year and a very high evaporative demand around 1500mm/year [START_REF] Duchemin | Monitoring wheat phenology and irrigation in Central Morocco: On the use of relationships between evapotranspiration, crops coefficients, leaf area index and remotely-sensed vegetation indices[END_REF], Chehbouni et al. 2007).

According to the regional public agency in charge of agricultural water management (ORMVAH), there are three dominant land classes that represent more than 80% of land surfaces: (1) orchards, most of it perennial (olive and citrus trees); (2) cereal crops, mainly wheat, to less extent barley; (3) non-cultivated areas. Additional land classes include forages (mainly alfalfa, colza and oat), vineyards, broad-leave orchards (apple, apricot and peach trees), and small vegetable crops.

[Insert Figure 1 about here] The High-Atlas mountain range experiences much higher precipitations and provides irrigation water to the plain [START_REF] Chaponniere | A combined high and low spatial resolution approach for mapping snow covered areas in the Atlas mountains[END_REF], Chehbouni et al. 2007). There are three types of irrigation systems: the modern network connected with dams, the traditional network, and pumping stations [START_REF] Duchemin | Agrometerological study of semi-arid areas: an experiment for analysing the potential of FORMOSAT-2 time series of images in the Marrakech plain[END_REF]). The main irrigated areas are supplied by dam water and managed by ORMVAH. They cover about 1200 km² with three distinct sub-regions (figure 1):

• The western NFIS sub-region, mainly cropped with orchards on fields of irregular size (~ 100 m² to ~ 10 ha);

• The central Haouz sub-region, mostly cropped with cereals, where the landscape appears rather uniform with relatively larger fields (3-4 ha);

• The eastern Tessaout sub-region, very patchy with a mixture of various annual crops and orchards cultivated on very small fields (100 to 1000 m²).

In order to evaluate land use maps, we use two sets of ground data collected during the 2002-2003 agricultural season. The first one is composed of 151 individual fields spread over the study area divided as following: 11 plots of orchard on bare soil, 80 plots of orchard on annual crop, 28 plots of non-cultivated areas and 32 plots of annual crop (see [START_REF] Simonneaux | Using high resolution image time series for crop classification and evapotranspiration estimate over an irrigated area in south Morocco[END_REF]). The second one exhaustively covers a 3 x 3 km² area within the Haouz sub-region during the 2002-2003 agricultural season (see [START_REF] Duchemin | Monitoring wheat phenology and irrigation in Central Morocco: On the use of relationships between evapotranspiration, crops coefficients, leaf area index and remotely-sensed vegetation indices[END_REF]. It is composed of 313 plots divided as following: 5 plots of orchard, 67 plots of non-cultivated areas and 241 plots of cereal crops (wheat and barley).

In order to study the space-time variability of land classes, we analyse data on dam irrigation water and precipitations. ORMVAH collects the annual amount of dam irrigation water supplied to the three sub-regions. As it is difficult to exactly know when and where irrigation occurs, we assume a uniform distribution: the amount of dam irrigation water is divided by the total area of each sub-region to provide average values in mm. Precipitations are collected from a network made of about 20 raingauges stations spread over the plain. There is a large seasonal variability of rainfall, both in terms of annual quantity and of seasonal distribution: accumulated values of 140 mm for the driest years (2000-2001 and 2004-2005) against 300 mm for the most humid years (2003-2004 and 2005-2006); early rainfall in 2003-2004 or delayed rainfall in 2001-2002.

Satellite data

High spatial resolution data are used to produce a reference land use map in order to evaluate classification maps obtained with MODIS data. We use a SPOT5 with a revisit time of approximately three weeks. These images were geometrically corrected using GPS ground control points and resampled to 30m. The radiometric processing (calibration and atmospheric correction) was performed using reflectance values recorded at field [START_REF] Duchemin | Monitoring wheat phenology and irrigation in Central Morocco: On the use of relationships between evapotranspiration, crops coefficients, leaf area index and remotely-sensed vegetation indices[END_REF], Simonneaux et al. 2007).

Terra-MODIS data are freely available from the NASA website (http://delenn.gsfc.nasa.gov/). We have downloaded 16-day composite images because they display geometric problems. All images are free of clouds. This is expected since the time step of compositing is rather long (16 days) and the cloudiness is low in the study area, around 30% [START_REF] Hadria | Monitoring of irrigated wheat in a semi-arid climate using crop modelling and remote sensing data: Impact of satellite revisit time frequency[END_REF].

Reference land use map (2002-2003 season)

The reference land use map is derived from high spatial resolution data on the common area between the Landsat images, the SPOT ones and the study area (about 1500 km², see figure 2). The classification identifies the three predominant land classes using a two-step procedure:

1) The orchards are depicted on the 2.5m panchromatic SPOT image using the "Olicount" software [START_REF] Simon | Computer-assisted recognition of Olive trees in digital imagery[END_REF]. The software operates with a set of input parameters that essentially define the morphology of trees (shape) and their radiometry (gray level). This first class groups all the areas where trees are detected, including case of intercropping (trees + wheat or trees + alfalfa) and the natural vegetation that may also grow between the trees or in the understory. The reference land use map is evaluated against the ground truth collected on individual fields (see §3.1). According to the confusion matrix (table 1), the overall accuracy, i.e. the number of well-classified pixels divided by the total number of pixels, is around 78%, with very low omission errors for the class orchard on bare soil (about 10%) and for the class annual crop (about 3%). Two types of confusion are detected: 1) between annual crop and orchard on annual understory, and 2) between bare soil and annual crop. The causes of these confusions were discussed in [START_REF] Simonneaux | Using high resolution image time series for crop classification and evapotranspiration estimate over an irrigated area in south Morocco[END_REF] and [START_REF] Benhadj | High spatial and temporal resolution FORMOSAT-2 images: first results and perspectives for land cover mapping of semi-arid areas (Marrakech/Al Haouz plain)[END_REF]. They are related to the disparities that exist for a same land class, which causes overlapping of signatures between the three land classes. For cereals, there is a large heterogeneity in cereal crop calendar as well as irrigation and fertilisation schedules. Non irrigated areas may include a wide range of vegetation type (colza, oat, grass). Finally, there are large variations of density and age in tree plantations, which may include an understory of vegetation cultivated as forage (wheat, grass, alfalfa…).

[Insert Table1 about here] The reference land use map is used for evaluating MODIS estimates for the 2002-2003 agricultural season at 1 km² scale. For this purpose, a co-registration between MODIS data and the reference land use map is done using an automatic correlation algorithm [START_REF] Benhadj | Land cover in semi arid area derived from NDVI images at high and low spatial resolution[END_REF]. Then the reference map is up-scaled at 1 km resolution by spatial averaging to obtain the fractions covered by orchards, bare soils and annual crops.

Linear unmixing of MODIS data

To predict the land use fractions of the three dominant land classes, the linear unmixing model is applied to MODIS multi-temporal NDVI images. The model calculates the NDVI of a mixed pixel as the sum of the NDVI values of the different land classes weighted by their corresponding fraction within the pixel (equation 1). We retrieve the typical NDVI time course of each land class (endmember) using the three-step procedure which is detailed below.

( ) ( ) ( ) t t NDVI t NDVI i j j ij i ε π + × = ∑ = 3 1 (1)
where i NDVI is the NDVI of MODIS mixed pixel i at the date t, ij π is the fraction of class j in pixel i, j NDVI is the endmember of class j (j = 1 to 3) and i ε is an error term of the pixel i.

Step 1. An unsupervised classification "k-means" [START_REF] Tou | Pattern Recognition Principles[END_REF]) is applied to MODIS multi-temporal images in order to group the pixels which have similar NDVI seasonal courses. The result is N mean NDVI profiles corresponding to N groups 1 of pixels. We set N to 20, which appears as a good compromise allowing a reasonable computing time cost while keeping a sufficient level of details to describe the NDVI space-time variability within the study area. Furthermore, the grouping of pixels with the same vegetation seasonality allows the reduction of local noise due to: (1) imperfect superimposition of MODIS data before temporal compositing, (2) inaccuracy in atmospheric correction and, (3) the variation in sun-target-sensor geometry between successive acquisitions.

Step 2. An iterative test is applied for all possible triplets of endmembers (three land classes) among the series of N mean NDVI profiles. The total number of iteration nb is

3 N C .
For each triplet, the land use fractions are retrieved for the remaining 17 (i.e. N-3) groups by minimizing the Root Mean Square Error (RMSE, equation 2) between the NDVI profiles observed by MODIS and those reconstructed from the endmembers.

( ) [ ] ∑ = × = T t i i t T RMSE 1 2 1 ε (2) With 0 ≥ ij π and 1 3 1 = ∑ = j ij π
Where T represents the number of MODIS data

Step 3. We calculate an error term ( k M , equation 3), which represents the ability of the triplet number k to explain the NDVI response for the 17 groups. Finally, the triplets are sorted according to this error term: the triplet for which k M is minimal is called triplet rank 1, the following is called triplet rank 2, etc.

( ) [ ] ( ) [ ] ∑∑ - = = × × - = 3 1 1 2 3 1 N i T t i k t T N M ε (3) With [ ] nb k , 1 ∈ and
( )! 3 ! 3 ! 3 - × = = N N C nb N
Once the endmembers are identified, they are assigned to the appropriate land use class and the surface covered by a class within a pixel (land use fraction) is retrieved by minimization (equation 2). This is applied pixel by pixel using land use fractions ranging from 0 to 1 and under the constraint that the sum of fractions is equal to 1.

We apply the algorithm using two different areas for the identification of endmembers.

The first one is the whole study area (figure 1). The second one is the reference area (figure 2) on which the reference land use map is available ( §3.3). In both case, the land use fractions maps are analysed at the scale of the whole area. MODIS estimates are evaluated against the reference land use map (see §3.3) and against the ground truth collected on the 3 x 3 km² area (see §3.1). In order to explain the difference between annual crop endmembers between the two investigated areas, we carry out a purity analysis. The pixels of each group resulting from the k-means classification are located in the reference land use map (figure 2) and their compositions are averaged.

Results and discussion

In this section, we successively present: a quantitative evaluation of the results obtained during the 2002-2003 agricultural season; a generalised analysis of inter-annual coherence and variability of the results through the 2000-2006 period; an error analysis with typical cases for which the results are not satisfactory.

Typical NDVI time series and endmembers (2002-2003 agricultural season)

The NDVI profiles of the 20 groups identified with K-means classification over the two areas of interest (whole and reference areas) can be discriminated through the combination of NDVI seasonal amplitude and average value (figure 3). It appears that the K-means method groups pixels according to the density of perennial vegetation (hierarchy of rather stable NDVI profiles with average values from 0.15 to 0.55) and according to the vegetation seasonality (contrast between high NDVI values during the agricultural season and low values in summer).

When looking at the endmembers (figure 3), it is noticeable that the algorithm tends to select the profiles that display extreme values and rejects intermediates ones.

Furthermore, the endmembers appear descriptive of the three dominant classes: the first one, with maximum NDVI values below 0.2, corresponds to the bare soil class; the second one, with NDVI always high (between 0.45 and 0.65), appears representative of a dense perennial vegetation (orchard class); the third one, with a large NDVI amplitude, can be associated to the class annual crop. The latter displays minimum values in November (at the sowing period), then a rapid increase to maximum values mid-March when cereal reaches full development, and a final decrease until June after total senescence of plants. This analysis makes easy to label each endmember.

[Insert Figure 3 about here] The case of annual crop is of particular interest since the endmembers are not the same for the two investigated areas (figure 3). In particular, there is a difference in the NDVI value at the beginning (September to November 2002, before day 90) and ending of the season (June to August 2003, after day 280). The level is around 0.25 for the endmember extracted on the whole area, while it is only 0.18 for the endmember extracted on the reference area. This last endmember appears to be more characteristics of annual crop, for which minimal NDVI values are close to those of bare soil (~0.15) out of the agricultural season.

In order to explain the difference between annual crop endmembers, their purity are analysed (table2). The endmembers display a high proportion of either bare soil or annual crop or orchard for the two areas comparing to the remaining 17 NDVI profiles non selected as endmembers. One exception is detected for the class annual crop when the whole area is considered (72% of annual crop in group 3 that is selected as endmember against 88% in group 20, see the main left column of table 2). The difference of endmembers purity between the two areas is small for the bare soil and orchards classes, but large for annual crop (purity of 88% for the reference area against 72% for the whole area, compare group 3 in the two main columns of table 2). This difference is due to significant presence of orchard in the annual crop endmember derived over the whole area (~27%, against ~9% for the reference area). When the whole area is used, the automatic extraction algorithm selects groups that include pixels of the Tessaout sub-region, where there is a mixture of olive orchards and annual crops cultivated on very small fields. In contrast, when the identification of endmembers is restricted to the reference area, the algorithm selects pixels in the irrigated Haouz subregion where fields are mainly cropped with cereals and of larger size. Therefore, this analysis demonstrates that: (1) our working hypothesis, i.e. pure pixels may exist at the spatial resolution of 250m, is valid; (2) the automatic extraction algorithm is able to identify the most pure areas; (3) there is an advantage to derive the endmembers on the reference area compared to the whole area.

[Insert table 2 about here]

Quantitative evaluation of land use map (2002-2003 agricultural season)

The comparison of land use fractions estimated with MODIS and the reference land use map (figure 4) shows the consistency of areas with low and high fractions between the two maps. This is true for the three land classes: high proportion of bare soil at South-West; high proportion of annual crop near High-Atlas foothills and on the Haouz and Tessaout irrigated areas in the eastern part; high proportion of orchard near the Tensift river at North and within the NFIS irrigated area at West. Average land use fractions derived from reference and estimated maps display an overall agreement (table 3), which denotes the global ability of the algorithm to describe the study area using three dominant land classes. However, the algorithm slightly underestimates the orchard and the annual crop fractions at the benefit of bare soil fractions when the whole area is considered. This underestimation is attenuated when the reference area is used to derive the endmembers.

[Insert Figure 4 about here]

[Insert Table 3 about here] The quantitative comparison of MODIS and the reference land use map (table 4 andfigure 5) shows that the two land use fractions always well correlate (R² around 0.8 with a minimal value of 0.68), and the efficiency is generally largely positive (>0.65). When the reference area is used to derive the endmembers, the method gives more accurate estimates of bare soil and orchard fractions (lower RMSE and bias, larger efficiency).

For both areas, the estimates of orchard fractions appear less accurate than for the two other classes (efficiency of 0.65-0.7 instead of 0.80). This is likely due to the fact that the orchard class is rather heterogeneous because trees are of different nature, age and spacing, with possible case of inter-cropping. In contrast, the endmember associated to this class is representative of dense perennial vegetation (mainly old olive and citrus tree with low spacing between crown). Despite this limitation, we consider that land use fractions are correctly estimated, though the study area is only described by three typical NDVI profiles.

[Insert Table 4 about here]

[Insert Figure 5 about here] Finally, the comparison of MODIS land use fractions and the ground truth available over the 9 km² area shows a global agreement of land use fractions for all classes (figure 6), with few orchards (less than 2% of the 9 km², see table 5). For the two others classes, we obtain accurate results, with R² larger than 0.85 and RMSE lower than 0.1.

The accuracy of estimates is improved when the endmembers are derived on the reference area (RMSE of 0.07 against 0.09 in figure 6). 'normal' years (2000-2001, 2002-2003, 2004-2005 and 2005-2006).

• For the 2003-2004 season, the NDVI display an early NDVI from 0.2 to 0.4 between November and December ("03-04 (rank 1)" profile in figure-bottom-right).

This pattern also appears coherent with the seasonal distribution of rainfall. Heavy rainfall at the very beginning of the season resulted in an early sowing or growth of natural vegetation. Here the analysis of other NDVI profiles allows to identify a substitute to represent the phenology of cereal crops. This endmember ("03-04 (rank 2)" profile in figure 8-bottom-right) is similar to the ones observed for the 'normal' years, and results in a low unmixing error (second rank in the minimisation process).

[insert Figure 8 about here]

Spatio-temporal variability of land use maps (2000-2006 period)

A visual examination of land use fractions maps (figure 9) shows that the algorithm always detects the same region with low or high proportion of each class. Orchard fractions appear especially stable during the six years, in coherence with the duration of tree plantations. On the contrary, there are some compensations in the fractions of the two other classes (bare soil and annual crop). In particular, there is a high proportion of bare soil and a low proportion of annual crop for the 2001-2002 agricultural season compared to others. These compensations are analysed on what follows.

[Insert Figure 9 about here] Land use statistics are calculated for the six years of study by averaging fractions over each of the three irrigated sub-regions (table 6). One can see that the proportion of orchard is quite stable, around 37% for NFIS, 18% for Haouz and 32% for Tessaout.

These values appear coherent with the qualitative knowledge of the study area ( §3.1).

Except for the 2001-2002 season, bare soil fractions are rather stable, between 50 and 56% for the NFIS sub-region, 35 and 46% for the Haouz and between 16 and 21% for the Tessaout. The variation of annual crop fractions around its average value is of the same order. The 2001-2002 season is very particular with an important reduction of annual crop fractions, by a factor 2.5 within NFIS (4% in 2001-2002 against 10% the other years) and Tessaout (20% against 45-50%) and a factor 5 within Haouz (8% against 40%).

[Insert Table 6 about here] The anomaly detected in annual crop fractions for the 2001-2002 agricultural season appears as an indicator of the water shortage experienced this year. We illustrate this for the Haouz sub-region, where the anomaly is of maximal amplitude (figure 10). The limitation of irrigation water during the driest year (annual average of 30 mm in 2001-2002 instead of 130 mm for the other years) results in a large decrease of annual crop fractions (by about 30%) and a large increase of non-cultivated areas (by about 30%).

The orchard fractions appear stable despite the shortage of irrigation water, consistent with the fact that orchards are irrigated in priority.

[Insert Figure 10 about here]

Error and limitation analysis

In order to identify the limitations of the approach, we calculate the relative error (RRMSE, equation 4) between MODIS observations ( ( ) obs t NDVI

) and the NDVI reconstructed from the linear combination of the endmembers associated to their land use fractions (

( ) sim t NDVI
). This criteria allows us to quantify the ability of the three endmembers to reproduce MODIS NDVI space-time patterns over the study area. Maps of RRMSE are computed for each season and averaged over the six seasons (figure 11).

(

)

100 1 1 2 × - = ∑ = obs T t sim obs t NDVI mean t NDVI t NDVI T RRMSE (4)
[Insert figure11 about here]

It can be seen that the MODIS NDVI time courses are generally well reproduced (figure 11). The histogram associated to the spatial variation of RRMSE displays a peak centred around a value of 10%, with 90% of pixels have a value of RRMSE lower to 20%. This confirms the efficiency of the algorithm to recover NDVI space time variations, but some anomalous pixels display high errors (RRMSE>40%). These pixels are mainly located in the NFIS irrigated sub-region at the western part of the study area.

There are two main cases where the capacity of the algorithm to fit MODIS observations is low:

• In case 1, the NDVI time course displays two peaks at the middle (January) and at the end (April) of the agricultural season; this indicates successive cropping of vegetables with a short growing period;

• In case 2, the NDVI time course displays an inverse pattern than the one of annual crop, with a large growing period between April and January; such pattern is consistent with the phenology of deciduous tree crops (apricot, apple, peach trees) and vineyards.

The two previous confusions concern a small part of the study area (0.2% with RRMSE>40%). Further investigations would be necessary to analyse the performance of the algorithm using more endmembers and more NDVI profiles as an input of the minimisation procedure (N>20 in equation 3). However, this may result in larger computation time and additional compensations/overlaps between land use classes.

Conclusion

In this study, we investigate the potential of time series of MODIS data (MOD13Q1 product, i.e. 16-day NDVI composite images by CVMVC algorithm, [START_REF] Huete | Overview of the radiometric and biophysical performance of the MODIS vegetation indices[END_REF] to monitor the land-use of the Tensift plain, a semi-arid region located in the surrounding of the Marrakech city. MODIS data offers a costless coverage of the Earth with a high temporal resolution, but its spatial resolution (250m) is large compared to the average field size in the study site. Thus, we develop an approach based on the linear unmixing of multi-temporal MODIS data. In this approach, the identification of endmembers -key point in linear unmixing -is performed on an annual basis following a two-step procedure: 1) pixels are grouped according to the vegetation seasonality;

(2) the set of groups that displays the best ability to explain all NDVI time courses are automatically extracted using a statistical analysis. Some advantages can be mentioned here. Firstly, there is no need of extra information such as a training set where the land use is known. Secondly, there are no substantial differences in the acquisition conditions between endmembers and the data that are unmixed. Thirdly, the regional conditions on which the vegetation growth (e.g. dry or humid year) are integrated to the endmembers.

This procedure provides a continuous description of the land use in term of fractions of three classes (orchard, annual crop, non-cultivated areas) and on an annual basis (September to August, i.e. the agricultural season). These three classes are the most important for agricultural water management because they are predominant and they corresponds to very different water needs. The use of these three broad categories also facilitate the analysis of the inter-annual variability of MODIS estimates of land use fractions as well as its evaluation against additional data sets (ground truth and high spatial resolution images).

The analysis of typical NDVI profiles firstly demonstrates that our working assumption, Maps of land use fractions are in coherence with the qualitative knowledge of the study area, in particular for the three main irrigated sub-regions (NFIS, Haouz and Tessaout).

Using both high spatial resolution data and ground truth, we quantify the error in land use fractions to around 0.1 at 1km spatial resolution (2002-2003 season). The analysis of land use maps derived for the six successive agricultural seasons (2000-2001 to 2005-2006) also confirms the performance of the approach. The orchard class is logically the most stable, with fractions around 37%, 18% and 32% for the NFIS, Haouz and Tessaout sub-regions, respectively. The compensations observed between the fractions of bare soil and annual crop show a high degree of space-time coherence with irrigation statistics. In particular, the algorithm retrieves a large reduction of annual crops after the severe drought that occurs at the beginning of the period of study. These results are promising in the perspective of the regional monitoring of water resources in the semi-arid Tensift/Marrakech plain.

Finally, the examination of some anomalous NDVI profiles, i.e. which are not well reproduced by the linear unmixing model, denotes the incapacity of the algorithm to describe the phenology of particular crop types (e.g. vineyards, vegetable crops).

Inclusion of other land use components would provide additional information and possibly more accurate results. Further tests should be performed to identify the optimal number of both the endmembers and the groups of pixels used as endmembers potential candidates. In this perspective, the availability of time series of images with both high spatial resolution and high temporal repetitivity (e.g. FORMOSAT-2, GMES-Sentinel, RapidEye or Venµs) would offer additional opportunities.

FIGURE CAPTIONS

Figure 1. Delimitation of the whole study area (in red) and its three main irrigated subregions -NFIS (in yellow), Haouz (in black) and Tessaout (in cyan) -on a Landsat7 image. The white square represents the coverage of Landsat and SPOT4/5 images. -2001, 2002-2003, 2004-2005 and 2005-2006 annual crop endmembers. In figure b (bottom), "03-04 (rank1)" and "03-04 (rank2)" correspond to the endmembers linked to the 1 st and the 2 nd ranks in the minimisation procedure, respectively. (1)

(2)

  The algorithm is applied to a six-year archive of MODIS NDVI to obtain maps of land use fractions on a yearly basis, from agricultural season2000-2001 to 2005-2006. The algorithm is applied on two different areas, the whole study area and a subpart of the study area where the landscape is rather regular and where more data are available for evaluation. The processing results in 12 land use maps (6 years x 2 training areas) in term of the fractions of the three predominant classes (orchard, non-cultivated areas and annual crop).

  panchromatic image at 2.5m resolution acquired the 23/07/2002 and 10 cloud-free Landsat/ETM7+ and SPOT4/5 images acquired during the 2002-2003 agricultural season. The Landsat/SPOT images were collected between 07/11/2002 and 20/06/2003

(

  MOD13Q1 product) from the 2000-2001 to the 2005-2006 agricultural seasons. These images contain atmospherically corrected reflectances and NDVI at 250m spatial resolution based on the Constrained View Maximum Value Composite algorithm (Huete et al. 2002). They were resampled at 270m (9x30m) spatial resolution using the cubic convolution technique, then subset to the Tensift-Marrakech plain. They were stacked into 6 multi-temporal NDVI images (from September 2000 to August 2001, September 2001 to August 2002 etc). A total of 141 images were processed and visually examinated in order to detect eventual anomalies. Most of images are of good quality excepted three images (18/02/2001, 23/04/2001 and 01/01/2003) that were eliminated

2)

  To discriminate the two remaining classes, NDVI maximum values are calculated from NDVI profiles derived from time series of SPOT and Landsat images. Pixels with a maximum NDVI below 0.4, which contain sparse vegetation, are assigned to the class bare soil. The remaining pixels are supposed to include irrigated areas and are assigned to the class annual crop. The threshold value (0.4) was calibrated to obtain a maximal global accuracy of the classification. [Insert Figure 2 about here] This processing leads to the partition of the area into three classes with about 20% of orchard, 50% of bare soil and 30% of annual crop. The land use map (figure 2) shows that: the bare soil class is predominant outside irrigated areas in western and southern parts of the region; the annual crop class is mainly depicted at the eastern part of the study area within Haouz and Tessaout irrigated sub-regions as well as downstream High-Atlas wadis; orchards are spread over the plain, with the maximal density in the western NFIS irrigated sub-region.

  whole area is used to retrieve the endmembers, the NDVI profiles display rather high value (>0.23) at the beginning and the end of the agricultural season for all years except2005-2006 (figure 8-top). A detailed investigation of the groups of pixels resulting from the k-mean classification shows that the annual crop endmember mainly include pixels of the Tessaout region for the first 5 years(2000)(2001)(2002)(2003)(2004)(2005), while it includes those of the Haouz region for the last year(2005)(2006). The selection of pixels in the Tessaout region results in a significant proportion of trees included in the annual crop class, as discussed in section §4.1 for the 2002-2003 season. This problem disappears when the endmembers are extracted on the reference area. In this case (figure 8-bottom), the seasonality of the annual crop endmember is generally consistent with the phenology of cereal crops (growing season from December to April, and NDVI values below 0.2 outside), but two exceptions can be noticed: • For the 2001-2002 season, the increase of NDVI is delayed and largely reduced (peak of NDVI around 0.4 after April, figure 8-bottom-left). This year is characterised by a shortage in irrigation water after the severe drought that occurs during the 1999-2001 period. In this case, the NDVI pattern matches the 2001-2002 seasonal distribution of rainfall, with most of precipitations recorded in March and April. Therefore, the 2001-2002 annual crop endmember appears not suitable for the retrieval of annual crop fractions. The analysis of other NDVI profiles for this year shows that no profile is representative of the phenology of cereal crop. As an alternative, we replace the 2001-2002 annual crop endmember by the average NDVI profile of the endmember identified on the four

  i.e. quite pure pixels exist at the spatial resolution of 250m, is valid. Secondly, the algorithm is able to identify the most pure areas associated to each of the three classes of interest. The NDVI profiles retained as endmembers match with phenological features of non-cultivated areas (flat profiles with low values on the bare soil class), dense perennial vegetation (flat profiles with rather high values on the orchard class) and cereals (largest NDVI seasonality on the annual crop class). Thirdly, the algorithm is robust since the endmembers generally slightly differ between years. The inter-annual stability of endmembers is particularly true for orchards and bare soils, while the endmembers associated to the annual crop class display a larger inter-annual variability, in relation with changes in water availability (dam irrigation water, seasonal amount and distribution of rainfall).

Figure 2 .

 2 FIGURE CAPTIONSFigure1. Delimitation of the whole study area (in red) and its three main irrigated subregions -NFIS (in yellow), Haouz (in black) and Tessaout (in cyan) -on a Landsat7 image. The white square represents the coverage of Landsat and SPOT4/5 images. Figure 2. Land use map derived from high spatial resolution data on the reference area (2002-2003 season, 30m spatial resolution). Figure 3. 2002-2003 NDVI profiles averaged over the 20 groups of pixels resulting from the k-means classification (gray lines) on the whole area (a) and on the reference area (b). Bold lines with symbols highlight the NDVI endmembers associated to orchard ( ), bare soil ( ), and annual crop ( ). The first day is September the 1 st , 2002.Figure 4. 2002-2003 land use fraction maps derived on each class from the reference land use map (left) and from linear unmixing of MODIS data with the endmembers extracted on the whole area (middle) and on the reference area (right).Figure 5. Estimated versus reference land use fractions (2002-2003 season, 1km spatial resolution): orchard (a), bare soil (b), annual crop (c). Estimates are provided by the linear unmixing model with the endmembers extracted on the whole area (1, at top) and on the reference area (2, at bottom). Black lines are X=Y lines; gray lines are regression lines. Figure 6. Estimated versus observed land use fractions (3 km x 3 km R3 irrigated area, 2002-2003 season, 1km spatial resolution). Estimates are provided by the linear unmixing model with the endmembers extracted on the whole area (a) and on the reference area (b). Black lines are X=Y lines.Figure7. Estimated endmembers through the six-year period of study(2000-2001 to 2005-2006 agricultural seasons) on orchard (a) and bare soil (b) classes. The endmembers are extracted on the whole area (top figures) and on the reference area (bottom figures). On all X-axis, the first day is 1 st September. Figure8. Same as Figure7for the annual crop class: (a)2000-2001 to 2002-2003 seasons, (b) 2003-2004 to 2005-2006 seasons. The endmembers are extracted on the whole area (top figures) and on the reference area (bottom figures). In figure a (bottom), the "4 year average" represents the average of the2000-2001, 2002-2003, 2004-2005 and 2005-2006 annual crop endmembers. In figure b (bottom), "03-04 (rank1)" and "03-04 (rank2)" correspond to the endmembers linked to the 1 st and the 2 nd ranks in the minimisation procedure, respectively. Figure9. Maps of land use fractions derived from linear unmixing of MODIS data for the six years of study(2000-2001 to 2005-2006 agricultural seasons): orchard (left), bare soil (middle) and annual crop (right).

  FIGURE CAPTIONSFigure1. Delimitation of the whole study area (in red) and its three main irrigated subregions -NFIS (in yellow), Haouz (in black) and Tessaout (in cyan) -on a Landsat7 image. The white square represents the coverage of Landsat and SPOT4/5 images. Figure 2. Land use map derived from high spatial resolution data on the reference area (2002-2003 season, 30m spatial resolution). Figure 3. 2002-2003 NDVI profiles averaged over the 20 groups of pixels resulting from the k-means classification (gray lines) on the whole area (a) and on the reference area (b). Bold lines with symbols highlight the NDVI endmembers associated to orchard ( ), bare soil ( ), and annual crop ( ). The first day is September the 1 st , 2002.Figure 4. 2002-2003 land use fraction maps derived on each class from the reference land use map (left) and from linear unmixing of MODIS data with the endmembers extracted on the whole area (middle) and on the reference area (right).Figure 5. Estimated versus reference land use fractions (2002-2003 season, 1km spatial resolution): orchard (a), bare soil (b), annual crop (c). Estimates are provided by the linear unmixing model with the endmembers extracted on the whole area (1, at top) and on the reference area (2, at bottom). Black lines are X=Y lines; gray lines are regression lines. Figure 6. Estimated versus observed land use fractions (3 km x 3 km R3 irrigated area, 2002-2003 season, 1km spatial resolution). Estimates are provided by the linear unmixing model with the endmembers extracted on the whole area (a) and on the reference area (b). Black lines are X=Y lines.Figure7. Estimated endmembers through the six-year period of study(2000-2001 to 2005-2006 agricultural seasons) on orchard (a) and bare soil (b) classes. The endmembers are extracted on the whole area (top figures) and on the reference area (bottom figures). On all X-axis, the first day is 1 st September. Figure8. Same as Figure7for the annual crop class: (a)2000-2001 to 2002-2003 seasons, (b) 2003-2004 to 2005-2006 seasons. The endmembers are extracted on the whole area (top figures) and on the reference area (bottom figures). In figure a (bottom), the "4 year average" represents the average of the2000-2001, 2002-2003, 2004-2005 and 2005-2006 annual crop endmembers. In figure b (bottom), "03-04 (rank1)" and "03-04 (rank2)" correspond to the endmembers linked to the 1 st and the 2 nd ranks in the minimisation procedure, respectively. Figure9. Maps of land use fractions derived from linear unmixing of MODIS data for the six years of study(2000-2001 to 2005-2006 agricultural seasons): orchard (left), bare soil (middle) and annual crop (right).

Figure 5 .Figure 6 .

 56 FIGURE CAPTIONSFigure1. Delimitation of the whole study area (in red) and its three main irrigated subregions -NFIS (in yellow), Haouz (in black) and Tessaout (in cyan) -on a Landsat7 image. The white square represents the coverage of Landsat and SPOT4/5 images. Figure 2. Land use map derived from high spatial resolution data on the reference area (2002-2003 season, 30m spatial resolution). Figure 3. 2002-2003 NDVI profiles averaged over the 20 groups of pixels resulting from the k-means classification (gray lines) on the whole area (a) and on the reference area (b). Bold lines with symbols highlight the NDVI endmembers associated to orchard ( ), bare soil ( ), and annual crop ( ). The first day is September the 1 st , 2002.Figure 4. 2002-2003 land use fraction maps derived on each class from the reference land use map (left) and from linear unmixing of MODIS data with the endmembers extracted on the whole area (middle) and on the reference area (right).Figure 5. Estimated versus reference land use fractions (2002-2003 season, 1km spatial resolution): orchard (a), bare soil (b), annual crop (c). Estimates are provided by the linear unmixing model with the endmembers extracted on the whole area (1, at top) and on the reference area (2, at bottom). Black lines are X=Y lines; gray lines are regression lines. Figure 6. Estimated versus observed land use fractions (3 km x 3 km R3 irrigated area, 2002-2003 season, 1km spatial resolution). Estimates are provided by the linear unmixing model with the endmembers extracted on the whole area (a) and on the reference area (b). Black lines are X=Y lines.Figure7. Estimated endmembers through the six-year period of study(2000-2001 to 2005-2006 agricultural seasons) on orchard (a) and bare soil (b) classes. The endmembers are extracted on the whole area (top figures) and on the reference area (bottom figures). On all X-axis, the first day is 1 st September. Figure8. Same as Figure7for the annual crop class: (a)2000-2001 to 2002-2003 seasons, (b) 2003-2004 to 2005-2006 seasons. The endmembers are extracted on the whole area (top figures) and on the reference area (bottom figures). In figure a (bottom), the "4 year average" represents the average of the2000-2001, 2002-2003, 2004-2005 and 2005-2006 annual crop endmembers. In figure b (bottom), "03-04 (rank1)" and "03-04 (rank2)" correspond to the endmembers linked to the 1 st and the 2 nd ranks in the minimisation procedure, respectively. Figure9. Maps of land use fractions derived from linear unmixing of MODIS data for the six years of study(2000-2001 to 2005-2006 agricultural seasons): orchard (left), bare soil (middle) and annual crop (right).

Figure 9 .

 9 FIGURE CAPTIONSFigure1. Delimitation of the whole study area (in red) and its three main irrigated subregions -NFIS (in yellow), Haouz (in black) and Tessaout (in cyan) -on a Landsat7 image. The white square represents the coverage of Landsat and SPOT4/5 images. Figure 2. Land use map derived from high spatial resolution data on the reference area (2002-2003 season, 30m spatial resolution). Figure 3. 2002-2003 NDVI profiles averaged over the 20 groups of pixels resulting from the k-means classification (gray lines) on the whole area (a) and on the reference area (b). Bold lines with symbols highlight the NDVI endmembers associated to orchard ( ), bare soil ( ), and annual crop ( ). The first day is September the 1 st , 2002.Figure 4. 2002-2003 land use fraction maps derived on each class from the reference land use map (left) and from linear unmixing of MODIS data with the endmembers extracted on the whole area (middle) and on the reference area (right).Figure 5. Estimated versus reference land use fractions (2002-2003 season, 1km spatial resolution): orchard (a), bare soil (b), annual crop (c). Estimates are provided by the linear unmixing model with the endmembers extracted on the whole area (1, at top) and on the reference area (2, at bottom). Black lines are X=Y lines; gray lines are regression lines. Figure 6. Estimated versus observed land use fractions (3 km x 3 km R3 irrigated area, 2002-2003 season, 1km spatial resolution). Estimates are provided by the linear unmixing model with the endmembers extracted on the whole area (a) and on the reference area (b). Black lines are X=Y lines.Figure7. Estimated endmembers through the six-year period of study(2000-2001 to 2005-2006 agricultural seasons) on orchard (a) and bare soil (b) classes. The endmembers are extracted on the whole area (top figures) and on the reference area (bottom figures). On all X-axis, the first day is 1 st September. Figure8. Same as Figure7for the annual crop class: (a)2000-2001 to 2002-2003 seasons, (b) 2003-2004 to 2005-2006 seasons. The endmembers are extracted on the whole area (top figures) and on the reference area (bottom figures). In figure a (bottom), the "4 year average" represents the average of the2000-2001, 2002-2003, 2004-2005 and 2005-2006 annual crop endmembers. In figure b (bottom), "03-04 (rank1)" and "03-04 (rank2)" correspond to the endmembers linked to the 1 st and the 2 nd ranks in the minimisation procedure, respectively. Figure9. Maps of land use fractions derived from linear unmixing of MODIS data for the six years of study(2000-2001 to 2005-2006 agricultural seasons): orchard (left), bare soil (middle) and annual crop (right).

Figure 10 .

 10 Figure 10. Estimated land use fractions averaged over Haouz irrigated sub-region for the six years of study (2000-2001 to 2005-2006 agricultural seasons), together with the annual average of irrigation. Figure 11. Left: map of the relative root mean square error (RRMSE) maps, averaged for the six years of study. Right: histogram associated to the spatial variation of RRMSE.

  Figure 6

  Figure 11 Table 1. Confusion matrix of the 2002-2003 reference land use map (in pixels)

					Field observations
				Orchard on annual understory	Orchard on bare soil	Bare soil	Annual crop	total	Commission error (%)
			Orchard	369	237	0	17	623 2.7
	Output	classification	Bare soil Annual crop	0 162	3 24	279 165	0 499	282 1.1 850 41.3
			total	531	264	444	516	1755
			Omission error (%)	30.5	10.2	37.2	3.3
					Overall Accuracy =77.6%

Table 2 .

 2 Reference land use fractions(%) averaged over the 20 groups of pixels resulting 1 from the k-means classification of 2002-2003 MODIS data; gray colors highlight the 2 composition of the groups selected as endmembers; numbers in bold indicates the 3 highest purity for each of the three classes of interest. 4 5

		Whole area			Reference area
	Group Orchard Bare soil	Annual crop	Orchard Bare soil	Annual crop
	70	3.4	26.6	71.1	3.7	25.2
	1.8	97.4	0.8	1.4	98.2	0.4
	26.6	1.3	72.1	8.8	2.9	88.3
	3.9	91.9	4.3	19.4	16.2	64.4
	57.2	12.5	30.3	54.0	9.2	36.8
	27.3	6.4	66.3	3.1	88.6	8.3
	3.8	73.5	22.7	29.4	57.6	13.0
	27.4	59.0	13.7	40.2	20.9	38.9
	50.1	19.8	30.1	4.3	93.1	2.5
	26.6	24.1	49.3	55.3	25.5	19.3
	55.4	26.5	18.1	3.5	10.7	85.8
	16.9	51.4	31.7	53.5	3.4	43.1
	41.4	41.2	17.4	65.5	12.3	22.2
	43.8	5.4	50.7	2.9	67.0	30.1
	64.0	9.1	26.9	24.8	41.5	33.7
	4.0	31.8	64.2	16.9	76.1	7.0
	15.5	77.2	7.3	33.7	4.4	62.0
	7.3	22.4	70.2	5.6	39.8	54.6
	52.9	1.8	45.4	15.4	65.0	19.5
	6.1	5.7	88.2	41.9	41.9	16.1
	6

Table 4 .

 4 Statistical variables calculated between the estimated and the reference land use fractions(2002)(2003) season, 1km spatial resolution); estimates are provided by the linear unmixing model applied with the endmembers extracted on the whole area (left part) and on the reference area (right part).

			Whole area			Reference area
	Land class	R²	RMSE EFF	Bias	R²	RMSE	EFF	Bias
	Orchard	0.69	0.11	0.65	0.04	0.71	0.10	0.70	-0.01
	Bare soil	0.90	0.12	0.82	-0.07	0.90	0.10	0.88	-0.02
	Annual crop	0.81	0.11	0.79	0.03	0.82	0.10	0.80	0.03

The term 'groups' is used to refer the classes identified by the K-means method in order to avoid confusion with those derived from MODIS data after unmixing.

Acknowledgements

This research was part of SUDMED program. In addition to IRD, this research has been partially supported by the PLEIADeS project of the European Commission (Contract GOCE 037095) and by the « Programme d'Action Intégrée du Comité Mixte Interuniversitaire Franco-Marocain » (PAI Volubilis 06/148). We are indebted to the staff and the directors of ORMVAH and ABHT for their technical assistance in the collection of irrigation and rainfall data. We thank Albert Olioso (UMR CSE, Avignon, France) and Hervé Cardot (UMR CESAER, Dijon, France) for their scientific help during the course of this study. French CNRS-PNTS ('Programme National de Télédétection Spatiale') and CNES-ISIS ('Incitation à l'utilisation Scientifique des Images SPOT') programs are acknowledged for their support.

[Insert Table 5 about

Generalised analysis of endmembers (2000-2006 period)

The algorithm is applied to the 2000-2006 period using successively each MODIS multi-temporal NDVI images. The endmembers expected for the orchard and the bare soil classes are always selected (figure 7a and7b, respectively), the first ones with rather high NDVI values (>0.4) and low seasonal amplitudes (~ 0.2), the second ones with the lowest values (six-year maximum of 0.22).

[Insert Figure 7 about here] For the bare soil and orchard classes, there is a general stability of the endmembers from one year to the other (figure 7). In contrast, the NDVI profiles with the highest amplitudes (annual crop endmembers, figure 8) display a higher variability. When the Figure 10 Table 3. Reference and estimated land use fractions (%) averaged over the reference area (2002)(2003)


Didier Arquès
email: arques@univ-mlv.fr

Jean-François Béraud
email: beraud@univ-mlv.fr

Rooted maps on orientable surfaces, Riccati's equation and continued fractions

Keywords: rooted map, Dyck's equation, Riccati's equation, continued fraction

We present a new approach in the study of rooted maps without regard to genus. We prove the existence of a new type of equation for the generating series of these maps enumerated with respect to edges and vertices. This is Riccati's equation. It seems to be the first time that such a differential equation appears in the enumeration of rooted maps. Solving this equation leads to different closed forms of the studied generating series. The most interesting consequence is a development of this generating function in a very nice continued fraction leading to a new equation generalizing the well known Dyck equation for rooted planar trees. In a second part, we also obtain a differential equation for the generating series of rooted trees regardless of the genus, with respect to edges. This also leads to a continued fraction for the generating series of rooted genus independent trees and to an unexpected relation between both previous generating series of trees and rooted maps.

Introduction

Recent history of maps began with W. T. Tutte in 1962. Three principal approaches characterize this combinatorial domain: 1. The bijective approach was the first used (see [START_REF] Tutte | A census of planar maps[END_REF]). Its principle is to construct a one-toone correspondence between the family of the studied maps and another family of objects for which are known efficient enumeration methods. This approach was developed by R. Cori, B. Vauquelin and D. Arquès (cf. [START_REF] Arquès | Énumération et codage des cartes et hypercartes planaires pointées, calcul formel sur les fractions multicontinues[END_REF][START_REF] Arquès | Les hypercartes planaires sont des arbres très bien étiquetés[END_REF][START_REF] Cori | Planar maps are well labelled trees[END_REF]). 2. The algebraic approach is more recent and was developed by D. M. Jackson and T. I.

Visentin in [START_REF] Jackson | A character theoric approach to embeddings of rooted maps in an orientable surface of given genus[END_REF]. A rooted map can be seen as a pair of permutations acting transitively on the set of half-edges. This approach applies algebraic combinatorics methods of the symmetric group to these permutations to derive a closed form for the exponential generating series of maps regardless of genus. One of the most interesting results in [START_REF] Jackson | A character theoric approach to embeddings of rooted maps in an orientable surface of given genus[END_REF] is the formula: for the generating function Me(y,z) enumerating orientable rooted maps regardless of genus with respect to vertices (y) and edges (z). While theoretically possible to go from these equations to those obtained topologically (see 3.), it seems difficult to do so.

3. The topological approach consists in applying to rooted maps a topological operation interpreted in terms of a functional equation for the generating series of the studied maps. Many articles based on this kind of approach can be found in the literature [1-5, 8-12, 18, 20].

Using a topological approach, W. Tutte obtained by contracting the root edge, the first equation (and an explicit enumeration, [START_REF] Tutte | On the enumeration of planar maps[END_REF]) for the generating series M 0 (v,z) of rooted planar maps with respect to edges (z) and degree of the root vertex (v): ( )

M v z v zM v z vz vM v z M z v 0 2 0 2 0 0 1 1 1 ( , ) ( , ) ( , ) , = + + - -
This equation can be generalized to any genus g and solved (cf. [START_REF] Giorgetti | Énumération des cartes pointées de genre quelconque en fonction des nombres de sommets et de faces, Institut Gaspard Monge[END_REF][START_REF] Bender | The number of rooted maps on an orientable surface[END_REF]). It has also been generalized for the maps on locally orientable surfaces [START_REF] Béraud | Énumération des cartes pointées sur la bouteille de Klein[END_REF][START_REF] Bender | The asymptotic number of rooted maps on a surface[END_REF][START_REF] Bender | The enumeration of maps on the torus and the projective plane[END_REF]. Unfortunately, this set of equations for the generating series M g (v,z) (g≥0) contains new unknown generating series for maps with multiple root vertices and cannot be gathered to derive a unique equation for the generating series M(v,z) of genus independent rooted maps.

The aim of this paper is to obtain such an equation. In the second section of this paper we quickly recall the basic notions of map theory. In Section 3 we derive a new functional equation for the generating series of rooted map without regard to genus, with respect to edges and vertices. This equation is obtained by a new interpretation of the classical operation consisting in removing the root edge of the map. We also obtain a functional equation for the generating series of rooted trees, with respect to edges. In Section 4 we enumerate rooted maps and rooted trees. We express the series as continued fractions, that leads to two unexpected relations, among them a generalization for genus independent maps of the well known Dyck equation for rooted planar trees.

Remark 1. Continued fractions never appeared in the literature to express the generating series of maps. However, in [START_REF] Arquès | Énumération et codage des cartes et hypercartes planaires pointées, calcul formel sur les fractions multicontinues[END_REF][START_REF] Arquès | Les hypercartes planaires sont des arbres très bien étiquetés[END_REF], using the bijective approach, D. Arquès introduced multicontinued fractions to express the generating series of rooted planar maps and hypermaps.

Definitions

For the convenience of the reader we recall quickly some definitions used in the following (for more details about combinatorial maps refer for example to [START_REF] Cori | Un code pour les graphes planaires et ses applications[END_REF]).

Topological map

A topological map C on an orientable surface Σ of 3 is a partition of Σ in three finite sets of cells:

i. The set of the vertices of C, that is a finite set of points ;

ii. The set of the edges of C, that is a finite set of simple open Jordan arcs, disjoint in pairs, whose extremities are vertices ; iii. The set of the faces of C. Each face is homeomorphic to an open disc, and its border is a union of vertices and edges. The genus of the map C is the genus of the surface Σ. A cell is called incident to another cell if one of them is in the border of the other. An isthmus is an edge incident on both sides to the same face.

Combinatorial map

We call an oriented edge of the map a half-edge. We denote by B the set of all half-edges of the map. To each half-edge, is associated in an evident way its initial vertex, its final vertex, and the underlying edge.

α (respectively σ) is the permutation on B associating each half-edge h with its opposite half-edge (resp. the first half-edge met by turning around the initial vertex of h in the positive way chosen on the orientable surface). α is a fixed point free involution. The cycles of α (resp. σ) represent the edges (resp. vertices) of the map.

σ is the permutation σοα on B. The cycles of σ represent the oriented borders of the faces of the map.

In the following, a vertex (resp. edge, face) will be, depending on the context, either the topological object defined at 2.1, or the cycle of σ (resp. α, σ ) according to the previous definitions.

The triplet (B, σ, α) is called the combinatorial definition of the associated topological map C.

Rooted map

A map is called a rooted map if a half-edge h is chosen. The half-edge h is called the root half-edge of the map, and its initial vertex is called the root vertex of the map. We call external face (or root face), the face σ * ( ~) h generated by the root half-edge h . The planar map with only a vertex and no edge is also regarded as rooted, even though it contains no half-edge.

Two rooted maps with the same genus are isomorphic if there exists a homeomorphism of the associated surface, preserving its orientation, mapping the vertices, edges, faces and the root half-edge of the first map respectively on those of the second one.

An isomorphic class of rooted maps of genus g will simply be called a rooted map. Our goal is to enumerate these equivalence classes of rooted maps independently of their genus.

A differential equation for the series of rooted maps

We present here the first topological equation for the generating series of orientable rooted maps regardless of the genus, with respect to vertices and edges.

An equation for maps

We denote by M(y,z) the generating series of orientable rooted maps of any genus (we simply call them rooted maps in the following), where the exponent of y (resp. z) refers to the number of vertices (resp. edges) in the map.

Theorem 1. The generating series M(y,z) of rooted maps is the solution of the following Riccati's differential equation:

[ ]

M y z y zM y z zM y z z z M y z ( , ) ( , ) ( , ) ( , ) = + + + 2 2 2 ∂ ∂ . (1)
Proof. Let C denote a rooted map with root half-edge h . The proof is based on the topological operation of deleting the root half-edge h as introduced by W. Tutte [START_REF] Tutte | On the enumeration of planar maps[END_REF] in the study of planar maps. Four cases are possible. The first two terms (in the right part of the equation) look like the first two terms of Tutte's equation recalled in the introduction, and are obtained in the same way (but generalized to maps enumerated without regard to genus). The last two terms are of a new type.

First case. If C is the rooted planar map reduced to a vertex, the contribution in Eq. ( 1) is y 1 z 0 =y.

Second case. If the edge supporting

h is an isthmus whose deletion disconnects the map into two maps C 1 and C 2 , the contribution of this case in Eq. ( 1) is zM y z ( , ) 2 . The contribution of the first map C 1 is M(y 1 ,z 1 ). The contribution of C 2 is M(y 2 ,z 2 ). During the reconstruction of the map C, all the edges add up (z 1 =z 2 =z), and we do not forget the added root isthmus (multiplication by z). The vertices add up too (y 1 =y 2 =y).

Third case. If

h is a planar loop and is the oriented border of the root face, we remove the loop and root the resulting map C 1 with σ( h ) (Fig. 1). C 1 is a general rooted map to which we add a loop during the inverse step of reconstruction. So the contribution of this case in ( 1) is zM y z ( , ) Fourth case. Here we group together all the cases which have not been studied above. Two sub-cases may be encountered: Subcase 1. The root half-edge h is not an isthmus and is not a planar loop border of the root face (Fig. 2) ; Subcase 2. The root half-edge h is an isthmus whose removal does not disconnect the map C (Fig. 3).

¢¡ £ ¥¤ ¦ § ¦ § © © !
a) Removal of the root half-edge h . We remove the root half-edge h . It gives the map C 1 (cf. points 1 and 2). We root C 1 with both root half-edges h h ~) to be able to reconstruct C from C 1 . We obtain a map C 1 with two root half-edges, possibly equal.

b) Addition of a root half-edge in a double rooted map: the inverse operation. Let C 1 be a map with two root half-edges h 1 and h 2 (with initial vertices v 1 and v 2 respectively). We reconstruct the rooted map C from C 1 by adding an edge from v 1 to v 2 , in the angular sectors a 1 and a 2 defined by h 1 and h 2 (more precisely for k=1,2 we consider the sector between σ -1 ( ) h k and h k ) and we root the half-edge oriented from v 1 to v 2 .

If the angular sectors a 1 and a 2 are in the same face of C 1 , the root half-edge h of C is embedded in this face, splitting it in two new faces of C: h is not an isthmus (Fig. 2). It could be a planar loop, if a 1 =a 2 , but in this case the root face is "outside" the loop. If a 1 and a 2 are not in the same face of C 1 , the root half-edge h of C is embedded on an added handle that collects both faces generated by h 1 and h 2 into one face: h is an isthmus of C whose removal does not disconnect the map (Fig. 3).

To summarize we must choose a half-edge anywhere (eventually equal to the first root halfedge) in any map. To convey this operation on the generating series, we first choose an edge by applying the operator z.∂/∂z to the series M(y,z). There are two possible half-edges associated to the chosen underlying edge, so we multiply by 2, and the contribution of the added edge is z. Thus the contribution of this case in ( 1) is [ ]

2 2 z z M y z ∂ ∂ ( , )
This concludes the proof of Theorem 1.

" # $ " ¢% # % '& ( ') 10 ¥2 # $ # 43 & 0 ¥2 # $ 5 % 5 3
Fig. 2: (Subcase 1) The second root half-edge h 2 of C 1 is in the face generated by h 1 : the root half-edge h of the map C is not an isthmus.

h $ ~h $ 1 C C 1 h 2 a 1 a 2
Fig. 3: (Subcase 2) The second root half-edge h 2 of C 1 is not in the face generated by h 1 : the root half-edge h of the map C is an isthmus whose removal does not disconnect the map.

Corollary 1. The generating series M(z) of rooted maps with respect to the number of edges is solution of the differential equation:

[ ]

M z zM z zM z z z M z ( ) ( ) ( ) ( ) = + + + 1 2 2 2 ∂ ∂ . (2)
Proof. Immediate from (1) by substituting y by 1.

An equation for trees

Let us recall that a tree (of any genus) is a map with only one face. By duality it exists a one-to-one correspondence between rooted trees and rooted monopoles (rooted maps with only one vertex). We denote by T(z) the generating series of rooted trees (or rooted monopoles), where the exponent of z refers to the number of edges in the tree.

Theorem 2. The generating series T(z) of rooted trees-or rooted monopoles-is solution of the following differential equation:

[ ] T(z zT(z z z T(z ) ) ) = + + 1 2 2 ∂ ∂ . (3)
Proof. The previous proof of Theorem 1 applies to monopoles. Let C denote a rooted monopole with root half-edge h . It is easy to see that we can apply the topological operation of removal of the root half-edge in the same way as for Eq. ( 1). Only the second case, with a deconnecting root isthmus, is impossible here.

Remark 2. Th. 2 can also be obtained by noting that the generating function for rooted trees is the coefficient of y 1 in M(y,z), i.e.

T z M y z y

y ( ) ( , ) =       =0 .

Corollary 3.

(1) The number of rooted trees with n edges is:

( )! ! 2 2 n n n . (10)
(2) The number of rooted maps with n edges is:

1 2 1 2 1 1 1 1 0 0 1 1 1 1 n i j j j i k k n k k i n k k i i + = + + + = + > = - ∏ ∑ ∑ + + ( ) ( )! ! ... ,..., . (11)
Proof.

(1) Assume that T is the series t z (2) From ( 8), we express the generating series M(z) with respect to T(z), and we obtain [START_REF] Bender | The enumeration of maps on the torus and the projective plane[END_REF] with little algebra.

Remark 6. The number of rooted trees with n edges is equivalent to the number of permutations σ with one cycle over 2n elements (with two chosen ones for the root edge) divided by 2 n-1

(n-1)!. We can also note that this is the very classical enumeration formula for fixed-point-free involutions, namely the odd factorial.

Here we present enumerating tables. In the first one, we compute the first terms of the generating series of rooted maps M(y,z) with respect to vertices and edges (see Table 1, where we group together maps, vertically with the same number of vertices, and horizontally with the same number of edges). In the second table, we compute the first terms of the series of rooted trees T(z) and maps M(z) (see Table 2). Note that continued fractions (4) and ( 9) are used for the computation (faster than the computation from explicit formulae [START_REF] Bender | The number of rooted maps on an orientable surface[END_REF] 2: The number of rooted trees T(z) and rooted maps M(z), regardless of genus, with respect to edges.

Conclusion

In the first part we proved the first functional equation for the generating series of rooted maps regardless of the genus, with respect to vertices and edges (Th. 1). We also derived an analogous equation for the series of rooted trees (Th. 2). On that occasion we proved that we could obtain a functional equation for the generating series of a family of rooted maps without the help of a variable whose exponent counts the valency-the length of the cycle-of the root face (or the root vertex, by duality).

In the last part, by solving these equations we presented the first expression, in terms of continued fractions, of the generating series of rooted maps regardless of genus, with respect to vertices and edges (Th. 3). We also gave a continued fraction form of the generating series of rooted trees regardless of genus, with respect to edges (Th. 4). Then we deduced the first explicit formula for the number of genus independent rooted maps with a given number of edges (Cor. 3).

We are now working on topological interpretations of the generalized Dyck Equation ( 6), and of the relation between trees and maps [START_REF] Arquès | Relations fonctionnelles et dénombrement des cartes pointées sur le tore[END_REF].

Fig. 1 :

 1 Fig. 1: The root half-edgeh is a planar loop, and the root face is "inside" the loop.

C 1 1

 1 = σ ( ~) (where σ C h 1 ( ~) is the first half-edge in C 1 among σ( ~) h and σ 2 ( ~) h ) and h h 2 = σ(

1 .

 1 This enumeration has been previously obtained through a very different way by T. Walsh and A. Lehman in[START_REF] Walsh | Counting rooted maps by genus. I[END_REF].

Acknowledgments

The authors would like to thank the referees for their useful comments and suggestions.

Enumeration of rooted maps and trees

In this section we solve Eq. ( 1) and (3) to obtain continued fraction forms of the generating series of rooted maps and trees regardless of genus (Th. 3 and4). Then we give explicit formulae for the enumeration of the number of rooted maps and trees with n edges (Cor. 3).

A continued fraction for rooted maps

Equation ( 1) is a Riccati differential equation. We present in Theorem 3 an iterative solution of Equation [START_REF] Béraud | Énumération des cartes pointées sur la bouteille de Klein[END_REF], which leads to a very nice continued fraction form of the generating series of rooted maps. Well known analytical methods may also solve such equations. In Remark 3, we recall the closed forms of such solutions.

Theorem 3. The generating series M(y,z) of rooted maps with respect to the number of vertices and edges is the continued fraction:

Proof. The proof is by recurrence. Let us first define the sequence (M k (y,z)) k ≥ 0 of series by the set of equations: 1) M 0 (y,z) is the desired generating function M(y,z) of rooted maps ; 2) For any integer k in , M k+1 (y,z) is obtained from M k (y,z) by the relation (α k )

Then for every k, M k (y,z) is a solution of the equation

This result is proved by recurrence:

• For k=0, this is Theorem 1 and Equation (1).

• Let k be a positive integer and let us suppose that M k (y,z) is a solution of (β k ):

Now we substitute M k (y,z) by its expression with respect to M k+1 (y,z) from (α k ). After little algebra one obtains the Equation (β k+1 ). The result is then proved by recurrence and the interpretation of this set of equations (α k ) k ≥ 0 gives the continued fraction form of M(y,z). In fact, (5)

This concludes the proof of Theorem 3.

where C is a "constant", depending on y, to be fixed, and ϕ 1 ( , ) y z is the confluent

, a particular solution of the second order differential equation (see [START_REF] Gradshtein | Table of integrals, series and products[END_REF]):

(we add the vertex map) with an algebraic combinatorics point of view (see [START_REF] Jackson | A character theoric approach to embeddings of rooted maps in an orientable surface of given genus[END_REF]). It can be verified that this formula is the solution of Eq. (4).

In Theorem 3 a new relation on maps appears: Corollary 2. The generating series M(y,z) of rooted maps with respect to vertices and edges is the solution of the following generalized Dyck equation:

M y z y zM y z M y z ( , )

Proof. Straightforward since we note that the continued fraction can be rewritten as

A continued fraction for rooted trees

We now give some similar results for rooted trees.

Theorem 4.

(1) Both generating series of rooted trees and rooted maps are linked by the relation:

(2) The generating series T(z) of rooted trees with respect to edges is the continued fraction:

Proof.

(1) From Formula ( 7) and Remark 2 we obtain

(2) By Theorem 3, applied with y=1, we deduce

.

This concludes the proof of Theorem 4. Using well known links between gaussian integrals and continued fraction (see for example [19, eq. 92.15]), we retrieve the result of Theorem 4 after some equivalent transformations of the continued fraction.

Explicit enumeration formulae

From previous expressions, we deduce explicit formulae enumerating rooted maps and trees with a given number of edges.