Keywords: opportunistic data, species distribution map, sampling effort, detection probability

With the internet, a massive amount of information on species abundance can be collected under citizen science programs. However, these data are often difficult to use directly in statistical inference, as their collection is generally opportunistic, and the distribution of the sampling effort is often not known. In this paper, we develop a general statistical framework to combine such "opportunistic data" with data collected using schemes characterized by a known sampling effort. Under some structural assumptions regarding the sampling effort and detectability, our approach allows to estimate the relative abundance of several species in different sites. It can be implemented through a simple generalized linear model. We illustrate the framework with typical bird datasets from the Aquitaine region, south-western France. We show that, under some assumptions, our approach provides estimates that are more precise than the ones obtained from the dataset with a known sampling effort alone. When the opportunistic data are abundant, the gain in precision may be considerable, especially for the rare species. We also show that estimates can be obtained even for species recorded only in the opportunistic scheme. Opportunistic data combined with a relatively small amount of data collected with a known effort may thus provide access to accurate and precise estimates of quantitative changes in relative abundance over space and/or time.

Introduction

How species abundance varies in space and time is a major issue both for basic (biogeography, macroecology) and applied (production of biodiversity state indicators) ecology. Professionals working on biodiversity thus spend considerable resources collecting data that are suitable for estimating this variation [START_REF] Yoccoz | Monitoring of biological diversity in space and time[END_REF]. Most of the scientific literature recommends the implementation of both a statistically valid sampling design and a standardized protocol for collecting such data (e.g. see [START_REF] Williams | Analysis and management of animal populations : modeling, and decision making[END_REF], for a review). Many methods have been developed to estimate species abundance in a defined location, e.g., using mark-recapture methods [START_REF] Seber | The estimation of animal abundance and related parameters[END_REF] or distance sampling approaches [START_REF] Buckland | Distance Sampling : Estimating Abundance of Biological Populations[END_REF]. However, these approaches require an intense sampling effort and are not always practical. Many authors have noted that most frequently, interest will not be in abundance itself, but either in the rate of population change, i.e., the ratio of abundance in the same location at two different time points, or in the relative abundance, i.e., the ratio of abundance at two separate locations [START_REF] Mackenzie | How should detection probability be incorporated into estimates of relative abundance ?[END_REF].

Relative abundance is frequently monitored with the help of simpler schemes. For instance, a set of sites is randomly sampled in the area of interest, and counts of organisms are organized on these sites using a given protocol. At a given location, the resulting count can be used as an index of the true abundance. Indeed, assuming constant detectability over space and time, the average number of animals counted per sampled site is proportional to the true abundance of the species in the area. Log-linear models can be used to represent this average number of animals detected per site as a function of space and/or time (and, possibly, other factors such as the habitat ; see for example [START_REF] Van Strien | Indexing european bird population trends using results of national monitoring schemes : a trial of a new method[END_REF], and thereby, to infer population trends. Thus, such programs have been implemented in many countries to monitor the changes in the abundance of several groups of species, such as birds (e.g., for the French Breeding Bird Survey, see [START_REF] Julliard | Common birds facing global changes : what makes a species at risk ?[END_REF] or butterflies (e.g., for the European Butterfly Monitoring Scheme, see [START_REF] Van Swaay | Butterfly monitoring in europe : methods, applications and perspectives[END_REF]. Estimates of relative abundance have also been commonly used for mapping the spatial distribution of several species (e.g., [START_REF] Gibbons | Mapping avian distributions : the evolution of bird atlases[END_REF].

In addition to such data characterized by a known sampling effort, a large amount of data can also be collected by non-standardized means, with no sampling design and no standardized protocol. In particular, the distribution of the observers and of their sampling effort is often unknown [START_REF] Dickinson | Citizen science as an ecological research tool : challenges and benefits[END_REF]. These so-called "opportunistic data" have always existed, and with the recent development of citizen science programs, we observe a massive increase in the collection of these data on a growing number of species (e.g., [START_REF] Dickinson | Citizen science as an ecological research tool : challenges and benefits[END_REF][START_REF] Hochachka | Dataintensive science applied to broad-scale citizen science[END_REF][START_REF] Dickinson | The current state of citizen science as a tool for ecological research and public engagement[END_REF]. Additionally, as the use of online databases facilitates the exchange and storage of data, such opportunistic data may now include millions of new observations per year that are collected in areas covering hundreds of thousands of square kilometres (e.g., the global biodiversity information facility, including more than 500 million records at the time of writing, see [START_REF] Yesson | How global is the global biodiversity information facility[END_REF].

The temporal and spatial distributions of the observations in such data reflect unknown distributions of both observational efforts and biodiversity. Thus, a report of a high number of individuals of a given species at a given location compared to other locations could be because the focus species is abundant at this location or because numerous observers were present at this location. Using such opportunistic data to estimate changes in the space and time of species abundance is therefore complex, since any modeling approach should include a submodel of the observation process [START_REF] Kéry | Trend estimation in populations with imperfect detection[END_REF][START_REF] Hochachka | Dataintensive science applied to broad-scale citizen science[END_REF] or an attempt to manipulate the data to remove the bias caused by unequal effort (see a discussion in [START_REF] Phillips | Sample selection bias and presence-only distribution models : implications for background and pseudo-absence data[END_REF].

As noted by [START_REF] Mackenzie | Improving inferences in population studies of rare species that are detected imperfectly[END_REF], "In some situations, it may be appropriate to share or borrow information about population parameters for rare species from multiple data sources. The general concept is that by combining the data, where appropriate, more accurate estimates of the parameters may be obtained." In this paper, we propose a general framework which enables to combine data with known observational effort (which we call "standardized" data) with "opportunistic" data with an unknown sampling effort. We focus on multi-species and multi-site data that correspond to the data typically collected in this context.

The purpose of this study is to estimate the relative abundance of the species at different sites (different locations and/or times). We base this estimation on two datasets recording the number of animals detected by observers for each species of a pool of species of interest and each spatial unit of a study area of interest : (i) one "standardized" dataset is collected under a program characterized by a known sampling effort, possibly varying among spatial units, (ii) one "opportunistic" dataset is characterized by a completely unknown sampling effort. We take into account the variation across species of their detectability, yet, as a first step, we assume that the observational bias towards some species are the same across the different sites. We show that, under this assumption, the information concerning both the distribution of the observational effort and the biodiversity can be efficiently retrieved from "opportunistic" data by combining them with standardized data. Moreover, we prove that such a combination returns more accurate estimates than when using the standardized data alone. Our statistical framework allowing this win-win combination can open numerous avenues for application. We used data on French birds, which are typical of existing data, to illustrate the numerous qualities of this framework. Note however that the work presented in this paper is a first step, and that further work will be required to fully account for varying observational bias towards some habitat types across the different sites.

During the reviewing process of this paper, we became aware of an independent and simultaneous work by [START_REF] Fithian | Bias correction in species distribution models : pooling survey and collection data for multiple species[END_REF] which develops similar ideas for combining multi-species and multi-sites data with thinned Poisson models.

Statistical modeling

We want to estimate the relative abundance (relative number of individuals) of I species in J sites. The "sites" j can either refer to different spatial sites, to different times, or to different combinations of sites and times. We suppose that we have access to K datasets indexed by k which gather counts for each species i at each site j. We have in mind a case where some datasets have been collected with some standardized protocol, while some others are of opportunistic nature.

Let X ijk be the count of individuals of the species i by the observers in the site j in the dataset k. In this paper, we propose to model the counts X ijk by

X ijk ∼ Poisson(N ij P ik E jk ), for i = 1, . . . , I, j = 1, . . . , J and k = 0, . . . , K -1, (1)
where N ij is the number of individuals (animals, plants, etc) of a species i at site j, and P ik , E jk are two parameters accounting for the bias induced by the observational processes. The parameter P ik reflects both the detectability of the species i (some species are more conspicuous than others, some are more easily trapped, etc.) and the detection/reporting rate of this species in the dataset k (the attention of the observers may systematically vary among species). The parameter E jk reflects the impact of the varying observational effort (including number and duration of visits, number of traps, etc.) and the varying observational conditions met during the counting sessions. In the next two sections, we explain the origin of our modeling, the hypotheses under which it is valid (see also the discussion Section 5), and we describe precisely the meaning of the two dimensionless parameters P ik and E jk . We refer to the Appendix A for a discussion on the link with models based on thinned Poisson processes. Before moving to these modeling issues, we point out that estimation can be easily carried out in the model (1), since it can be recast into a linear generalized model, see Section 2.4.

Count modeling

The count X ijk of individuals of the species i in the site j for the dataset k is assumed to gather the counts from all visits in the site j. We assume that an individual is only counted once during a single visit, yet it can be counted several times in any dataset due to the possible multiple visits to a site j for a dataset k. In particular, we may have X ijk larger than the number N ij of individuals of the species i in the site j. In the following, we neglect identification errors and false positives.

For an individual a ij of the species i in the site j and a visit v jk in the site j for the dataset k, we define the random variable Z aij v jk which equals 1 if the individual a ij has been seen and recorded during the visit v jk , and 0 otherwise. Assuming that there is no multiple count of an individual during a single visit, the count X ijk is then given by

X ijk = v jk ∈V jk Nij aij =1 Z aij v jk ,
where V jk is the set of all the visits v jk in the site j for the dataset k. In the following, we denote by p aij v jk = P Z aij v jk = 1 the probability for the individual a ij to be seen and recorded during the visit v jk .

If we assume that the random variables {Z aij v jk : a ij = 1, . . . , N ij and v jk ∈ V jk } are independent and that

v jk ∈V jk Nij aij =1 p 2 aij v jk is small compared to v jk ∈V jk Nij aij =1 p aij v jk ,
(which happens when p aij v jk is small), then, according to Le Cam Inequality (Le Cam, 1960), the count X ijk follows approximatively the Poisson distribution

X ijk ∼ Poisson v jk ∈V jk Nij aij =1 p aij v jk = Poisson N ij v jk ∈V jk piv jk , with piv jk = 1 N ij Nij aij =1 p aij v jk .
(2) The parameter piv jk corresponds to the average probability to detect and report during the visit v jk an individual of the species i which has been sampled at random in the site j. We observe that the mean of the Poisson distribution

N ij v jk ∈V jk piv jk = N ij O ijk
is the product of a first term N ij , which is the number of individuals of the species i present in the site j, by a second term O ijk , which is a nuisance term due to the observational process. We underline that the term O ijk can be larger than 1 when the number V jk of visits in the site j for the dataset k is large, since an individual can be counted several times during the V jk visits.

Main modeling assumption

The main hypothesis of our modeling (1) is that the observational parameter O ijk can be decomposed as

O ijk = P ik E jk . (3)
Let us give three examples where such a decomposition holds.

Example 1. (single habitat type) Assume that the ratios piv jk /p i v jk depend only on the species i and i and on the dataset k, so that piv jk /p i v jk = piv j k /p i v j k for all i, i , j, j , v jk , and v j k . This means that the detection/reporting probability piv jk of an individual of the species i during the visit v jk can be decomposed as piv jk = P ik q v jk ,

with P ik the mean detection/reporting probability of the species i during a visit for the dataset k and q v jk depending only on the visit v jk (not on the species i). The parameter q v jk represents the influence of the observational conditions during the visit v jk on the detection/reporting probability. The parameter q v jk is then a very complex function of the observational duration, the visibility conditions (weather conditions during the visit, vegetation met, etc.) and many other variables that affect the detection/reporting probability (number of traps, length of line transects, etc.). When the decomposition (4) holds, we have the decomposition (3) with E jk = v jk ∈V jk q v jk .

The decomposition (4) enforces that the detection/reporting probability piv jk does not depend on interactions between the species i and the visit v jk . This property is quite restrictive and it is not likely to be met when several habitat types are present within a site j. Actually, if two visits v jk and v j k take place in two different habitat types h jk and h j k then the ratios piv jk /p i v jk and piv j k /p i v j k are not likely to be equal for all i and i since some species may be specialized to the habitat type h jk and some others to the habitat type h j k . We can weaken the assumption (4) by allowing interactions iv jk between the species i and the visit v jk as long as they cancel on average on each site j piv jk = P ik q v jk + iv jk , with v jk ∈V jk iv jk 0.

(5)

When (5) holds, we again have the decomposition (3) with E jk = v jk ∈V jk q v jk . Such interactions iv jk can take into account heterogeneous observer attention bias toward the species i, but it does not allow for some systematic bias induced by heterogeneous habitat types. Actually, assume that the site j has two habitats h and h and the site j has only the habitat h . Then if the species i (respectively i ) is specialized to habitat type h (respectively h ) we will have either v j k iv j k < 0 or v jk i v jk < 0. So (5) cannot hold. The next two examples focus on the impact of heterogeneous habitat types.

Example 2. (known habitat types) Assume that for each count, we know in which habitat type it has occurred. Let us introduce the parameter k = (h, k) where h represents the habitat type h and k the dataset. For each dataset k, we can then pool together the counts occurring in the same site j and habitat type h. Let us denote by X ij(h,k) the counts of the species i in the site j, the habitat type h for the dataset k. We assume in the following that each visit occurs in a single habitat type : If not, we can artificially split a single visit in H different habitat types into H different visits, each occurring in a single habitat type.

Our main modeling assumption in this example is that the ratios piv j(h,k) /p i v j(h,k) depend only on the species i and i , the dataset k and the habitat type h. This means that for each i, i , j, j and k = (h, k) we have piv j k /p i v j k = piv

j k /p i v j k
for all visits v j k , v j k in the same dataset and the same habitat type. In this case, the probability piv j(h,k) can be decomposed as

piv j(h,k) = P i(h,k) q v j(h,k) , (6)
with P i(h,k) the mean detection/reporting probability of a typical individual of the species i during a visit in the habitat type h for the dataset k and q v j(h,k) not depending on i. We then have for k = (h, k)

O ij k = P i k E j k , with E j k = v j k ∈V j k q v j k and P i k = P i(h,k) defined by (6).
As above, we can allow some non-systematic heterogeneity by merely assuming that piv j(h,k)

= P i(h,k) q v j(h,k) + iv j(h,k) with v j(h,k) ∈V j(h,k) iv j(h,k) 0.
Example 3. (homogeneous habitat type proportions) We assume again that each visit v jk occurs in a single habitat type h(v jk ) (by artificially splitting non-homogeneous visits). Yet, we assume that this habitat type is not reported in the dataset. As in the second example, we also assume that the ratios piv jk /p i v jk depend only on the species i and i , the dataset k and the habitat type h(v jk ). Hence, the probability piv jk can be decomposed as

piv jk = P ih(v jk )k q v jk , (7)
with P ihk the mean detection/reporting probability of a typical individual of the species i during a visit in the habitat type h for the dataset k and q v jk not depending on i. Writing V jk (h) for the set of the visits v jk in the habitat type h we have

O ijk = H h=1 v jk ∈V jk (h) piv jk = H h=1 P ihk E jhk , with E jhk = v jk ∈V jk (h) q v jk .
The parameters E jhk are likely to depend on h since there can be some observational bias towards some habitat types. If we assume that the observational bias is the same for each site j, which means that E jhk /E j hk does not depend on h, we have the decomposition

E jhk = E jk Q hk , (8)
where Q hk reflects the observational bias towards the habitat type h in the dataset k. When the decompositions ( 7) and ( 8) hold, we have

O ijk = H h=1 P ihk Q hk E jk = P ik E jk , with P ik = H h=1 P ihk Q hk ,
so O ijk fulfills the decomposition (3). Again, as in the two first examples, we can weaken (7) by merely assuming that piv jk = P ih(v jk )k q v jk + iv jk , with

v jk ∈V jk iv jk 0.
Let us explore when the decompositions ( 7) and ( 8) are likely to hold. We first observe that the decomposition (7) will be met as long as we include in the definition of the "habitat type" h(v jk ) all the exogenous variables which induces an interaction between the species i and the visit v jk . The decomposition ( 8) is much more stringent. It requires that, for each dataset k, the observational bias towards some habitat types is the same across the different site j. It may not hold when the proportions on habitat types differ among the different sites. For example, if an habitat h is missing in a site j, then E jhk = 0, so (8) cannot hold if E j hk = 0 for another site j . An example where this property is more likely to be met is when the "sites" j correspond to the same spatial unit observed at different years j. In such a case, we can expect that the observational bias towards some habitat types remains stable years after years. When the observational bias towards some habitat types is not constant across the site, the decomposition (3) is not met in general. This case requires a substantial additional modeling that will be developed elsewhere.

Interpretation. Let us interpret more precisely the parameters P ik and E jk in the decomposition (3). Writing V jk for the number of visits in the site j for the dataset k, we first observe that

1 J J j=1 1 V jk v jk ∈V jk piv jk = 1 J J j=1 O ijk V jk = 1 J J j=1 P ik E jk V jk = P ik Ēk , with Ēk = J -1 J j=1 E jk /V jk .
We can always replace (P ik , E jk ) in the decomposition ( 3) by (P ik , E jk ) = (P ik Ēk , E jk / Ēk ). Applying this renormalization step and dropping the prime (for notational simplicity), we obtain

P ik = 1 J J j=1 1 V jk v jk ∈V jk piv jk , (9)
which means that P ik is the mean detection/reporting probability of a typical individual of the species i during a typical visit for the dataset k.

As explained in the three above examples, the parameter E jk in ( 3) is a complex function of the conditions met during the visits in the site j for the dataset k, including the observational effort. This parameter E jk can be (much) larger than 1 when the number V jk of visits in the site j for the dataset k is very large. We point out that we can have E jk very large even if O ijk is smaller than 1, when the probability P ik of detection/reporting of a typical individual of the species i is very small. In the remaining of the paper, we call observational intensity at the site j in the dataset k the parameter E jk .

Identifiability issues

In the following, we deal with two datasets. A first dataset labeled by k = 0, in which we suppose that the observational intensities E j0 are known up to a constant. Henceforth, we will call this dataset the standardized dataset. We also consider a second dataset labeled by k = 1, characterized by unknown observational intensities E j1 . We will refer to this dataset as the opportunistic dataset.

A single opportunistic dataset is not enough

We consider first the case where we have a single dataset, i.e. K = 1. For notational simplicity, we drop the index k in this paragraph. Our observations X ij then follows a Poisson distribution with intensity λ ij , where λ ij = N ij P i E j . We cannot recover the IJ + I + J parameters N ij , P i , and E j from the IJ intensities λ ij . Yet, if we are only interested by the relative abundances N ij /N ij with respect to a reference site, say j = 1, can we recover the I(J -1) ratios

{N ij /N i1 : j = 2, . . . , J, i = 1, . . . , I} from the IJ parameters λ ij ? Let us write λ ij = N ij P i E j with N ij = N ij P i E 1 , P i = 1 and E j = E j /E 1 .
The parameters N ij differ from the N ij by a multiplicative constant P i E 1 depending only on the species i. Therefore, we have

N ij /N ij = N ij / N ij ,
which means that the parameters N ij give access to the relative abundances N ij /N i1 of the species i. When the dataset has been collected with a known sampling design, the observational intensity in a given site E j is known up to an unknown constant, so that the ratios E j /E j are known and we can recover the N ij (and hence the relative abundances) from λ ij since the P i and E j are known. The situation is different with opportunistic datasets characterized by unknown ratios E j /E j . In this case, the E j are also unknown, so we cannot recover the N ij from the parameters λ ij . Hence, we do not have access to the relative abundance N ij /N i1 . As explained in the next paragraph, we need to combine different datasets.

Combining an opportunistic dataset with a standardized one

Let us now investigate the identifiability issues when we combine a standardized dataset (labeled by k = 0) with an opportunistic one (labeled by k = 1). In this case, we have 2IJ parameters λ ijk = N ij P ik E jk for IJ +2(I +J) parameters N ij , P ik and E jk . For IJ > 2(I +J), which typically holds for large J and I ≥ 3, we have more parameters λ ijk than parameters N ij , P ik and E jk . Nevertheless, as explained in the Appendix B, the model is not identifiable without J +I +1 additional identifiability conditions. As in Section 2.3.1, we introduce some renormalisation N ij , E jk of P ik of N ij , E jk and P ik , which enables us to easily express these identifiability conditions while preserving the identity

N ij E jk P ik = λ ijk = N ij E jk P ik .
In the following, we assume that the ratios {E jk /E j k : j = j } are known for the dataset k = 0 (standardized dataset), but not for the dataset k = 1 (opportunistic one). As above, we define E j0 = E j0 /E 10 (which is known) and P i1 = 1 for all i. We could have set P i0 = 1 instead of P i1 = 1, but the latter choice is more suited for handling species i monitored in the dataset k = 1 but not in the dataset k = 0, as we will show later. We must still set one more constraint. We choose P 10 = 1 for convenience. These I + J + 1 constraints combined with the identity N ij E jk P ik = λ ijk = N ij E jk P ik lead to the change of variables :

N ij = N ij P i1 E 10 P 10 P 11 , E jk = E jk E 10 × P 1k P 10 (10)
P ik = P ik P i1 × P 11 P 1k .
In terms of these new variables, we have the simple statistical model X ijk ∼ Poisson( N ij E jk P ik ) with E j0 = E j0 /E 10 for all j, P i1 = 1 for all i and P 10 = 1. These J + I + 1 quantities are known, and the resulting statistical model is identifiable.

Let us interpret these new quantities. The parameter N ij is proportional to the abundance N ij by an unknown factor P i1 E 10 P 10 /P 11 depending only on the species i. As in Section 2.3.1, these parameters give access to the relative abundance N ij /N i1 = N ij / N i1 of each species i in each site j. The parameters E j1 are equal, up to a constant, to the observational intensity E j1 ; therefore, they provide the relative observational intensities E j1 /E 11 for each site j in the dataset 1. Finally, P i0 is proportional to the ratio P i0 /P i1 by an unknown factor P 11 /P 10 , so we can compare the ratios P i0 /P i1 across the different species. The ratio P i0 /P i1 reflects the systematic difference of attention toward some species among the observers of the two schemes.

In addition, we emphasize that we can consider the case where some species i are not monitored in the dataset 0 but are recorded in the dataset 1. This case can be handled by merely adding the constraints P i0 = P i0 = 0 for the concerned species i.

Estimation via a Generalized Linear Model

We can estimate the parameters N ij , E jk and P ik by the maximum likelihood estimators ( N ij , E jk , P ik ) with the constraints E j0 = E j0 for all j, P i1 = 1 for all i and P 10 = 1. This estimation can be carried out with the help of a generalized linear model. Indeed, with the notations n ij = log( N ij ), e jk = log( E jk ) and p ik = log( P ik ), Model (1) can be recast as a classical generalized linear model from the Poisson family with a log link :

X ijk ∼ Poisson(λ ijk ), with log(λ ijk ) = n ij + e jk + p ik . (11)
Indeed, we only have to define e j0 = log E j0 as a known offset in the model, p i1 = 0 for all i, and fit the resulting model with any statistical package (see Supplementary materials).

Theoretical gain of combining two datasets

It is important to investigate whether the estimates of the relative abundance obtained by combining the dataset 1 with unknown observational intensity ratios E j1 /E j 1 to the dataset 0 with known observational intensity ratio E j0 /E j 0 improves upon the estimates obtained with the single dataset 0. In this section, we investigate this issue analytically. An improvement is expected simply by looking at the balance between the number of observations and the number of free parameters. With the dataset 0, we have IJ observations, and we want to estimate IJ free parameters ; whereas with the two datasets 0 and 1, we have 2IJ observations for IJ + J + I -1 free parameters. The balance between the number of observations and the number of free parameters is better in the second case. Below, we quantify the theoretical improvement more precisely by comparing the variance of the maximum-likelihood estimators in the two cases. Then, we show that dataset combination also allows to estimate relative abundance for species i not monitored in the dataset 0.

Variance reduction

For mathematical simplicity, we assume in the following that the ratios P i0 /P i1 are known for all i. In terms of the normalized variables, this means that the P i0 are known.

When we work with the single dataset 0, we can estimate N ij with the maximum likelihood estimator N 0 ij = X ij0 /( E j0 P i0 ). Let us investigate how the maximum likelihood estimator N ij associated with the model X ijk ∼ Poisson( N ij E jk P ik ) improves upon N 0 ij . We consider the case where the (unknown) observational intensities E j1 in the dataset 1 is much larger than the observational intensities E j0 in the dataset 0. Hence, we consider the asymptotic setting where E j1 goes to infinity. In the Appendix B, we show that the limit variance of

N ij when E j1 → ∞ is given by var( N ij ) Ej1→∞ → var( N 0 ij ) × P i0 N ij l P l0 N lj . (12)
In particular, the variance of the estimate is reduced by a factor

var( N ij ) var( N 0 ij ) Ej1→∞ ≈ P i0 N ij l P l0 N lj ,
when working with the two datasets instead of the sole dataset 0. This factor can be very small for rare species (N ij small), hardly detectable species (P i0 small), or when the number I of monitored species is large.

Let us explain the origin of this variance reduction in the simple case where the ratios P i0 /P i1 are the same for all the species i (which formally corresponds to P i0 = 1 for all i). In this case, we have a closed-form formula for N ij (see Formula (21) in the Appendix B)

N ij = X ij0 + X ij1 l (X ij0 + X lj1 ) × l X lj0 E j0 ,
which reveals the contribution of each dataset to the estimation of the (normalized) relative abundance. Actually, the estimator N ij is the product of two terms, where the first term mainly depends on the opportunistic dataset 1 when the observational intensities E j1 are large, whereas the second term only depends on the dataset 0

N ij Ej1→∞ ≈ X ij1 l X lj1 × l X lj0 E j0 .
Let us interpret these two terms. The first ratio on the right-hand side provides an estimation of the proportion N ij / l N lj of individuals in a site j that belong to a species i. This proportion is estimated by the ratio of the number X ij1 of individuals of the species i observed at site j in the opportunistic dataset to the total number l X lj1 of individuals observed at site j in the same data. When the observational intensities E j1 in the opportunistic dataset 1 is large, the ratio X ij1 / l X lj1 provides a very accurate estimation of the abundance proportion N ij / l N lj , and we have (see Formula ( 22) in the Appendix B)

N ij Ej1→∞ ≈ N ij l N lj × l X lj0 E j0 . ( 13
)
The second term in the right-hand side of ( 13) provides an estimation of the total (normalized) relative abundance l N lj at the site j. This total (normalized) abundance is estimated from the dataset 0 by dividing the total number l X lj0 of individuals counted at the site j in the dataset 0 by the (normalized) observational intensity E j0 . Let us now explain the reduction of variance observed in (12). The formula (13) shows that we estimate N ij by first estimating the total (normalized) relative abundance l N lj with the dataset 0 and then renormalize this estimation with the ratio N ij / l N lj which has been accurately estimated with the dataset 1. The reduction of variance observed in (12) then results from the use of the whole counts l X lj0 at site j in the dataset 0 for estimating N ij instead of the sole counts X ij0 of the species i at site j.

Species not monitored in the scheme characterized by a known sampling observational intensity

As already mentioned, combining the two datasets also allows to estimate N ij for some species i that are not monitored in the dataset 0, but are monitored in the opportunistic dataset 1. This situation formally corresponds to the case where P i0 = 0. For E j1 → ∞, the limit variance of the estimator N ij is (see Formula (25) in the Appendix B)

var( N ij ) Ej1→∞ ∼ N 2 ij l P l0 N lj E j0 .
Because the species i is not monitored in dataset 0, the (normalized) relative abundance N ij cannot be estimated with the sole dataset 0. Thus, there is an obvious improvement to be made by using our estimation scheme that combines the two datasets. To reveal the power of our approach, let us compare the variance var( N ij ) of our relative abundance estimator with the variance of the imaginary estimator N 0,imaginary ij based on an imaginary dataset 0 where the species i would have been monitored with some (imaginary) detection/reporting probability P imaginary i0 . The variance of the maximum likelihood estimator N 0,imaginary ij of N ij with this imaginary dataset 0 would be

N ij /( E j0 P imaginary i0 ) so that var( N ij ) Ej1→∞ ∼ var( N 0,imaginary ij ) × P imaginary i0 N ij l P l0 N lj .
In particular, the estimation provided by N ij can significantly outperform the imaginary estimation we would have obtained with the sole imaginary dataset 0 (where the species i would have been monitored).

Moreover, if we compare the estimator N ij with the imaginary estimator N imaginary ij based on both the imaginary dataset k = 0 and the dataset k = 1, we observe that the ratio of their variance

var( N ij ) var( N imaginary ij ) = P imaginary i0 N ij + l P l0 N lj l P l0 N lj
remains close to one when P imaginary i0 N ij l P l0 N lj . This means that with our estimation scheme, there is not much difference between the estimation based on a dataset collected with known observational intensities where a species i is rare and the estimation based on a dataset collected with known observational intensities where a species i is not monitored. In other words, there is no instability on the estimation of the relative abundance of a species when it is not present in the dataset collected with known observational intensities.

Illustration

Datasets

In this section, we investigate on some datasets the predictive power of our modeling approach. We estimated the relative abundance of 34 bird species in the non-urban habitat of 63 sites in the Aquitaine region (South West of France). We fitted our model with an opportunistic dataset and a dataset collected with known observational intensity. We then assessed the predictive power of our approach with the help of an independent dataset collected with known observational intensity in the same area, hereafter referred as "validation dataset". We therefore illustrate the ability of our approach to provide better predictions of species relative abundance than other approaches based on either of the two datasets alone.

We first describe the opportunistic dataset. We used the recent online database developed by the Ligue de Protection des Oiseaux (LPO, Bird Life representative in France, largest French bird watcher NGO, with regional delegations). This online system was launched successively by the different regional LPO groups, and we acquired data from one of the first groups to start, Aquitaine, South-Western France, with data collection starting in 2007 (www.fauneaquitaine.org). Any citizen who can identify bird species can register on this website and record any bird observation s/he wishes, noting the species, date, and location (to the nearest 500 m). Hundreds of observers thus record hundreds of thousands observations. We typically ignore why these observations were made, e.g., the motivation of the observer, the reason for choosing to report these observations over others, whether they report all the species they have seen at a given place and time, the underlying observational intensity, etc. We selected all such opportunistic records between April and mid-June 2008-2011. For each record, we considered the number of animals detected by the observer. Data were pooled over years, because we will focus here only on spatial variation in relative abundance. Over 115 000 species records detected in a non-urban habitat were considered in this study (see Fig. 1B).

We then describe the dataset collected with known observational intensity, used for the fit of the model. We used the data from the ACT monitoring plan jointly carried out by the French National Game and Wildlife Agency (ONCFS, Office National de la Chasse et de la Faune Sauvage), the national hunter association (FNC, Fédération Nationale des Chasseurs) and the French departemental hunters associations (FDC, Fédérations Départementales des Chasseurs). The main objective of the ACT survey was to monitor the breeding populations of several migratory bird species in France [START_REF] Boutin | Breeding bird monitoring in France : the ACT survey[END_REF] ; ACT stands for Alaudidae, Columbidae, Turdidae, which were the main bird clades of interest for this monitoring, though this program also monitors several Corvidae species (see table 1 for the list of species of interest for our study). Thus, only a fraction of the species recorded by the LPO program was also studied by the ACT survey. The Aquitaine region was discretized into 64 quadrat cells, and in each cell, a 4km long route was randomly sampled in the non-urban habitat of the cell (see Fig. 1(A)). Each route included 5 points separated by exactly 1 km. Each route was traveled twice between April and mid-June, and every point was visited for exactly 10 minutes within 4 hours after sunrise in appropriate weather conditions. Every bird heard or seen was recorded, and for each point and each species, the maximum count among the two visits was retained. The observers were professionals from the technical staff of either the ONCFS or the hunters associations. Note that due to organization constraints, some listening points in a site were not necessarily counted every year. Between 2008 and 2011, over 9 500 birds were counted.

Finally, we describe the validation dataset, used to assess the predictive power of our model. We used the data from the STOC program (Suivi temporel des oiseaux communs), a French breeding bird survey carried out by the French museum of natural history (MNHN, Museum National d'Histoire Naturelle) for the same region and the same years. The STOC survey [START_REF] Jiguet | French citizens monitoring ordinary birds provide tools for conservation and ecological sciences[END_REF] is based on a stratified random sampling, with each volunteer observer being assigned a 2 × 2 km square randomly chosen within 10 km of his house. The observer then homogeneously distributed 10 points within the square. Each point was visited twice between April and mid-June (before and after May 8th, with at least 4 weeks between visits) for exactly 5 minutes within 4 hours after sunrise in appropriate weather conditions (no rain or strong winds). Every bird heard or seen was recorded, and for each point and each species, the maximum count among the two visits was retained. These counts were then summed for a given square, year and species. Between 2008 and 2011, 251 listening points belonging to 29 such squares have been surveyed in non-urban habitat (to allow the comparison with the other datasets, we removed the listening points located in urban habitat), most of them for several years, and over 15241 birds were detected by the observers.

Our aim was to test our model ability to provide a better prediction of the spatial variation in species relative abundance than any model based on either of the two datasets alone. The "sites" of our model were the 63 quadrat cells defined for the ACT survey ; we removed the quadrat cell containing the metropolitan area of Bordeaux (a large town with a population of > 1 million inhabitants), where the sampling process in the opportunistic dataset could not be supposed to be the same as in the other areas (see Fig. 1(A)). We focused on I = 34 bird species (see Table 1). Note that the smaller number of species monitored in the ACT survey allowed to demonstrate the ability of our approach to estimate the relative abundance of species monitored only in opportunistic dataset. For both the ACT survey and the STOC survey, the observational intensity in the site j was measured as the number of points-years sampled in the quadrat cell j during the period 2008-2011. We used the validation STOC dataset to assess the predictive power of our modeling approach. Only 24 sites contained at least one STOC listening point (Fig. 1(C)), so that this assessment was restricted to these sites.

Comparison of the predictive power

Let X ijk be the number of animals of the species i detected in the site j in the dataset k. Let k = a denote the dataset with known observational intensity collected by the ACT survey ; let k = denote the opportunistic dataset collected by the LPO ; finally, let k = s denote the validation dataset collected by the STOC survey. We compared different statistical approaches to estimate the relative abundances of the species in the sites.

Let N m ij be the relative abundance estimated for the species i in the site j with the statistical approach m. We estimated the relative abundance of each species i in each site j with the following approaches : where π k j denotes the number of listening points of the site j sampled in the dataset k, and S j denotes the area of the site j (determined by intersecting each ACT quadrat with the Aquitaine region). For the LPO dataset k = , we had to account for the site-specific unknown intensity. We estimated this intensity with two proxies that are commonly used in such cases. First, we assumed that observational intensity was spatially uniform so that it varied only with quadrat cell area S j (the resulting approach is labeled 1). Another proxy considered that the observational intensity within a site was proportional to the total number of records across the sites (pooled over all species ; the resulting approach is labeled 2).

N a ij = X ija /π a j (14) N s ij = X ijs /π s j (15) N 1 ij = X ij /S j (16) N 2 ij = X ij / i X ij (17)
Finally, we fitted the model described in the previous sections, using the ACT dataset a as the dataset collected with known observational intensity (k = 0), and the LPO dataset as the opportunistic dataset (k = 1). Note that we supposed a quasi-Poisson distribution, to account for moderate overdispersion in our dataset. Thus, we could estimate the value of N +a ij with our approach.

The relative abundance is the absolute abundance multiplied by an unknown constant, and this constant may vary among approaches. Therefore, to allow the comparison between the various approaches, we standardized the relative abundance estimates in the following way :

N m ij = N m ij j N m ij
We want to investigate whether the estimates obtained by our model are closer or not to the true densities than any of the estimates that could be obtained from the individual datasets. We used the value N s ij estimated with the validation STOC dataset as the value of reference. We assessed the predictive power of each approach m by calculating, for each species, the Pearson correlation coefficient between the standardized relative abundance N m ij estimated with the method m and the standardized relative abundance N s ij estimated with the validation dataset. We summarized this power by calculating the median and interquartile range (IQR) of these coefficients over the different species of interest. Although the relative abundance estimates were calculated on the complete dataset, these results were presented by separating the species monitored in the ACT survey, and the species not monitored in this survey. This allowed to evaluate the ability of our approach to estimate the relative abundance of species not monitored in the standardized dataset.

We also investigated the stability of our statistical approach when the standardized dataset is small. We therefore assessed this stability by replacing our big standardized ACT dataset a by a much smaller dataset a . We subsampled the dataset a : for each site, we randomly sampled only one listening point in every site, and we considered the bird counts of only one randomly sampled year for every point. Thus, we artificially divided the observational intensity by 18 in average in this dataset : the complete ACT dataset a stored the bird counts carried out in 1107 listening points-years, whereas the reduced dataset a stored the bird counts carried out in only 63 listening-points-years (one in every site). We also estimated the standardized relative abundance N a ij = X a ij / j X a ij with this reduced dataset. Finally, we estimated the relative abundance N +a ij by combining this reduced standardized dataset with the opportunistic dataset according to our model. We also assessed the predictive power of these two approaches by comparing the estimates with the reference values obtained with the STOC dataset.

The online supplementary material contains the data and the code for the R software (R Core Team, 2013) that will allow the reader to reproduce our calculations.

Results

We fitted our model on the LPO and ACT datasets. There was only a small amount of overdispersion in our data (the coefficient of overdispersion was equal to 1.22) ; the examination of the residuals did not reveal any problematic pattern and the quality of the fit was satisfying. We observe in table 2 that the predictive power was larger for our statistical approach than for all other approaches, whether based on the dataset a or alone.

The predictive power of our statistical approach did not decrease much when model was fit on the smaller standardized dataset a , despite the fact that the observational intensity in this dataset was divided by about 20. In particular, the predictive power of our approach with a reduced dataset remained larger than the predictive power of the other approaches. We observe a strong positive correlation bet-Table 2 -Predictive capabilities of the various possible approaches to estimate the relative abundance of 34 bird species in 63 sites in the Aquitaine region. For each possible estimation approach m, we present the median (calculated over the species) of the Pearson's correlation coefficient between the relative abundance N m ij estimated by the approach m and the relative abundance N s ij estimated by the "reference" STOC approach. In parentheses, we present the interquartile range of this coefficient. These quantities are calculated for the set of species only monitored in the ACT survey and for the set of species not monitored in this survey.

Ratio Species only in ACT Species not monitored in ACT N a+ ween the estimates N +a ij obtained with the full standardized dataset and the estimates N +a ij obtained with the reduced standardized dataset (median Pearson's R = 0.84, IQR = 0.81 -0.90). This illustrates clearly the gain of precision obtained by combining the small standardized dataset with a large amount of opportunistic data, which we demonstrated in section 3.1. The very fine-grained distribution of observations contained in the opportunistic dataset can more efficiently predict site-specific variation in relative abundance than can the standardized dataset.

We investigated the ability of our method to estimate the relative abundance of species not monitored in the ACT survey. Note that the between-site variance of the log relative abundance estimated with our method N a+ ij was larger in average for the species monitored in the ACT survey (median = 2.41, IQR = 1.1 -104) than for the species not monitored in this survey (median = 1.15, IQR = 1.04 -1.34), which resulted in smaller Pearson's coefficient for the latter species (Tab. 2). Our approach performed better than the approaches based on the dataset a or 1 alone. The predictive power of our approach and the approach 2 were similar. Actually, the log observational intensity estimated in a site by our approach for the LPO dataset was strongly correlated with the logarithm of the total number of birds detected in this site (Pearson's R = 0.85), which supports to some extent the common practice of biologists to use the total number of birds detected in a place as a measure of the observational intensity.

Discussion

Overview

We propose a general approach to estimate relative abundances of multiple species on multiple "sites" (corresponding to different times and/or locations) by combining one or several datasets collected according to some standardized protocol with one or several datasets of opportunistic nature. The estimation is performed with the generalized linear model ( 11). This modeling relies on several assumptions, including : (i) the datasets have the same spatiotemporal extent, (ii) the individuals of the monitored species do not cluster into large groups, (iii) either the habitat types are known or the observational bias towards some habitat types are the same across the different sites. In particular, the third hypothesis is quite restrictive and handling cases where it is not met requires significant additional modeling.

We have demonstrated both theoretically (under the assumption that the model is well-specified) and numerically on some datasets, that combining opportunistic data with standardized surveys produces more reliable estimates of the relative abundances than either dataset alone. In particular, we observe an improvement in our example Section 4 even if the above hypothesis (iii) is probably violated. We have also shown that combining opportunistic data with standardized data allows for estimating relative abundance for species which are not monitored in the standardized dataset.

Our approach for combining opportunistic data with survey data is quite general : It requires to be extended in order to overcome the current limitations (see the discussion in the next section) and to be adapted to the specifics of each case study. Yet, we highlight two already promising applications of our framework. First, we emphasize that our framework can be readily used to estimate temporal changes. In such cases, the "sites" j correspond to different times j and E j1 represent the parameters describing the unknown observational intensity at time j for the opportunistic dataset. For temporal variation, biased attention for some habitats in the opportunistic dataset will meet the hypothesis (iii) as long as this biased attention is constant over time. As explained in Section 2.2, such biases will be entirely captured in the estimation of the P ik . For example, the accuracy of bird population trends for France will be considerably improved by the addition of opportunistic data to the current Breeding Bird Surveys.

Another very interesting feature of our framework is its ability to estimate the relative abundance of very rare species, even if these species are not monitored with a scheme with known sampling effort. This has important practical implications. For example, [START_REF] Guisan | Using niche-based models to improve the sampling of rare species[END_REF] noted "in a sample of 550 plots surveyed in a random-stratified way based on the elevation, slope, and aspect of the plot during two consecutive summers in the Swiss Alps (704.2 km 2 ), not one occurrence of the rare and endangered plant species Eryngium alpinum L. was recorded. This was despite the species being easily detectable if present and independent records of the species existing in the area within similar vegetation types." Our framework would be very useful in this context. In particular, if a citizen science program collects opportunistic data on this species along with some other more common species, then the relative abundance of the rare species can be estimated by combining these opportunistic data with standardized surveys monitoring the same common species.

Limitations and extensions

We derived from our analysis Section 2.1 a model based on the Poisson distribution. In practice, we may observe some overdispersion in the data. Causes of overdispersion include clustering of individuals, spatial auto-correlation, identification errors, etc. It is then wise to account for overdispersion in the modeling (see Section 4).

The main assumption in our modeling (1) is that the observational bias O ijk can be decomposed into O ijk = P ik E jk . As explained in Section 2.2, this mainly amounts to assume that the habitat types are known or the observational bias towards some habitat types are the same across the different sites. This assumption will not be met in many cases and we can expect a significant improvement by taking habitat types heterogeneity into account. This issue requires a significant additional modeling and it will be developed elsewhere.

In our estimation framework, we did not take into account any variable affecting the distribution of the relative abundance in the different sites. However, it is well-known that there might be a spatial (if the "sites" are spatial units) or temporal (if the "sites" are time units) autocorrelation in the densities. For example, it is frequent that if the abundance of a given species is high in a given spatial unit, it will also be high in neighboring units. Moreover, spatial units with a similar environmental composition will often be characterized by similar abundances. Explicitly accounting for these patterns in the estimation process could lead to an increased accuracy of the estimation (by reducing the effective number of parameters). This could be done by modeling the relative abundances N ij as a function of environmental variables, or as a function of spatial effects (e.g. using conditional autoregression effects in a hierarchical model, see [START_REF] Banerjee | Hierarchical modeling and analysis for spatial data[END_REF]. Alternatively, it is possible to maximize a regularized log-likelihood, i.e. to maximize for example : log L -

I i=1 J j=1 J m=1 νπ jm ( N ij -N im ) 2
where L is the likelihood of the model, π jm is a measure of "environmental and spatial proximities" between the unit j and the unit m, and ν is a positive parameter that determines the strength of the penalty. The proximities could be of any sort (e.g. taking the value 1 if the two spatial units are neighbours, and 0 otherwise ; inverse Euclidean distances between the units in the space defined by the environmental variables, etc.). This kind of regularization would reduce the number of effective parameters in the model and thereby increase the accuracy of the estimation (for example, see [START_REF] Malbasa | Spatially regularized logistic regression for disease mapping on large moving populations[END_REF].

Our statistical approach relies on the assumption that the measurement errors (identification errors, false positive) were negligible. This is a common assumption in this type of study, although recent studies seem to indicate that (i) even a small number of false positives can lead to biases in estimates [START_REF] Royle | Generalized site occupancy models allowing for false positive and false negative errors[END_REF], and (ii) even highly trained professionals may be subject to such errors (e.g. [START_REF] Mcclintock | Experimental investigation of observation error in anuran call surveys[END_REF]. As a solution to this problem, [START_REF] Miller | Improving occupancy estimation when two types of observational error occur : non-detection and species misidentification[END_REF] proposed to combine data collected using different approaches characterized by different probabilities of identification errors (e.g. hear counts vs. visual counts). This approach has not yet been thoroughly tested though, especially in the context of (relative) abundance estimation. Taking into account measurement errors in our framework, e.g. by integrating the approach of [START_REF] Miller | Improving occupancy estimation when two types of observational error occur : non-detection and species misidentification[END_REF], still requires further study.

The detectability of a given species is not necessarily constant across sites j in the standardized dataset, as documented in the literature [START_REF] Link | Estimation of population trajectories from count data[END_REF][START_REF] Mackenzie | How should detection probability be incorporated into estimates of relative abundance ?[END_REF]. This unaccounted variation of detection probability will result into an unaccounted variation of the observational intensity. Because the knowledge of this intensity plays a crucial role in the fit of the model, such errors may bias the estimates if this variation of detection probability is structured according to some exogenous variables (e.g. habitat types). Many statistical frameworks based on a particular sampling design have been suggested to estimate detectability, such as using mixture models based on repeated counts [START_REF] Royle | N-mixture models for estimating population size from spatially replicated counts[END_REF]. Further work is required to adapt such methods to our proposed framework.

A Link with thinned-Poisson processes

In Section 2.1, we described a first modeling of the count data X ijk leading to our model (1). In this appendix, we explain how the model (1) can also be motivated by another point of view relying on the inhomogeneous point process (IPP, see [START_REF] Cressie | Statistics for Spatial Data Wiley[END_REF]. Indeed, IPPs have recently been shown to be a central approach to model species distribution in ecology. Aarts et al. (2012) have shown the close connections existing between IPPs and resource selection functions, a commonly used approach to model habitat selection by the wildlife [START_REF] Boyce | Relating populations to habitats using resource selection functions[END_REF]. Moreover, IPPs have also been shown to generalize other statistical approaches commonly used to model species distribution, such as the MaxEnt approach [START_REF] Renner | Equivalence of MAXENT and Poisson Point Process Models for Species Distribution Modeling in Ecology[END_REF] or the classical logistic regression [START_REF] Fithian | Finite-sample equivalence in statistical models for presence-only data[END_REF]. We compare the IPP with our approach in this section.

The framework of IPPs suppose that the individuals of the species i are distributed on a domain D according to a Poisson point process with intensity λ i (s). If we assume that the individual at location s is detected and recorded in the dataset k with probability b ik (s), then the individuals of the species i recorded in the dataset k are distributed according to a Poisson point process with intensity λ i (s)b ik (s). The multiplication of λ i (s) with b ik (s) results in a "thinning" of the IPP ; for this reason, the resulting point process is sometimes called thinned-Poisson process (e.g. [START_REF] Fithian | Bias correction in species distribution models : pooling survey and collection data for multiple species[END_REF]. Note that in the context of IPPs, each individual is supposed to be counted at most once in each dataset (undercounting). On the contrary, in Section 2.1, we allowed multiple counts of a single individual during the multiple visits in a site, which makes our development more sensible for studies characterized by a strong observational intensity (which is generally the case of citizen science data).

However, even with this difference, our model (1) can be motivated in the context of IPPs. We can adopt different points of view for estimating relative abundances with this modeling based on IPPs. A first point of view is to introduce a model for the abundance intensities λ i (s) and the probabilities b ik (s) and then estimate these quantities accordingly. Such a point of view has been successfully developed in a simultaneous and independent work by [START_REF] Fithian | Bias correction in species distribution models : pooling survey and collection data for multiple species[END_REF] : They model the abundances intensities by λ i (s) = e αi+β T i x(s) with x(s) some observed environmental variables, the probabilities by b i1 (s) = e γi+δ T z(s) with z(s) some other observed environmental variable and b i0 (s) = 1 at locations where survey data are available and b i0 (s) = 0, else. The abundance intensities are then estimated by λ i (s) = e αi+ β T i x(s) , with α i and β i some penalized maximum likelihood estimators of α i and β i .

An alternative point of view, which corresponds to the point of view developed in this paper, is not to try to infer the intensities λ i (s) for each s, but instead, to work at the scale of a whole site S j ⊂ D and infer the mean abundance Λ ij = Sj λ i (s) ds of the species i on S j . An important feature is that we do not model the abundance intensities λ i (s) and the probabilities b ik (s) in terms of some observed environmental variables, but rather simply assume some structural properties on these functions. In particular, the mean abundance Λ ij in the site j is not assumed to be completely driven by some observed environmental variables.

Let us explain how the model (1) can arise in such a context. Let us denote by d ij (s) = λ i (s)/Λ ij the probability density distribution describing the probability for a given individual of the species i in the site j to be located in s ∈ S j . The number X ijk of individuals the species i counted in the site j in the dataset k is then distributed according to The decomposition does not give a decomposition O ijk = P ik E jk in general. Yet, such a decomposition arises in the three scenarii described below (which are the counterparts of the three examples described in Section 2.2).

Example 1 : sites with homogeneous habitat type. Assume that the species intensity ratios λ i (s)/λ i (s) depend on the species i, i and the site j, but not on the location s ∈ S j . Such a property is likely to be met if the site j has an homogeneous habitat type. In this case, we have λ i (s)/λ i (s) = Λ ij /Λ i j and hence λ i (s) = Λ ij g(s) for all i and s ∈ S j . Then, we have O ijk = P ik E jk with P ik = p ik and E jk = Sj g(s)φ k (s) ds.

Example 2 : observations with known habitat type. In this example, we assume that for each observation we know in which habitat type h(s) it has occurred (in particular, it will be the case if we know the location s of each observation). Exactly as in the Example 2 in Section 2.2, we define k as the couple k = (h, k). Assume that the density distribution d ij (s) depends on the species i only through the habitat h(s) of s : For any i, i and s, s ∈ S j such that h(s) = h(s ) we have d ij (s)/d i j (s) = d ij (s )/d i j (s ). In this case, we have a decomposition d ij (s) = α ih(s) g(s) for all s ∈ S j . Let us denote by S jh = {s ∈ S j : h(s) = h} the portion of the site S j with habitat type h. For any i, j and k = (h, k), the counts X ij(h,k) of individuals of the species i in the habitat h in the site j for the dataset k is distributed according to

X ij(h,k) ∼ Poisson Λ ij P i(h,k) E j(h,k)
with P i(h,k) = α ih p i(h,k) and E j(h,k) = S jh g(s)φ (h,k) (s) ds.

We then have the decomposition O ij k = P i k E j k with k = (h, k). We emphasize that in this case the probability p i(h,k) appearing in the decomposition ( 18) is allowed to depend on the habitat type h (the bias towards some species may differ depending on the habitat type).

Example 3 : homogeneous distribution of habitat types. We do not assume anymore that the habitat type h(s) for each observation is known. We assume again that we have the decomposition d ij (s) = α ih(s) g(s) for all s ∈ S j , hence When the ratios |S jh |/|S j | do not depend on j, the above ratio depends on h and k only, so (19) holds. This case corresponds to sites S j all having a similar distribution of habitat types. This property will be met if the sites S j correspond to the same location at different times j.

O ijk =

B Mathematical proofs B.1 Identifiability conditions

With the notations n ij = log(N ij ), e jk = log(E jk ) and p ik = log(P ik ), the model (1) described in our paper can be recast as a classical generalized linear model X ijk ∼ Poisson(λ ijk ), with log(λ ijk ) = n ij + e jk + p ik .

The kernel of the design matrix associated with this linear regression has a dimension equal to I + J + 1. Therefore, we need I + J + 1 constraints to ensure the identifiability of the model.

B.2 Properties of the estimators

The negative log-likelihood of the parameters ( N ij , E jk , P ik ) is

L = i∈I j∈J k∈{0,1}
N ij E jk P ik -X ijk log( N ij E jk P ik ) + log(X ijk !) where the parameters { E j0 , j ∈ J} and { P i0 , i ∈ I} are known, P 10 = 1 and P i1 = 1 for all i ∈ I.

To keep the mathematical analysis of the maximum likelihood estimators comprehensible, we focus below on the case where the P i0 are known. The maximum likelihood estimators of N ij and E j1 are then the solutions of

N ij = X ij0 + X ij1 P i0 E j0 + E j1 and E j1 = X #j1 N #j , (20)
where X #jk = i X ijk and N #j = i N ij .

We first treat the simplest case where the P i0 are all equal.

Figure 1 -

 1 Figure 1 -The datasets used to illustrate our statistical framework. The location of the Aquitaine region in France is displayed in the insert. (A) distribution of the ACT listening points in the region ; (B) distribution of the LPO records (opportunistic dataset) in the region ; (C) distribution of the STOC listening points in the region. The grey quadrat cells are used as the "sites" in our analysis (they measure ≈ 30 × 20 km). Note that the quadrat cell containing the Bordeaux metropolitan area (indicated by an asterisk in (A)) has been removed from the dataset.

X

  ijk ∼ Poisson (Λ ij O ijk ) with O ijk = Sj d ij (s)b ik (s) ds .Let us describe some scenarii, where the observational bias O ijk can be decomposed as O ijk = P ik E jk , leading to the model (1).In the three examples below, we will assume that the detection/reporting probability b ik (s) can be decomposed inb ik (s) = p ik φ k (s)(18) with φ k (s) not depending on i. This means that the detection/reporting bias b ik (s)/b 1k (s) = p ik /p 1k towards the species i in the dataset k is independent of the location s (in other words the functions b 1k (s), . . . , b Ik (s) are proportional one to the others). When this property is met we have the decomposition O ijk = p ik Sj d ij (s)φ k (s) ds.

Table 1 -

 1 List of the 34 bird species under study. The 13 species monitored only by the ACT survey are indicated by an asterisk. All species were surveyed by the STOC and the LPO program.

	Latin name	species
	Aegithalos caudatus	Long-Tailed Tit
	Alauda arvensis *	Eurasian Skylark
	Alectoris rufa *	Red-Legged Partridge
	Carduelis carduelis	European Goldfinch
	Carduelis chloris	European Greenfinch
	Certhia brachydactyla	Short-Toed Treecreeper
	Columba palumbus *	Common Wood Pigeon
	Coturnix coturnix *	Common Quail
	Cuculus canorus	Common Cuckoo
	Dendrocopos major	Great Spotted Woodpecker
	Erithacus rubecula	European Robin
	Fringilla coelebs	Common Chaffinch
	Garrulus glandarius *	Eurasian Jay
	Hippolais polyglotta	Melodious Warbler
	Lullula arborea *	Woodlark
	Luscinia megarhynchos Common Nightingale
	Milvus migrans	Black Kite
	Cyanistes caeruleus	Eurasian Blue Tit
	Parus major	Great Tit
	Passer domesticus	House Sparrow
	Phasianus colchicus *	Common Pheasant
	Phoenicurus ochruros	Black Redstar
	Phylloscopus collybita	Common Chiffchaff
	Pica pica *	Eurasian Magpie
	Pica viridis	Eurasian Green Woodpecker
	Sitta europaea	Eurasian Nuthatch
	Streptopelia decaocto *	Eurasian Collared Dove
	Streptopelia turtur *	European Turtle Dove
	Sylvia atricapilla	Eurasian Blackcap
	Troglodytes troglodytes	Eurasian Wren
	Turdus merula *	common Blackbird
	Turdus philomelos *	Song Thrush
	Turdus viscivorus *	Mistle Thrush
	Upupa epops	Eurasian Hoopoe

  If we assume in addition thatS jh g(s)φ k (s) ds = Q hk Sj g(s)φ k (s) ds,(19)thenO ijk = P ik E jk with P ik = p ik h α ih Q hk and E jk =Let us investigate when the decomposition (19) can be met. Assume first that φ k (s) = β kh(s) γ k (s) where γ k (s) reflects local fluctuations independent of the habitat type. The function g(s)γ k (s) then represents small scale fluctuations and we expect to haveS g(s)γ k (s) ds ≈ q k |S|,for S large enough. It would be the case for example if g(s)γ k (s) was the outcome of a stationary process.We then have S jh g(s)φ k (s) ds Sj g(s)φ k (s) ds≈ β hk q k |S jh |

	g(s)φ k (s) ds.
	Sj

p ik h α ih S jh g(s)φ k (s) ds. h β hk q k |S jh | .

Acknowledgements

We warmly thank the members of the CiSStats group for stimulating and fruitful discussions on opportunistic data and related statistical issues. We also thank Laurent Couzy and Ondine Filippi-Codaccioni for facilitating access to the LPO-Aquitaine database. Many thanks are also due to the coordinators of the ACT survey at the French wildlife management organization (ONCFS), the French national hunters association (FNC) and the departmental associations (FDC), for allowing us to use this dataset in our study. This work was partially supported by the Fondation Mathématiques Jacques Hadamard through the grant no ANR-10-CAMP-0151-02 in the "Programme des Investissements d'Avenir", by the Labex LMH, by the Mastodon program from CNRS, by the CiSStats program from INRA and by the Chaire de Modélisation Mathématiques et Biodiversité from VEOLIA-Ecole Polytechnique-MNHN.

B.2.1 Case of constant ratios P i0 /P i1

We consider in this paragraph the case where P i0 = P 10 for all i ∈ I. This corresponds to the case where for all the species i, the detection/reporting probability ratios P i0 /P i1 are the same and equal to P 10 /P 11 . We derive from (20)

and inserting this expression in the formula for E j1 we obtain E j1 = E j0 X #j1 /X #j0 . As a consequence, we obtain the closed-form expression for N ij

According to the strong law of large numbers for Poisson processes, we have

and

.

If we estimate N ij with the sole "known-effort" data X ij0 , the maximum likelihood estimator is given by N 0 ij = X ij0 / E j0 and its variance equals var(

B.2.2 Case of arbitrary ratios P i0 /P i1

We no longer assume that the P i0 are all equal. In this case, we have no closed-form formula for N ij but we can compute a first-order expansion of N ij in terms of the inverse of X #j1 .

The first step is to check that N #j is upper-bounded independently of the X ij1 . When P i0 > 0 for all i (which means that the same species are monitored in the datasets 0 and 1), we have from ( 20)

Summing these inequalities we obtain the upper-bound

which does not depend on X ij1 . The case where P i0 = 0 for some i can be treated similarly : splitting apart the indices in I 0 = {i ∈ I : P i0 = 0} and those out of I 0 , we get from ( 20)

This inequality is equivalent to

In the asymptotic E j1 → ∞ we obtain the asymptotic upper-bound

. Now that we have checked that N #j is (asymptotically) upper-bounded independently of the X ij1 , we can write a first-order expansion of the formula ( 20)

Summing these expansions over i ∈ I and simplifying the expression gives

.

Plugging this formula in (24) gives

where the last limit follows again from the law of large numbers for Poisson processes. Computing the asymptotic variance when E j1 → ∞, we find after simplification

As in the previous case, we can compare this variance to the variance of the maximum likelihood estimator N 0 ij = X ij0 /( P i0


Ayoub Tamim
email: ayoubtamim@gmail.com

Hussein Yahia

Khalid Daoudi

Khalid Minaoui

Abderrahman Atillah

Driss Aboutajdine

Mohammed Faouzi Smiej

Mohammed V-Agdal

Detection of Moroccan Coastal Upwelling Fronts in SST Images using the Microcanonical Multiscale Formalism

Nonlinear signal processing using the Microcanonical Multiscale Formalism (MMF) is used to the problem of detecting and extracting the upwelling fronts in coastal region of Morocco using Sea Surface Temperature (SST) satellite images. The algorithm makes use of the Singularity Exponents (SE), computed in a microcanonical framework, to detect and analyze the critical transitions in oceanographic satellite data. The objective of the proposed study is to develop a helpful preprocessor that transforms SST images into clean and simple line drawing of upwelling fronts as an input to a subsequent step in the analysis of SST images of the ocean. The method is validated by an oceanographer and it is shown to be superior to that of an automatic algorithm commonly used to locate edges in satellite oceanographic images. The proposed approach is applied over a collection of 92 SST images, covering the southern Moroccan Atlantic coast of the years 2006 and 2007. The results indicate that the approach is promising and reliable for a wide variety of oceanographic conditions.

Introduction

The sea surface temperature (SST) images obtained from Advanced Very High Resolution Radiometer (AVHRR) sensor have allowed the oceanographers to study and monitor the mesoscale features in the oceans [START_REF] Holyer | Edge detection applied to satellite imagery of the oceans[END_REF]. On SST images, mesoscale structures are delimited by thermal fronts that occur between two water masses of different temperatures; they may be often associated to a convergence zone which is shown to contain loci of high biological activity [START_REF] Bakun | Fronts and eddies as key structures in the habitat of marine fish larvae : Opportunity, adaptative response and competitive advantage[END_REF]. As a consequence, the knowledge of thermal fronts is an important factor for the fish management and aggregation [START_REF] Chassot | Satellite remote sensing for an ecosystem approach to fisheries management[END_REF].

In coastal upwelling systems, the upwelling fronts and filaments are perfect examples of such thermal structures, which are, on top of SST data, evidenced by the presence of low temperature values and high chlorophyll concentration in onshore water.

A variety of techniques and algorithms have addressed the problem of automated detection of upwelling fronts in SST images. The straightforward and the most popular methods used to distinguish these structures, are the conventional horizontal gradients computed from classical convolution operators [START_REF] Belkin | An algorithm for oceanic front detection in chlorophyll and SST satellite imagery[END_REF], based on the fact that upwelling fronts are characterized by a very complex horizontal alteration in the upper ocean, generating strong thermal gradients. In this context, [START_REF] Oram | Gradient-based edge detection and feature classification of sea-surface images of the southern california bight[END_REF] developed an edge detection algorithm specifically for detecting features in sea surface satellite images. The latter is started with the application of gradient based edge detection using the Canny edge detector with an automatic selection of the upper and lower thresholds values from the cumulative global histogram calculated for each image. Then series of steps (edge following, edge validation and edge merging) are used to accurately and objectively detect the upwelling fronts in satellite images of the ocean. Another class of methods deals directly with the gray-level histogram based on the idea that edges in SST images are usually regions of separation between two regions of constant temperatures [START_REF] Nieto | Mesoscale frontal structures in the canary upwelling system: New front and filament detection algorithms applied to spatial and temporal patterns[END_REF]. Therefore an edge is considered to be present if the algorithm detects two significantly different populations.

Another method that combines the histogram based separation and structure growing process has also been applied to satellite images [START_REF] Marcello | Automatic tool for the precise detection of upwelling and filaments in remote sensing imagery[END_REF], where the bimodality in an SST image is used to detect the first approximation of the upwelling fringe separating the cold waters near the coast and warmer offshore waters followed by the structure growing process in order to attain the fine detail of the upwelling structures.

The foremost shortcoming encountered in the preceding approaches, is that they have been constructed for use other than the oceanographic context, and consequently they often work poorly when used to detect the mesoscale structures in ocean satellite images [START_REF] Cayula | Edge detection algorithm for SST images[END_REF]. The conventional methods used in image processing often rely on simple hypotheses, compared to the physical of the phenomena encountered. In particular, an SST image is a measure of a thermodynamical variable, the temperature of the ocean's microscopic upper layer in a turbulent flow (the surface of the ocean); the existence of turbulence implies the formation of extremely complex structures that conventional (linear) approaches cannot apprehend, mainly because of the absence of localization in the Fourier transform. As a matter of fact, the intermittency and the existence of multiscale organization present in Fully Developed Turbulence (FDT) [START_REF] Frisch | Turbulence : The legacy of A. N. Kolmogorov[END_REF] make the upwelling region a complex system relevant to nonlinear signal processing.

In [START_REF] Isern-Fontanet | Microcanonical multifractal formalism: Application to the estimation of ocean surface velocities[END_REF], the authors use the Microcanonical Multifractal Formalism (MMF) and the Most Singular Manifold (MSM) in order to reconstruct the stream function from a specific reconstruction formula defined by a universal kernel propagator applied to satellite data. The resulting stream function is then used to derive mesoscale ocean dynamics. While the authors make use of the singularity exponents and the MSM to determine mesoscale ocean dynamics, the work proposed in our study focuses on the thermal fronts detection from the knowledge of the singularity exponents. Additionally, the thermal fronts detected are validated by a professional oceanographer using a database of 92 SST images of Moroccan atlantic coast. The present paper is part of a larger effort to test the feasibility and the relevance of the Microcanonical Multifractal Formalism to the Moroccan coastal upwelling. In particular, we show that edge pixels computed in such a complex systems framework provides superior validated results in the case of oceanographic data.

The flowchart of the algorithm is presented in Fig. 1, and the remainder of the paper is organized as follows. Section 2 presents the database and the geographic area of study used in this work. Section 3 describes the method used for the automatic detection and extraction of upwelling fronts in SST images. Section 4 reports the results conducted to evaluate the proposed approach. The conclusion and the future work are reported in Section 5.

Study Area and Data

In this study, we make use of a database of 92 AVHRR SST images covering the southern part of Moroccan Atlantic coast (20 • 0 -27 • 60 N and 13 • 35 -23 • 30 W) during the years 2006 and 2007. This region, belonging to the Canary upwelling system, is characterized by strong and persistent thermal structures [START_REF] Atillah | Produits opérationnels d'océanographie spatiale pour le suivi et l'analyse du phénomène d'upwelling Marocain[END_REF], including upwelling fronts which separate the cold coastal waters and warmer offshore oceanic waters, and also include the filaments which, represent the extensions of upwelling fronts in open ocean.

The SST images used in this study are received and processed at the Royal Centre of Remote Sensing (CRTS) of Morocco, including geometric, radiometric and atmospheric corrections and also include the cloud and land masks. The cloud overlay is generated using the multispectral radiance measurements [START_REF] Stowe | Scientific basis and initial evaluation of the CLAVR-1 global clear/cloud classification algorithm for the advanced very high resolution radiometer[END_REF] tested and evaluated over the AVHRR satellite images. The size of each SST image is 770 × 990 pixels, with spatial resolution of 1.1×1.1 km. The temperatures, in Celsius degrees, are estimated using the non-linear SST algorithm (NLSST) [START_REF] Walton | The development and operational application of nonlinear algorithms for the measurement of sea surface temperatures with the noaa polar-orbiting environmental satellites[END_REF], which greatly reduces any dependence of the gamma parameter upon the satellite measurements and has been shown to be more accurate than the linear algorithm, Multi-Channel SST (MCSST) [START_REF] Mcclain | Comparative performance of avhrr-based multichannel sea surface temperatures[END_REF].

Due to the a priori knowledge of our region of interest and the fact that the coastal upwelling is largely conditioned by the width and the direction of the continental shelf [START_REF] Atillah | Produits opérationnels d'océanographie spatiale pour le suivi et l'analyse du phénomène d'upwelling Marocain[END_REF], we explore the whole database and retain for our analysis only SST images that have at least 50% of valid information (no clouds or missing pixels) in the region of interest (ROI) between the coast and 200 km offshore [START_REF] Nieto | Mesoscale frontal structures in the canary upwelling system: New front and filament detection algorithms applied to spatial and temporal patterns[END_REF]. Additionally, this ROI is chosen to be parallel to the coastline direction (NE -S W), in order to efficiently explore the effect of the width of Moroccan continental shelf [START_REF] Atillah | Produits opérationnels d'océanographie spatiale pour le suivi et l'analyse du phénomène d'upwelling Marocain[END_REF][START_REF] Metzner | Interpretation of dynamical processes in the upwelling area off northwest africa using simulated and remotely sensed data[END_REF]. The dates of the 92 SST images selected from the years 2006 and 2007 are detailed in Table 1.

Fig. 2 shows two SST images selected from our benchmark of 92 SST images, illustrating two scenarios of upwelling variability: well-defined upwelling in terms of relatively strong thermal fronts and with continuity along the coast (Fig. 2(a)), and upwelling with weak thermal gradients and where the cloud contamination affects the continuity of upwelling along the coast (Fig. 2(b)). These two images will be used throughout this paper to illustrate the validity of our proposed algorithm.

Fronts detection methodology

Fronts detection

Transitions of the graylevel values inside complex natural images, such that SST images, are related to the existence of a hierarchy of multiscale structures characteristic of FDT [START_REF] Frisch | Turbulence : The legacy of A. N. Kolmogorov[END_REF]. Singularities in SST images corresponding to sets of pixels with sharp and strong variations in intensities values, unlock the multiscale organization (Turiel et al., 2008b). In this work we use a novel formalism called the Microcanonical Multifractal Formalism (MMF) for the analysis of these singularities, based on the strength of variations around any given pixel, which yields to a precise description of upwelling fronts in SST images.

For an SST image I defined by its graylevel function I( x) and its norm of gradient ∇I ( x) at each pixel x, we introduce a measure µ defined by its density with respect to the Lebesgue measure, which corresponds to summing the norm of the gradient:

dµ( x) = d x ∇I ( x), (1)
So that for a ball B r ( x) of radius r centered around the pixel x, the measure µ is given by:

µ(B r ( x)) = B r ( x) d x ∇I ( x) (2)
We work with a normalized luminance intensities :

I( x) = I( x) -I 0 , I 0 = Ω I( x)d x, (3)
(Ω : image's domain) such that, the luminance function I( x) has zero average. The measure µ defined above is a multifractal measure, and for any pixel x over a set of fine scales r the following equality holds :

µ(B r ( x)) = α( x)r d+h( x) + o(r d+h( x) ) (r → 0) ( 4
)
where d = 2 indicates the dimension of space (d = 2 for an image), and the α( x) is a signal-dependent amplitude prefactor. The exponent h( x) defined in equation ( 4) is the singularity exponent (SE) at x [START_REF] Turiel | Numerical methods for the estimation of multifractal singularity spectra on sampled data : A comparative study[END_REF]; it quantifies the strength of the SST image at the pixel x in a statistically coherent way across the scales.

It is worth mentioning that the MMF is not relying on "thermalized" values associated to large ensembles averages and ergodic hypotheses; it has been shown that such "canonical" approaches show difficulties in the processing of empirical data [START_REF] Turiel | Numerical methods for the estimation of multifractal singularity spectra on sampled data : A comparative study[END_REF]; instead, in a microcanonical analogy, the behaviour of the measure is evaluated locally and geometrically at every given pixel location, without relying on thermalized averages [START_REF] Maji | Edges, transitions and criticality[END_REF], it has been shown that such an approach is computationally more efficient and, most importantly, it gives access to the geometric localization of transition fronts.

A central concern in MMF is the proper determination, at high numerical precision, of the singularity exponents h( x). In this context, an approach presented in (Turiel et al., 2008a) [START_REF] Pont | Singularity analysis of digital signals through the evaluation of their Unpredictable Point Manifold[END_REF] provides numerically stable computation of the SE at each pixel x, based on the wavelet projection of the measure:

h( x) = log(τ ψ µ( x, r 0 ))/ τ ψ µ(., r 0 ) log r 0 + o( 1 log r 0 ) (5)
where τ ψ µ( x, r 0 ) is a projection of measure µ evaluated at point x and scale r 0 with wavelet ψ (not necessarily admissible, as positive mother wavelets display interesting computational properties) and τ ψ µ(., r 0 ) is the average value of the wavelet projection over the whole signal, with r 0 chosen to reduce the relative amplitude of the correction term o( 1 log r 0 ). The scale r 0 is chosen such that a discretized image of size N × M corresponds to size 1, so that one chooses r 0 = 1 √ N×M . It is expected that the set of pixels having the strongest transitions, i.e. the least values of singularity exponents, also known as the Most Singular Manifold (MSM) contains most of a turbulent's signal information and provides indications in the natural image about the pixels where sharp and sudden local variations take place. Moreover, from a complex systems approach, the MSM is made of most unpredictable points. The MSM is defined as:

F ∞ = { x : h( x) = h ∞ = min(h( x))} (6)
Visual inspection of the MSM pixels in natural images reveals the presence of the contours and edges. Fig. 3(a) shows the set of MSM pixels, corresponding to the smallest values of singularities exponents calculated over the original image in Fig. 2(a).

Fronts extraction

The presence of cloud contamination in the input image (pixels with gray color on the left side of Fig. 2 serious effects on the results of edge pixel detection. Indeed, the MSM points detected in Fig. 3(a) include the MSM pixels of upwelling fronts separating cold and warmer water, and the MSM pixels of cloud fronts separating the cloud and water masses.

To avoid the misclassification of the cloud frontal boundaries as upwelling fronts, we used the cloud mask cited in [START_REF] Stowe | Scientific basis and initial evaluation of the CLAVR-1 global clear/cloud classification algorithm for the advanced very high resolution radiometer[END_REF], such that for each MSM pixel in the Fig. 3(a), if the 8-neighboring pixels contain at least one cloud pixel, the MSM of the given pixel is discarded. Fig. 3(b) shows the MSM map after removing all the MSM pixels adjacent to clouds in Fig. 3(a).

In order to examine each individual front in Fig. 3(b), it is necessary to first extract and link independent MSM pixels into continuous contours. The algorithm implemented for this purpose scans the MSM edge map (Fig. 3(b)) until a MSM pixel is found. That pixel (MSM(i,j)) becomes the 1st pixel of the nth contour. The algorithm then adds the adjacent pixels by examining the 8-connected neighbors of pixel MSM(i,j). The contour is terminated when none of the 8-connected MSM pixels neighbors are found.

In an effort to discard short, possibly insignificant fronts, a minimum length criteria is imposed on each extracted front: each edge not exceeding a minimum segment length (10 pixels) is discarded [START_REF] Oram | Gradient-based edge detection and feature classification of sea-surface images of the southern california bight[END_REF]. Fig. 3(c) shows the upwelling fronts that are linked and extracted from the MSM pixels map in Fig. 3(b).

Fig. 4(a) and (b) illustrate the results of upwelling fronts detected with the proposed algorithm, in which the fronts are automatically overlaid on the top of the respectively original SST images in Fig. 2(a) and (b).

Results analysis on the whole database

In order to validate our algorithm in terms of over-detection or missing-detection of upwelling fronts, we compared our results to the algorithm defined in [START_REF] Oram | Gradient-based edge detection and feature classification of sea-surface images of the southern california bight[END_REF]. The latter extended the standard Canny edge detector [START_REF] Canny | A computational approach to edge detection[END_REF] by using an automatic selection of the hysteresis thresholds values from the global cumulative histogram in each image. The main advantages of the [START_REF] Oram | Gradient-based edge detection and feature classification of sea-surface images of the southern california bight[END_REF] algorithm are: 1) the use of the Canny edge detector, which is a well-known and wellestablished gradient-based edge detector; and 2) the fact that the Canny edge detector is used as comparative method with a recent developed algorithm [START_REF] Nieto | Mesoscale frontal structures in the canary upwelling system: New front and filament detection algorithms applied to spatial and temporal patterns[END_REF]. Fig. 5 shows the results of the [START_REF] Oram | Gradient-based edge detection and feature classification of sea-surface images of the southern california bight[END_REF] method (black contours) overlaid on the original images in Fig. 2.

It is necessary to mention that, the in-situ measurements, e.g., buoys data, can give us a sounder support to the results obtained by the proposed approach. However, these measures are very complex to establish and would be extremely challenging for a use in coastal zone context, because of lack of these measures for all the SST satellite images used in this study and for every specific region of the study area. In this sense, it is a well-known fact that validation in the case of ocean data is often and truly performed by assessing the results by professional oceanographers [START_REF] Nascimento | Automated computational delimitation of SST upwelling areas using fuzzy clustering[END_REF][START_REF] Tamim | A simple and efficient approach for coarse segmentation of moroccan coastal upwelling[END_REF][START_REF] Sousa | Identification of upwelling areas on sea surface temperature images using fuzzy clustering[END_REF], which has been done in the case of this work.

Throughout this evaluation, we used 5 grades: "Bad", "Poor", "Acceptable", "Good" and "Excellent". "Bad" is assigned when the upwelling fronts are not well detected by the algorithms (MSM or [START_REF] Oram | Gradient-based edge detection and feature classification of sea-surface images of the southern california bight[END_REF]), and "Excellent" is assigned when all the thermal fronts are well identified.

Regarding the results of the evaluation presented in Fig. 6, we can conclude that the MSM method outperforms the [START_REF] Oram | Gradient-based edge detection and feature classification of sea-surface images of the southern california bight[END_REF] method which is based on the Canny edge detector with automatic selection of the hysteresis thresholds values. Indeed, for the grade "Bad", the value of 2% is reached by our proposed method while for the [START_REF] Oram | Gradient-based edge detection and feature classification of sea-surface images of the southern california bight[END_REF] the value of 9% is attained. For the grade "Good" a value of 48% is achieved using the MSM method and 40% with [START_REF] Oram | Gradient-based edge detection and feature classification of sea-surface images of the southern california bight[END_REF]. Overall, 54% and 42% are reached by the two grades "Good" and "Excellent" together respectively for MSM and [START_REF] Oram | Gradient-based edge detection and feature classification of sea-surface images of the southern california bight[END_REF] methods.

After this evaluation of the two fronts detection methods, over this representative database, we can conclude that the proposed algorithm using the MSM method has provided a satisfactory results and shown to be superior to that of the method used in [START_REF] Oram | Gradient-based edge detection and feature classification of sea-surface images of the southern california bight[END_REF]. In fact, the MSM method detects almost all the thermal fronts in SST images and has the advantage of remove all the detected contours in the cloud contamination regions. However, the method used in [START_REF] Oram | Gradient-based edge detection and feature classification of sea-surface images of the southern california bight[END_REF] does not detect all the real fronts (Fig. 6) and has the limitation to misclassify the cloud contamination regions as frontal regions (Fig. 5(a) and (b)), namely in the area near shoreline and islands. More importantly, the fronts detection quality obtained by this method suggests that it can be used to obtain a first approximation of the upwelling fronts and serve as a basis for a further and subsequent steps in the analysis of the SST images of the ocean.

Conclusion and future prospects

The singularity exponents computed in a microcanonical formulation show to be an effective approach for the identification and detection of upwelling fronts and filaments in the Moroccan coastal upwelling. The implemented methodology starts by detecting the mesoscale structures in the SST images using the MSM method, followed by an extraction process in order to link the individual edge pixels. The work reported in this paper presents the first and the basis step in the analysis of the SST images. The next step will include the automated classifying of the thermal structures, delineated previously, in order to indicate their association with some various oceanographic structures, such as fronts, filaments and eddies characterizing the variability of the Moroccan surface ocean. The developed algorithms have been successfully applied over a database of 92 AVHRR daily SST images, covering the southern part of Moroccan Atlantic coast.

Fig. 1 :

 1 Fig. 1: Simplified block diagram of the automatic algorithm used to detect the thermal upwelling structures in SST images.

Fig. 2 :Fig. 3 :

 23 Fig. 2: SST images obtained respectively on (a) 2007-01-02 and (b) 2006-01-06, showing different upwelling scenarios.

  Fig. 4: (a) Fronts shown in Fig. 3(c) overlaid on the original image in Fig. 2(a); (b) Fronts detected in the same way (MSM method) overlaid on the original SST image in Fig. 2(b).

Fig. 5 :

 5 Fig. 5: Results of applying the algorithm Oram et al. (2008) respectively to the original images in Fig. 2(a) and (b).

Fig. 6 :

 6 Fig. 6: Qualitative evaluation made by the oceanographer of fronts detection results produced by MSM and Oram et al. (2008) methods applied over the 92 SST images.

Table 1 :

 1 Dates of the 92 SST images used in this study, covering the years 2006 and 2007.

	Years Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
	2006 06, 14,	03, 16,	05, 13,	02, 07,	04, 15,	01, 07,	07, 19,	11, 19,	10, 20,	01, 10,	02, 14,	06, 15,
	24, 28	25	15, 25	22, 25	22	17, 23,	27, 31	21, 31	25	19, 24	24, 26	26, 29
						30
	2007 02, 14,	09, 14,	02, 12,	04, 14,	08, 11,	08, 16,	07, 15,	03, 05,	04, 11,	03, 12,	01, 09,	02, 03,
	20, 30	22	14, 22	22, 27	19, 30	19, 28	22	16, 25	14, 24	21, 29	17	17, 26

Acknowledgement

This work is funded by the French-Moroccan research program Volubilis (MA/11/256) and the project n • MPI 12/2010.