Reliable Crop Identification with Satellite Imagery in the Context of Common Agriculture Policy Subsidy Control

Agricultural subsidies in the context of the Common Agricultural Policy (CAP) represent over 40% of the EU’s yearly budget. To ensure that funds are properly spent, farmers are controlled by National Control and Paying Agencies (NCPA) using tools, such as computer-assisted photo interpretation (CAPI), which aims at identifying crops via remotely-sensed imagery. CAPI is time consuming and requires a large team of skilled photo interpreters. The objective of this study was to develop a reliable control system to partially replace CAPI for crop identification, with the overreaching goal of reducing control costs and completion time. Validated control data provided by the Portuguese Control and Paying Agency and an atmospherically-corrected Landsat ETM+ time series were used to perform parcel-based crop classification, leading to an accuracy of only 68% due to high similarity between crops’ spectral signatures. To address this problem, we propose an automatic control system (ACS) that couples crop classification to a reliability requirement. This allows the decision-maker to set a reliability level, which restricts automatic crop identification to parcels that are classified with high certainty. While higher reliability levels reduce the risk of misclassifications, lower levels increase the proportion of automatic control decisions (ACP). With a reliability level of 80%, more than half of the parcels in our study area are automatically identified with an overall accuracy of 84%. In particular, this allows automatically controlling over 85% of all parcels classified as maize, rice, wheat or vineyard. Remote Sens. 2015, 7 9326


Introduction
In this section, the problem we address is described, and the main motivations and goals of this research are discussed.The reader can find a list of the most important acronyms at the end of the paper.

Use of Remote Sensing for CAP Subsidy Control
The Common Agricultural Policy (CAP) is a system of European Union (EU) agricultural subsidies and programs that is very significant in financial terms, representing over 40% of the EU's budget, equivalent to e58 billion in 2011 [1].To ensure that CAP funds are spent appropriately, Member State Authorities have to comply with legal management and control mechanisms [2].Toward that end, European Commission (EC)'s Joint Research Centre (JRC) provides technical support to Member States.To date, Control with Remote Sensing (CwRS), Digital Land Parcel Identification System (LPIS) and parcel area measurement using Global Navigation Satellite Systems (GNSS) devices have become the keystones of the efficient administration and control of CAP subsidies [3].Each Member State is responsible for subsidy administration and control, which are carried out by a National Control and Paying Agency (NCPA).
In order to obtain area-based financial support, farmers are required to submit an application to their NCPA early in the year, where they declare the precise location of all of their agricultural parcels, as well as the crop type.The National Agency is responsible for controlling at least 5% of those declarations and penalizing farmers who submit incorrect information by performing so-called On-The-Spot (OTS) checks.For area-based subsidies, an agricultural parcel must be controlled at two different levels: both the declared crop and area must be correct [4].The EC in turn controls the NCPAs.When discrepancies between the control result and the reality are found, a Member State is penalized and has to return to the EU part of the subsidies that were distributed to farmers.
The complex process of subsidy control requires computational tools: NCPAs rely on Integrated Administration and Control System (IACS)s, which includes a LPIS.The main functions of those spatial databases are localization, identification and quantification of agricultural land via detailed geospatial data, in order to facilitate the distribution of CAP subsidies [5].
CwRS, the goal of which is to perform OTS checks "in the office" as much as possible, has become an official method for NCPAs to carry out part or all of their OTS checks of EU farms [2].It uses remotely sensed data, which are acquired by the EU and distributed to the NCPAs by JRC.The CwRS program enabled in 2013 the control of over 400,000 farmers' area-based subsidy applications and represented 70% of the EU OTS checks [3].NCPAs have been relying on Computer Assisted Photo-Interpretation (CAPI) of High Resolution (HR) and Very High Resolution (VHR) imagery time series to perform CwRS.However, CAPI is time consuming; it depends on the skills of the photo interpreters and requires a large team.Therefore, there is an obvious need for the development of cost-effective methods to automatize crop identification, which are essential for control in the context of area-based agricultural subsidies, in order to lower costs, speed up work and improve reliability compared to CAPI [4,6].However, the accuracy of automatic classification in an operational context is critical and must be taken into account during the development of such automatic identification methodologies.

Parcel-Based Crop Identification
Multitemporal and multispectral remote sensing imagery has been widely used for crop identification in past years [7,8] since time series of satellite images are believed to be a cost-effective data source to assess land cover, such as agricultural crops over large areas [9].The basis for separating one crop from another is the supposition that each crop species has a specific spectral signature in a time series of multispectral images.However, major limitations on crop identification with satellite imagery, like the similarity of the plant reflectance of different crops, parcel-to-parcel variability of the plant reflectance of the same crops and changes in the pattern of individual crop phenology, may occur [8].
Pixel-based classification approaches, where pixels are classified individually regardless of their spatial aggregation, generally lead to poor results [8,10].To overcome this problem, object-based techniques have been increasingly used in remotely-sensed image analysis.Object-Based Image Analysis (OBIA) typically relies on segmentation algorithms to define geographical objects that can be classified; this is known as "object-based classification" [10,11].However, when a spatial database, such as the LPIS, is available, segmentation algorithms might not be necessary, since agricultural parcels extracted from LPIS can in principle be used instead of objects resulting from image segmentation; this is what we call "parcel-based classification".The operational applicability of parcel-based approaches is conditioned by the geographic information available in the national LPIS.As described in [5], the LPIS Conceptual Model (LCM) requires the geometry of the reference parcel to be stored.However, and depending on the type of reference parcel used by each Member State, it might not represent a single crop.As discussed in [5], the LCM class Reference Parcel is associated with a land cover type, which includes, among others, arable land, irrigated rice or grassland, but information about specific crops, like wheat or barley, are stored in a separate class (Agricultural Parcel) with no geometry.To complicate matters, the reference parcel might include very small non-eligible features.Nonetheless, the location of individual agriculture parcels is included in LPIS, in the form of the farmer's sketch, as required by EU regulation [12].When this information is available in digital format, it can be used for parcel-based crop identification.Moreover, some classes of interest for aid application correspond to land cover types (e.g., vineyard) and are therefore represented geometrically in LPIS.
Several studies have been conducted on automatic crop identification in the context of CAP subsidies using object and parcel-based image classification.Matikainen et al. [13] developed a method for automatic change detection in the Finnish LPIS.They used segmentation of aerial orthoimagery inside of each of the LPIS parcels.Each of the resulting objects were classified with a decision tree and then compared to detect changes.Their results suggest that real changes in LPIS can be detected relatively well, but a large number of false detections occur.Blaes et al. [4] developed a parcel-based strategy for crop identification in the operational context of subsidy controls.Those authors used farm plot outlines digitized from farmers' declarations and determined the best combination of acquisition dates and image sensors, including radar sensors, that maximize the detection of incorrect farmer's declarations.Oesterle and Hahn [14] also used segmentation of VHR aerial orthoimagery to investigate the automatic update of the German LPIS.

Reliable Crop Identification
Given the large amounts of funds involved in the CAP and the risk of penalization by the EC, it is crucial that automatic crop identification is as precise as possible in the context of operational subsidy control.We believe that the conventional "one size fits all" approach, where the crop classification is applied to all parcels, is not appropriate to develop a control system that produces reliable results.Instead, the automatic approach should be able to identify the largest possible set of agriculture parcels that can be automatically and correctly classified as determined by a reliability level set by the decision-maker.The remaining parcels, in turn, should be subjected to visual analysis in the form of CAPI or farm visits, since this is the only form to guarantee the quality of the decision process.

Objectives of the Study
The goal of this study was to develop a simple and cost-effective Automatic Control System (ACS) to automatize the control process and reduce OTS check costs and completion time.The proposed solution consists of the automatic classification of as many agricultural parcels as possible, eliminating the need of CAPI for controlling those parcels.Our ACS comprises a reliability level, chosen by the decision-maker, which allows our approach to focus on parcels that can be classified automatically with high certainty.For those parcels, we performed a parcel-based classification of multispectral and multitemporal land cover signatures retrieved from remotely-sensed imagery, taking advantage of the object-based nature of agricultural parcels.We used field-validated OTS check data from the Portuguese NCPA (Instituto de Financiamento da Agricultura e Pescas (IFAP)) and an atmospherically-corrected multispectral Landsat 7 Enhanced Thematic Mapper (ETM+) time series to develop and train our methodology for twelve major land cover classes in the Portuguese Ribatejo agricultural landscape.
Although this study responds to a practical problem raised by the Portuguese NCPA, our approach is conceptual and could be in principle applied to different areas, crops and remote sensing imagery.Our case study happens to be situated in an area with a Mediterranean climate, where cloud-free images are typically available with high temporal resolution during a large part of crops' development cycles (spring and summer).However, our ACS could be used as long as a dense time series of images is available for the area of interest.In particular, SAR data could be combined with optical imagery for cloud-persistent regions, as in [4].

Study Area Description
The study area is located within the Portuguese province Ribatejo, northeast of Lisbon, including mainly parts of the Lisbon and the Santarém districts (Figure 1).It extends over an area of about 6390 km 2 and is situated between longitudes 9 • 6 0 W and 8 • 9 36 W and latitudes 38 • 43 47 N and 39 • 27 36 N (Datum WGS84).Agriculture is the main activity in Ribatejo with a wide range of different commodities, containing some of Portugal's richest agricultural land.This agricultural area is situated within the Tagus River basin, which plays an important role for both the agricultural activity and climate of the region.The Tagus River is also responsible for the relatively flat relief of the whole study area, with a height ranging between zero and 200 m above sea level.The study area is characterized by a Mediterranean climate with long hot and dry summers and moderate, rainy winters.On average, annual mean temperature is between 15 • C and 16 • C and rainfall ranges from 600 to 800 mm/year, with many crops being irrigated during the dry summer growing season.The very diverse cropping pattern found in the Ribatejo province was the main reason that this area was selected for our study, along with the good availability of field-validated OTS control data from IFAP.

Remote Sensing Data
Table 1 shows the five optimal periods used for multitemporal image acquisition by IFAP to perform CAPI [15].IFAP usually requests one HR image per period to JRC.These periods were chosen in order to optimally follow the annual agricultural cycle of both winter and summer crops cultivated in Portugal.Images for this study were also screened for atmospheric effects and clouds.
Based on these criteria, six Landsat-7 ETM+ multispectral images over the study area (WRS-2 Path 204, Row 33) from a period between November 2004 and August 2005 were acquired through USGS EarthExplorer (Table 1).The imagery consisted of atmospherically-corrected surface reflectance data (climate data records), generated from the LEDAPS tool.LEDAPS uses MODIS atmospheric correction routines and complex 6S radiative transfer models in order to generate surface reflectance [16].The six visible and short-wave infrared bands (Bands 1 to 5 and 7) with a spatial resolution of 30 m were used.The bands from all multispectral images in the temporal sequence (six bands for each image) were stacked, forming a single multi-date image with 36 bands.Landsat data was primarily chosen due to its cost-effectiveness and good temporal coverage.In addition, Landsat ETM+ imagery is known to be appropriate for detailed large-area crop mapping [17].We were limited to Landsat 7 ETM+ imagery for this study, since no other Landsat images were available at EarthExplorer for the study area and period.All acquired images were affected by the Scan Line Corrector (SLC)-off issue, which results in the loss of approximately 22% of the normal scene area for Landsat ETM+ acquisitions after 31 May 2003.However, the SLC-off issue has no impact on the quality of the remaining valid pixels [18].The study area covered only a relatively small portion of the scenes, including their central part, and was therefore not highly affected.More precisely, only 12.6% of the pixels within the analyzed crop parcels had "fill-values" in the images provided by USGS.Those pixels were ignored in further processing steps.

Agricultural Parcels and Crop Classes
Field validated data from agricultural parcels controlled in 2005 in the context of CAP were used in this study to develop and train our ACS.IFAP produced this data by OTS checks through CAPI, followed by farm visits when necessary.Since the Portuguese LPIS in 2005 represented single crop parcels (Figure 7 in [12]), which were refined through field validation, boundaries of single crops were available for this study.The year 2005 was chosen due to the relatively high control rate in Ribatejo in that year, allowing us to analyze a broad range of different crops, as mentioned before.After 2005, CAP's Single Payment Scheme, which decoupled subsidies from the production of specific crops, was progressively adopted by Portuguese farmers, with the consequence of reducing the diversity of declared crops.
A total of 32,062 agricultural parcels lying within the study area were initially provided by IFAP.Only pixels lying entirely within a parcel were considered in order to avoid mixed pixels at the boundary.
A subset of the initial parcels was then selected for further analysis based on two criteria.Firstly, all parcels had to contain at least one entire pixel.Secondly, only parcels containing crop classes accounting for approximately 95% of the total parcel area were included in the subset.In detail, 12 classes (Table 2) covered 94.4% of the total parcel area of 111,963 ha, with the remaining 5.6% being covered by 48 classes.This resulted in a total of 11,852 parcels used for this work, representing 105,702 ha of the study area.Of those parcels, 68.5% were not at all affected by the SLC-off problem of Landsat ETM+.Note that, strictly speaking, our 11,852 parcels do not reflect with certainty the most cultivated crops in the region, but the most controlled crops in 2005 in the study area.

Methods
In simple terms, the proposed ACS is a classifier coupled with a reliability requirement.In Section 3.1, we describe the techniques used for variable selection.The ACS itself is divided into three steps: classification, which is done by an SVM classifier using 10-fold cross-validation (Section 3.2); calibration, which is at the core of our approach and is described in Section 3.3; and application in an operational context (Section 3.4).Finally, in Section 3.5, we explain how we assess the proposed ACS. Figure 2 shows the three processing steps, as well as their data inputs and outputs.All data processing steps were performed using the freely-available software environment R, Version 3.0.2[19].Step 3: Application Step 2: Calibration

Variable Selection
Each single parcel is described by the set of pixel reflectances in its interior, and each pixel is observed in several spectral bands and dates.Therefore, the dimensionality is very high, which means that feature selection is required.This is done at two levels.Firstly, we discuss how the set of pixel measurements within the parcel can be reduced.Then, we discuss how to select the best set of combinations of spectral bands and acquisition dates from the 36 possible combinations (6 bands and 6 dates).Dimensionality reduction is a common preprocessing step in supervised classification, since it prevents possible redundancy in the datasets and overfitting [20].Moreover, dimensionality reduction may permit reducing the amount of data (fewer acquisition dates) and computational effort.
To address the first issue of within-parcel measurements, we compared the variability between parcels with the variability within parcels for each combination of band and date.Our goal was to investigate if the variability between parcels was much higher than within parcels, which is expected for parcel-based classification where parcels are homogeneous.Since our population includes several crops, we decomposed the total variability, for all pixels within the 11,852 parcels, as in a two-way nested ANOVA, where crops are the main factor and parcels are the nested factor.This allows us to isolate the crop effect and to compute the quotient: between-parcel mean square within-parcel mean square .
which measures only the effect of the parcels on the overall reflectance variability.A high F -value suggests that single parcels are homogeneous and that variability among pixel reflectances is mostly due to differences between parcels.Moreover, RMSE = within-parcel mean square (2) measures the standard deviation within parcels.A low RMSE means that the set of all pixel reflectances within each parcel can be replaced, without much information loss, by its simple average.We stress that this analysis is simply descriptive, since parcels were not selected according to a sampling design.
To address the second level of dimensionality reduction, we performed a variable selection based on Principal Component Analysis (PCA) over the initial 36 variables (X 1 ).Our goal was to select the best set of original variables (the best pairs of bands and dates).However, standard PCA returns a set of orthogonal linear combinations of the variables.To obtain a set of original variables from PCA, we applied an iterative technique based on Jolliffe [21], which allows us to exclude in each step the least important variable, until linear dependencies within the dataset start to vanish.As a result, a subset X 2 of variables was generated.Since PCA only explores linear combinations of variables, it does not guarantee that a non-linear classifier will always perform better over X 2 than X 1 .Therefore, we also compared the accuracy of classification over X 1 and X 2 .

Classification Step
Parcel-based crop classification was carried out using the SVM classifier with X 1 and X 2 as training data.SVM have been successfully used in remote sensing applications [22,23].This classifier is particularly attractive due to its ability to successfully handle small training datasets and for being less susceptible to problems of overfitting than other methods [24,25].The Radial Basis Function (RBF) was used as the kernel function, because it has fewer parameter values to define and has been found at least as robust as other kernel types for remote sensing applications [24].The SVM parameters C and γ were set by performing an optimum parameter search using 10-fold cross-validation.
In the classification process, a classifier assigns an object to a class based on known input variables describing the object, x.This is the classification decision.We also consider the posterior probabilities derived from SVM according to [26].Those are the probabilities P (ω j |x) of the true class being ω j given x [27].For example, P ("wheat"|x) denotes the estimated posterior probability of a parcel being a wheat crop when its signature is x.For the i-th parcel, p i denotes the posterior probability for the most likely class.Therefore, {p 1 , . . ., p k } are the posterior probabilities of the k classification decisions, with k being the number of parcels in the training dataset.
An important concern in remote sensing applications is to quantify the agreement between the performed classification and ground truth data by performing accuracy assessment [6].Accuracy assessment was undertaken using the error matrix approach [28].Robust estimates of the error matrix, Overall Accuracy (OA), Producer's Accuracy (PA) and User's Accuracy (UA) for each dataset were obtained using 10-fold cross-validation.Note that PA is related to the commonly-used omission error (error of exclusion), defined as 1−PA, while UA is related to commission error (error of inclusion), which is 1 − UA.In summary, the classification step features the following inputs and outputs: Input: dataset X; Outputs: classifier, error matrix, maximum posterior probabilities {p 1 , . . ., p k }.

Calibration Step
As outlined before, it is crucial to take into account the accuracy of automatic crop classification in the operational context of CAP subsidy controls, in order to build a system that produces reliable results.The rationale behind our approach is that attempting to automatically classify all parcels will inevitably lead to poor classification accuracies, at least for some classes.Our ACS allows us to control the commission errors of crop classification by using a reliability level λ.The reliability level is a user-defined lower bound on user's accuracy for all classes and can be described as an "overall UA".Some classification results are more reliable than others, as revealed by their respective posterior probabilities P (ω j |x).The calibration of the control system is done by excluding parcels that were classified with low reliability, i.e., with low posterior probabilities, until the user's accuracy in each class matches the user-defined reliability level.As a result, UA ≥ λ is guaranteed over all classes for the selected subset of parcels.
In practice, after λ is chosen, the calibration step determines a minimum posterior probability for each class, {q 1 , . . ., q c }, above which a classification is considered reliable (with c being the number of classes).Generally speaking, the higher the value of q j , the more difficult it is to make a reliable decision regarding class j.In Section 4, we discuss why UA, and not PA, is the correct criterion for the calibration step.The calibration step requires and produces the following inputs and outputs: Inputs: reliability level λ, error matrix, maximum posterior probabilities {p 1 , . . ., p k }; Output: required posterior probabilities {q 1 , . . ., q c }.

Application Step
The ACS is designed to be used by a control agency in an operational context in the application step once it was successfully calibrated.It is applied to each new parcel with some signature x derived from an atmospherically-corrected multispectral imagery time series and the parcel's location.
In the application step, a given parcel is classified using the trained SVM classifier.The ACS then accepts or rejects that classification decision, with the chosen reliability level, according to the following decision rule: where P (ω j |x) is the posterior probability and q j is the lower bound of the posterior probability for class j defined in the calibration step.Accepting a classification decision means that the corresponding parcel can be controlled automatically by the system, while rejecting it means that the parcel does not meet the desired reliability.Rejected parcels need to be subjected to CAPI.The following summarizes the application step for a given parcel: Inputs: classifier, required posterior probabilities {q 1 , . . ., q c }, signature x; Output: accept or reject classification decision for parcel with signature x.

Automatic Classification Proportion
For the purpose of assessing the ACS calibrated with a given reliability level, accuracy statistics OA, PA and UA were estimated only for the parcels that could be classified automatically.Moreover, to measure the degree of applicability of our ACS, i.e., the proportion of parcels that can be reliably classified, we defined the Automatic Classification Proportion (ACP) for each class as follows: No. parcels classified in class j with P (ω j |x) ≥ q j No. parcels classified in class j the overall ACP is computed by simply dividing the total number of automatically-classified parcels by the total number of parcels.

Variable Selection
As described in Section 3, we performed a variance decomposition to investigate if pixel reflectances within each parcel and for each combination of ETM+ bands and acquisition dates could be replaced by the average reflectances within the parcels.Towards that end, we determined F parcels in Equation (1) for each band and date combination.The F -values were always high (between 25.8 and 62.8) for the combinations we used for classification, as shown in Table 3, which indicates that variability within parcels is at least 25-times lower than variability between parcels.In addition, RMSE values show that the overall standard deviation of pixel reflectance values with respect to the parcels' average range between 0.011 and 0.037 units of reflectance, i.e., just 1.1% to 3.7%, for the combination of bands and dates we used for classification.Those results suggest that replacing the whole set of reflectance values within a parcel by the parcel's average reflectance can be done without losing relevant information.Table 3. F parcels (left columns) and RMSE (right columns) statistics for 12 selected ETM+ band/date combinations for dimensionality reduction.F parcels compares between-parcel with within-parcel variability after removing the crop effect, while RMSE measures the within-parcels standard deviation.To select the best combinations of ETM+ bands and acquisition dates, we performed the PCA-based analysis described in Section 3. As a result, we obtained twelve combinations X 2 out of the 36 original combinations X 1 .The subset X 2 is also described in Table 3. From a remote sensing perspective, the results of the selection of the best combinations of spectral bands and acquisition dates are not surprising, since the Near Infrared (NIR) spectral region (Band 4) was always selected, and the "red" band was chosen for three distinct dates.

SVM Classification
Classification was performed SVM with X 1 using C = 1.5 and γ = 0.4 and with X 2 using C = 2.5 and γ = 0.4.This indicates that with reduced dimensionality, we need to increase the penalty for misclassification errors in order to achieve good classification results.Analyzing the effect of dimensionality reduction, we found that the classifications yielded OA values of 68.4% and 68.11% for X 1 and X 2 , respectively.Given these almost identical results and the aforementioned advantages of dimensionality reduction, we considered that X 2 was the best option for the next steps, in particular calibration and application.
Generally speaking, OA results were rather unsatisfactory with both datasets.This confirms that attempting to classify all parcels leads to poor results.Specifically, our results show that classification OA was significantly lower than accuracies reported by other authors also performing object-based crop classification.Peña et al. [23] yielded 88% accuracy of SVM classification of nine different crops.Duro et al. [22] reported an OA of 94% with six different classes, also using the SVM classifier.Conrad et al. [7] classified six crops with an OA of 80.1%.Castillejo-González et al. [6] reported an OA of 90.7 with 10 land cover classes, mostly crops.Our classification accuracy is lower due to the fact that some of our classes are spectrally very similar, as can be seen in the spectral signatures in Figure 3.The reason for this behavior is that some classes represent very similar land cover, such as forage crops, fallow, poor grassland and non used area, and are therefore very difficult to discriminate.Cereal crops (wheat, barley and oat) also show a similar spectral pattern.However, those classes are the ones that the Portuguese NCPA currently identifies by CAPI, and therefore, the real problem at hand is intrinsically difficult.As we will see, our ACS will still lead to high accuracies.Error bars are standard deviations.The y-axis shows parcel reflectance, and the x-axis shows acquisition dates.Some classes, such as permanent grassland, poor grassland and non used area, reveal similar signatures and are therefore difficult to discriminate.Other classes, such as maize, rice and vineyard, show quite unique behaviors and are thus expected to be easily classified.

Effect of the Reliability Level on ACP
In this section, we analyze the effect of the reliability level on the proportion of agriculture parcels that can be classified automatically.This information is crucial to make a the reliability level to adopt for automatic classification.The ACS was calibrated with λ ranging from 50% to 100%, with a 5% step, using the classification results with X 2 .ACP and overall ACP were estimated for each reliability level.
The relationship between reliability level and overall ACP is shown in Figure 4.As expected, overall ACP decreases as the reliability level increases.This happens because higher reliability levels result in higher q j values, which according to the system's decision rule will inevitably lead to less accepted parcel classification decisions.With a low reliability level of 50%, almost all parcels can be classified automatically (97.6%).However, this is not acceptable from the perspective of subsidy control.On the other hand, if a 100% accurate classification is demanded, only 5.5% of all parcels can be classified in an automatic way.This would in fact be more accurate than CAPI, since photo interpretation is known not to be in general 100% reliable [28].An interesting reliability level is 80%, which leads to classification errors below 20% and classifies more than 50% of all parcels in an automatic way.The relationship between reliability level and ACP for each individual land cover class is described in Figure 5.The general trend visible in Figure 4 is obviously also present here: higher λ values lead to lower ACP values.In particular, more than 85% of maize (MAI), rice (RIC), wheat (WHE) and vineyard (VYA) parcels can be controlled using λ = 80%.

Accuracy Assessment of the ACS
For the purpose of describing and assessing ACS results for all classes, we considered a reliability level of 80%.Table 4 compares error matrices and accuracy statistics for crop classification of all parcels and classification of parcels with λ = 80%.The transition from the first to the second error matrix provides a good idea of how the calibration step works.The parcels that were removed from the first error matrix are parcels with posterior probability below the required q j determined by the calibration step.A summary of q j and ACP values for λ = 80% is provided in Table 5.
The ACS allowed improving the classification OA from 68.1% to 84.1%.However, this improvement comes with the cost of reducing the proportion of parcels that can be automatically controlled to only 55.4%.The crop rice yielded the best classification results.Other authors have also found this crop to be the best performing, attributing this efficiency to the fact that this crop grows in flooded fields, which are very distinguishable due to the effect of water in the NIR and SWIR spectral regions [10].maize also presents high accuracy and 100% of automatically-classified parcels (Table 5).permanent grassland (PGL) shows a very high PA of 98.5%, but a UA equal to the minimum required value of 80%, which is a result of the difficulty to discriminate this class from other classes, as revealed by the error matrix.Our best performing crops maize, rice and vineyard are among the most accurate classes reported by other studies, such as the work of Peña Barragán et al. [10], showing that our results are consistent with the literature.poor grassland (POG) and non used area (NUA) have completely accurate classifications in terms of UA, but a more careful analysis of these results shows that only one and two parcels were classified as poor grassland and non used area, respectively.Furthermore, PA values of those classes are near-zero, 0.8% and 2.2%, respectively, revealing that the classifier is not able to deal with those classes in a satisfactory manner.Results from Table 5 further sustain the conclusion that poor grassland and non used area are very hard to identify, as suggested by the small proportion of automatically-classified parcels (ACP values of respectively 2.3% and 1.5%).The situation in fallow is even worse, with no parcel being automatically assigned to this class.Based on this discussion, we recommend that no automatic decisions should be done in classes fallow, poor grassland and non used area.They actually represent very similar land cover types on the ground, constituted by an ill-defined mixture of different vegetation types or even bare soil.
Figure 6 shows both automatically-controlled and rejected parcels in the study area.This kind of map can be useful for NCPAs to carry out field work in the future, since it shows areas where more rejections that may require farm visits occur.As expected from the information in Tables 2 and 5, PGL occupies a vast amount of the automatically-classified area.Figure 6.Agricultural parcels that were automatically identified with λ = 80% (solid filled parcels) and parcels that were rejected by the ACS (outlined parcels).The Tagus River and its tributaries are represented in the background.
We also investigated how overall ACP varies when fewer dates are used.Towards that end, we ran the classification and calibration steps with a reliability level λ = 80% excluding each one of the six dates at the time.The overall ACP decreased from 55.4% to 51.7% when November was excluded and to only 25.1% when July was excluded, with other date removals leading to values close to 50%.This is consistent with the fact that maize and rice have a very distinct SWIR (Band 5) response in July, as shown in Figure 4. We repeated this procedure to determine what was the second least important date.We concluded that excluding May would lead to an additional reduction of the overall ACP to 47.9%.Since both May and June belong to the same optimal period used by IFAP (see Table 1), this suggests that excluding one of those dates results in a reduced loss of performance.
It is important to emphasize the usage of UA, rather than PA, to define the reliability level λ.The goal of ACS is to replace, at least partially, CAPI in assigning a crop to each parcel.To measure how well ACS replaces CAPI, one can use the probability that it returns the correct class for each classified parcel, which is precisely what UA estimates.For example, even if just one out of 130 poor grassland parcels are correctly classified, the confidence in ACS to replace CAPI for that class is still high, since no parcels classified as poor grassland by the ACS were misclassified.Since the remaining 129 parcels were classified in the permanent grassland class, the concern is if ACS is reliable enough to replace CAPI for that class.This can be addressed by increasing λ, with the drawback of reducing ACP.Although PA is important to evaluate a classifier, our ACS does not use it explicitly.However, and since high PA values tend to correspond to high ACP values (see Tables 4 and 5), ACP gives a complementary measure to UA of the performance of the classifier.

An Application Example
A practical application example is provided to illustrate how ACS works in the operational context.Classification results with two different reliability levels are discussed in this example: 80% and 95%. Figure 7 shows false color composites of three selected parcels with known land cover.The corresponding classification results can be found in Table 6.6).
Table 6.Automatic classification results for the parcels represented in Figure 7. q j values and corresponding results are with λ = 80% and λ = 95%.Parcel (a) shows a homogeneous maize cropping with a very high spectral response.This led to a clear classification of the parcel as maize with a high posterior probability P (ω j |x) of 0.95.This result was accepted for both reliability levels (λ = 80%, λ = 95%), and the parcel was classified automatically, since P (ω j |x) was higher than both q j (80%) and q j (95%).Parcel (b) is covered by permanent grassland and shows a rather heterogeneous crop pattern.In the classification process, it was assigned to the incorrect class forage crops with posterior probability 0.647.However, the ACS correctly rejected that classification and decided not to classify the parcel in an automatic way, for both reliability levels.For 80% reliability, it was a close rejection, given that P (ω j |x) almost matched the minimum of 0.686.As for 95%, the classification is clearly rejected, which makes sense, since demanding a higher reliability in the classification decisions increases the threshold to accept classifications, ultimately resulting in more rejected parcels.This parcel exemplifies the usefulness of the ACS, rejecting parcels that cannot be classified with the required certainty.The third example is Parcel (c), which corresponds to rice, revealing a very heterogeneous behavior.This parcel was also assigned to the wrong class by the classifier, in this case to maize, with an associated probability of 0.389.The low q j for maize allows for misclassifications when parcels are classified as maize, even with relatively low posterior probabilities, which happened in this case.The classification was wrongly accepted by the system with a reliability level of 80% and correctly rejected with a reliability of 95%.This example clearly shows the positive effect a higher reliability level has on the accuracy of an automatic classification decision, with the drawback of automatically classifying a smaller proportion of the total number of parcels.

Conclusions
Every year, Common Agriculture Policy area-based subsidies are controlled using remote sensing imagery.This process, known by EU authorities as control with remote sensing, is done over each national Land Parcel Identification System.We used Portuguese official land parcels and an atmospherically-corrected Landsat 7 ETM+ time series to automatize crop identification for subsidy control.However, automatic classification of all parcels, as done in Blaes et al. [4], Matikainen et al. [13], Oesterle and Hahn [14], leads in our case to an overall accuracy of only 68%, which is explained by the similarities of the spectral characteristics of the crops relevant to agriculture subsidies.To address this problem, we coupled crop classification to a reliability requirement.This allows the decision-maker to set the reliability of the classifier, restricting automatic crop identification to parcels that are classified with high certainty.In this paper, we quantify the accuracy of the proposed approach and analyze the trade-off between the reliability level and the proportion of parcels that can be automatically controlled.Moreover, we estimate the proportion of land parcels, for 12 common agricultural occupations, that can be automatically classified.
When the reliability level increases from 50% to 100%, the overall proportion of parcels that are classified automatically decreases from 97.6% to 5.5%.In particular, our results suggest that with a reliability level of 80%, traditional photo interpretation of remote sensing imagery can be replaced by automatic crop identification for more than half of the parcels, with this proportion raising above 80% for maize, rice, wheat or vineyard.Furthermore, maize and rice were found to be easily identified, allowing us to control 100% of parcels with these classes in an automatic way.On the contrary, no parcels containing fallow, poor grassland and non-used area should be controlled due to extreme confusion between these classes.In our experiments, overall accuracy raised from 68% to 84% for that same 80% reliability level.Therefore, we have shown that automatic classification of agricultural land parcels can be reliably performed even when crops are difficult to discriminate.

Figure 1 .
Figure 1.(A) Study area and Landsat 7 scene location over southern Portugal (1: Lisbon District, 2: Santarém District).(B) ETM+ true color composite of the whole study area at 8 July 2005, with data loss due to the SLC-off issue being clearly noticeable.Yellow areas are the 11,852 parcels used in this work, covering approximately 1057 km 2 of the study area.

Figure 2 .
Figure 2. Data processing workflow for the proposed Automatic Control System (ACS) (ACS).Gray boxes indicate inputs/outputs, and white boxes are processing steps.

Figure 4 .
Figure 4. Overall proportion of agriculture parcels that can be classified automatically with different classification reliability level values.At a low value of λ = 50%, almost all parcels are classified in an automatic way.On the other hand, a completely accurate classification only applies to 5.5% of classified parcels.A trade-off between reliability and the number of automatic classifications with respective time and money savings has to be undertaken.

Figure 7 .
Figure 7. False color composite (R = NIR, G = red, B = green) on 8 July 2005 of three selected parcels with different known classes for exemplification of the system's application step (see Table6).

Table 1 .
[15]mal periods for image acquisition according to Carmona[15]used by IFAP and acquisition dates of Landsat ETM+ images selected for this study (date formats: dd/mm and dd/mm/yyyy).

Table 2 .
Summary of the used crop classes.More than half of the area is covered by permanent grassland.Note the great variability in parcel size found in the data (measured as the standard deviation of the parcel area).
Proportion of agriculture parcels that can be classified automatically with different classification reliability level values in each class.Note the very distinct behavior of the different classes.Many classes show a near 100% ACP for λ = 60% and 70%, while other classes show low classification proportions, even at low reliability levels.

Table 4 .
Error matrices and accuracy statistics UA, PA and OA for crop classification: all parcels (top) and parcels with reliability level λ = 80% (bottom).Note that the chosen reliability level is reflected in the UA ≥ 80% yielded in each class in the second error matrix.

Table 5 .
Estimates of q j and ACP, both with a reliability level λ = 80%.