AGRICULTURALLY CONSISTENT MAPPING OF SMALLHOLDER FARMING SYSTEMS USING REMOTE SENSING AND SPATIAL MODELLING

Smallholder agriculture provides 90 % of primary food production in developing countries. Its mapping is thus a key element for national food security. Remote sensing is widely used for crop mapping, but it is underperforming for smallholder agriculture due to several constraints like small field size, fragmented landscape, highly variable cropping practices or cloudy conditions. In this study, we developed an original approach combining remote sensing and spatial modelling to improve crop type mapping in complex agricultural landscapes. The spatial dynamics are modelled using Ocelet, a domain-specific language based on interaction graphs. The method combines high spatial resolution satellite imagery (Spot 6/7, to characterize the landscape structure through image segmentation), high revisit frequency time series (Sentinel 2, Landsat 8, to monitor the land dynamic processes), and spatiotemporal rules (STrules, to express the strategies and practices of local farmers). The method includes three steps. First, each crop type is defined by a set of general STrules from which a model-based map of crop distribution probability is obtained. Second, a preliminary crop type map is produced using satellite image processing based on a combined Random Forest (RF) and Object Based Image Analysis (OBIA) classification scheme, after which each geographical object is labelled with the class membership probabilities. Finally, the STrules are applied in the model to identify objects with classes locally incompatible with known farming constraints and strategies. The result is a map of the spatial distribution of crop type mapping errors (omission or commission) that are subsequently corrected through the joint use of spatiotemporal rules and RF class membership probabilities. Combining remote sensing and spatial modelling thus provides a viable way to better characterize and monitor complex agricultural systems.


INTRODUCTION
The ESA (European Spatial Agency) Sentinel 2 (S2) mission provides free imagery with an unprecedented combination of high spatial (10-60m), temporal (5 days) and spectral (13 bands) resolutions, having promising applications in land use/cover characterization and crop monitoring (Drusch et al., 2012). A study carried out on 12 contrasted agricultural sites around the world showed that image time series acquired at a decametric resolution are adapted to field crop characterization for the majority of intensive agricultural systems . However, methods usually used in northern countries are not efficient in the context of smallholder tropical agriculture. In such complex systems, the small size of fields and the landscape fragmentation are significant constraints to crop characterization using remote sensing. To overcome these limitations, studies proposed OBIA (Objects Based Image Analysis) approaches to improve crop identification (Peña-Barragán et al. 2011), or combined with high-resolution time series, very high resolution imagery and auxiliary spatialized data (i.e. Digital Elevation Model; Lebourgeois et al., 2017). However, limitations remain for a crop characterization in smallholder agriculture areas.
To improve current approaches, several studies combine remote sensing with other disciplines. For example, the number of papers where remote sensing is coupled with neural networks or deep learning is increasing. They integrate satellite data, and environmental or socio-economic data (called hereafter context data) as inputs into the network, in order to extract essential variables for each class. However, these researches have mainly focused on cultivated areas with relatively large fields, a large * Corresponding author number of ground truth points, and have not, or only slightly, been applied to smallholder agricultural areas.
Only a few studies focused on the combination of remote sensing and expert rules or spatial modelling. Bailly et al. (2018) highlight that taking into account the spatio-temporal structures of fields improves small objects classification. Two main approaches can be used for modelling agro-systems. The first one consists in using expert knowledge to create crop rotation models (e.g. Dury et al. 2012). To do so, it is necessary to understand the economic, social and environmental stakes of these agroecosystems and therefore to rely on experts in these disciplines. The second consists in using machine learning methods (such as deep learning and neural network) to free from the limitations of expert knowledge. The direct integration of context data into the classifier allows the latter to learn automatically from the ground truth data. Nevertheless, this approach is particularly suitable to address one issue at a time, such as for crop rotation in Bailly et al. (2018) and with a large and diversified dataset.
To go beyond traditional approaches in remote sensing, this research aims at investigating the complementarity of remote sensing and spatio-temporal modelling to study: • How can spatial modelling help converting remote sensing data into more reliable thematic products? • How can remote sensing provide evidence of the landscape processes expressed as formalized knowledge in models?
Such an approach is expected to be useful for a better understanding of complex fragmented and heterogeneous landscapes, such as those shaped by smallholder tropical agriculture in developing countries. To illustrate our method, we describe its application for the characterization of smallholder agriculture in Madagascar Highlands. This agricultural landscape is an interesting case study because of its high complexity that led to mapping difficulties when remote sensing approaches were used alone, as shown in Inglada et al., (2015) and Lebourgeois et al. (2017).
The approach is divided into three main parts: the first one concerns remote sensing data processing (producing a first land use map), the second one focuses on the spatio-temporal modelling (based on the implementation of spatio-temporal rules expressing farmers' strategies and practices), and the last one on the combination of remote sensing and spatio-temporal modelling (for obtaining a more agriculturally consistent map).

Study Area
Our study area is 20 km x 30 km in size and is part of the JECAM network (Joint Experiment for Crop Assessment and Monitoring, www.jecam.org). It is located, close to Antsirabe, the administrative center of Vakinankaratra (coordinates: 19°31'56" to 19°51'32"S,46°57'24" to 47°13'14"E) in Madagascar Highlands. This area is characterized by complex heterogeneous and predominantly agricultural landscapes. Within the cropland, the complexity is mainly due to the very low average area cultivated per household (ranging from 0.1 to 0.5 ha, according to Sourisseau et al. (2016)). There is also a large variety of crop types with different crop calendars, and the practice of intercropping is widespread. Moreover, the landscape is fragmented by the presence of natural vegetation between the fields that is phenologically synchronized with the crops. The primary growing season occurs during the rainy season, from November to April, constraining the acquisition of optical cloudfree satellite images. The crops are diversified, but rice and maize dominate. They are cultivated under rainfed (in slopes and plateaus) or irrigated conditions (mainly lowland rice paddies).

Data
For this study, we used several types of data: (i) both high spatial resolution and high revisit frequency time series satellite data, (ii) ground truth data, (iii) agro-environmental and socio-economic knowledge data on the farmers' strategies (for the spatiotemporal rules definition), and (iv) auxiliary data (for the rules implementation).

Remote sensing data and pre-processing:
The remote sensing dataset corresponds to (i) one very high spatial resolution SPOT 6 image, and (ii) a high spatial resolution (10m-30m) time series of 37 Sentinel 2 and 9 Landsat 8 images, for the period between October 1 st , 2017 and April 29 th , 2018, corresponding to the main growing season.
The SPOT 6 image was acquired on February 14 th 2017 around the peak of the growing season. The image was ortho-rectified, converted to Top Of Atmosphere reflectance and pansharpened using the open-source Orfeo ToolBox (OTB; Grizonnet et al., 2019) to obtain multispectral images at 1.5 m spatial resolution. Sentinel-2 image time series (Donadieu, L'Helguen, 2016) were acquired from THEIA website (https://theia.cnes.fr). These images are provided in surface reflectance, and with a cloud mask, thanks to the MAJA (MACCS ATCOR Joint Algorithm) processing chain (Hagolle et al., 2017). Then, additional preprocessing steps were performed using the MORINGA chain (Gaetano et al., 2019). The functions used in MORINGA include: a) the automatic resampling (with bicubic interpolation) of S2 20 m bands to 10 m resolution, b) a multitemporal gap-filling (Inglada et al., 2017), and c) a co-registration of S2 and L8 time series on the SPOT 6 image used as a reference.

2.2.2
In situ data on land use: A ground truth database on actual land use was built during field surveys around the peak growing season (beginning of February to end of March 2018). GPS waypoints were recorded following an opportunistic sampling approach. Waypoints were chosen according to their accessibility (but some parts of the area could not be reached due to heavy rains affecting the road conditions), and the size of the object (at least 20 * 20 m, so that they can be identified on satellite images). We recorded land uses according to a modified version of the JECAM nomenclature (see http://jecam.org/wpcontent/uploads/2018/10/JECAM_Guidelines_for_Field_Data_ Collection_v1_0.pdf) to better take into account the specificities of tropical agriculture. In the case of intercropping, we recorded the main and the second crop in the field. Once the waypoints were acquired, the boundaries of each field were digitized on the SPOT 6 image, and the class labels (and other attributes) were joined to the polygon database. Additional non-cropland polygons were also added from photo-interpretation of the SPOT 6 image for easily recognizable classes such as built-up surfaces, water bodies, wetlands, mineral soils, or natural forests. At the end, 2170 polygons were obtained (1759 crop and 511 non-crop).

Knowledge on farmer strategies and practices:
Our approach aimed at integrating agronomic, environmental or socio-economic expert knowledge that express farmer strategies and practices. To do so, three sources of information were considered: (i) documents (research papers, books, reports, technical documents, territorial analysis, etc.) describing the agricultural functioning of the region, (ii) interviews of a set of experts with different backgrounds (agronomy, geography or socio-economy), (iii) semi-directive interviews of farmers, local officials or persons involved in rural development. These last interviews were conducted at the Fokontany (municipality/village) scale using a questionnaire. The questions of the surveys concerned the local agricultural practices (crop calendars, cropping constraints and spatial distribution), and the socio-economic context that has an impact on farmer strategies.

Auxiliary Data:
These include existing spatial information such as altitude, slope and crop calendars. The full list of auxiliary data used is given in Table. The topographic position index (TPI) was used to generate the local and global topography using a multi-scale grid from the digital elevation model (Weiss, 2001) . The hydrological network was modelled with GRASS GIS (GRASS Development Team, 2017) to obtain an active network for the rainy season.

METHODOLOGY
The overall approach involves remote sensing data processing, spatio-temporal modelling and their combination (Figure 1). It is divided into five steps: (i) A preliminary RS-based crop type map is produced using satellite image processing based on a combined OBIA and Random Forest (RF) classification scheme. Each object is labelled with its corresponding class membership probabilities, i.e. the probability that the object pertains to each of the classes (hereafter referred to as CMVPsat, for class membership vector of probabilities based on satellite image processing). (ii) Each crop type is defined by a set of general spatio-temporal rules (STrules) that are formalized from the analysis of knowledge gathered on farmer strategies and practices.
(iii) The STrules are implemented in a modelling platform to produce a model-based map of crop distribution probability (CMVPrule for STrules-based class membership probabilities) of each object. (iv) Possible misclassified objects are detected following specific rules. (v) The membership probability of the selected objects is recalculated (to obtain a corrected class membership vector probabilities, CMVPcorr) to produce a corrected and more agriculturally consistent land use map (Figure 1).

Remote sensing
The preliminary land use map was produced from remote sensing. This involved the following steps: (i) segmentation of the very high spatial resolution (VHSR) image into homogeneous objects, (ii) extraction of spectral and textural features from satellite images and other products (DEM), and (iii) training and validation of a Random Forest classifier which was then applied to the whole study area. For this study, all steps were performed using the MORINGA processing chain (Gaetano et al. 2019), which uses OTB and Gdal functions driven by Python scripts. Figure 1: Framework of the method. The processes are presented in red color. CMVPsat corresponds to class membership vector probabilities obtained after random forest (RF) classification of satellite images. CMVPrules refers to the same vector obtained after spatio-temporal rules (STrules) implementation in the model. Finally, CMVPcorr presents the redefinition of the class membership vector probabilities from the two previous ones using fusion methods.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-3/W11, 2020 PECORA 21/ISRSE 38 Joint Meeting, 6-11 October 2019, Baltimore, Maryland, USA 3.1.1 OBIA: According to the OBIA approach (Blaschke et al., 2014), the segmentation of the VHSR image aimed at producing objects delineating the boundaries of agricultural fields, or at least groups of homogeneous fields to characterize the landscape structure. This segmentation determines the landscape structure of the study area, i.e. the boundaries of the "entities" that will be used during the complete workflow. The Baatz and Schäpe segmentation algorithm (Baatz and Schäpe, 2000) implemented in OTB was used with the following parameters: 0.7 for the homogeneity criteria, 0.5 for the spectral homogeneity and the threshold was fixed at 200.

Feature extraction:
Satellite features were extracted for each entity of the ground truth database to build a learning database. This step involved the extraction of the following features: (i) reflectance from the Sentinel-2 and Landsat-8 time series (471 features), (ii) spectral indices from the same time series (267 features), (iii) textural indices from the Spot 6 panchromatic image (6 features), and (iv) altitude and slope from the digital elevation model (2 features). Those features were chosen according to a previous study by Lebourgeois et al. (2017). Concerning the spectral indices, we computed the most common ones: NDVI (Rouse et al., 1974), NDWI (Gao, 1996), MNDWI (Xu, 2006), MNDVI (Jurgens, 1997), Brightness Index, and RNDVI (Jorge et al., 2019) for each image of Sentinel-2 and Landsat 8. Haralick textures (Haralick et al. 1973) were computed using mean and correlation for internal textures and contrast as edge texture. Two windows were chosen (7 pixels and 14 pixels) to traduce the diversity of fields organisation???. SRTM-30m ("SRTM V3.0, 1 arcsec") (NASA JPL, 2014) was used to extract altitude and to compute slope.

Classification:
In the final step of the remote sensing approach, the objective was to obtain a land use map and to extract and a class membership vector of probabilities for each object (CMVPsat). This vector is useful for understanding the probabilities distribution between classes, and in the following steps, for redefining class membership probabilities for selected objects. The Random forest (RF) algorithm (Breiman, 2001) is widely used in remote sensing (e.g. Lebourgeois et al., 2017;Debats et al., 2016) and has been shown to outperform other classifiers in benchmark studies (e.g. Fernández-Delgado et al., 2014;Inglada et al., 2015). The RF classifier was used to produce a land-use classification based on the learning database at the end of the growing season. The main output, the class membership vector of probabilities (CMVPsat), was extracted for each object.

Spatio-temporal modelling
The model was built with the Ocelet Modelling Platform (OMP, www.ocelet.org). OMP allows the spatial modelling of a system and their dynamics (Degenne and Lo Seen, 2016) using both raster and vector data. Ocelet is based on the concept of interaction graphs to describe how entities composing a system (e.g. fields, rivers, weather, municipality) are inter-connected. Each entity is characterized by a number of properties that can evolve through time with the application of interaction functions.

Formalization and selection of STrules:
The knowledge of farmer strategies and practices gathered during literature review and interviews was formalized so that crop type can be described by a set of general spatio-temporal rules (STrules). The spatial rules can concern distances or intersection with specific entities (such as roads, villages, markets and hydrographic network) or zones (such as soil map units, altitude and slope classes). The temporal rules can be linked to the cropping calendars or rainfall regime during the growing season, or the previous crop rotation on the same field. The spatiotemporal rules for each crop type were selected according to the following criteria: (i) the rule should be representative of the whole study area, (ii) the rule should be frequently mentioned by the persons interviewed or in the documents, (iii) the rule must be consistent with agronomic knowledge and observable in the field, and (iv) the data must be available for the rule to be applied. A crop type can be defined by exclusive rules (e.g. maximum distance from a point) or optional rules (e.g. generally in fertile soil). The assignment of an entity to a given class must comply with all exclusive rules and part (or none) of the optional rules.

Implementation:
In the model, the landscape structure of the system is based on the segmentation produced during step 3.1.1, which provide the boundaries of the entities. For each crop, rules are then implemented in the OMP at the entity level to produce a map giving the probability of presence of the corresponding crop in the study area. From this step onwards, the application of all the rules will define a class membership vector of probabilities for each entity, based on formalized knowledge on farmer strategies (referred to as CMVPrule for STrules-based membership probabilities).

Redefinition of class membership probabilities
Traditional evaluation methods used in remote sensing to validate a classification are based on the computation of indices such as overall accuracy, Kappa, F-score (Labatut, Cherifi, 2012) or the quantity and allocation disagreement (Pontius, Millones, 2011). These methods can be derived from confusion matrices, giving the number of well-classified and misclassified ground truth occurrences per class. These methods allow evaluating the results of a given classifier, but only ground truth information. Other entities are not evaluated. This step proposes a spatial evaluation of the remote sensing classification with respect to farmer strategies (expressed through STrules in the model), and to redefine class membership probabilities to obtain more consistent agricultural maps.
When comparing the CMVPrule to CMVPsat for the whole study area, a map can be generated highlighting areas or entities where the rules disagree, and for which corrections are required. For a given entity, in order to decide which of the two probability vectors is more reliable, other information such as the number of satellite acquisition for an object during the growing season, the set of expert rules concerned, or the confidence value obtained during the random forest vote are used. Then, to redefine class membership, various fusion methods can be used such as fuzzy logic, average between CMVPsat and CMVPrule, or rules weighting.

RESULTS
In the following sections, we give examples of STrules used and present three main outputs of the approach. Based on the application of STrules for each entity in the model, probabilities of presence of each crop type were mapped. This map was combined with the map produced during the remote sensing step to create a disagreement map. The final output of the method was a new land use map with an associated class membership probability for each object.

Spatio-temporal rules
Literature review and interviews allowed determining a set of STrules for each crop. The type of dominant rules can differ between crops. For example, some crops are strongly dependent on economic rules (e.g. good connection to markets) while other crops are more constrained to a given biophysical environment or linked to a specific crop calendar.
Examples of different exclusive and optional rules are given below for three crop classes of the study area: Market gardening: Exclusive rules: R1. Less than 30 minutes by truck or less than 15 km from markets (the presence of market gardening is strongly linked to the presence of markets). R2. Distance from villages less than 150 m (to avoid steals because market gardening is a high-value crop).
Optional rules: R3. The field is usually in fertile soil (sandy soil, recent alluvial soil). R4. The field could be irrigated or close to a watercourse. R5. The slope can be between 0% and 7%.

Irrigated rice:
Exclusive rules: R1. The field must be in an irrigated lowland or in the riverbed. R2. The parcel must be close to a watercourse. R3. The crop calendar spreads from 11/21 to 12/20 (for transplantation phase) and harvest occurs from 04/20 to 05/10. Optional rules: R4. The field is frequently in a hydromorphic soil.
Optional rules: R3. The crop is frequently in intercropping with rainfed rice or soya (to be implemented, this rule thus relies on the probabilities of presence of rainfed rice and soya) R4. Maize is frequently present in the livestock sector, which corresponds to slopes or plateaus.

Map of probabilities
The maps of probabilities are the result of STrules Implementation in the model. Figure 2 shows the irrigated rice probability map in the study area. For each crop, the same map of probability of presence can be computed based on STrules as defined during step 3.3.1. For irrigated rice, this map highlights high probabilities of presence around watersheds and lowlands. In the center of the map, the Ambohibary lowland mainly consists in rice and market gardening crops. On plateaus and slopes, the probability of the presence of irrigated rice is low as expressed by STrules.
Maps of probabilities allow to analyze and understand the possible spatial distribution of the crop in the landscape, and to identify the main underlying explanatory factors of this repartition. Figure 2. Map of probability for irrigated rice from spatiotemporal modelling (STrules implementation).

Disagreement map
A comparison between CMVPsat and CMVPrule (see 3.4) allowed producing a disagreement map. This map highlights the entities where remote sensing and STrules seem to be in disagreement. These disagreements mean that either the expert rules or the classifier are not accurate for the concerned entities .
Assuming that the STrules have been verified by experts and were found to be consistent, this map thus allows the spatialization of remote sensing image classification errors.
The spatialization of the classification errors can also help to identify production basins or specific zones where remote sensing has more difficulties to discriminate between crop types. These maps support the analysis of the possible factors limiting the success of remote sensing approach in such context. In our study case, the disagreement map presented in Figure 3 shows remote sensing classification errors for maize, irrigated rice, rainfed rice and market gardening.. Majority of irrigated rice is well classified with a predominance in lowlands, but a part of market gardening and other "minor" crops (e.g. cassava, taro, etc.) in lowlands is misclassified. For these minor crops, the weakness of the learning set, leading to an overlearning of the RF classification algorithm can explain the disagreement . In the North-East part of the map, an important area is misclassified. This inaccessible zone had no ground truth data (used in the learning of the RF classification algorithm) to represent its land use variability, and few auxiliary data to support STrules. Consequently, few entities reached an agreement in this zone. A photo-interpretation of this place showed that, except near watercourses, a majority of natural vegetation and rainfed crops were present. In the South-East of the map, another large disagreement area is identified. In this place, ground truth and auxiliary data were present in sufficient quantity. Disagreements are related to the important misclassification of market gardening and orchards with other rainfed crops. An analysis of CMVPsat showed that all these crop types had a very close membership probability. . From this first result, an analysis of the type and characteristics of the entities where a disagreement is found (the type of crop assigned according to CMVPsat and CMVPrule, and values of probability of presence of each crop for the same entity) can be performed to improve the understanding of processes and the characterization of such complex system.

Corrected map
The main output of the method is a new modified land use map, more consistent with agricultural strategies and practices of the region. An example of a corrected map, produced after the redefinition of probability values by averaging CMVPsat and CMVPrule for each entity in disagreement, is presented in Figure  4. In Madagascar highlands, a first comparison with ground truth data showed that this method allowed better discrimination of rainfed crop types. This crop type is challenging to characterize due to the very small field sizes and the high level of fragmentation of the rainfed cropland area, as shown previously in Lebourgeois et al. (2017). Figure 4. Zoom on land use maps obtained using the remote sensing approach, before (a) and after (b) the redefinition of CMVP values of each entity for which rule-based and remote sensing-based class assignment differs for main crops.

CONCLUSION AND OUTLOOKS
This work proposes a new framework in which remote sensing and expert knowledge on farmer practices and strategies are used in a spatio-temporal modelling approach to produce more agriculturally consistent crop type maps in complex agricultural systems. The method has the originality of combining remote sensing and modelling with knowledge from human and social sciences that reflect farmers' behaviour. Such approach can help to improve (i) the understanding of processes driving land use distribution in complex agricultural areas such as in smallholder farming systems and (ii) the accuracy of crop type maps in regions where traditional remote sensing approaches present limitations.
This approach is also being currently tested in the same study area with a new set of satellite images with higher spatial and temporal resolutions (Pléiades instead of SPOT6/7 for the very high spatial coverage, and addition of Venμs data in complement of the Landsat-8 / Sentinel-2 time series), building on the same set of STrules. An analysis of the model sensitivity to the integration of each spatio-temporal rule will also be performed. Beyond its application to the complex smallholder agricultural landscape of Madagascar Highlands, the approach proposes a generic framework that explicitly brings rich thematic expert knowledge into the analysis of satellite time series data. Such a framework can certainly be adapted to other regions of interest or thematic applications (e.g. biodiversity). However, in a broader perspective, this method can be considered as a first attempt to combine knowledge-driven and data-driven approaches for extracting useful information from one of the most abundant sources of geographical data that are earth observation satellites.