Integration of remote sensing data into national statistical office sampling designs for agriculture

Ambrosio, Luis; Iglesias, Luis; Marín, Carmen; Deffense, Nicolas

doi:10.3233/SJI-220116

Integration of remote sensing data into national statistical office sampling designs for agriculture

Article type: Research Article

Authors: Ambrosio, Luis^{a; *} | Iglesias, Luis^a | Marín, Carmen^a | Deffense, Nicolas^b

Affiliations: [a] Universidad Politécnica de Madrid, Madrid, Spain | [b] Université Catholique de Louvain, Ottignies-Louvain-la-Neuve, Belgium

Correspondence: [*] Corresponding author: Luis Ambrosio, Universidad Politécnica de Madrid, Madrid, Spain. E-mail: [email protected].

Keywords: Agricultural statistics, point sampling, categorical data, multinomial logit models, cost efficiency

DOI: 10.3233/SJI-220116

Journal: Statistical Journal of the IAOS, vol. 39, no. 2, pp. 473-489, 2023

Published: 30 May 2023

Get PDF

Abstract

The integration of remote sensing data in agricultural statistics is a research topic with a long history. The research focus is on using statistical models to link ground and remote sensing data such that the resulting estimators are design-consistent.

A design-consistent estimator assisted by linear models is well established in the literature. However, it requires enough geographic information about the boundaries of agricultural parcels to develop a simple sample of areas. Many countries use complex samples based on non-georeferenced list frames of households or farms and reduce to point data the georeferenced information required for linking ground and remote sensing data.

Data on crop acreage observed at a point are necessarily categorical because a point is dimensionless. Little work has been done on the integration of categorical ground data within complex list samples using remote sensing data. Our focus was on using multinomial logit models for this integration. Special attention was paid to evaluate the cost efficiency of remote sensing data.

1.Introduction

Recently, the European Space Agency (ESA) launched the Sentinels for Agricultural Statistics(Sen4Stat) project. The aim was to facilitate the use of remote sensing information at the National Statistical Offices (NSOs) supporting agricultural statistics. Other international programs to improve agricultural statistics were begun around the same time. An additional Sen4Stat target is synergy with other initiatives, including the Living Standards Measurement Study – Integrated Survey on Agriculture (LSMS-ISA), which was launched in 2008 by the World Bank, and the Global Strategy to Improve Agricultural and Rural Statistics (GSARS), which was endorsed in 2010 by the United Nations Statistical Commission and implemented by the Food and Agriculture Organization (FAO) of the United Nations. An FAO output is the Agricultural Integrated Survey (AGRIS) method [1], together with the 50 × 2030 Initiative, which builds on work of the LSMS-ISA and AGRIS programs [2]. Four of five partner countries in Sen4Stat were also partners in other initiatives, i.e., Ecuador and Senegal of AGRIS, Malawi of LSMS-ISA, and Tanzania of GSARS. Accordingly, through existing collaborations and initiatives, the FAO coordinated the provision of ground data from respective agencies in Ecuador, Senegal and Tanzania, while ground data collection from Malawi was coordinated by the World Bank.

Remote sensing (RS) furnishes a very useful auxiliary data source to improve the cost efficiency of design-based inference, and the integration of ground and remote sensing data is a research topic with a long history [3, 4, 5, 6, 7]. Statistical models are used for this integration, in such a way that the estimator assisted by these models is design-consistent [8].

The focus has been on linear models, assuming that the sample design is simple. In countries where geographic information of agricultural parcel boundaries is available, it is possible to use an individual or cluster of parcels as the sampling unit and to design a simple sample of areas (polygons) using area frames. A linear model is suitable for integrating ground and remote sensing data on crop acreage or yield because these are continuous variables. However, many countries use complex samples based on non-georeferenced list frames (of households or farms) and reduce to point data the georeferenced information required for linking ground and remote sensing data. Crop data observed at a point are necessarily categorical because a point is dimensionless, and multinomial logit is an appropriate statistical model to integrate remote sensing and categorical ground data [9, 10].

Little work has been done on the integration of categorical ground data in complex list samples of households or farms with remote sensing data. In this paper, we focus on using multinomial logit models for this integration. We give details on statistical models to be used for crop acreage, yield and production estimation, as well as on the sample design used by five Sen4Stat partner countries that are representative of the large types of sample designs proposed in the literature. These are list, area, and point sampling [11].

Considering the literature on the applications of remote sensing to agricultural statistics, as well as expectations of the NSO partner countries for these applications, we elaborated four use cases: 1. Cost efficiency; 2. Granularity; 3. Multi-seasonal estimation; 4. Optimizing the sample design. A prototype for integrating ground and remote sensing data in each use case was developed for a large set of sampling designs, including those currently used by the partner countries. We include results of the application of these prototypes to data provided by the NSOs, for cost efficiency and granularity. Finally, we outline prototypes for multi-seasonal estimation and optimizing the sample design.

2.Methods

The approach to elaborate official statistics is well established in the literature. It is design-based in the sense that the inference is based on the sampling distribution generated by the probabilistic sampling scheme designed to select the sample [12, 13]. This sampling scheme associates with each population unit, i=1,2,…,N, a known and positive inclusion probability πi, which play a key role in design-based inference. Usually, the sample size is sufficiently large to attain design consistency (the design-based distribution of the estimator is tightly concentrated around the true population value) at national and regional levels, and design-consistent estimators are used. Estimate accuracy is also based on the estimator’s design-based distribution.

Our focus was on improving the cost-efficiency of the currently used sampling designs for elaborating official statistics on crop acreage, using RS as auxiliary data. The approach to using auxiliary data in design-based inference is well stablished in the literature [8], i.e., the auxiliary extra-sample information is integrated in the sampling design using statistical models, without loss of design-consistency. In minor administrative areas (municipalities) the sample size is small or nil, so that design consistency is meaningless. Thus, a model-based approach is used to get estimates, whose accuracy strongly depends on the auxiliary information reliability.

Model suitability for ground and remote sensing data integration depends on the nature of ground data. For continuous ground data a linear regression model is used, but for categorical and counts ground data, a generalized linear model is more suitable [9, 10]. Both types of models allow remote sensing data of any character, continuous or discrete. The efficiency of RS as auxiliary data for crop acreage has already been extensively demonstrated in the scientific literature, but only for continuous ground data and with the assumption that the sample design is simple. To date, little work has been done on the integration of categorical ground data in complex list samples of households or farms with remote sensing data. In this paper, we focus on using multinomial logit models for this integration.

2.1Sample designs

The design of the sampling scheme for use in sample selection is key in design-based inference because it associates with each population unit i=1,2,…,N a known and positive inclusion probability πi, which is used to define the estimators’ statistical distribution. This distribution is used to evaluate the estimators’ characteristics, including design consistency and accuracy.

A sampling frame is required to select a probabilistic sample and we considered the two main types of sampling frames used in practice, area frame and list frame. The former is based on maps and/or satellite images and the latter on population censuses. Sampling units are either polygons or points in area frames and households or farms in list frames. Usually, the former is georeferenced, but not the latter [11].

The area frame has important advantages over list frames, including versatility, reduction of non-sampling errors such as coverage errors, and longevity. However, it has some inconvenience (sensitivity to outliers and sub representation of rare or minor activities) [14]. To overcome these inconveniences, a dual frame composed of an area and list frame is recommended [15]. In Section 4, we give details on the dual frame used in Senegal to integrate ground and RS data.

Using area frames, it is possible to define sampling units (polygons) of the same size and simple samples with equal inclusion probabilities. In contrast, using list frames, the sample selection scheme entails two or more stages, the sampling units are usually of distinct size, and the probabilistic scheme assigns unequal (size-proportional) inclusion probabilities. Therefore, the sample is not simple but complex.

The integration of ground and RS data will be illustrated for the set of sample designs most often used in practice for collecting ground data [11]. With area frames, systematic sampling with one or more random starts and stratified random sampling are most often used as the sample schemes when geographic information on parcel boundaries is available. If the parcel boundaries are unknown, the sampling unit is usually a point, and we will illustrate how to integrate a sample of points selected from a list frame with RS data.

Agricultural data are frequently collected using national household surveys, for which the sample is selected from a list frame (a population census) in two or more stages, with unequal inclusion probabilities. The sampling unit is a household and geographic information at parcel level is rarely available. For RS data integration, georeferenced information is required at a minimum. To have at least one georeferenced point by agricultural parcel within the household sample, an additional sampling stage is used to select the parcel centroid.

2.1.1Area sampling

Three of the five Sen4Stat partner countries use area samples, namely, Spain, Ecuador, and Tanzania. In Spain, the area frame is the national topographic map. The sampling unit is a square segment of side 700 m (49 hectares). The sample design is non-stratified systematic with three random starts. The sample of segments for crop acreage estimation is selected in a single stage with inclusion probability πi=n/N, where n and N are the number of sampling units in the sample and population, respectively.

In Ecuador, the sampling frame is a land-use map, stratified into four strata according to the percentage of cultivated area. The sampling unit is a segment with geometric boundaries, whose size changes with strata: 9 hectares (ha) in stratum 1a, 36 ha in stratum 1b, 144 ha in stratum 2, and 576 ha in stratum 3. The inclusion probability is equal within strata, πh⁢i=nh/Nh, where nh and Nh are the number of sampling units in the sample and in the population of stratum h, respectively.

In Tanzania the sampling frame used to be a list frame based on an agricultural census. Recently, in the framework of the GSARS program, the Government of Tanzania in collaboration with the FAO, United States Department of Agriculture (USDA) and African Development Bank, designed an area sample [16]. The sampling frame is a map of the country stratified by land use and the sampling unit is a point. Point sampling has a long history [17, 18]. In agricultural statistics, it is used for selecting a sample of farms when a sampling frame is not available; only an accurate map is required [11]. The sample of farms is used for data on a large set of items, generally by direct interview of farmers.

Some of the aforementioned quantities, including crop acreage and yield and natural resources, can be observed directly on the ground, bypassing the need to contact farmers for data collection. France is an example of countries using point sampling for crop acreage estimation, using data collected directly on the ground [19]. Furthermore, the USA is an example of countries using point sampling for a national inventory of natural resources [20, 21].

The sampling design is basically the same in the two aforementioned countries, i.e., a two-stage sampling in which the primary sampling unit (PSU) is a square segment and the secondary sampling unit (SSU) is a point. In the USA National Resources Inventory survey, the segments are stratified, and the side of the typical segment is 800 meters (1/2 mile), i.e., 64 hectares, with some smaller (200 meters) and some larger (1600 meters). The primary sampling rates are generally 2% to 6% of the land area. There are about 300,000 PSUs and 844,000 points, and the second-stage sample size is between 1 and 3 points per PSU.

In the French TER-UTI, the sampling design is non-stratified, one-stage, and systematic. The national territory is divided into square blocks with side 12000 meters. Each block is divided into four square sections of side 6000 meters. In each section, a systematic sample of 36 points is selected, aligned in both row and column directions, separated by a distance of 300 meters, forming a 6 × 6 point grid.

This sampling scheme can be seen as two-stage, so the inclusion probability is of the form πi1⁢i2=πi1⁢πi2|i1, where πi1 and πi2|i1 are the first and second stage inclusion probabilities, respectively.

2.1.2List sampling

Two of the five Sen4Stat partner countries use list sampling, Malawi and Senegal. In Malawi, the ground data were collected in the framework of the LSMS-ISA program. The IHS5 sampling frame is based on the 2018 Malawi Population and Housing Census. It is stratified into rural and urban strata and each of the 27 districts were considered a separate sub-stratum, part of the main rural stratum. A two-stage sampling scheme was used to select the sample from each district. The PSU is the Enumeration Area (EA) and the first-stage sample is systematic without replacement, with probabilities proportional to the number of households in the PSUs. The SSU is a household and a systematic sample of 16 households was used to choose the second-stage sample, using equal probabilities. To integrate with RS data, a georeferenced point is selected in each parcel of the sample [22].

In Senegal, the sampling frame is based on the 2013 population census [23]. The sampling unit is a household and the sample was designed within the AGRIS program framework. This is a two-stage sampling scheme in which the PSU is the EA, and the first-stage sample is selected with replacement and probabilities proportional to the number of households in the PSUs. The SSU is a farm (households with agricultural activities) and the second-stage sample is selected without replacement and with equal probabilities among households with agricultural activities.

To integrate the data in the list sample of farms with RS data, a point is selected in each parcel within the two-stage sample. The inclusion probability of this point is πi1⁢i2⁢i3⁢i4=πi1⁢πi2|i1, where πi1 and πi2|i1 are the first and second stage inclusion probabilities, respectively. This is so because if the SSU i1⁢i2 is included in the sample, then the set of parcels i1⁢i2⁢i3 are also included, πi3|i1i2=1; the georeferenced point i1⁢i2⁢i3⁢i4 is not selected at random but is the parcel centroid πi4|i1i2i3=1.

2.2Statistical models

Statistical models are used to integrate the auxiliary extra-sample information into the sampling design without loss of design-consistency [8]. We used linear models if field data were continuous and multinomial models if those data were categorical [9, 24].

2.2.1Linear models

We consider a population of N units, say agricultural parcels. If geographic information on parcel boundaries is available, it is possible to generate the required remote sensing (RS) data at unit level. We then consider a linear model yi=xi⁢β+εi relating ground data in the it⁢h unit yi, with the values of a set of RS variables associated with this same unit, xi=𝑚𝑖𝑠𝑠𝑖𝑛𝑔r⁢o⁢w1⩽l⩽L⁢(xl⁢i). The parameter vector β=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽l⩽L⁢(βl) is unknown but can be estimated. εi represents unobservable zero-mean random perturbations that, conditionally on xi, account for the variability of yi about its expected value E⁢yi=xi⁢β.

National-level estimators

The survey variable total in the population is yN=∑i=1Nyi and the mean is y¯N=1N⁢∑i=1Nyi. If xi contains 1 (say x1⁢i= 1), then yN=xN⁢BN. Here, xN=∑i=1Nxi, and BN=(XNT⁢XN)-1⁢XNT⁢yN designates the vector of population regression coefficients BN=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽l⩽L⁢(Bl), where XN=𝑚𝑖𝑠𝑠𝑖𝑛𝑔r⁢o⁢w1⩽i⩽N⁢𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽l⩽L⁢(xl⁢i) and yN=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽i⩽N⁢(yi). Since xi is known for every population unit i=1,2,…,N, then xN and XN are known. However, yN is unknown and so is BN.

To estimate BN, we use the data collected in the sample of size n units, selected from the population with inclusion probabilities {πi;i=1,2,…,N}. The design-based estimator B^π=(∑i=1nxiT⁢xiπi)-1⁢∑i=1nxiT⁢yiπi is design-consistent for BN and as a result, the synthetic or projective estimator y^N=xN⁢B^π is design-consistent for yN. This same estimator can be written in the form y^N=yπ+(xN-xπ)⁢B^π, called a generalized regression estimator, where yπ=∑i=1nyiπi and xπ=∑i=1nxiπi. The projective estimator can be expressed as weighted sum y^N=∑i=1nwi⁢yi with weight wi=giπi, where gi=1+(∑i=1Nxi-∑i=1nxiπi)⁢(∑i=1nxiT⁢xiπi)-1⁢xiT [12].

The error of y^N as an estimator of yN is y^N-yN=xN⁢(B^π-BN) and the sampling variance is V(y^N-yN)=xNV(B^π-BN)xNT, where V⁢(B^π-BN)=(∑i=1NxiT⁢xi)-1⁢V⁢∑i=1nxiT⁢εi⁢Nπi⁢(∑i=1NxiT⁢xi)-1 and εi⁢N=yi-xi⁢BN. A designconsistent estimator of the sampling variance is V^⁢(y^N-yN)=V⁢(∑i=1nε^i⁢Nπi), where V(.) is the design-based variance and ε^i⁢N=yi-xi⁢B^π.

The asymptotic distribution of y^N is [V(y^N-yN)]-1/2(y^N-yN)→N(0,1) and can be used for constructing confidence intervals for yN.

A design-consistent estimator for the mean is y¯^N=1N⁢y^N. Its sampling variance is V⁢(y¯^N-y¯N)=1N2⁢V⁢(y^N-yN) and can be estimated using V^(y¯^N-y¯N)=1N2V^(y^N-yN). Its asymptotic distribution is the same as that of y^N.

Domain-level estimators

A domain is a part of the population, for instance, a region R or province within a country. Let NR=∑i=1NIi be the domain size, where Ii=1 if unit i belongs to the domain, and Ii= 0 otherwise. The survey variable total in the domain is yNR=∑i=1NRyi and the mean is y¯NR=1NR⁢∑i=1NRyi. To estimate yNR, we use the sample s of size n selected from the population with inclusion probabilities {πi;i=1,2,…,N} and the estimator y^NR=xNR⁢B^π+NRN^R⁢∑i=1nRyi-xi⁢B^ππi. Here, nR=∑i=1nIi is the number of units in the sample belonging to the study domain and N^R=∑i=1nIiπi is an estimator of NR. If NR is unknown, then we use the estimator y^NR=xNR⁢B^π+∑i=1nRyi-xi⁢B^ππi. The sampling variance is V⁢(y^NR-yNR)=NR2⁢V⁢1N^R⁢(∑i=1nRεi⁢Nπi). A design-based estimator of V⁢(y^NR-yNR) is V^⁢(y^NR-yNR)=NR2⁢V⁢1N^R⁢(∑i=1nRε^i⁢Nπi).

Model-based small area estimation

The sample was designed to achieve the required estimate accuracy at national level. However, reliable estimates in small administrative areas such as a municipality are also required without increasing sample size n. In a minor administrative area, n will always be smaller than at national level and, as a result, so will be the estimators’ accuracy. A small area is a part of the population of which, owing to the small sample size, the design-based estimator is not sufficiently accurate for most uses and the design consistency requirement is meaningless.

For estimates at municipality level, we used a model-based estimator to “borrow strength” from related small areas to obtain accurate estimates for a given small area. Let {(yd⁢i,xd⁢i);i=1,2,…,nd;d=1,2,…,D} be the available dataset, where yd⁢i represents ground data in unit i of the small area d, xd⁢i is the vector of RS data, nd is sample size in the small area d, and D is the number of small areas, so that n=∑d=1Dnd.

We consider the unit-level linear mixed model yd⁢i=xd⁢i⁢β+ud+εd⁢i, where (ud,εd⁢i) are zero-mean independent random variables of variances (σu2,σe2). We assume the same regression parameters β and same variance components (σu2,σe2), through small areas.

The model-based estimator of the total survey variable in a small area yNd=∑i=1Ndyd⁢i is y^Nd=Nd[(1-γ^d)⁢x¯Nd⁢β^+γ^d⁢(y¯nd+(x¯Nd-x¯nd)⁢β^)], where β^=(XT⁢V^-1⁢X)-1⁢XT⁢V^-1⁢y is the generalized least-square estimator of β, with y=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽d⩽D⁢𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽i⩽nd(yd⁢i), X=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽d⩽D⁢𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽i⩽nd⁢(xd⁢i), V^-1=𝑚𝑖𝑠𝑠𝑖𝑛𝑔d⁢i⁢a⁢g1⩽d⩽D⁢(V^d-1) and V^d-1=1σ^e2⁢Ind-γ^dnd⁢σ^e2⁢1nd⁢1ndT, where γ^d=σ^u2σ^u2+σ^e2/nd. Here, σ^e2=eT⁢en-D-1 and σ^u2=uT⁢u-(n-2)⁢σ^e2n∗, with n∗=∑d=1Dnd⁢[1-nd⁢x¯d⁢(XT⁢X)-1⁢x¯dT], are unbiased estimators of the variance components, where eT⁢e is the sum of squared residuals in the model fitted by ordinary least squares, taking as fixed the small-area effect ud. uT⁢u is the sum of squared residuals in the model fitted by ordinary least squares, with ud= 0. x¯Nd and x¯nd are the population and sample means of the vector xd⁢i, respectively.

An unbiased estimator of the total mean-square error estimator 𝑀𝑆𝐸⁢(y^Nd) is [25]:

M^⁢𝑆𝐸⁢(y^Nd)=Nd2⁢(1-ndNd)2

[h1⁢d⁢(σ^u2,σ^e2)+h2⁢d⁢(σ^u2,σ^e2)+2⁢h3⁢d⁢(σ^u2,σ^e2)],

where:

h1⁢d⁢(σ^u2,σ^e2)=γ^d⁢(σ^e2nd)+(1-ndNd)2

(Nd-nd)Nd2

h2⁢d⁢(σ^v2,σ^e2)=σ^e2⁢(X¯Nd-nd-γ^d⁢x¯d)⁢A-1

(X¯Nd-nd-γ^d⁢x¯d)T

h3⁢d(σ^u2,σ^e2)=1nd21(σ^u2+σ^e2nd)3[(σ^e2)2Vσ^u2+

(σ^u2)2Vσ^e2-2σ^e2σ^u2𝐶𝑜𝑣(σ^e2,σ^u2)],

where

A=∑d=1D(∑i=1ndxd⁢iT⁢xd⁢i-γ^d⁢nd⁢x¯ndT⁢x¯nd)

V(σ^u2)=2n∗2[1n-D-1(D-1)(n-2)

(σ^e2)2+2n∗σ^e2σ^u2+n∗∗(σ^u2)2]

V⁢(σ^e2)=2⁢(σ^e2)2n-D-1

𝐶𝑜𝑣⁢(σ^e2,σ^u2)=-1n∗⁢(D-1)⁢V⁢σ^e2

Here, n∗∗=∑d=1Dnd2⁢(1-nd⁢x¯nd⁢A1-1⁢x¯ndT)+t⁢r⁢(A1-1⁢∑d=1Dnd2⁢x¯ndT⁢x¯nd)2. Because A1=∑d=1D∑i=1ndxd⁢iT⁢xd⁢i, n∗∗ may be simplified to n∗∗=∑d=1Dnd2⁢[1-x¯nd⁢(XT⁢X)-1⁢x¯ndT]=n∗-n+∑d=1Dnd2.

2.2.2Multinomial logit models

If geographic information on parcel boundaries is unavailable, then it is not possible to generate RS data at parcel level. In this case, we use a point as the unit of observation. Because a point is dimensionless, observed crop data are necessarily categorical. For integrating categorical ground and RS data, multinomial logit is a suitable model [9, 10, 24].

We consider the survey vector yi=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K⁢(yi⁢k), where yi⁢k= 1 if crop k covers pixel i and yi⁢k= 0 otherwise. K is the total number of crops so that the constraint ∑k=1Kyi⁢k=1 holds. We focus on high-resolution satellite images and the population size N is the number of pixels in the study area A=a⁢N, where a is the area of the piece of land represented by a pixel.

National-level estimators

The population total of the survey vector is yN=∑i=1Nyi=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K⁢(∑i=1Nyi⁢k)=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K⁢(yN⁢k), where yN⁢k is the number of pixels covered by crop k. It is assumed that a pixel is covered by only one crop (or that a class of mixed pixels is included), so that the area covered by k is Ak=a⁢yN⁢k and the population area is A=∑k=1KAk. We want to estimate the total yN or mean y¯N=1N⁢∑i=1Nyi=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K⁢(1N⁢∑i=1Nyi⁢k)=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K⁢(yN⁢kN), where yN⁢kN is the proportion of pixels covered by k.

We assume that the survey vector yi=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K⁢(yi⁢k) follows a multinomial distribution 𝑀𝑁⁢(1,μi), where μi=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K⁢(μi⁢k) and μi⁢k is the probability that unit i is of category k, i.e., the probability of yi⁢k= 1 and yi⁢k′= 0; ∀k′≠k, with the constraint ∑k=1Kμi⁢k= 1. We estimate yN using an estimator of μi⁢k based on a sample of pixels {(yi,xi⁢s);i=1,2,…,n} selected with inclusion probability πi, where xi=𝑚𝑖𝑠𝑠𝑖𝑛𝑔r⁢o⁢w1⩽l⩽L⁢(xi⁢l) is the vector of RS data.

To link μi⁢k with RS data, we consider a multinomial logit model μi⁢k=exp⁡(xi⁢βk)∑k=1Kexp⁡(xi⁢βk), where βk=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽l⩽L⁢(βk⁢l) is an unknown parameter vector. To estimate the probability μi⁢k that the cover of pixel i is crop k=1,2,…,K, it is sufficient to acquire estimates of βk for k=1,2,…,K-1 [9]. The probability of category K is μi⁢K=1-∑k=1K-1μi⁢k and can be estimated using the estimates B^π⁢k for k=1,2,…,K-1.

The design-based parameter estimator B^π=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K-1(B^π⁢k)=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K-1⁢𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽l⩽L⁢(B^π⁢k⁢l) is design-consistent for the maximum likelihood estimator of β=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K-1(βk) based on the population data {(yi,xi);i=1,2,…,N}, considered a simple random sample of the multinomial model. This can be found iteratively using B^π(m+1)=B^π(m)+[∑i=1nH⁢(yi,β)πi|β=B^π(m)]-1∑i=1nb⁢(yi,β)πi|β=B^π(m) [13], where b⁢(yi,β)=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K-1𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽l⩽L⁢[(yi⁢k-μi⁢k)⁢xi⁢l] and H(yi,β)=[𝑚𝑖𝑠𝑠𝑖𝑛𝑔col1⩽k⩽K-1𝑚𝑖𝑠𝑠𝑖𝑛𝑔row1⩽k′⩽K-1(δk⁢k′μi⁢k-μi⁢kμi⁢k′)]⊗xiTxi, with δk⁢k′=1 if k=k′ and δk⁢k′= 0 otherwise A design-consistent estimator of yN=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K⁢(yN⁢k) is y^N=∑i=1Nμ^i+NN^p⁢∑i=1nyi-μ^iπi=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K-1⁢(y^N⁢k), where μ^i=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K-1⁢(μ^i⁢k), μ^i⁢k=exp⁡(xi⁢B^π⁢k)1+∑k=1K-1exp⁡(xi⁢B^π⁢k) for k=1,2,…,K-1, μ^i⁢K=1-∑k=1K-1μ^i⁢k=11+∑k=1K-1exp⁡(xi⁢B^π⁢k), N^p=∑i=1n1πi, and y^N⁢k=∑i=1Nμ^k⁢i+NN^p⁢∑i=1nyk⁢i-μ^k⁢iπi is the estimator of the total number of pixels covered by crop k for k=1,2,…,K-1; for category K it is y^N⁢K=N-∑k=1K-1y^N⁢k. The sampling covariance matrix of y^N is given approximately by V⁢y^N=N2⁢V⁢1N^p⁢∑i=1nyi-μiπi=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K-1⁢𝑚𝑖𝑠𝑠𝑖𝑛𝑔r⁢o⁢w1⩽k′⩽K-1⁢(N2⁢𝐶𝑜𝑣⁢(y^r⁢k,y^r⁢k′)), where y^r⁢k=1N^p⁢∑i=1nyk⁢i-μk⁢iπi is a function of N^p⁢k=∑i=1nyk⁢iπi, the estimators of the total number of sampling units of category k for k=1,2,…,K (because N^p=∑k=1KN^p⁢k), and of the estimator of the total residuals r^k⁢N=∑i=1nyk⁢i-μ^k⁢iπi of category k. We focus on the sampling variance V⁢y^N⁢k for k=1,2,…,K-1, which are the elements along the diagonal of V⁢y^N. Let G^k=𝑚𝑖𝑠𝑠𝑖𝑛𝑔r⁢o⁢w⁢(gj)1⩽j⩽K+1=[r^k⁢N⁢N^p⁢1⁢N^p⁢2⁢…⁢N^p⁢K] be a row vector whose K+1 components are the estimators on which y^r⁢k depends. The y^N⁢k sampling variance is given approximately by V⁢y^N⁢k=𝑚𝑖𝑠𝑠𝑖𝑛𝑔r⁢o⁢w1⩽j⩽K+1⁢(∂⁡y^r⁢k∂⁡gj)⁢V⁢G^k⁢𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽j⩽K+1⁢(∂⁡y^r⁢k∂⁡gj), where V⁢G^k=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽j⩽K+1⁢𝑚𝑖𝑠𝑠𝑖𝑛𝑔r⁢o⁢w1⩽j′⩽K+1⁢(𝐶𝑜𝑣⁢(gj,gj′)) is the design-based covariance matrix of G^k. The sampling covariance matrix y^N is estimated replacing μi by μ^i in V⁢y^N [13].

Domain-level estimators

The survey variable total in the domain R is yNR=∑i=1NRyi=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K⁢(∑i=1NRyi⁢k)=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K⁢(yNR⁢k), where yNR⁢k is the number of pixels in R covered by crop k. To estimate yNR, we use the sample s of size n selected from the population with inclusion probabilities {πi;i=1,2,…,N} and the estimator y^NR=∑i=1NRμ^i+NRN^R⁢p⁢∑i=1nRyi-μ^iπi=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K⁢(y^NR⁢k). Here, nR=∑i=1nIi is the number of units in the sample belonging to the study domain, and N^R⁢p=∑i=1nRIiπi. The sampling variance of yNR is given approximately by V⁢y^NR=NR2⁢V⁢1N^R⁢p⁢∑i=1nRyi-μiπi=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽k⩽K-1⁢𝑚𝑖𝑠𝑠𝑖𝑛𝑔r⁢o⁢w1⩽k′⩽K-1⁢(NR2⁢𝐶𝑜𝑣⁢(y^𝑅𝑟𝑘,y^R⁢r⁢k′)). Let y^𝑅𝑟𝑘=1N^𝑅𝑝⁢∑i=1nRyk⁢i-μk⁢iπi and G^𝑅𝑘=𝑟𝑜𝑤⁢(gR⁢j)1⩽j⩽K+1=[r^𝑘𝑁R N^NR⁢p⁢1N^NR⁢p⁢2…N^NR⁢p⁢K]. The y^NR⁢k sampling variance is given approximately by V⁢y^NR⁢k=𝑚𝑖𝑠𝑠𝑖𝑛𝑔r⁢o⁢w1⩽j⩽K+1(∂⁡y^𝑅𝑟𝑘∂⁡gR⁢j)⁢V⁢G^𝑅𝑘⁢𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽j⩽K+1⁢(∂⁡y^𝑅𝑟𝑘∂⁡gR⁢j), where V⁢G^R⁢k=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽j⩽K+1⁢𝑚𝑖𝑠𝑠𝑖𝑛𝑔r⁢o⁢w1⩽j′⩽K+1⁢(𝐶𝑜𝑣⁢(gR⁢j,gR⁢j′)). The sampling covariance matrix y^NR is estimated replacing μi by μ^i in V⁢y^NR.

Works in the literature [26, 27, 28] for small area estimation based on multinomial mixed models follow a penalized quasi-likelihood approach. As pointed out by McCulloch and Searle [24], these methods are not completely satisfactory in practice. Those authors recommended instead a linearization of the non-linear multinomial models and used linear mixed models. Additional research is required to achieve a completely satisfactory solution to this problem

2.2.3Cost efficiency

Cost efficiency is the usual criterion for comparing a set of sampling strategies developed for estimating the same characteristic in the same population. Let C𝐺𝐷 be the cost and V𝐺𝐷 the sampling error of the current sampling strategy using only ground data, and let C𝐺𝐷+𝑅𝑆 and V𝐺𝐷+𝑅𝑆 be the cost and sampling error, respectively, of the strategy integrating ground and RS data. The cost efficiency is C𝐺𝐷⁢V𝐺𝐷 for the former and C𝐺𝐷+𝑅𝑆⁢V𝐺𝐷+𝑅𝑆 for the latter.

Although the Sentinel images and software required for the integration of ground and RS data are provided by ESA for free, there are costs for the NSOs such as the time of experts and commercial cloud services required for storage and computations. Because estimation of these other costs is a difficult task, we assume in the following that C𝐺𝐷+𝑅𝑆≃C𝐺𝐷, so the comparison criterion reduces to the efficiency.

The relative efficiency of RS data with respect to ground data is 𝑅𝐸𝑅𝑆=V𝐺𝐷⁢(V𝐺𝐷+𝑅𝑆)-1. If 𝑅𝐸𝑅𝑆>1, then, using RS data, the current sampling error V𝐺𝐷 may be reduced to V𝐺𝐷-V𝐺𝐷+𝑅𝑆=𝑅𝑆𝑒𝑓𝑓𝑒𝑐𝑡⁢V𝐺𝐷 without increasing the survey cost. Here, 𝑅𝑆𝑒𝑓𝑓𝑒𝑐𝑡=1-𝑅𝐸𝑅𝑆-1 is the effect of RS data on the sampling efficiency. In other words, using RS data, the current sample size n can be reduced without loss of accuracy in a quantity equal to n-n𝑅𝑆=𝑅𝑆𝑒𝑓𝑓𝑒𝑐𝑡⁢n. Thus, using RS data, the ground sample size can be n𝑅𝑆=(1-𝑅𝑆𝑒𝑓𝑓𝑒𝑐𝑡)⁢n=𝑅𝐸𝑅𝑆-1⁢n.

We evaluate the effect of RS data on the cost efficiency using either linear or multinomial models. Spain and Ecuador are partner countries using area samples selected in only one stage, with equal inclusion probabilities, and we consider linear models for RS data integration. In Spain, the sampling frame is not stratified and the inclusion probability is πi=n/N. In Ecuador, the sampling frame is stratified and the inclusion probability is πh⁢i=nh/Nh for every sampling unit of the same stratum i=1,2,…,Nh, with h=1,2,…,H.

The sampling variance using only ground data is VG=V⁢y^N=V⁢∑i=1nyiπi=V⁢∑i=1nyin/N=N2⁢1/n⁢∑i=1nyi=N2⁢(1-n/N)⁢1/n⁢(n-1)⁢∑i=1n(yi-y¯)2 in Spain, and VG=V⁢y^N=V⁢∑h=1H∑i=1nhyh⁢iπh⁢i=V⁢∑h=1H∑i=1nhyh⁢inh/Nh=∑h=1NhNh2⁢1/nh⁢∑i=1nhyh⁢i=∑h=1HNh2⁢(1-nh/Nh)⁢1/nh⁢(nh-1)⁢∑i=1nh(yh⁢i-y¯h)2 in Ecuador.

Using RS data, the sampling variance is V𝐺𝐷+𝑅𝑆=V⁢y^N=V⁢∑i=1nyi-xi⁢B^n/N=N2⁢V⁢1/n⁢∑i=1n(yi-xi⁢B^)=N2⁢(1-n/N)⁢1/n⁢(n-2)⁢∑i=1n(yi-xi⁢B^)2 in Spain, where B^=(∑i=1nxiT⁢xi)-1⁢∑i=1nxiT⁢yi, and V𝐺𝐷+𝑅𝑆=Vy^N=∑h=1HNh2(1-nh/Nh)[1/nh(nh-2)]∑i=1nh(yh⁢i-xh⁢iB^π)2 in Ecuador, where B^π=(∑h=1H∑i=1nhxh⁢iT⁢xh⁢iπh⁢i)-1⁢∑h=1H∑i=1nhxh⁢iT⁢yh⁢iπh⁢i.

Using RS data, the current sampling error V𝐺𝐷 can be reduced to a quantity V𝐺𝐷-V𝐺𝐷+𝑅𝑆=𝑅𝑆𝑒𝑓𝑓𝑒𝑐𝑡⁢V𝐺𝐷. In Spain, 𝑅𝐸𝑅𝑆≃∑i=1n(yi-y¯)2⁢(∑i=1n(yi-xi⁢B^)2)-1 and 𝑅𝑆𝑒𝑓𝑓𝑒𝑐𝑡=1-𝑅𝐸𝑅𝑆-1. In Ecuador, 𝑅𝑆𝑒𝑓𝑓𝑒𝑐𝑡=∑h=1Hch⁢nh⁢V𝐺ℎ⁢𝑅𝑆effect,h⁢(∑h=1Hch⁢nh⁢V𝐺ℎ)-1, where 𝑅𝑆effect,h=1-𝑅𝐸RS,h-1, 𝑅𝐸RS,h≃∑i=1nh(yh⁢i-y¯h)2(∑i=1nh(yh⁢i-xh⁢i⁢B^π)2)-1, and ch is the cost per sampling unit in stratum h.

Senegal, Malawi, and Tanzania use point sampling with unequal inclusion probabilities {πi;i=1,2,…, N}, and we consider multinomial models for RS data integration. The sampling variance of the y𝑁𝑘 estimator using only ground data is approximately VG=N2⁢V⁢1N^p⁢∑i=1nyk⁢iπi. Using ground and RS data, the sampling error is VG+𝑅𝑆=N2⁢V⁢1N^p⁢∑i=1nyk⁢i-μk⁢iπi. The relative efficiency of RS data for estimating the total number of pixels covered by crop k is 𝑅𝐸𝑅𝑆𝑘=V⁢1N^p⁢∑i=1nyk⁢iπi×(V⁢1N^p⁢∑i=1nyk⁢i-μk⁢iπi)-1.

RS cost-efficiency relies on the RS contribution to reduce the sampling variance of the residuals yk⁢i-μk⁢i with respect to the sampling variance of the ground data yk⁢i. The sampling variance reduction depends on the reliability of the RS data to estimate the yk⁢i expected values, 𝐸𝑦k⁢i=μk⁢i. Consequently, the accuracy and timeliness of the xi datasets required to estimate μi are key to achieve a significant gain of cost-efficiency.

2.2.4Other applications

As mentioned in the Introduction, “multi-seasonal estimation” and “optimizing the sample design” are applications of RS data useful for agricultural statistics. We developed prototypes for these two applications but we cannot test them here because the required data are unavailable. Thus, we are limited to outlining the developed prototypes.

Multi-seasonal estimation

We consider the annual agricultural cycle divided into, say, four seasons. We use the subscript tm to refer to season m=1,2,3,4 of year t=1,2,…. We want to estimate the total yNtm of the survey variable in m of t and its annual aggregate yNt.

2.2.4.1.1. The sample

We examine a two-phase sampling strategy. The first-phase sample is the NSO sample already in place and used for agricultural statistics. We use the subscript 1 to refer to this sample, for instance, n1 denotes its size. We select a secondphase sample specific to each season from among the sampling units within n1. We choose a splitpanel or supplementary panel sampling design consisting of a panel component and specific component [29].

For the first season, the panel component is a sample of n2⁢p sampling units selected from n1, and the specific sample for that season is a sample of n2⁢et1 sampling units selected from among those included in n1 but not included in n2⁢p. In the second season, the panel component is the same as in season 1, n2⁢p, and the specific sample is a number n2⁢et2 of sampling units chosen from among those included in n1 but not included in either n2⁢p or n2⁢et1. For the third season, the panel component is the same as in seasons 1 and 2 (n2⁢p) and the specific sample is a number n2⁢et3 of sampling units selected from among those included in n1 but not included in n2⁢p, n2⁢et1, or n2⁢et2. In the fourth season, the panel component is the same as in seasons 1, 2 and 3, and the specific sample is a number n2⁢et4 of sampling units chosen from among those included in n1 but not included in n2⁢p, n2⁢et1, n2⁢et2, or n2⁢et3.

2.2.4.1.2. Composite estimators

A composite estimator is a function of single estimators, each defined separately for panel and specific samples. We consider the linear model ytm⁢i=xtm⁢i⁢βtm+εtm⁢i, referring to season tm and both panel and specific samples. Model parameters are estimated using the weigthed estimator and data from the second-phase sample n2⁢s⁢tm (with s=p for the panel sample and s=e for the specific sample), with B^π⁢2⁢stm=(∑i=1n2⁢jtmxstm⁢iT⁢xstm⁢iπ2⁢stm⁢i)-1⁢∑i=1n2⁢jtmxstm⁢iT⁢ystm⁢iπ2⁢stm⁢i. The seasonal survey variable total is estimated separately, using data from the panel and specific samples and projective estimator Y^N⁢st⁢m=xNtm⁢B^π⁢2⁢stm. A design-consistent estimator of the sampling error variance is V^⁢(y^N⁢stm-yN⁢stm)=V⁢∑i=1n2⁢stmε^i⁢N⁢stmπ2⁢stm⁢i, where ε^i⁢2⁢N⁢st⁢m=y2⁢stm⁢i-x2⁢stm⁢i⁢B^π⁢2⁢stm.

Let yNt=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽m⩽4⁢(yNtm) be the (4×1) vector of survey variable totals in the four seasons of year t, and let y^Nt=𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢l1⩽m⩽4⁢𝑚𝑖𝑠𝑠𝑖𝑛𝑔c⁢o⁢lj=e,p⁢(y^N⁢jt⁢m) be the (8×1) vector of total estimators based on the panel and specific samples. We consider the model y^Nt=𝑍𝑦Nt+eNt, where eNt=y^Nt-yNt is the vector of sampling errors and Z=I4⊗12 is an (8×4) matrix indicator of the panel or specific sample.

We assume that panel and specific samples in the same season are independent and, as a result, in the error covariance matrix 𝑉𝑒Nt the covariances between y^𝑁𝑝tm and y^N⁢etm′ are nil. A composite estimator of yNt is y^𝐶𝑁t=(ZT⁢(𝑉𝑒Nt)-1⁢Z)-1⁢ZT⁢(𝑉𝑒Nt)-1⁢y^Nt and its covariance matrix is V⁢y^𝐶𝑁t=(ZT⁢(𝑉𝑒Nt)-1⁢Z)-1. The composite estimator of the annual total yNt=14T⁢yNt is y^𝐶𝑁t=14T⁢y^𝐶𝑁t and its variance is V⁢y^𝐶𝑁t=14T⁢V⁢y^𝐶𝑁t⁢14.

To estimate 𝑉𝑒Nt, we use the empirical correlation and an autoregressive model yi⁢tm=Y¯tm+ei⁢tm, where Y¯tm=1Ntm⁢∑i=1Ntmyi⁢tm and ei⁢tm are modeled using ei⁢tm=ui⁢tm+εi⁢tm. Here, (ui⁢tm,εi⁢tm) designates the zeromean independent random variable 𝐶𝑜𝑣⁢(ui⁢tm,εi⁢tm)= 0 and the variance of εi⁢tm is σε2. ui⁢tm=ρ⁢ui,tm-1+ηi⁢tm is an AR(1) stationary process, where ηi⁢tm is a random zeromean perturbation term whose variance is σ2 [30].

As a result, we have the linear mixed model yi⁢tm=Y¯tm+ui⁢tm+εi⁢tm, whose temporal random component ui⁢tm=ρ⁢ui,tm-1+ηi⁢tm is autocorrelated and 𝐶𝑜𝑣⁢(ηi⁢tm,εi⁢tm)= 0. In this last model, yi⁢tm represents the observed data in sampling unit i in season tm and the right term in the equation is interpreted as follows: Y¯tm+ui⁢tm is the true value of the survey variable and the sum of the population mean Y¯tm and ui⁢tm. The latter is a specific component associated with sampling unit i in season tm and represents deviation from Y¯tm Finally, εi⁢tm is the measurement error of the aforementioned true value.

The autocovariance function of the process ui⁢tm is 𝐶𝑜𝑣⁢(ui⁢tm,ui⁢tm′)=Cu⁢(s), where s=|m′-m| and Cu⁢(s)=σ2⁢ρs⁢(1-ρ2)-1. The variance of ui⁢tm is 𝑉𝑢i⁢t=Cu⁢(0)=σ2⁢(1-ρ2)-1 and the autocorrelation function is 𝐶𝑜𝑟𝑟⁢(ui⁢tm,ui⁢tm′)=Cu⁢(s)/V⁢ui⁢tm=ρs. The autocovariance function of the process ei⁢t is 𝐶𝑜𝑣⁢(ei⁢tm,ei⁢tm′)=𝐶𝑜𝑣⁢(ui⁢tm,ui⁢tm′)=Cu⁢(s) and the variance is 𝑉𝑒i⁢tm=𝑉𝑢i⁢tm+V⁢εi⁢tm=σ2⁢(1-ρ2)-1+σε2. The autocorrelation function for the process ei⁢tm is ρe⁢(s)=𝐶𝑜𝑟𝑟⁢(ei⁢tm,ei⁢tm′)=Cu⁢(s)/V⁢ei⁢tm=I(s=0)+ρs⁢(1+(1-ρ2)⁢κ)-1I(s≠0), where κ=σε2/σ2 and I(.) is an indicator variable whose value is 1 if the argument is true and zero otherwise.

The ratio v=σε2/𝑉𝑒i⁢tm=κ⁢(1-ρ2)⁢(1+κ⁢(1-ρ2))-1 is a measure of the weight of the measurement error relative to the total error. The parameters ρ and κ can be estimated through fitting by nonlinear minimum least squares the theoretical correlation model ρe⁢(s)=I(s=0)+ρs⁢(1+(1-ρ2)⁢κ)-1⁢I(s≠0) to the empirical correlations in 𝑉𝑒Nt.

Optimizing the sampling design

The sample design is optimized to find the design variable values that minimize the sampling variance subject to cost (budget) constraints. The sampling variance depends on the spatial correlation structure of the survey variable and we follow a superpopulation approach to identify this structure. We limit ourselves to area samples and assume that the population values of the survey variable yi constitute a sample generated according to a second-order stationary random process with the following characteristics [31]: the mean is 𝐸𝑦i=μ, variance is 𝑉𝑦i=σ2, and the covariance 𝐶𝑜𝑣⁢(yi,yi′)=σ2⁢ρy⁢(𝑑𝑖𝑠𝑡⁢(si,si′)) between two elementary observations (yi,yi′) at the points of coordinates si and si′ is positive, decreasing when the distance between these points d=𝑑𝑖𝑠𝑡⁢(si,si′) increases.

To assess the correlation structure ρy⁢(d), theoretical variogram and correlogram models have been proposed in the literature. Two frequently used correlogram models are the exponential, ρ(u,v|a,τ)=(1-τ)e-d/a, and the spherical ρ(u,v|a,τ)=(1-τ)e-d/a[1-32⁢da+d32⁢a3] if d⩽a and ρ(u,v|a,τ)= 0 if d>a. Here u is the number of sampling units between si and si′ in the row direction, v the number of sampling units between si and si′ in the column direction, and d=u2+v2. The model parameters are the range rate a and ratio τ=τ0/(τ+0τd). Here τ0 is the nugget effect, i.e., the variation at or near the origin (independent of distance), τd is the partial sill (a function of distance between sampling points) and (τ0+τd) is the sill, i.e. the maximum variation far from the origin.

It was demonstrated by Ambrosio et al. [32] how this approach should be followed to design systematic samples using a land-use map. Here, we propose instead to use a crop-type map [33] to estimate the parameters of the correlogram model. The estimated correlogram can be used to evaluate the anticipated variance [34] as a function of design variables for the set of sampling strategies used in practice [35, 36].

2.2.4.2.1. Anticipated variance

We illustrate the proposed approach considering segments of size n0 and a simple random sample of segments of size n. The variance of the sampling error is V⁢y^N=N2⁢(1-n/N)⁢S2/n, where S2=1/(N-1)⁢∑i=1N(yi-y¯N) is the population variance. The anticipated variance is the model-based expected value of this design-based sampling variance 𝐸𝑉⁢y^N=N2⁢(1-n/N)⁢𝐸𝑆2/n, where 𝐸𝑆2=σ2Ψ(N,n0|a,τ) and Ψ(N,n0|a,τ)=n0(𝑁𝑛0-1)(N-1)-1[1-Φ(N,n0|a,τ)]-𝑁𝑛0(n0-1)(N-1)-1[1-Φ(n0|a,τ)]. Here, Φ(N,n0|a,τ) is the average correlation between pairs of observations over the C𝑁𝑛02 pairs in the population. Φ(n0|a,τ) is the average correlation between pairs of observations over the Cn02 pairs in a segment.

2.2.4.2.2. The optimization problem

The design variables are the segment and sample sizes, and we find optimal values of these variables by solving the following optimization problem: min{(n,n0)}𝐴𝑉Y^=min{(n,n0)}Nh2(1-n/N)σ2/nΨ(N,n0|a,τ), subject to C0+∑h=1LCc⁢n+∑h=1LCw⁢n⁢n0+∑h=1LCk⁢A⁢n⩽C. Here, Cc is the cost of adding a segment to the sample, excluding travel cost but including positioning cost (travel to the first segment visited from the interviewer home base and then back to that base from the last segment visited during the data-collection trip), and Cw is the observing cost, including the cost of locating the segment

The solution to this problem is the optimum segment size n0 and optimum sample size of segments n. In addition to the budget, this optimum solution is conditioned to the correlogram model parameters (a,τ).

3.Data

We used field data collected by the NSO of Spain to illustrate our approach to the integration of continuous ground data in the sampling design using linear models, and field data collected by the NSO of Senegal to illustrate our approach to the integration of categorical ground data into the sampling design using multinomial models.

The two RS products traditionally used as auxiliary data for crop acreage and yield estimation are pixel classification by crop type (for the former) and a set of vegetation indices (for the latter). For training pixel classification models, ground georeferenced data of crop type is required. In many countries, data from the agricultural sector are collected using national household or farm surveys, and no geographic information at parcel level is available.

Table 1

Barley acreage estimates in a 100 × 100 km area of Castilla y León, 2018

Data	Acreage	95% Confidence interval (hectares)		Coeff. of	Relative efficiency
		Uncertainty
	(hectare)	Limits	Amplitude	variation (%)*	of RS data
Ground	236165.4	Lw: 215951.7	40427.2	4.37	–
		Up: 256379.0
Ground + RS	228550.1	Lw: 219699.8	17700.5	1.98	5.2
		Up: 237400.3

*Quotient between root square of the sampling variance and estimate.

Table 2a

Barley yield estimates in a 100 × 100 km area of Castilla y León, 2018. NDVI

Data	Yield	95% Confidence interval (kg/hectare)		Coeff. of	Relative efficiency
		Uncertainty
	(kg/hectare)	Limits	Amplitude	variation (%)	of RS data
Ground	4213.6	Lw: 4093.9	239.2	1.45	–
		Up: 4333.2
Ground + RS	4155.6	Lw: 4033.2	244.9	1.50	0.95
		Up: 4278.1

Table 2b

Barley yield estimates in a 100 × 100 km area of Castilla y León, 2019. LAI

Data	Yield	95% Confidence interval (kg/hectare)		Coeff. of	Relative efficiency
		Uncertainty
	(kg/hectare)	Limits	Amplitude	variation (%)	of RS data
Ground	2352.0	Lw: 2227.4	249.1	2.70	–
		Up: 2476.6
Ground + RS	2327.8	Lw: 2215.1	225.4	2.47	1.22
		Up: 2440.6

Azzari et al. [22] compared several ways of generating the data needed to train pixel classification models when no such parcel-level information is available. Focusing on integrating ground survey data of maize in Malawi and Ethiopia from national household surveys using the Sentinel-2 satellite, the authors evaluated the accuracy of pixel classification producing georeferenced data at parcel level, ranging from the full parcel boundary to only one point (centroid).

The authors concluded that collecting full-parcel boundary data or GPS coordinates of the polygon defined by complete parcel corner points yields the best-quality information for model training; however, the use of mid-sized sample (3000–4000 parcels) plot centroids could perform similarly to full plot boundaries. The authors did not consider the statistical model required to achieve the integration of ground data with Sentinel-2 to improve the accuracy of crop acreage estimates, obtaining design-based consistency.

In Sen4Stat, the focus was on developing an open-source system that permits any user to generate from Sentinel 1 SCL Sentinel-2 L1C and/or L2A images the RS data required to improve agricultural statistics. In addition to vegetation indices, using cloudmask and cloud-free mosaics for optical sensors, and for transforming data from SAR sensors, the system allows the choice of several techniques of supervised classification (including random forest) for pixel classification and crop-type map generation.

Raw data in spectral bands can be directly integrated with ground data using multinomial models. Indeed, as shown by Hogland et al. [10], using spectral-band digital numbers in conjunction with multinomial models is a pixel classification method.

4.Results and discussion

4.1Spain

4.1.1Crop acreage estimates at national level

We consider the area covered by an image of size 100 × 100 km in Castilla y León (Spain), with X coordinates (in meters): (300000,400000) and Y coordinates (4600000,4700000). The sampling unit is a square segment of side 700 meters and sample size in the area is 419 segments.

Ground data are observed at parcel level. However, we aggregated data on crop acreage at segment level for computations, so that yi represents ground data on acreage of the study crop in segment i. Further, the RS data are xi=[1⁢xi], where xi is the number of pixels classified as belonging to the study crop in segment i. Model parameter estimates are B^1⁢π= 0.046 for the independent term and B^2⁢π= 0.903 for the angular coefficient of xi. Results for barley data observed in 2018 are in Table 1.

These results show that using RS data, estimator accuracy improved considerably; the amplitude of the confidence interval decreased and the estimation error was reduced by half. The relative efficiency was high; using RS data, the current sample size could be reduced to less than one fifth without loss of accuracy. This is so because RS data are reliable for crop acreage. These results are design-based, which implies that the results change with the sample design. Thus, the reduced sample should be chosen using the same design as the NSO is currently using.

4.1.2Crop yield estimates

Ground data on yield are observed at parcel level, with yi⁢j representing these data for the study crop in parcel j of segment i. The RS data are denoted xi⁢j=𝑚𝑖𝑠𝑠𝑖𝑛𝑔r⁢o⁢w1⩽l⩽L⁢(xi⁢j⁢l), a row vector with 1 in the first position and a set of vegetation indices in the remaining positions. The latter include the normalized difference vegetation index (NDVI) and leaf area index (LAI) (sum of LAI simulated by a Savitsky-Golay interpolation, fitting all LAI observations in the growing season) during (i) sprouting, (ii) flowering, and (iii) ripeness, plus (iv) maximum value of the LAI S-G interpolation and yield simulated by the simple algorithm for yield (SAFY) crop-growth model. Results for barley data observed in 2018 are in Table 2a for NDVI only and in Table 2b for LAI and yield simulated by SAFY.

The relationship between crop yield and the vegetation indices is statistically significant. However, their correlation is not as strong as for crop acreage and as a result, using RS data, yield estimator accuracy improves little. Using RS data, the estimation error is of the same order of magnitude as using only ground data. The relative efficiency is nearly 1; using RS data, the current sample size could be reduced little without loss of accuracy. This is because RS data reliability for crop yield is not as great as for crop acreage. Additional research is required to make the production estimates more reliable through an improved estimate for yield from RS; physically based models are more reliable and can be integrated with RS data.

4.1.3Crop production estimates

We estimate crop production as the product of the crop acreage and yield estimates. Results are in Table 3.

Table 3

Barley production estimates in a 100 × 100 km area of Castilla y León, 2018

Data	Production	95% confidence interval (tons: 1000 kg)		Coeff. of
		Uncertainty
	(tons) (1000 kg)	Limits	Amplitude	variation (%)
Ground	987841	Lw: 903640	168402	4.35
		Up: 1072042
Ground + RS	965733	Lw: 915194	101078	2.67
		Up: 1016272

Thanks to improvement in the crop acreage estimator using RS data, the production estimator accuracy increased considerably; the estimation error decreased by half, even if the RS data failed to improve the yield estimate.

4.1.4Crop acreage estimates at province level

To illustrate provincial estimation using areal sampling, we consider barley acreage estimates at the provincial level within the 100 × 100 km area of Castilla y León (Spain). Results are in Table 4.

Table 4

Barley acreage estimates at provincial level

Province	Using only ground data		Using ground & RS data		Relative efficiency
	Acreage (has.)	Error (CV%)	Acreage (has.)	Error (CV%)	of RS data
León	6853.2	24.15	6834.5	16.20	2.3
Palencia	88602.0	7.36	90535.3	3.33	4.7
Valladolid	128209.5	5.57	119707.4	2.66	5.1
Zamora	12324.2	17.71	10948.2	8.16	6.4
Total area	235989.1	4.37	228028.5	1.98	5.2

The estimate accuracy at provincial level is, as expected, less than at the national level, but the RS relative efficiency is of the same order of magnitude. In provinces where the sampling error using only ground data is very high (León), the RS contribution is qualitative in the sense that it allows the estimates to be labeled as official statistics by reducing the coefficient of variation below the standard limits (20%) considered acceptable in official statistics.

4.1.5Crop acreage estimates at the municipality level based on linear mixed models

To illustrate small-area estimation using areal sampling, we consider barley acreage estimates at municipality level in that part of Zamora province within the 100 × 100 km area of Castilla y León (Spain). Results are in Table 5.

Table 5

Barley acreage estimates at municipality level

Municipality	Sample	Acreage
	size	Hectares	Coeff. of
	(segments)		variation (%)
Belver de los Montes	1	212.96	29.1
Castroverde	3	2914.22	8.0
Pinilla de Toro	4	963.30	10.0
Quintanilla del Monte	1	466.65	20.3
Toro	1	615.91	14.0
Vezdemarbán	3	1358.22	12.6
Villalpando	2	560.05	39.1
Villamayor de Campos	1	1056.23	11.1
Villanueva del Campo	1	784.03	13.2
Villar de Fallaves	1	844.16	11.0
Villardondiego	1	516.40	11.5
Villavendimio	1	656.07	10.4
Total Zamora	20	10948.20	8.2

Considering that the sample size at municipality level is small or null, the accuracy of municipality estimates is good, thanks to the RS contribution. In most municipalities, the estimates could be labeled as official statistics because the coefficient of variation is smaller than the standard limit (20%), even when the sample size is only one sampling unit

4.2Senegal

4.2.1Crop acreage estimates in Nioro Department

To demonstrate crop acreage estimation using point sampling and multinomial models, we consider the Nioro Department of Senegal, in the Kaolack Region. The ground data were observed in a sample of 345 points pixels selected using a list frame. The source of RS data is a crop-type map. In fact, we used a dual frame, since we complemented the NSO list frame with an area frame based on satellite images, in which agricultural areas were distinguished from non-agricultural areas and were stratified into crop types.

In this case, the NSO list frame alone is not sufficient to integrate field and RS data This is because the number of pixels in the agricultural areas is required for expanding the sample estimates to the entire population. This expansion factor is provided by the area frame.

We focus on the two main crops (millet and groundnut) observed in the field sample. The remaining crops (mainly maize) are included in a third category (K=3) called other, whose probability estimate is one minus the estimate’s probability of millet and groundnut. The RS data are coded according to the Earth Observation (EO) crop type into which pixels are classified: xi=[1 0 0 0] for any pixel i in the EO class maize, xi=[0 1 0 0] for any pixel i in the EO class millet, xi=[0 0 1 0] for any pixel i in the EO class groundnut, and xi=[0 0 0 1] for any pixel i in the EO class other crops. Model parameter estimates are in Table 6.

Table 6

Model parameter estimates (B^π⁢𝑚𝑖𝑙𝑙𝑒𝑡 and B^π⁢𝑔𝑟𝑜𝑢𝑛𝑑𝑛𝑢𝑡)

Crop	EO crop type map
	Maize	Millet	Groundnut	Other crops
Millet	-0.010655	1.47958334	0.88931346	8.76504093
Groundnut	-0.659204	-0.59871700	2.89873948	1.21898514

Table 7

Crop acreage estimates, Nioro (Senegal)

		Uncertainty
Crop type	Hectare	Standard	Coefficient of	Limits of 95% confidence interval
		error	variation (%)	Lower	Upper	Amplitude
Millet	89215	3661.1	4.11	81978.8	96330.4	14351.5
Groundnut	78815	2923.9	3.71	73089.1	84550.9	11461.8

Table 8

Efficiency of RS data for crop acreage estimation Nioro (Senegal)

	Standard errors of proportion estimators		Relative efficiency
Crop type	Using only ground data	Using ground & RS data	of RS data
Millet	3.37	1.90	3.13
Groundnut	3.34	1.52	4.80

Table 9

Confusion matrix based on the sample

Ground data		Number of pixels in the EO class (%)*
Crop	Number of points/pixels in the sample	Maize	Millet	Groundnut	Other
Maize	48	9 (18.8)	26 (54.2)	12 (25.0)	1 (2.0)
Millet	134	10 (7.4)	97 (72.4)	23 (17.2)	4 (3.0)
Groundnut	163	5 (3.1)	14 (8.6)	143 (87.7)	1 (0.6)
Total	345	24	137	178	6

*% Percentage of pixels correctly classified.

Crop acreage estimates of millet and groundnut based on ground and RS data are in Table 7.

We evaluate the RS data efficiency 𝑅𝐸𝑅𝑆𝑘 for the acreage estimation of millet and groundnut. Results are in Table 8.

The effect of integrating RS data in the ground sample data was a reduction in the sampling error of millet and groundnut and, as a result, in the confidence interval of these two crops, without loss of design-based consistency. In other words, using RS data, the cost of estimating millet and groundnut acreage could be reduced to less than a third of the current cost, without loss of accuracy. These results are design-based, so it is understood that the reduced sample size should be selected using the currently sampling design used by the NSO.

4.2.2Estimation directly based on pixel classification

There is consensus in the official statistics community on using methods providing design-consistent estimates of the population characteristics under study and measures of estimator uncertainty, such as sampling error, coefficient of variation, or confidence intervals. The method proposed in this paper for integrating RS data in the NSO sampling design agrees with this consensus.

Table 10

Probability that the cover of a pixel of the EO class j is the crop k

EO class (j)	Millet (k=1)	Groundnut (k=2)	Maize & Other_crops (k=3)
Maize (j=1)	0.395	0.206	0.399
Millet (j=2)	0.739	0.093	0.168
Groundnut (j=3)	0.113	0.841	0.046
Other (j=4)	0.997	0.001	0.002

Table 11

Estimates number of pixels, N^k

EO class		Estimates of the number of pixels covered by crop k in each EO crop type class, N^k⁢j=Nj×μ^k⁢j
Crop type (j)	Number of pixels	Millet	Groundnut	Maize & Other
	Nj	(k=1)	(k=2)	(k=3)
Maize (j=1)	1032628	407587	213088	411953
Millet (j=2)	9470572	7000329	876040	1594203
Groundnut (j=3)	8076448	910538	6791735	374175
Other (j=4)	604853	603038	605	1210
Estimates of the total number of pixels by crop N^k=∑j=1JN^k⁢j	19184501	8921492	7881468	2381541

Table 12

Crop acreage estimates at the district (arrondissement) level, Nioro (Senegal)

Arrondissement	Millet		Groundnut
	Acreage (has.)	Error (CV%)	Acreage (has.)	Error (CV%)
Medina Sabakh	20067.2	8.6	19765.3	7.3
Paoskoto	38316.0	5.3	35018.2	4.0
Wack Ngouna	30831.7	11.9	24030.7	10.7
Total Nioro	89215,0	4.1	78815,0	3.7

Although the proposed approach allows for any form of RS data, including raw reflectance data in the form of pixel digital numbers, we focus on the use of croptype maps as auxiliary information. In this context, a natural question that arises is: Why not directly use the number of pixels in each EO croptype class to estimate crop acreage instead of using it as auxiliary data in a statistical model?

We follow the suggestion of a referee to clarify this question, comparing the results of our design-based approach (shown in Table 7) with estimates directly based on the croptype map. The comparison is necessarily limited to the crop acreage estimates because the usual algorithms used for croptype map generation do not provide measures of uncertainty comparable to those of Table 7.

Using a random forest classifier, four EO crop types were considered in the croptype map of Nioro: maize, millet, groundnut, and others. The remaining pixels were classified in a non-agricultural landuse class. The number of pixels in each EO croptype class is in Table 11.

As seen in Table 9, there is uncertainty in the pixels classification and any estimate based on it are subject to this uncertainty. Most pixels of millet in the sample (72.4%) are correctly classified in the EO class millet, but a non-negligible percentage of them were confounded with groundnut (17.2%), maize (7.4%), and others (3.0%). For groundnut, the percentage of pixels that were correctly classified is higher (87.7%) than for millet, but the confusion with millet (8.6%) and maize (3.1%) is non-negligible.

In the proposed approach, the multinomial logit model is used for estimating the probability μk⁢j that the crop covering a pixel classified in the EO class j is actually crop k. The model parameter estimates are in Table 6 and the probability estimates are in Table 10.

To estimate the number of pixels Nk⁢j in EO class j that are actually covered by crop k, we multiplied the total number Nj of pixels in EO class j by the estimator of the aforementioned probability μ^k⁢j: N^k⁢j=Nj×μ^k⁢j. The estimator of the total number of pixels covered by crop k in Nioro, N^k=∑j=1JN^k⁢j is the sum of the estimators in each EO class. This estimator is design-consistent and the estimates based on it are in Table 11.

Both the number of pixels in EO crop type millet (9470572, i.e., 94705.72 hectares, as a pixel represents 100 m2) and in EO crop type groundnut (8076448, i.e., 80764.48 hectares) are larger than the design-based number estimates: 8921492 (89214.92 hectares) for millet and 7881468 (78814.68 hectares) for groundnut. For millet, the difference was 5490.8 hectares (6.2%) and for groundnut 1949.8 hectares (2.5%). For the former, the difference is greater than the design-based standard error (3661.10 hectares) and the coefficient of variation (4.11%). For the latter the difference is smaller than the design-based standard error (2923.94 hectares) and the coefficient of variation (3.71%). These differences are not statistically significant, since the estimates directly based on the EO croptype classes are within the design-based confidence limits, namely, [81978.88, 96330.4] hectares for millet and [73089.15, 84550.98] hectares for groundnut.

However, the comparison must focus not on the results, which are always uncertain, but on the methods. The methods make the major difference; whereas the design-based estimators are in the mainstream of official statistics, the estimators based directly on croptype maps are not. The proposed approach provides designconsistent estimates together with the usual measures of uncertainty (sampling error, coefficient of variation, and confidence intervals). The two main objections to estimators based directly on the croptype map are that (i) they are not design-consistent and (ii) the usual uncertainty measures are not available.

4.2.3Crop acreage estimates at district level in Nioro Department

To illustrate crop acreage estimation at district level using point sampling and multinomial models, we considered the three districts (arrondissements) of Nioro, Medina Sabakh, Paoskoto, and Wack Ngouna. The estimates are in Table 12.

The estimate accuracy at district (arrondissement) level is, as expected, lower than at the national level. However, the sampling error is within the standard limit (coefficient of variation < 20%) for labeling as official statistics.

The differences between the number of pixels in the EO crop type millet and the design-based estimate of the number of pixels covered by millet are smaller than the sampling error in Medina Sabakh (0.7%) and Paoskoto (4.1), and larger than the design-based estimate in Wack Ngouna (12.4%). The differences between the number of pixels in the EO crop type groundnut and the design-based estimate of the number of pixels covered by groundnut are smaller than sampling error in the three districts: Medina Sabakh (3.4%), Paoskoto (2.5%), and Wack Ngouna (1.2%).

5.Concluding remarks

A result well established in the literature is that RS data improve the cost-efficiency of design-consistent crop acreage estimators based on linear models and simple area samples of georeferenced polygons (segments). The strong RS reliability for crop acreage estimation facilitates a notable improvement of crop production estimates, even if RS data are not sufficiently reliable for yields.

However, many countries use complex list samples of non-georeferenced households or farms for agricultural statistics. We have demonstrated in this paper that with the use of multinomial models, a point georeferenced by parcel in these complex samples is sufficient to enhance the cost efficiency of design-consistent, national-level crop acreage estimators using RS data.

The open-source Sen4Stat system is being designed in such a way that any NSO can use it to improve the cost efficiency of the procedure currently used to obtain crop acreage and production estimates, at both national and minor administrative area levels (province and municipality). The only inputs required by the system are the observed data in the sample currently designed by the NSO for field data collection. The system allows for both a segment simple sample, based on area frames, and a point complex sample, based on list frames of households/farms. In both cases, the system generates from Sentinel images a crop-type map and integrates it with ground data, using linear models in the segment case and multinomial models in the point case.

We have extended the method developed in this paper to two additional RS applications, multi-season estimates and optimization of the sample design. Additional work is required to obtain the data required to test these two additional RS applications.

Acknowledgments

This research was carried out in the framework of a consortium integrating the Université Catholique de Louvain (UCLouvain) as contractor and the Universidad Politécnica de Madrid (UPM), CS ROMANIA S.A. (C-S RO) and Systèmes d’Information à Référence Spatial (SIRS/CLS) as sub-contractors. The national statistical offices of Ecuador, Malawi, Senegal, Spain, and Tanzania participated in the work as partners and provided the required ground data. SIRS/CLS elaborated the ground databases, C-S RO developed the system for image classification, and UPM developed statistical prototypes for the integration of ground and RS data. UCLouvain coordinated the project and contributed mainly to the specification of products and services, system design RS database elaboration for integration with the ground database, and the full-scale demonstration plan.

We thank the referees and the subject editor for their very helpful comments and suggestions.

Funding

This work was supported by the European Space Agency in the Sen4Stat project framework [Contract No. 4000127181/19/I-NS]. This agency contributed to the study design, together with a steering committee made up of members of international organizations such as the World Bank and FAO, with the aim of achieving synergies between this project and those launched by these organizations to improve agricultural statistics.

References

[1]	Fonteneau F, Delincé J. Surveying Farms in the 21st Century. In: FAO editor. Handbook on the Agricultural Integrated Survey (AGRIS); Global Strategy to Improve Agricultural and Rural Statistics. Rome: FAO; (2018) ; pp. 1-6.
[2]	Zezza A, Gourlay S, Molini V. Closing the data gap in agriculture through sustainable investment in the data value chain: Realizing the vision of the 50x2030 Initiative. Stat J IAOS. (2022) ; 38: : 57-62. doi: 10.3233/SJI-220933.
[3]	Hanuschak GA, Allen RD, Wigton WH. Integration of Landsat data into the crop estimation program of USDA’s Statistical Reporting Service (1972-1982). Invited paper at: 8th International Symposium on Machine Processing of Remote Sensed Data. 1982 July 7–9; West Lafayette, Indiana: Purdue University. Available from: https//www.lars.purdue.edu/home/references/sym_1982/sym_1982.html.
[4]	Allen JD. A look at the Remote Sensing Applications Program of the National Agricultural Statistics Service. J Off Stat. (1990) ; 6: (4): 393-409.
[5]	Ambrosio L, Alonso R, Villa A. Estimación de superficies cultivadas por muestreo de áreas y teledetección. Precisión relativa. Estad Esp. (1993) ; 35: (132): 91-103.
[6]	Delincé J. Cost-Effectiveness of Remote Sensing for Agricultural Statistics in Developing and Emerging Economies. In: FAO editor. GSARS Technical Report GO09-2015. Rome: FAO; 2015. doi: 10.13140/RG.2.2.25985.45927.
[7]	Gallego FJ. Remote sensing and land cover area estimation. Int J Remote Sens. (2004) ; 25: (15): 3019-47. doi: 10.1080/01431160310001619607.
[8]	Firth D, Bennett KE. Robust models on probability sampling. J. R. Stat. Soc. Series B Stat. Methodol. (1998) ; 60: :3-21. doi: 10.1111/1467-9868.00105.
[9]	Agresti A. Categorical data analysis. New York: Wiley; (2002) ; 710 p.
[10]	Hogland J, Billor N, Anderson N. Comparison of standard maximum likelihood classification and polytomous logistic regression used in remote sensing. Eur J Remote Sens. (2013) ; 46: (1): 623-40. doi: 10.5721/EuJRS20134637.
[11]	FAO. Handbook on Master Sampling Frames for Agricultural Statistics. Frame Development, Sample Design and Estimation. Improving Agricultural and Rural Statistics. Global Strategy. Rome: FAO; (2015) Dec. 170 p.
[12]	Särndal CE, Swensson B, Wretman, J. Model Assisted Survey Sampling. New York: Springer; (1992) ; 695 p.
[13]	Fuller WA. Sampling statistics. New York: Wiley; (2009) ; 472 p. doi: 10.1002/9780470523551.
[14]	Cotter J, Nealon J. Area frame design for agricultural surveys. Research and Applications Division. National Agricultural Statistics Service. United States Department of Agriculture. Washington, D.C., (1987) .
[15]	FAO. Multiple frame agricultural surveys. Volume 1: Current surveys based on area and list sampling methods. Volume 2: Agricultural survey programmes based on area frame or dual frame (area and list) sample designs. Statistical development series 7 and 10. Rome: FAO; 1996, 1998.
[16]	Ministry of Agriculture, Livestock and Fisheries. 2014/15 Annual Agricultural Sample Survey Report. Zanzibar (TZ); (2016) ; 72 p. Available from: http//hdl.handle.net/20.500.12018/2905.
[17]	Yates F. Sampling Methods for Census and Surveys. London: Griffin. (1949) ; 318 p.
[18]	Jinguji I. Dot sampling method for area estimation. In: Crop Monitoring for Improved Food Security. Srivastava M.K. editor, Bangkok (TH): FAO and ADB; (2015) RAP PUBLICATION 2014/28: pp. 27-48.
[19]	Gay C, Porchier JC. Land cover and land use classification using TER-UTI. In: Holland TE and van den Broecke MPR, editors. Proceedings of Agricultural Statistics 2000, an International Conference on Agricultural Statistics 1998 Mar 18-20: pp. 93-201. Available from: https//www.isi-web.org/isi.cbs.nl/iamamember/Books/agric2000/page-193.pdf.
[20]	Nusser SM, Goebel JJ. The National Resources Inventory: a long-term multi-resource monitoring programme. Environ Ecol Stat. (1997) ; 4: (3): 181-204. doi: 10.1023/A1018574412308.
[21]	Nusser SM, Goebel JJ, Fuller WA. Design and estimation for investigating the dynamics of natural resources. Ecol Appl. (1998) ; 8: (2): 234-45. doi: 10.1890/1051-0761(1998)008[0234:DAEFIT]2.0.CO;2.
[22]	Azzari G, Jain S, Jeffries G, Kilic T, Murray S. Understanding the requirements for surveys to support satellite-based crop type mapping: evidence from Sub-Saharan Africa. Remote Sens. (2021) ; 13: (23), 4749. doi: 10.3390/rs13234749.
[23]	Direction de l’Analyse de la Prévision et des Statistiques Agricoles. Ministère de l’Agriculture et de l’Equipement Rural. Senegal – Enquete Agricole Annuelle. Methodologie et plan de sondage de l’enquete agricole. Senegal. (2017) . http//anads.ansd.sn/index.php/catalog/218/download/1846.
[24]	MDcCulloch CE, Searle SR. Generalized, linear, and mixed models. New York: Wiley; (2001) ; 325 p.
[25]	Ambrosio L, Iglesias L. Land cover estimation in small areas using ground surveys and remote sensing. Remote Sens. Environ. (2000) ; 74: : 240-8. doi: 10.1016/S0034-4257(00)00114-0.
[26]	Saei A, Chambers R. Small area estimation under linear and generalized linear mixed models with time and area effects. Methodology Working Paper M03/15. Southampton (UK): Southampton Statistical Sciences Research Institute, (2003) : 31 p. http//eprints.soton.ac.uk/id/eprint/8165.
[27]	Molina I, Saei A, Lombardia MJ. Small area estimates of labour force participation under a multinomial logit mixed model. J. R. Stat. Soc. Ser. A Stat. Soc. (2007) ; 170: (4): 975-1000. doi: 10.1111/j.1467-985X.2007.00493.x.
[28]	López-Vizcaíno E, Lombardía MJ, Morales D. Small area estimation of labour force indicators under a multinomial model with correlated time and area effects. J. R. Stat. Soc. Ser. A Stat. Soc. (2015) ; 178: (3): 535-65. doi: 10.1111/rssa.12085.
[29]	Fuller WA. Environmental surveys over time. J Agric Biol Environ Stat. (1999) ; 4: (4): 331-45. doi: 10.2307/1400493.
[30]	Fuller WA, Breidt FJ. Estimation for supplemented panels. The Indian Journal of Statistics, Serie B. (1999) ; 61: (1): 58-70. https//www.jstor.org/stable/25053068.
[31]	Cochran WG. Relative accuracy of systematic and stratified random samples for a certain class of populations. Annals of Mathematical Statistics. (1946) ; 17: (2): 164-77.
[32]	Ambrosio L, Iglesias L, Marin C. Systematic sample design for the estimation of spatial means. Environmetrics. (2003) ; 14: (1): 45-61. doi: 10.1002/env.564.
[33]	Defourny P, Bontemps S, Bellemans N, Cara C, Dedieu G, Guzzonato E, el al. Near real-time agriculture monitoring at national scale at parcel resolution: performance assessment of the Sen2-Agri automated system in various cropping systems around the world. Remote Sens. Environ. (2019) ; 221: : 551-68. doi: 10.1016/j.rse.2018.11.007.
[34]	Fuller WA, Isaki CT. Survey design under superpolutation models. In: Current Topics in Survey Sampling. Krewski D, Rao JNK and Platek R, editors, New York: Academic Press; (1981) ; pp. 199-226.
[35]	Ambrosio L, Iglesias L. Identifying the most appropriate sampling frame for specific landscape types. Technical Report Series. GO-01-2014. FAO. https//www.fao.org/3/ca6436en/ca6436en.pdf.
[36]	Ambrosio L, Iglesias L, Marín C. A model-assisted approach to identify a cost-efficient spatial sampling strategy. Paper presented at: Eighth International Conference on Agricultural Statistics (ICAS-VIII); (2019) Nov 18-21; New Delhi (IN).