Knauer, K., Gessner, U., Fensholt, R., et Kuenzer, C. (2016). An ESTARFM fusion framework for the generation of large-scale time series in cloud-prone and heterogeneous landscapes. *Remote Sensing, 8*(5). http://doi.org/10.3390/rs8050425

L

Labatut, V., et Cherifi, H. (2012). Accuracy measures for the comparison of classifiers. In *The 5th International Conference on Information Technology*. Amman, Jordanie.

LADA. (2008). *Mapping Land Use Systems at global and regional scales for Land Degradation Assessment Analysis*. Nachtergaele F. et Petri M. LADA Technical report n.8, version 1.1. Rome, Italy.

Lam, N. S. N. (2008). Fractals and Scale in Environmental Assessment and Monitoring. In *Scale and Geographic Inquiry: Nature, Society, and Method* (pp. 23–40). http://doi.org/10.1002/9780470999141.ch2

Lambert, M. J., Waldner, F., et Defourny, P. (2016). Cropland mapping over Sahelian and Sudanian agrosystems: A Knowledge-based approach using PROBA-V time series at 100-m. *Remote Sensing, 8*(3). http://doi.org/10.3390/rs8030232

Lambin, E. F., Geist, H. J., et Lepers, E. (2003). Dynamics of land-use and land-cover change in tropical regions. *Annual Review of Environment and Resources, 28*(1), 205–241. http://doi.org/10.1146/annurev.energy.28.050302.105459

Lambin, E. F., Geist, H. J., et Rindfuss, R. R. (2006). Introduction: Local Processes with Global Impacts. *Land-Use and Land-Cover Change, Local Processes and Global Impact*, (Turner), 2-8. http://doi.org/10.1007/3-540-32202-7_1

Lassalle, P., Inglada, J., Michel, J., Grizonnet, M., et Malik, J. (2015). A Scalable Tile-Based Framework for Region-Merging Segmentation. *IEEE Transactions on Geoscience and Remote Sensing, 53*(10), 5473–5485. http://doi.org/10.1109/TGRS.2015.2422848

Leenhardt, D., Angevin, F., Biarnès, A., Colbach, N., et Mignolet, C. (2010). Describing and locating cropping systems on a regional scale. A review. *Agronomy for Sustainable Development*. http://doi.org/10.1051/agro/2009002

Lelong, C., Alexandre, C., et Dupuy, S. (2014). Discrimination of tropical agroforestry systems in very high resolution satellite imagery using object-based hierarchical classification: A case-study in Cameroon. *South-Eastern European Journal of Earth Observation and Geomatics*.

Leutner, B., et Horning, N. (2017). RStoolbox: tools for remote sensing data analysis. R Package version 0.1.8.

Ligges, U., Short, T., et Kienzle, P. (2015). signal: Signal processing. R package version 0.7.6. Disponible à l'adresse : https://cran.r-project.org/package=signal
Lima, A. A. C., Oliveira, F. N. S., et de Aquino, A. R. L. (2000). *Solos e aptidão agrícola das terras do Estado do Tocantins. Embrapa Agroindústria Tropical*. *Documentos*, **31**. Fortaleza.

Lloyd, S. P. (1982). Least Squares Quantization in PCM. *IEEE Transactions on Information Theory*, **28**(2), 129–137. http://doi.org/10.1109/TIT.1982.1056489

Lobell, D. B., et Field, C. B. (2007). Global scale climate–crop yield relationships and the impacts of recent warming. *Environmental Research Letters*, **2**(1), 7pp. http://doi.org/10.1088/1748-9326/2/1/014002

Loveland, T. R., et DeFries, R. S. (2004). Observing and monitoring land use and land cover change. *Ecosystems And Land Use Change*, **153**, 231–246.

M

MacQueen, J. (1967). Some methods for classification and analysis of multivariate observations. In *Proceedings of the fifth Berkeley symposium on mathematical statistics and probability* (Vol. 1, pp. 281–297).

Malingreau, J. P., Bartholomé, E., et Barisano, E. (1987). Surveillance de la production agricole en Afrique de l'Ouest. Nécessité d'une intégration de différentes plates-formes satellitaires. In *Proceedings of the Symposium SPOT 1, utilisation des images, bilans, résultats* (pp. 353–370). Paris.

Massey, R., Sankey, T. T., Congalton, R. G., Yadav, K., Thenkabail, P. S., Ozdogan, M., et Sánchez Meador, A. J. (2017). MODIS phenology-derived, multi-year distribution of conterminous U.S. crop types. *Remote Sensing of Environment*, **198**, 490–503.
http://doi.org/10.1016/j.rse.2017.06.033

Mathys, E., et Gardner, A. (2009). USAID Office of Food for Peace Burkina Faso Food Security Country Framework FY 2010 – FY 2014. Washington, D.C: Food and Nutrition Technical Assistance II Project (FANTA-2), Academy for Educational Development (AED).

Matson, P. A., Parton, W. J., Power, A. G., et Swift, M. J. (1997). Agricultural Intensification and Ecosystem Properties. *Science*, **277**(5325), 504–509.
http://doi.org/10.1126/science.277.5325.504

Maus, V., Câmara, G., Cartaxo, R., Sanchez, A., Ramos, F. M., et De Queiroz, G. R. (2016). A Time-Weighted Dynamic Time Warping Method for Land-Use and Land-Cover Mapping. *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, **9**(8), 3729–3739. http://doi.org/10.1109/JSTARS.2016.2517118

McNairn, H., Champagne, C., Shang, J., Holmstrom, D., et Reichert, G. (2009). Integration of optical and Synthetic Aperture Radar (SAR) imagery for delivering operational annual crop inventories. *ISPRS Journal of Photogrammetry and Remote Sensing*, **64**(5), 434–449. http://doi.org/10.1016/j.isprsjprs.2008.07.006
McNemar, Q. (1947). Note on the sampling error of the difference between correlated proportions or percentages. *Psychometrika*, *12*(2), 153-157.

Mingwei, Z., Qingbo, Z., Zhongxin, C., Jia, L., Yong, Z., et Chongfa, C. (2008). Crop discrimination in Northern China with double cropping systems using Fourier analysis of time-series MODIS data. *International Journal of Applied Earth Observation and Geoinformation*, *10*(4), 476-485.
http://doi.org/10.1016/j.jag.2007.11.002

Ministère de l’éducation nationale et ministère de l’enseignement supérieur et de la recherche. Arrêté du 28 novembre 1980 relatif à la terminologie de la télédétection arospatiale. Journal Officiel numéro complémentaire du 11 décembre 1980, p. 10848 (1980).

Monfreda, C., Ramankutty, N., et Foley, J. A. (2008). Farming the planet: 2. Geographic distribution of crop areas, yields, physiological types, and net primary production in the year 2000. *Global Biogeochemical Cycles*, *22*(1). http://doi.org/10.1029/2007GB002947

Mücher, S., De Simone, L., Kramer, H., de Wit, A., Roupioz, L., Hazeu, G., ... Cormont, A. (2016). A new Global Agro-Environmental Stratification (GAES). *Wageningen Environmental Research. Report 2761*, 70. http://doi.org/10.18174/400815

Mücher, S., De Simone, L., Kramer, H., de Wit, A., Roupioz, L., Hazeu, G., ... Cormont, A. (2016). A new Global Agro-Environmental Stratification (GAES). http://doi.org/10.18174/400815

Müller, C., Bondeau, A., Popp, A., Waha, K., et Fader, M. (2010). *Climate Change Impacts on Agricultural Yields*. Washington D.C. © World Bank: License: CC BY 3.0 IGO.

N

Nachtergaele, F., et Petri, M. (2013). *Mapping Land Use Systems at Global and Regional Scales for Land Degradation Assessment Analysis*. Lada.

NASA LP DAAC. (s.d.-a). MOD09A1 data product. NASA EOSDIS Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center.

NASA LP DAAC. (s.d.-b). MOD13Q1 version 5 data product. NASA EOSDIS Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center. Disponible à l'adresse : https://lpdaac.usgs.gov/

Nir, D. (1987). Regional geography considered from the systems' approach. *Geoforum*, *18*(2), 187-202. http://doi.org/10.1016/0016-7185(87)90005-4
O

O'Neill, R. V, DeAngelis, D. L., Waide, J. B., et Allen, T. F. H. (1986). *Hierarchical Concept of Ecosystems. Monographs in Population Biology Series*. Vol.23. Consulté à l'adresse : http://search.barnesandnoble.com/booksearch/isbninquiry.asp?r=1etean=0691084378

OCHA ROWCA. (2015). Burkina Faso: Settlements. United Nations Office for the Coordination of Humanitarian Affairs (OCHA) Regional Office for West and Central Africa (ROWCA). Disponible à l'adresse : https://data.humdata.org/dataset/settlements

OCHA ROWCA. (2017). Second Administrative Level Boundaries Project (SALB), Institut Géographique du Burkina Faso (IGB). United Nations Office for the Coordination of Humanitarian Affairs (OCHA) Regional Office for West and Central Africa (ROWCA). Disponible à l'adresse : https://data.humdata.org/dataset/burkina-faso-administrative-boundaries

OSS. (2015). *Burkina Faso : Atlas des cartes d'occupation du sol. Projet amélioration de la résilience des populations sahéliennes aux mutations environnementales - REPSAHEL*. Observatoire du Sahara et du Sahel (OSS).

Ozdogan, M., et Woodcock, C. E. (2006). Resolution dependent errors in remote sensing of cultivated areas. *Remote Sensing of Environment*, 103(2), 203–217.
http://doi.org/10.1016/j.rse.2006.04.004

P

Pacheco, A. M., McNairn, H., et Merzouki, A. (2010). Evaluating TerraSAR-X for the identification of tillage occurrence over an agricultural area in Canada. *Remote Sensing for Agriculture, Ecosystems, and Hydrology Xii*, 7824, 1-7. http://doi.org/78240p\r10.1117/12.868218

Pacheco, A., et McNairn, H. (2010). Evaluating multispectral remote sensing and spectral unmixing analysis for crop residue mapping. *Remote Sensing of Environment*, 114(10), 2219-2228.
http://doi.org/10.1016/j.rse.2010.04.024

Panigrahy, S., Manjunath, K. R., et Ray, S. S. (2005). Deriving cropping system performance indices using remote sensing data and GIS. *International Journal of Remote Sensing*, 26(12), 2595-2606. http://doi.org/10.1080/01431160500114698

Panigrahy, S., Ray, S. S., Manjunath, K. R., Pandey, P. S., Sharma, S. K., Sood, A., ... Parihar, J. S. (2011). A Spatial Database of Cropping System and its Characteristics to Aid Climate Change Impact Assessment Studies. *Journal of the Indian Society of Remote Sensing*, 39(3), 355-364. http://doi.org/10.1007/s12524-011-0093-3

Parmentier, S. (2014). *Scaling-up agroecological approaches: what, why and how*. Oxfam-Solidarity, Brussels.
Pearson, K. (1901). Principal components analysis. *The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science*, **6**(2), 559.

Pedroso da Silva, A. R., et Almeida, M. G. (2007). O agronegócio e o Estado do Tocantins: o atual estágio de consolidação. *Caminhos de Geografia*, **8**(21).

Plet, F. (2003). La geographie rurale française : quelques jalons. *Sociétés Contemporaines*, **49-50**(1), 85-106. http://doi.org/10.3917/soco.049.0085

Plexida, S. G., Sfougaris, A. I., Ispikoudis, I. P., et Papanastasis, V. P. (2014). Selecting landscape metrics as indicators of spatial heterogeneity-Acomparison among Greek landscapes. *International Journal of Applied Earth Observation and Geoinformation*, **26**(1), 26-35. http://doi.org/10.1016/j.jag.2013.05.001

Plourde, J. D., Pijanowski, B. C., et Pekin, B. K. (2013). Evidence for increased monoculture cropping in the Central United States. *Agriculture, Ecosystems and Environment*, **165**, 50-59.
http://doi.org/10.1016/j.agee.2012.11.011

Pontius, R. G., et Millones, M. (2011). Death to Kappa: birth of quantity disagreement and allocation disagreement for accuracy assessment. *International Journal of Remote Sensing*, 32(15), 4407-4429. http://doi.org/10.1080/01431161.2011.552923

Q

Qiu, B., Zhong, M., Tang, Z., et Wang, C. (2014). A new methodology to map double-cropping croplands based on continuous wavelet transform. *International Journal of Applied Earth Observation and Geoinformation*, **26**(1), 97-104. http://doi.org/10.1016/j.jag.2013.05.016

R

R Core Team. (2016). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

Ramoino, F., Pera, F., et Arino, O. (2018). “S2 prototype LC map at 20m of Africa 2016” *Users Feedback Compendium* [6th February 2018]. Rome, Italy. Consulté à l'adresse :
http://due.esrin.esa.int/files/S2_prototype_LC_map_at_20m_of_Africa_2016-Users_Feedback_Compendium-6-Feb-2018.pdf

Ribeiro, P. F., Santos, J. L., Santana, J., Reino, L., Leitão, P. J., Beja, P., et Moreira, F. (2016). Landscape makers and landscape takers: links between farming systems and landscape patterns along an intensification gradient. *Landscape Ecology*, **31**(4), 791-803.
http://doi.org/10.1007/s10980-015-0287-0

Richards, J. a, et Jia, X. (2006). *Remote Sensing Digital Image Analysis: An Introduction*. New York. http://doi.org/10.1007/978-3-642-30062-2
Rizzo, D., Marraccini, E., Lardon, S., Rapey, H., Debolini, M., Benoît, M., et Thenail, C. (2013). Farming systems designing landscapes: Land management units at the interface between agronomy and geography. *Geografisk Tidsskrift*, **113**(2), 71–86.
http://doi.org/10.1080/00167223.2013.849391

Robbez-Masson, J. M., Foltête, J. C., Cabello, L., et Flitti, M. (1999). Prise en compte du contexte spatial dans l'instrumentation de la notion de paysage. Application à une segmentation géographique assistée. *Revue Internationale de Géomatique*, **9**, 173-195.

Rodrigues de Melo, K. D., Polastrine, A., Amaral, R. M., Grando Lopes, A., Sayão Lobato, C. L. D., et Rosanova, C. (2010). Caracterização da bovinocultura leiteira nos municípios de Taguatinga, Combinado e Lavandeira, região sudeste do Estado do Tocantins, Brasil. In *XX CONGRESSO BRASILEIRO DE ZOOTECNIA* (p. 3). Palmas, Tocantins.

Rosenfield, G. H., et Fitzpatrick-Lins, K. (1986). A coefficient of agreement as a measure of thematic classification accuracy. *Photogrammetric Engineering and Remote Sensing*, **52**(2), 223-227.

Rouse, J. W., Haas, R. H., Schell, J. A., et Deering, D. W. (1973). Monitoring Vegetation Systems in the Great Okains with ERTS. *Third Earth Resources Technology Satellite-1 Symposium*, **1**, 325-333. http://doi.org/10/citeulike-article-id:12009708

Rudorff, B., Risso, J., Aguiar, D., Gonçalves, F., Salgado, M., Perrut, J., ... Virtuoso, M. (2015). *Análise Geoespacial da Dinâmica das Culturas Anuais no Bioma Cerrado: 2000 a 2014*.

Rufino, M. C., Atzberger, C., Baldi, G., Butterbach-Bahl, K., Rosenstock, T. S., et Stern, D. (2016). Targeting landscapes to identify mitigation options in smallholder agriculture. In *Methods for Measuring Greenhouse Gas Balances and Evaluating Mitigation Options in Smallholder Agriculture* (pp. 15–36). http://doi.org/10.1007/978-3-319-29794-1_2

S

Sahajpal, R., Zhang, X., Izaurralde, R. C., Gelfand, I., et Hurtt, G. C. (2014). Identifying representative crop rotation patterns and grassland loss in the US Western Corn Belt. *Computers and Electronics in Agriculture*, **108**, 173–182.
http://doi.org/10.1016/j.compag.2014.08.005

Salmon, J. M., Friedl, M. A., Frolking, S., Wisser, D., et Douglas, E. M. (2015). Global rain-fed, irrigated, and paddy croplands: A new high resolution map derived from remote sensing, crop inventories and climate data. *International Journal of Applied Earth Observation and Geoinformation*, **38**, 321–334. http://doi.org/10.1016/j.jag.2015.01.014

Sauer, S., et Leite, S. P. (2012). Agrarian structure, foreign investment in land, and land prices in Brazil. *Journal of Peasant Studies*, **39**(3-4), 873–898.
http://doi.org/10.1080/03066150.2012.686492
Savitzky, A., et Golay, M. J. E. (1964). Smoothing and Differentiation of Data by Simplified Least Squares Procedures. *Analytical Chemistry*, **36**(8), 1627–1639.
http://doi.org/10.1021/ac60214a047

Sayre, N. F. (2005). Ecological and geographical scale: parallels and potential for integration. *Progress in Human Geography*, **29**(3), 276–290. http://doi.org/10.1191/0309132505ph546oa

Sebillotte, M. (1982). Les systèmes de culture. Réflexion sur l’intérêt et l’emploi de cette notion à partir de l’expérience acquise en région de grande culture. In *Séminaire du département d’agronomie de l’INRA* (pp. 63–80). INRA, Vichy.

Sebillotte, M. (1990). *Système de culture, un concept opératoire pour les agronomes*. (L. Combe et D. Picard, Eds.). INRA, Paris, 165-196.

See, L., Fritz, S., You, L., Ramankutty, N., Herrero, M., Justice, C., ... Obersteiner, M. (2015). Improved global cropland data as an essential ingredient for food security. *Global Food Security*. http://doi.org/10.1016/j.gfs.2014.10.004

SEPLAN. (2012). Base de Dados Geográficos do Tocantins - Atualização 2012. Palmas: Secretaria do Planejamento (SEPLAN), Superintendência de Planejamento e Gestão Central de Políticas Públicas, Diretoria de Zoneamento Ecológico-Econômico. Disponible à l'adresse : http://seplan.to.gov.br/zoneamento/bases-vetoriais/base-de-dados-geograficos-do-tocantins-atualizacao-2012/

Siebert, S., Hoogeveen, J., et Frenken, K. (2006). Irrigation in Africa, Europe and Latin America - Update of the Digital Global Map of Irrigation Areas to Version 4. *Frankfurt Hydrology Paper*, (5), 135p.

Siebert, S., Kummu, M., Porkka, M., Döll, P., Ramankutty, N., et Scanlon, B. (2015). A global data set of the extent of irrigated land from 1900 to 2005. *Hydrology and Earth System Sciences*, **19**(3), 1521–1545. http://doi.org/10.5194/hess-19-1521-2015

Silva, L. A. G. C. (2007). *Biomas presentes no Estado do Tocantins*. Consultoria Legislativa da Câmara dos Deputados.

Skakun, S., Kussul, N., Shelestov, A. Y., Lavreniuk, M., et Kussul, O. (2016). Efficiency Assessment of Multitemporal C-Band Radarsat-2 Intensity and Landsat-8 Surface Reflectance Satellite Imagery for Crop Classification in Ukraine. *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, **9**(8), 3712–3719.
http://doi.org/10.1109/JSTARS.2015.2454297

Slonim, N., Aharoni, E., et Crammer, K. (2013). Hartigan's k-means versus Lloyd's k-means - Is it time for a change? In *IJCAI International Joint Conference on Artificial Intelligence* (pp. 1677–1684).

Solano, R., Didan, K., Jacobson, A., et Huete, A. (2010). *MODIS Vegetation Index User's Guide (MOD13 Series)*.
Soudani, K., le Maire, G., Dufrêne, E., François, C., Delpierre, N., Ulrich, E., et Cecchini, S. (2008). Evaluation of the onset of green-up in temperate deciduous broadleaf forests derived from Moderate Resolution Imaging Spectroradiometer (MODIS) data. *Remote Sensing of Environment*, 112(5), 2643-2655. http://doi.org/10.1016/j.rse.2007.12.004

Spangler, K. R., Lynch, A. H., et Spera, S. A. (2017). Precipitation Drivers of Cropping Frequency in the Brazilian Cerrado: Evidence and Implications for Decision-Making. *Weather, Climate, and Society*, 9(2), 201-213. http://doi.org/10.1175/WCAS-D-16-0024.1

Spera, S. (2017). Agricultural Intensification Can Preserve the Brazilian Cerrado: Applying Lessons From Mato Grosso and Goiás to Brazil’s Last Agricultural Frontier. *Tropical Conservation Science*, 10, 194008291772066. http://doi.org/10.1177/1940082917720662

Stomph, T. J., Fresco, L. O., et van Keulen, H. (1994). Land use system evaluation: Concepts and methodology. *Agricultural Systems*, 44(3), 243-255. http://doi.org/10.1016/0308-521X(94)90222-2

Strahler, A. H., Boschetti, L., Foody, G. M., Friedl, M. A., Hansen, M. C., Herold, M., ... Woodcock, C. E. (2006). *Global Land Cover Validation: Recommendations for Evaluation and Accuracy Assessment of Global Land Cover Maps*. Scientific and Technical Research series (Vol. EUR 22156).

Strahler, A. H., Woodcock, C. E., et Smith, J. A. (1986). On the nature of models in remote sensing. *Remote Sensing of Environment*, 20(2), 121-139. http://doi.org/10.1016/0034-4257(86)90018-0

T

Thiombiano, A., et Kampmann, D. (2010). *Atlas de la Biodiversité de L’Afrique de L’Ouest, Tome II: Burkina Faso*. (A. Thiombiano et D. Kampmann, Eds.). Ouagadougou et Frankfurt/Main: Goethe-Universität Frankfurt am Main, Geowissenschaften/Geographie.

Trimble Germany GmbH. (2014). eCognition Developer© 9.0 User Guide. Munich, Germany.

U

United Nations. Department of Economic and Social Affairs. Population Division. (2017). *World Population Prospects: The 2017 Revision, Key Findings and Advance Tables* (No. ESA/P/WP/248).

USGS. (s.d.). Landsat Level-2 Science Products. U.S. Geological Survey Earth Resources Observation and Science Center (USGS/EROS), Science Processing Architecture (ESPA) On Demand Interface. Disponible à l'adresse : https://espa.cr.usgs.gov/
USGS. (2016). *Landsat 8 (L8) Data Users Handbook*. Sioux Falls, South Dakota: U.S. Geological Survey (USGS) Landsat Project Science Office at the Earth Resources Observation and Science (EROS) Center. Consulté à l'adresse : https://landsat.usgs.gov/sites/default/files/documents/Landsat8DataUsersHandbook.pdf

V

Van Wart, J., van Bussel, L. G. J., Wolf, J., Licker, R., Grassini, P., Nelson, A., ... Cassman, K. G. (2013). Use of agro-climatic zones to upscale simulated crop yield potential. *Field Crops Research, 143*, 44–55. http://doi.org/10.1016/j.fcr.2012.11.023

Veldkamp, A., Kok, K., De Koning, G. H. J., Schoorl, J. M., Sonneveld, M. P. W., et Verburg, P. H. (2001). Multi-scale system approaches in agronomic research at the landscape level. *Soil and Tillage Research, 58*(3–4), 129–140. http://doi.org/10.1016/S0167-1987(00)00163-X

Vermote, E. F., El Saleous, N. Z., et Justice, C. O. (2002). Atmospheric correction of MODIS data in the visible to middle infrared: First results. *Remote Sensing of Environment, 83*(1–2), 97-111. http://doi.org/10.1016/S0034-4257(02)00089-5

Vintrou, E., Desbrosse, A., Bégué, A., Traoré, S., Baron, C., et Lo Seen, D. (2012). Crop area mapping in West Africa using landscape stratification of MODIS time series and comparison with existing global land products. *International Journal of Applied Earth Observation and Geoinformation, 14*(1), 83–93. http://doi.org/10.1016/j.jag.2011.06.010

Vintrou, E., Soumaré, M., Bernard, S., Bégué, A., Baron, C., et Lo Seen, D. (2012). Mapping Fragmented Agricultural Systems in the Sudano-Sahelian Environments of Africa Using Random Forest and Ensemble Metrics of Coarse Resolution MODIS Imagery. *Photogrammetric Engineering et Remote Sensing, 78*(8), 839–848. http://doi.org/10.14358/PERS.78.8.839

Viña, A., Peters, A. J., et Ji, L. (2003). Use of multispectral ikoson imagery for discriminating between conventional and conservation agricultural tillage practices. *Photogrammetric Engineering and Remote Sensing*.

W

Waldner, F., Lambert, M. J., Li, W., Weiss, M., Demarez, V., Morin, D., ... Defourny, P. (2015). Land cover and crop type classification along the season based on biophysical variables retrieved from multi-sensor high-resolution time series. *Remote Sensing, 7*(8), 10400–10424. http://doi.org/10.3390/rs70810400

Wang, T., Kou, X., Xiong, Y., Mou, P., Wu, J., et Ge, J. (2010). Temporal and spatial patterns of NDVI and their relationship to precipitation in the Loess Plateau of China. *International Journal of Remote Sensing, 31*(7), 1943–1958. http://doi.org/10.1080/01431160902929263
Wardlow, B. D., et Egbert, S. L. (2008). Large-area crop mapping using time-series MODIS 250m NDVI data: An assessment for the U.S. Central Great Plains. *Remote Sensing of Environment*, **112**(3), 1096–1116. http://doi.org/10.1016/j.rse.2007.07.019

Warrens, M. J. (2015). Properties of the quantity disagreement and the allocation disagreement. *International Journal of Remote Sensing*, **36**(5), 1439–1446. http://doi.org/10.1080/01431161.2015.1011794

Wiens, J. A. (1989). Spatial Scaling in Ecology. *Functional Ecology*, **3**(4), 385. http://doi.org/10.2307/2389612

Willer, H., et Lernoud, J. (2017). *The World of Organic Agriculture. Statistics and Emerging Trends 2017*. Research Institute of Organic Agriculture (FiBL), Frick, and IFOAM - Organics International. Bonn.

Woodcock, C. E., et Strahler, A. H. (1987). The factor of scale in remote sensing. *Remote Sensing of Environment*, **21**(3), 311–332. http://doi.org/10.1016/0034-4257(87)90015-0

World Bank. (2010). *World Development Report 2010 : Development and Climate Change*. Washington D.C. © World Bank: License: CC BY 3.0 IGO.

World Bank. (2017). World Development Indicators (Agriculture, value added, % of GDP). Disponible à l'adresse : https://data.worldbank.org/

World Bank. (2018). Population estimates and projections (Population growth, annual %). Disponible à l'adresse : http://databank.worldbank.org/

X

Xiong, J., Thenkabail, P. S., Gumma, M. K., Teluguntla, P., Poehnelt, J., Congalton, R. G., ... Thau, D. (2017). Automated cropland mapping of continental Africa using Google Earth Engine cloud computing. *ISPRS Journal of Photogrammetry and Remote Sensing*, **126**, 225–244.

Xiong, J., Thenkabail, P. S., Tilton, J. C., Gumma, M. K., Teluguntla, P., Oliphant, A., ... Gorelick, N. (2017). Nominal 30-m cropland extent map of continental Africa by integrating pixel-based and object-based algorithms using Sentinel-2 and Landsat-8 data on google earth engine. *Remote Sensing*, **9**(10). http://doi.org/10.3390/rs9101065

Y

You, L., et Wood, S. (2006). An entropy approach to spatial disaggregation of agricultural production. *Agricultural Systems*, **90**(1–3), 329–347. http://doi.org/10.1016/j.agsy.2006.01.008

You, L., Wood, S., Wood-Sichra, U., et Wu, W. (2014). Generating global crop distribution maps: From census to grid. *Agricultural Systems*, **127**, 53–60. http://doi.org/10.1016/j.agsy.2014.01.002
Young, S. S., et Anyamba, A. (1999). Comparison of NOAA/NASA PAL and NOAA GVI data for vegetation change studies over China. *Photogrammetric Engineering and Remote Sensing, 65*, 679–696.

Z

Zhong, L., Yu, L., Li, X., Hu, L., et Gong, P. (2016). Rapid corn and soybean mapping in US Corn Belt and neighboring areas. *Scientific Reports*, 6. http://doi.org/10.1038/srep36240

Zonneveld, I. S. (1989). The land unit - A fundamental concept in landscape ecology, and its applications. *Landscape Ecology*, 3(2), 67–86. http://doi.org/10.1007/BF00131171
Annexe - Publications

Bellón, B., Bégué, A., Lo Seen, D., de Almeida, C.A., et Simões, M. (2017). "A Remote Sensing Approach for Regional-Scale Mapping of Agricultural Land-Use Systems Based on NDVI Time Series". *Remote Sensing*, 9(6), 600. http://doi.org/10.3390/rs9060600

Bellón, B., Bégué, A., Lo Seen, D., Lebourgeois, V., Evangelista, B. A., Simões, M., et Demonte Ferraz, R.P. (2018). "Improved regional-scale Brazilian cropping systems' mapping based on a semi-automatic object-based clustering approach". *International Journal of Applied Earth Observations and Geoinformation*, 68C, 127-138. http://doi.org/10.1016/j.jag.2018.01.019
<figure><img src="image_70.png" /></figure>

Article

# A Remote Sensing Approach for Regional-Scale Mapping of Agricultural Land-Use Systems Based on NDVI Time Series

Beatriz Bellón<sup>1,*</sup>, Agnès Bégué<sup>1</sup>, Danny Lo Seen<sup>1</sup>, Claudio Aparecido de Almeida<sup>2</sup>
and Margareth Simões<sup>3,4</sup>

Academic Editors: Raul Zurita-Milla, James Campbell and Prasad S. Thenkabail

Received: 2 May 2017; Accepted: 10 June 2017; Published: 13 June 2017

**Abstract:** In response to the need for generic remote sensing tools to support large-scale agricultural monitoring, we present a new approach for regional-scale mapping of agricultural land-use systems (ALUS) based on object-based Normalized Difference Vegetation Index (NDVI) time series analysis. The approach consists of two main steps. First, to obtain relatively homogeneous land units in terms of phenological patterns, a principal component analysis (PCA) is applied to an annual MODIS NDVI time series, and an automatic segmentation is performed on the resulting high-order principal component images. Second, the resulting land units are classified into the crop agriculture domain or the livestock domain based on their land-cover characteristics. The crop agriculture domain land units are further classified into different cropping systems based on the correspondence of their NDVI temporal profiles with the phenological patterns associated with the cropping systems of the study area. A map of the main ALUS of the Brazilian state of Tocantins was produced for the 2013–2014 growing season with the new approach, and a significant coherence was observed between the spatial distribution of the cropping systems in the final ALUS map and in a reference map extracted from the official agricultural statistics of the Brazilian Institute of Geography and Statistics (IBGE). This study shows the potential of remote sensing techniques to provide valuable baseline spatial information for supporting agricultural monitoring and for large-scale land-use systems analysis.

**Keywords:** geographic object-based image analysis (GEOBIA); Moderate Resolution Imaging Spectroradiometer (MODIS); principal components analysis (PCA); cropping systems; Stratification

## 1. Introduction

Ensuring food security through sustainable agricultural development sets a double challenge to agricultural systems: Food production must substantially increase to satisfy the demand of a continuously growing population while improving the stewardship of natural resources and minimizing environmental impacts [1]. This transition towards sustainability is further challenged by the consequences of climate change and an increased competition for land, generating an urgent need

<footnote n="1"><sup>1</sup> Cirad, UMR TETIS (Land, Environment, Remote Sensing and Spatial Information), Maison de la Télédétection, Rue Jean-François Breton, 34090 Montpellier, France; agnes.begue@teledetection.fr (A.B.); loseen@teledetection.fr (D.L.S.)</footnote>

<footnote n="2"><sup>2</sup> National Institute for Space Research (INPE), Av. dos Astronautas, 1758, São José dos Campos, SP 12227-010, Brazil; claudio.almeida@inpe.br</footnote>

<footnote n="3"><sup>3</sup> Embrapa Solos, Rua Jardim Botânico, 1024, Rio de Janeiro, RJ 22460-000, Brazil; margareth.simoes@embrapa.br</footnote>

<footnote n="4"><sup>4</sup> Department of Computer Engineering, Rio de Janeiro State University (UERJ/FEN/DESC/PPGMA), Rua São Francisco Xavier, 524, 5031 D, Maracanã, Rio de Janeiro, RJ 20550-900, Brazil</footnote>

<footnote>* Correspondence: beatriz.bellon@teledetection.fr; Tel.: +33-467-558-645</footnote>
for repetitive spatial information to help monitor the dynamics of the agricultural land-use systems at the regional and global scales.

Land use is defined as “the sequence of operations carried out with the purpose to obtain goods and services from the land” [2]. Land-use systems (LUS) can thus be considered as coupled human–environment systems that can be characterized by two main aspects: land resources and management practices [3]. LUS mapping originally consisted in the delineation of relatively homogeneous areas of land (referred to as land units), which are directly linked to a specific type of land use [4,5]. Regional-scale agricultural land-use systems (ALUS) maps typically represent the major agricultural zones for a given region, each agricultural zone (land unit) being directly linked to a particular agricultural system. The land units of a regional-scale ALUS map are naturally composed of multiple agricultural fields which are spatially “concentrated” and representative of a particular agricultural system, but also contain a mosaic of other land-cover types, e.g., patches of natural vegetation among the agricultural fields, due to their extent. The plotting of the land units’ boundaries is, however, challenged by limited access to spatially explicit and detailed land-use information over large areas, which may account for the general lack of LUS maps.

The few existing large-area LUS mapping approaches, such as the one developed by FAO [2], are usually based on a subjective choice of socio-economic and environmental variables. These are derived from heterogeneous sources which are empirically categorized and combined to map the different land-use systems. The resulting LUS maps are subject to error propagation due to the disparity of spatial resolution, production date, and quality of the original data. In addition, access to the data is not always guaranteed, which limits the reproducibility of the LUS maps in time and across regions. As a result, the description and location of the different agricultural land-use systems remain highly unclear for most world regions.

Earth observation satellite systems can significantly contribute to LUS mapping since they provide timely and detailed land-use information over large areas due to their synoptic coverage and high revisiting frequency. In particular, remote sensing-derived time series of vegetation indices allow monitoring of the phenology (seasonal vegetation variation) and the intra-seasonal variations of the agricultural land cover (i.e., cropland, pastures, and rangelands), from which the agricultural land use (including some management practices) can be determined [6–9]. Furthermore, the geographic object-based image analysis (GEOBIA) [10] segmentation techniques, which automatically delineate homogeneous objects at multiple scales, seem particularly adapted for the extraction of land units. Bisquert et al. [11] introduced an innovative approach based on GEOBIA of spectral and textural variables derived from coarse-resolution vegetation index images that produced a spatial segmentation of the French territory in broad-scale land units. The resulting land units proved to be consistent in terms of land cover [12], but remained unclassified, and the variables used are not descriptive enough for effective delineation of agricultural land-use systems since they do not account for the intra-seasonal variations that characterize agricultural land-use practices.

The aim of this work is to provide a simple and reproducible approach for regional-scale mapping of ALUS based on GEOBIA and vegetation index time series analysis. The presented approach involves the extraction of temporal variations of the vegetation cover in relation to agricultural land use from a vegetation index time series, by performing a principal components analysis and selecting the principal component images that capture the seasonal and intra-seasonal variability of the vegetation index time series. Segmentation and classification of the selected principal components images follow, finally producing an ALUS map. The approach was applied to an annual Moderate Resolution Imaging Spectroradiometer Normalized Difference Vegetation Index (MODIS NDVI) dataset and tested on a new agricultural expansion region of Brazil. The spatial distribution of the cropping systems in the final ALUS map was evaluated with the annual crop estimates of the Municipal Agricultural Production database (PAM) of the Brazilian Institute of Geography and Statistics (IBGE).
## 2. Materials

### 2.1. Study Area

The study area corresponds to the Brazilian State Tocantins, a region that spans 50°W to 45°W and 5°S to 13°S, covering 277,621 km² between the Amazon and the Cerrado biomes (Figure 1). The geographical location between two biomes provides this region with a rich diversity of environmental conditions, which have led to a diversification of human activities and thus a variety of agricultural landscapes [13].

<figure><img src="image_71.png" /><figcaption>**Figure 1.** The study area corresponds to the Tocantins state (TO, red-striped), located in the center-north of Brazil, which is part of the MATOPIBA region (color-coded orange). This region is the latest agricultural frontier of Brazil, composed of parts of the states of Maranhão (MA), Piauí (PI), Bahia (BA), and the full extent of Tocantins.</figcaption></figure>

While livestock production on extensive rangelands still dominates most of the land outside of the protected area's limits, this region has experienced rapid agricultural expansion within the last decade. This trend is mainly characterized by an expansion of the soybean-cultivated area, which is mostly produced for export (with an increase of 804,800 hectares in the last 15 years [14]). Sequential cropping systems (mainly soybean-based double-cropping systems, where soybean is preceded or succeeded by a cereal) are widespread, since the cultural season benefits from a 7-month rainy season (from October to April) and agriculture is mostly mechanized, with field sizes around 100 hectares.

### 2.2. TerraClass Land-Cover Map

Among the efforts to produce datasets which will effectively help monitor the recent land-use transitions in Brazil, the National Institute for Space Research (INPE) (in partnership with Brazilian Agricultural Research Corporation (EMBRAPA)) has contributed to the TerraClass project [15]. TerraClass products include land-cover maps of the Amazon biome for every two years since 2004, together with a first version of the Cerrado biome that was recently released for the year 2013 [16,17].

For this study, we combined the results of the latest TerraClass products (TerraClass Amazônia 2014 and TerraClass Cerrado 2013) to generate a land-cover map at 30 m spatial resolution covering the
full extent of the Tocantins state. The combined product is used to extract an agricultural land-cover map, composed of the land-cover classes' pasture and rangelands, and annual agriculture (Figure 2).

<figure><img src="image_72.png" /><figcaption>**Figure 2.** Agricultural land-cover map extracted from the combination of the TerraClass Amazônia 2014 map (covering the North East of the State of Tocantins) and the TerraClass Cerrado 2013 map (covering the rest of the State of Tocantins) [16,17].</figcaption></figure>

## 2.3. Agricultural Statistics

The official agricultural statistics at the Brazilian municipality administrative level were acquired from the Municipal Agricultural Production database (PAM) of the Brazilian Institute of Geography and Statistics (IBGE) [18]. This database provides annual estimates of planted and harvested area, average yield, average income, and production of the different crops since 1973.

For our study, we collected the harvested area estimates of the 2014 annual crops of the municipalities of Tocantins. The collected data were processed and used as the reference for the spatial distribution of the main cropping systems of Tocantins, against which our results were evaluated.

## 2.4. MODIS NDVI Time Series

Although the NDVI is known to be affected by background, aerosol effects, and saturation in high biomass regions [19], it is highly sensitive to small increases in the amount of photosynthetic vegetation [20] and is therefore well adapted to capture the vegetation dynamics of the agricultural land-use systems.

A time series of 23 16-day composited NDVI images from October (start of the growing season) 2013 to October 2014 (of the MOD13Q1 v5 product at a 250 m spatial resolution) is used in this
study. The dataset was retrieved from the NASA Land Processes Distributed Active Archive Center (LP DAAC) through the processing R package MODISStp Version 1.3.2 [21]. Two MOD13Q1 tiles (h13v09 and h13v10) are needed for covering the entire state of Tocantins.

### 3. Methods

The methods applied in this study follow a typical GEOBIA approach [10] through automatic image partition into meaningful image objects (the land units), and subsequent classification of the resulting objects by an assessment of their characteristics to generate new geographic information (the agricultural land-use systems (ALUS) map). The methodological framework is presented in Figure 3 and described in the following three sub-sections: Pre-processing of the NDVI time series (Section 3.1); Extraction of the land units (Section 3.2); Classification of the land units (Section 3.3). A final sub-section presents the unsupervised evaluation approach that was used to evaluate the cropping systems' classification results of the ALUS map (see Section 3.4. Unsupervised Evaluation of the Classification Results).

<figure><img src="image_73.png" /><figcaption>**Figure 3.** Representation of the methodological framework used in this study.</figcaption></figure>

#### 3.1. Pre-Processing of the NDVI Time Series

The compositing process of the MOD13Q1 product uses the Constrained View angle-Maximum Value Composite (CV-MVC) algorithm, which involves selecting (pixel-wise) the observation with the highest NDVI value and the smallest view angle out of all the observations made during a 16-day cycle [22]. However, remaining effects such as cloudiness or bi-directional effects are still present in the dataset.

In order to generate a noise-reduced time series, an adapted version of the function-fitting smoothing method introduced by Chen et al. [23] was applied. This noise-reduction algorithm is based on a linear interpolation of cloudy pixels, followed by the use of the Savitzky-Golay (SG) filter through an iteration process until the upper envelope of the original NDVI profile is best approached [23,24].

The cloud flag data supplied with the MOD13Q1 product (Pixel Reliability layer) showed systematic overestimations of clouds in the cropland areas, which are particularly evident during the harvest stages in the crop calendar. Consequently, only the pixels followed by an unnatural NDVI increase (greater than 0.4 between two composited images) were considered as cloudy dates and linearly interpolated. In order to reduce the noise of other residual effects such as consecutive cloudiness, the iterative SG filter was applied. After several tests, the parameters *d* (degree of the smoothing polynomial) and *m* (half-width of the smoothing window) of the SG filter (for the long-term change curve and for the fitting iteration) were set to *d* = 2, *m* = 2, and *d* = 6, *m* = 4, respectively. This
smoothing method effectively improved the quality of the time series, reducing signal noise while preserving the detail of the NDVI temporal profiles (Figure 4).

To avoid the temporal discontinuities generated by the compositing process (which affect the shape of the NDVI profiles), the values of the smoothed time series were linearly interpolated. This was done using their actual acquisition dates provided in the "Composite Day of Year (CDOY)" layer of the MOD13Q1 product, to place them at 16-day regular intervals.

<figure><img src="image_74.png" /><figcaption>**Figure 4.** Sample Normalized Difference Vegetation Index (NDVI) time series profile of (a) the raw data; (b) the filtered and smoothed data; and (c) the final temporally-corrected data.</figcaption></figure>

## 3.2. Extraction of the Land Units

The extraction of the land units is developed in two steps which are described in the following sub-sections: Principal Components Analysis (Section 3.2.1) and Delineation of the Land Units (Section 3.2.2).

### 3.2.1. Principal Components Analysis

A principal components analysis (PCA) transformation was applied to the corrected NDVI time series. Originally formulated by Pearson (1901), the PCA technique consists of a linear orthogonal transformation (or rotation) of the original variables into a set of uncorrelated variables known as the principal components (PC). When applied to remote sensing time series, PCA computes new values for each pixel of the original time series' images to produce a set of PC images (also known as eigenchannels). The resulting PC images are uncorrelated with one another and ordered according to the amount of variance retained within the total variance of the original time series dataset (from most to least). Low-order PC images (PC1, PC2) concentrate most of the variance and therefore inevitably capture the information redundancy among the original images, while the higher-order PC images capture less redundant information (associated with specific localized change events).

PCA has proven to be particularly effective (when applied to NDVI time series) at identifying seasonal changes of the land cover and revealing their relationship to different factors—mainly climate variability and human activities factors [25–30]. The type of information captured by each PC image will vary according to the characteristics of the input NDVI time series (temporal resolution and length of the time series, spatial resolution and extent) and the PCA method used (either unstandardized if the PC images are calculated using the covariance matrix, or standardized if the correlation matrix is used instead). However, some common outcomes are observed when the area under analysis is large enough to include different types of land cover: The first PC image captures the major element of variability in the NDVI time series (which for large areas corresponds to the spatial variability of the NDVI values, and is correlated to the integrated NDVI over all seasons), whilst the other PC images
capture the seasonal and intra-seasonal variability of NDVI, each related to a particular factor or a combination of factors [25–30].

We expect the phenological patterns of different agricultural land-use systems (configured by seasonal development of the vegetation and intra-seasonal variability and short-term fluctuations induced by the agricultural practices), to be detectable in these latter PC images. Therefore, after performing a PCA on the corrected NDVI time series of 23 images, we retained the last 22 PC images (the high-order PC images) out of the 23 resulting PC images for further analysis, which contain 27% of the total variance. The first PC image was deliberately discarded after verification that the temporal variability of NDVI was not significantly captured, and that (as expected) the 73% of the total variance captured by the first PC image is mostly related to the spatial variability of the NDVI values, as shown in other studies (e.g., [27,28]).

The PCA was performed by means of the rasterPCA tool in the processing R package RStoolbox Version 0.1.8 [31] and was based on the covariance matrix (unstandardized PCA). Unstandardized PCA is preferred since sudden changes in the magnitude of the NDVI related to land use (such as a rapid NDVI decrease due to crop harvest) will contribute significantly to the development of the new component images. Conversely, normalizing the input images through standardized PCA would weaken the contribution of these sudden change events by forcing all the original input images to have equal weight in the derivation of the PC images.

### 3.2.2. Delineation of the Land Units

The segmentation of the high-order PC images was performed using the iterative region-growing algorithm based on the Baatz and Schäpe criterion, implemented in the Multiresolution Segmentation tool of the object-orientated image analysis software eCognition® Developer 9.0 [32]. This algorithm starts with each pixel of the image forming an image object. At each iteration, neighboring objects are merged based on their similar characteristics described by a homogeneity criterion, until a user-defined threshold (known as the scale parameter) is reached. The multiresolution segmentation algorithm ensures a regular growth of the image objects over the whole scene, and so the resulting objects have a similar size, which is determined by the scale parameter threshold. In this way, a low scale parameter will lead to many objects of small size and, conversely, a high scale parameter will allow the growth of few large objects. The influence of the radiometry versus the shape in the definition of the homogeneity criterion is determined by the weight assigned to the color and shape parameters [33].

The objective of the segmentation being the delineation of land units with relatively homogeneous phenological patterns, the shape of the objects in the definition of the homogeneity criterion becomes irrelevant and, therefore, we assign the total weight to the color by setting the color parameter to 1 and the shape parameter to 0. In practice, setting the color parameter to 1 means that the relative homogeneity of the resulting image objects (land units) will be entirely based on the pixel values of the PC images and their implicit information related to the temporal variability of the vegetation cover. Besides, each of the 22 PC images is made to contribute equally to the segmentation by setting the weight of each input PC image to 1.

A scale (according to Benz et al. [34]) is defined as the level of aggregation and abstraction at which an object can be clearly described. Since our geographic objects of interest (the land units) do not relate to an inherent scale but can rather be described at different scales of analysis, we consider that there is no optimal segmentation scale and that different scale parameters should be tested in each particular case to approach adequately sized, thematically meaningful land units. We assume that the multiresolution segmentation algorithm based on homogeneity criterion delineates relatively homogeneous objects regardless of the scale parameter used. However, some simple guiding principles were followed in this study to choose among different scale parameters aiming at meaningful land units (in terms of relatively homogeneous phenological patterns): Objects should be larger than the minimum mapping unit of an agricultural land-use system (i.e., the field) and there should be a sensible
trade-off between the number of objects for further analysis and their degree of heterogeneity. After testing several segmentation levels, the scale parameter threshold was empirically set to 850 (unitless).

### 3.3. Classification of the Land Units

The classification of the land units is developed in two steps which are described in the following sub-sections: Rule-based Classification (Section 3.3.1) and Classification through Phenological Pattern Analysis (Section 3.3.2).

#### 3.3.1. Rule-Based Classification

The agricultural land-cover map extracted from the TerraClass product was used for the classification of the land units belonging to the agricultural domain (either to the crop agriculture domain or to the livestock domain) into four major ALUS: Dominant crop agriculture system, coexistence (crop-livestock) system, semi-intensive livestock system, and intensive livestock system. This preliminary classification of the ALUS was performed following the decision tree ruleset developed by Almeida et al. [35] for the localization of agricultural production systems in the Brazilian State Rondônia at the municipal level, presented in Figure 5. This ruleset seems well adapted for the classification of the ALUS in the Tocantins state, owing to its closeness to the Rondônia state (in terms of land-cover and land-use systems' characteristics).

<figure><img src="image_75.png" /><figcaption>**Figure 5.** Decision tree ruleset used for the classification of the land units into four major agricultural land-use systems adapted from Almeida et al. [35].</figcaption></figure>

#### 3.3.2. Classification through Phenological Pattern Analysis

The land units belonging to the crop agriculture domain were further classified into different cropping systems based on the analysis of their phenological patterns. The mean temporal NDVI profile of the cropland cover inside each land unit was extracted from the MODIS NDVI time series using the TerraClass annual agriculture class as a mask. The general NDVI variations over the growing season were analyzed through visual interpretation of the shape of the resulting mean temporal NDVI profiles, which revealed the principal crop development stages of the cropland cover inside each land unit. The NDVI profiles with visually equivalent shapes were grouped together and were
finally classified into different cropping systems by correspondence with the crop calendar of the main
cropping systems of Tocantins.

3.4. *Unsupervised Evaluation of the Classification Results*

Since there are no available thematically-equivalent maps with which to validate the accuracy
of the classification of the ALUS, an unsupervised evaluation was carried out with a reference map
extracted from the official agricultural statistics (area estimates of the annual crops in the municipalities
of Tocantins for the year 2014). The evaluation was only carried out for the cropping systems' results
since no official statistics were available on the pasture and rangelands with which to evaluate the
livestock systems' classification results.

In order to produce a synthetic reference map out of the detailed agricultural statistics dataset
that could be compared to the cropping systems in the ALUS map, we first calculated the proportion
of area occupied by each crop type in the municipalities with more than 2000 hectares of cultivated
area. The predominant crop types (accounting for most of the cultivated area) were differentiated from
minor crops and were used to characterize the main cropping system related to each municipality.

The ALUS map was finally evaluated based on the spatial agreement between the classified
cropping systems and the cropping systems in the reference map.

4. Results

4.1. Land Units

The segmentation of the high-order principal component images produced a total of 90 land
units (Figure 6), with an average size of 3090 km², 44 km² for the smallest one, and 15,069 km² for the
largest one.

4.2. Agricultural Land-Use Systems

The decision tree ruleset based on the agricultural land-cover map extracted from the TerraClass
product identified a total of 44 land units belonging to the agricultural domain (114,114 km²) among
the 90 automatically delineated land units. Out of the 44 land units, 21 were classed as dominant crop
agriculture systems (11,193 km²), 15 as semi-intensive livestock systems (82,576 km²), 7 as intensive
livestock systems (18,329 km²) and 1 was classed as a coexistence system (2016 km²) (the latter
corresponding to an area where both a crop agriculture system and a livestock system are present).

Three distinct phenological patterns were identified from the analysis of the shape of the temporal
NDVI profiles of the 22 land units belonging to the crop agriculture domain. The plots of the 22
temporal NDVI profiles are presented in Figure 7 in three groups of phenological patterns.

In the first group a major single crop cycle prevails, with a peak NDVI period around February
(from day 33 to day 65) (Figure 7). Since only minor increases of NDVI appear before and after the
main growing cycle, we can consider that there is no secondary meaningful development cycle of the
vegetation cover, and that, therefore, the land units belonging to the first group are mainly related to a
single cropping system. Since the soybean crop is the dominant crop of Tocantins, the first group of
land units seems distinctive of the soybean single cropping system.

The second and the third groups of land units both show a characteristic double NDVI peak
representing separate crop cycles (Figure 7). Although the first cycle of both groups has a common
peak NDVI period around January (from day 1 to day 33), these two groups are easily dissociable
by their temporal NDVI profile shapes, for their second cycle appears with more than a month
difference between the two groups. The phenological pattern of the second group is distinctive of the
soybean–maize double cropping system, which is characterized by a summer soybean crop during the
rainy season followed by a winter maize crop. This phenological pattern exhibits low NDVI values
around Mars (from day 65 to day 81) as a result of the soybean harvest, after which the winter maize is
usually immediately sown and subsequently harvested during June and July (shown by the low NDVI values after day 161).

<figure><img src="image_76.png" /><figcaption>**Figure 6.** Land units' boundaries over a color composition of the first three principal component (PC) images (RGB PC2, PC3, PC4) of the 22 PC images used as the segmentation variables.</figcaption></figure>

The phenological pattern of the third group is distinctive of the rice–soybean double cropping system, which is the characteristic cropping system of the lowlands located in the southwest of Tocantins in the Formoso river basin. This particular cropping system is characterized by a flood-irrigated rice crop during the rainy season followed by a sub-irrigated soybean crop from May to September (mostly for seed production). The low NDVI values around May (from day 113 to day 145), followed by an increase until the second peak NDVI period is reached around July (from day 193 to day 225—shown in Figure 7), evidences the late second soybean crop cycle, which is unique to the rice–soybean double cropping system.
<figure><img src="image_77.png" /><figcaption>**Figure 7.** The temporal NDVI profiles' plots of the 22 crop agriculture domain land units. The plots show the mean (solid black curve) and the standard deviation (shaded dark grey curve) of the NDVI pixel values inside the annual agriculture TerraClass class for each land unit against the day of the year (DOY). A representative land-unit profile for each of the three phenological pattern groups is shown in detail at the left and the profiles of the other land units belonging to each group are represented assembled at their right side. The "monthly" (32-day) periods (at which the NDVI profiles peak) are represented by shaded light grey vertical bands, and are consistent between the plots of each group.</figcaption></figure>

The crop agriculture domain land units were therefore sub-classified into the three main cropping systems of Tocantins according to the correspondence of their phenological patterns with the characteristic cropping pattern of each cropping system. The map of the main ALUS of Tocantins with the final classification is presented in Figure 8.

As shown in Figure 8, the livestock systems occupy a larger area (about eight times larger) than the cropping systems, and are mainly distributed along a longitudinal transect located in the center-west of the state. The cropping systems are concentrated in specific locations that are sparsely distributed over most of the state—except for the northern third of the state, where the livestock systems are predominant. The livestock systems and the cropping systems are mostly separated, except for the only coexistence area located in the southern region of the livestock-dominant longitudinal transect.
<figure><img src="image_78.png" /><figcaption>**Figure 8.** Map of the main agricultural land-use systems of the Tocantins state in the 2013–2014 growing season.</figcaption></figure>

### 4.3. Unsupervised Evaluation of the Classification Results

The agricultural statistics at the municipal level revealed three predominant crops accounting for more than 94% of the total cultivated area: Soybean (75%), rice (11.1%), and winter maize (8.6%). All of the municipalities with more than 2000 hectares of cultivated area present a high proportion of the soybean crop, with 29–100% of their cultivated area dedicated to its production (which confirms the dominance of this crop type in Tocantins). However, three distinct classes of municipalities were observed. A first class (composed of 15 municipalities)—where the cultivated area is almost exclusively covered by soybean—was related to the soybean single cropping system. A second class (composed of 7 municipalities)—where the proportion of winter maize is considerable—was related to the soybean–maize double cropping system. Finally, a third class (composed of 3 municipalities)—where the proportion of rice is nearly equivalent to the proportion of soybean—was related to the rice–soybean double cropping system. The average proportion of the different crop types for each class of municipalities is presented in Figure 9.
<figure><img src="image_79.png" /><figcaption>**Figure 9.** The average proportion of the crops of the municipalities belonging to the first class (**left**), the second class (**center**), and the third class (**right**), related to the soybean single cropping system, to the soybean–maize double cropping system, and to the rice–soybean double cropping system, respectively.</figcaption></figure>

The resulting classification of the municipalities is represented on a map together with the classification of the cropping systems in the ALUS remote sensing-based map (Figure 10).

<figure><img src="image_80.png" /><figcaption>**Figure 10.** The spatial distribution of the three main cropping systems (a) in the reference map at the municipal level, extracted from the agricultural statistics, and (b) in the agricultural land-use systems remote sensing-based map. The boundaries of the municipalities are included in both maps to facilitate their comparison.</figcaption></figure>

Although the reference map and the ALUS map are produced at different scales, an apparent coherence between both maps can be recognized for the spatial distribution of the soybean single cropping system and the rice–soybean double cropping system (Figure 10). The spatial distribution of the soybean–maize double cropping system shows general coherence between both maps but seems
overestimated locally for some regions, covering some municipalities which are classified as being
dominated by the soybean single cropping system in the reference map.

**5. Discussion**

The study area was partitioned into relatively homogeneous land units that were linked to a specific type of agricultural land use, resulting in a regional level agricultural land-use systems' map. Two livestock systems could be differentiated by the proportion of pasture and rangelands in the land units, and the three main cropping systems of the study area could be differentiated through the proportion of annual cropland in the land units and through their specific phenological patterns (Figure 7); this allowed not only to characterize the land units in terms of the crop type but also in terms of the cropping pattern.

The results showed a general agreement between the location of the three main cropping systems in the final agricultural land-use systems map and the reference map, except for some localized overestimations of the presence of the soybean–maize double cropping system over some of the reference soybean single cropping system municipalities. This overestimation might be due to the non-inclusion of the non-commercial summer grass cover crops such as millet, sorghum, or brachiaria in the agricultural statistics and, therefore, in the reference map. Cover crops are sometimes sown after the winter soybean harvest instead of maize, as part of the widely adopted no-tillage system in the study area. These cover crops might contribute to the second peak of NDVI that is apparent in the NDVI temporal profiles of the land units, which partially intersect municipalities where only the soybean surface area has been officially censused.

The reference map that was used for the evaluation of the classification results shows that land-use
systems' maps can be produced from the agricultural statistics. However, the thematic resolution
of this type of map is limited by the data available in the statistics, thus not accounting for certain
systems (e.g., the livestock systems). Furthermore, the spatial accuracy is greatly limited by the data
aggregation level (e.g., the municipal level)—the administrative boundaries not being representative
of the actual spatial distribution of the land-use systems. The presented remote sensing approach
produces a land-use systems' map which overcomes these limitations, by automatically delineating
boundaries which are representative of the actual physical differences between zones related to the
temporal variations of the vegetation cover.

While the ability of principal components analysis to uncover significant intra-annual change
events out of vegetation index time series has been widely proven (e.g., [25–30]), studies mainly focus
on the low-order principal component (PC) images—generally the first three—which capture most
of the variance. We decided to include the higher-order PC images in our analysis, since they are
highly sensitive to subtle changes in the data and therefore capture change events which are localized
both spatially and in time, such as the intra-seasonal variations and rapid changes induced by the
agricultural practices. Overall, the results confirm that the high-order PC images do capture the
information related to the temporal variability of the vegetation cover, which was the essential portion
of information of the NDVI time series (represented by 27% of the total variance in the time series)
needed for the segmentation of the study site in land units with relatively homogeneous phenological
patterns. Furthermore, the phenological patterns shown in Figure 7 reflected the intra-seasonal
variations which were linked to cropping practices such as the approximate planting and harvesting
dates of different crops.

The automatic extraction of relatively homogeneous land units in terms of phenological patterns,
through a principal components analysis and an automatic segmentation, is a simple and innovative
approach that can be reproduced in other regions. This approach overcomes the potential error
propagation of the expert-based land units' border delineation and the subjective choice and use
of variables in the traditional land-use systems' mapping approaches. Furthermore, the presented
approach is scale-independent and therefore allows multi-scale analysis, which might help improve
the understanding of the land-use systems by characterizing them at different organizational levels.
In addition, the decision tree ruleset used for the classification of four major agricultural land-use
systems (presented in the Figure 5—adapted from Almeida, et al. [35]), is conceived with simple
threshold rules based on the percentage cover of agricultural land, which can be easily adapted
or directly applied, to other agricultural regions. This classification is, however, dependent on a
land-cover map. Even if land-cover maps are now produced and available for many regions, their
accuracy will condition the agricultural land-use systems' classification results, and they may be rarely
updated. Consequently, the rule-based classification step would largely benefit from the inclusion of
an automatic remote sensing-based classification method to produce timely and accurate regional-scale
agricultural land-cover maps at the field level. Some studies have shown that partitioning the study
region into sub-regions improves the accuracy of regional-scale land-cover classifications [36,37].
Hence, the presented land units' extraction approach could therefore also be used as a preliminary step
for automatic remote sensing-based classifications, producing regional-scale land-cover and land-use
maps at the field level.

**6. Conclusions**

Remote sensing-derived time series of NDVI represent an important source of spectral and
temporal information related to agricultural land use and, when combined with GEOBIA techniques
in the present study, allowed the discrimination of different land-use systems (in particular, cropping
systems) at a regional scale. Remote sensing-based approaches, therefore, present a significant potential
for agricultural land-use system mapping over large areas and benefit from the spatial and temporal
continuity of the satellite data and technological advances such as the new high temporal and spatial
resolution Sentinel-2 satellite systems of the Copernicus program.

The presented approach can be potentially reproduced in other regions with minimal adaptation
to specific context, thus contributing to the development of generic and simple tools for the location
and characterization of the agricultural land-use systems across regions. This baseline-detailed spatial
information is highly valuable for further large-scale agronomic and environmental assessments.
If reproduced periodically, this type of approach can help with long-term analysis of the dynamics of
the major agricultural land-use systems, which can help monitor large-scale land-use change processes
such as cropland expansion or intensification (being particularly useful in changing regions such as
the MATOPIBA region, for effective land-use planning leading to sustainable use of the resources).

**Acknowledgments:** Beatriz Bellón received a fellowship from CIRAD (French Agricultural Research Centre for International Development) and from the FP7 SIGMA (Stimulating Innovation for Global Monitoring of Agriculture) project. This work was also supported by the French Space Agency (CNES) through the funding of the TOSCA AGRIZONE project and by the CAPES/COFECUB program through the funding of the GeoABC Project (Methodologies and technological innovation for satellite monitoring of low carbon agriculture in support to Brazil's ABC Plan), project No. 845/15. The authors would like to gratefully acknowledge the members of the NTSA (Thematic Unit on Agricultural Systems) of EMBRAPA (Brazilian Agricultural Research Corporation) Fisheries and Aquaculture for their valuable contribution to this work.

**Author Contributions:** Beatriz Bellón, Agnès Bégué and Danny Lo Seen conceived and designed the experiments; Beatriz Bellón processed the data and performed the proposed methods; Beatriz Bellón, Agnès Bégué and Danny Lo Seen analyzed the results; Beatriz Bellón wrote the paper; and Agnès Bégué, Danny Lo Seen, Claudio Aparecido De Almeida and Margareth Simões supervised the research and contributed to the editing and review of the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

**References**

1. FAO. How to feed the world in 2050. In Proceedings of the Expert Meeting on How to Feed the World in 2050, Rome, Italy, 24–26 June 2009.

2. Nachtergaele, F.; Petri, M. *Mapping Land Use Systems at Global and Regional Scales for Land Degradation Assessment Analysis*; FAO: Rome, Italy, 2013.
3. Bégué, A.; Arvor, D.; Lelong, C.; Vintrou, E.; Simoes, M. Agricultural systems studies using remote sensing. In *Land Resources Monitoring, Modeling, and Mapping with Remote Sensing*; CRC Press: Boca Raton, FL, USA, 2015; pp. 113–130.

4. Driessen, P.M.; Konijn, N.T. *Land-Use Systems Analysis*; Wageningen Agricultural University: Wageningen, The Netherlands, 1992.

5. FAO. *Guidelines for Land-Use Planning*; FAO: Rome, Italy, 1993.

6. Arvor, D.; Jonathan, M.; Meirelles, M.S.P.; Dubreuil, V.; Durieux, L. Classification of MODIS EVI time series for crop mapping in the state of Mato Grosso, Brazil. *Int. J. Remote Sens.* **2011**, *32*, 7847–7871. [CrossRef]

7. Cheema, M.J.M.; Bastiaanssen, W.G.M. Land use and land cover classification in the irrigated Indus Basin using growth phenology information from satellite data to support water management analysis. *Agr. Water Manag.* **2010**, *97*, 1541–1552. [CrossRef]

8. Kiptala, J.K.; Mohamed, Y.; Mul, M.L.; Cheema, M.J.M.; Van Der Zaag, P. Land use and land cover classification using phenological variability from MODIS vegetation in the Upper Pangani River Basin, Eastern Africa. *Phys. Chem. Earth* **2013**, *66*, 112–122. [CrossRef]

9. Wardlow, B.D.; Egbert, S.L. Large-area crop mapping using time-series MODIS 250m NDVI data: An assessment for the U.S. Central Great Plains. *Remote Sens. Environ.* **2008**, *112*, 1096–1116. [CrossRef]

10. Hay, G.J.; Castilla, G. Geographic Object-Based Image Analysis (GEOBIA): A new name for a new discipline. In *Object-Based Image Analysis*; Blaschke, T., Lang, S., Hay, G.J., Eds.; Springer: Berlin/Heidelberg, Germany, 2008; pp. 75–89.

11. Bisquert, M.; Bégué, A.; Deshayes, M. Object-based delineation of homogeneous landscape units at regional scale based on MODIS time series. *Int. J. Appl. Earth Obs.* **2015**, *37*, 72–82. [CrossRef]

12. Bisquert, M.; Bégué, A.; Deshayes, M.; Ducrot, D. Environmental evaluation of MODIS-derived land units. *GIsci. Remote Sens.* **2016**, *54*, 64–77. [CrossRef]

13. Silva, L.A.G.C. *Biomas Presentes no Estado do Tocantins*; Consultoria Legislativa da Câmara dos Deputados: Brasília, Brazil, 2007.

14. CONAB. Soja: Série Histórica de Produtividade (Safras 2000/01 a 2015/16). Available online: http://www.conab.gov.br/ (accessed on 15 April 2017).

15. Almeida, C.A.; Coutinho, A.C.; Esquerdo, J.C.D.M.; Adami, M.; Venturieri, A.; Diniz, C.G.; Dessay, N.; Durieux, L.; Gomes, A.R. High spatial resolution land use and land cover mapping of the Brazilian Legal Amazon in 2008 using Landsat-5/TM and MODIS data. *Acta Amaz.* **2016**, *46*, 291–302. [CrossRef]

16. INPE. Projeto Terraclass Cerrado: Mapeamento do uso e Cobertura Vegetal do Cerrado. Available online: http://www.dpi.inpe.br/tccerrado/ (accessed on 22 March 2017).

17. INPE. Projeto Terraclass Amazônia: Mapeamento do uso e Cobertura da Terra na Amazônia Legal Brasileira. Available online: http://www.inpe.br/cra/projetos_pesquisas/dados_terraclass.php (accessed on 22 March 2017).

18. IBGE, Instituto Brasileiro de Geografia e Estatística. Produção Agrícola Municipal (PAM) 2014. Available online: http://www.ibge.gov.br/home/estatistica/economia/pam/2014/ (accessed on 15 April 2017).

19. Huete, A.; Didan, K.; Miura, T.; Rodriguez, E.P.; Gao, X.; Ferreira, L.G. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. *Remote Sens. Environ.* **2002**, *83*, 195–213. [CrossRef]

20. Soudani, K.; le Maire, G.; Dufrêne, E.; François, C.; Delpierre, N.; Ulrich, E.; Cecchini, S. Evaluation of the onset of green-up in temperate deciduous broadleaf forests derived from Moderate Resolution Imaging Spectroradiometer (MODIS) data. *Remote Sens. Environ.* **2008**, *112*, 2643–2655. [CrossRef]

21. Busetto, L.; Ranghetti, L. MODISStp: An R Package for Preprocessing of MODIS Time Series. R Package Version 1.3.2. 2017. Available online: https://cran.rstudio.com/web/packages/MODISStp/index.html (accessed on 15 April 2017).

22. Solano, R.; Didan, K.; Jacobson, A.; Huete, A. *MODIS Vegetation Index User's Guide (MOD13 Series)*; Vegetation Index and Phenology Lab, University of Arizona: Tucson, AZ, USA, 2010.

23. Chen, J.; Jönsson, P.; Tamura, M.; Gu, Z.; Matsushita, B.; Eklundh, L. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter. *Remote Sens. Environ.* **2004**, *91*, 332–344. [CrossRef]

24. Savitzky, A.; Golay, M.J.E. Smoothing and differentiation of data by simplified least squares procedures. *Anal. Chem.* **1964**, *36*, 1627–1639. [CrossRef]
25. Eastman, J.R.; Fulk, M. Long sequence time series evaluation using standardized principal components. *Photogramm. Eng. Remote Sens.* **1993**, *59*, 991–996.

26. Gurgel, H.C.; Ferreira, N.J. Annual and interannual variability of NDVI in Brazil and its connections with climate. *Int. J. Remote Sens.* **2003**, *24*, 3595–3609. [CrossRef]

27. Hall-Beyer, M. Comparison of single-year and multiyear NDVI time series principal components in cold temperate biomes. *IEEE Trans. Geosci. Remote Sens.* **2003**, *41*, 2568–2574. [CrossRef]

28. Hirosawa, Y.; Marsh, S.E.; Kliman, D.H. Application of standardized principal component analysis of land-cover characterization using multitemporal AVHRR data. *Remote Sens. Environ.* **1996**, *58*, 267–281. [CrossRef]

29. Wang, T.; Kou, X.; Xiong, Y.; Mou, P.; Wu, J.; Ge, J. Temporal and spatial patterns of NDVI and their relationship to precipitation in the Loess Plateau of China. *Int. J. Remote Sens.* **2010**, *31*, 1943–1958. [CrossRef]

30. Young, S.S.; Anyamba, A. Comparison of NOAA/NASA PAL and NOAA GVI data for vegetation change studies over China. *Photogramm. Eng. Remote Sens.* **1999**, *65*, 679–696.

31. Leutner, B.; Horning, N. Rstoolbox: Tools for Remote Sensing Data Analysis. R Package Version 0.1.8. 2017. Available online: https://cran.r-project.org/web/packages/RStoolbox/index.html (accessed on 15 April 2017).

32. Trimble. *Ecognition Developer*© 9.0 User Guide; Trimble Germany GmbH: Munich, Germany, 2014.

33. Baatz, M.; Schäpe, A. Multi-resolution segmentation: An Optimization Approach for High Quality Multi-Scale Image Segmentation. In *Angewandte Geographische Informations-Verarbeitung, XII*; Wichmann Verlag: Karlsruhe, Germany, 2000; pp. 12–23.

34. Benz, U.C.; Hofmann, P.; Willhauck, G.; Lingenfelder, I.; Heynen, M. Multi-resolution, object-oriented fuzzy analysis of remote sensing data for GIS-ready information. *ISPRS J. Photogramm.* **2004**, *58*, 239–258. [CrossRef]

35. Almeida, C.A.; Mourão, M.; Dessay, N.; Lacques, A.-E.; Monteiro, A.; Durieux, L.; Venturieri, A.; Seyler, F. Typologies and spatialization of agricultural production systems in Rondônia, Brazil: Linking land use, socioeconomic and territorial configuration. *Land* **2016**, *5*, 18. [CrossRef]

36. Cai, H.; Zhang, S.; Bu, K.; Yang, J.; Chang, L. Integrating geographical data and phenological characteristics derived from MODIS data for improving land cover mapping. *J. Geogr. Sci.* **2011**, *21*, 705–718. [CrossRef]

37. Cano, E.; Denux, J.-P.; Bisquert, M.; Hubert-Moy, L.; Chéret, V. Improved forest-cover mapping based on MODIS time series and landscape stratification. *Int. J. Remote Sens.* **2017**, *38*, 1865–1888. [CrossRef]

<figure><img src="image_81.png" /></figure>

© 2017 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
<figure><img src="image_82.png" /></figure>

# Improved regional-scale Brazilian cropping systems' mapping based on a semi-automatic object-based clustering approach

<figure><img src="image_83.png" /></figure>

Check for
updates

Beatriz Bellón<sup>a,b,*</sup>, Agnès Bégué<sup>a,b</sup>, Danny Lo Seen<sup>a,b</sup>, Valentine Lebourgeois<sup>a,b</sup>, Balbino Antônio Evangelista<sup>c</sup>, Margareth Simões<sup>d,e</sup>, Rodrigo Peçanha Demonte Ferraz<sup>d</sup>

## ARTICLE INFO

### Keywords:

Cropping systems

Geographic object-based image analysis
(GEOBIA)

Stratification

Clustering

Moderate resolution imaging

spectroradiometer (MODIS)

Landsat

## ABSTRACT

Cropping systems' maps at fine scale over large areas provide key information for further agricultural production and environmental impact assessments, and thus represent a valuable tool for effective land-use planning. There is, therefore, a growing interest in mapping cropping systems in an operational manner over large areas, and remote sensing approaches based on vegetation index time series analysis have proven to be an efficient tool. However, supervised pixel-based approaches are commonly adopted, requiring resource consuming field campaigns to gather training data. In this paper, we present a new object-based unsupervised classification approach tested on an annual MODIS 16-day composite Normalized Difference Vegetation Index time series and a Landsat 8 mosaic of the State of Tocantins, Brazil, for the 2014-2015 growing season. Two variants of the approach are compared: an hyperclustering approach, and a landscape-clustering approach involving a previous stratification of the study area into landscape units on which the clustering is then performed. The main cropping systems of Tocantins, characterized by the crop types and cropping patterns, were efficiently mapped with the landscape-clustering approach. Results show that stratification prior to clustering significantly improves the classification accuracies for underrepresented and sparsely distributed cropping systems. This study illustrates the potential of unsupervised classification for large area cropping systems' mapping and contributes to the development of generic tools for supporting large-scale agricultural monitoring across regions.

## 1. Introduction

Agricultural systems are constantly evolving to meet the increasing global demand for food production (FAO, 2009). Both the cropland expansion and the intensification of management practices contribute to meet such demand, increasing the concern on the environmental impact associated with these land-use changes (Tilman et al., 2011). In view of these large-scale dynamics, there is a growing need for regularly updated spatial information on the cropping systems at regional scales.

Cropping systems result of a combination of human and environmental factors, which naturally vary through space. A given cropping system thus presents local specificities, leading to a high spatial variability over large areas. In addition, when considering large areas, the cropping systems are embedded in a complex landscape mosaic composed of a high diversity of land use and land cover (LULC) types.

Therefore, large-area mapping of cropping systems and their management practices, remains a challenging task, and as a result, their description and location remain rather unclear for most regions (Leenhardt et al., 2010).

Recent studies have however shown the usefulness of earth observation (EO) data to successfully map cropping systems, and some management practices over large areas (Atzberger, 2013; Bégué et al., 2015). In the last decade, various of these studies have successfully used vegetation index (VI) time series, mostly MODIS-derived, to map, at a regional scale, intensification practices such as the adoption of sequential cropping (e.g. Arvor et al., 2011; Brown et al., 2013; Cheema and Bastiaanssen, 2010; Guan et al., 2016; Gumma et al., 2016; Mingwei et al., 2008; Wardlow and Egbert, 2008).

Most of these studies, however, rely entirely on reference data, usually ground observations, to train supervised classifications (e.g.

<footnote><sup>a</sup> CIRAD, UMR TETIS, F-34398 Montpellier, France</footnote>

<footnote><sup>b</sup> TETIS, Univ Montpellier, AgroParisTech, CIRAD, CNRS, IRSTEA, Montpellier, France</footnote>

<footnote><sup>c</sup> Embrapa Pesca e Aquicultura, Prolongamento da Avenida NS 10, cruzamento com a Avenida LO 18, sentido Norte, loteamento Água Fria, 77008-900 Palmas, TO, Brazil</footnote>

<footnote><sup>d</sup> Embrapa Solos, Rua Jardim Botânico, 1024, 22460-000, Rio de Janeiro, RJ, Brazil</footnote>

<footnote><sup>e</sup> Rio de Janeiro State University (UERJ), Rua São Francisco Xavier, 524, Maracanã, 20550-900, Rio de Janeiro, RJ, Brazil</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: beatriz.bellon@cirad.fr (B. Bellón), agnes.begue@cirad.fr (A. Bégué), danny.lo_seen@cirad.fr (D. Lo Seen), valentine.lebourgeois@cirad.fr (V. Lebourgeois), balbino.evangelista@embrapa.br (B.A. Evangelista), margareth.simoes@embrapa.br (M. Simões), rodrigo.demonte@embrapa.br (R.P. Demonte Ferraz).</footnote>
Arvor et al., 2011; Brown et al., 2013; Mingwei et al., 2008). Gathering reference data twice or more during a growing season to train classifications of sequential cropping systems represents a time and resource consuming task. This is especially the case for large areas since it is difficult to extensively cover the study area to guarantee representativeness of the intra-class spectral variability linked to the high spatial variability (Franklin and Wulder, 2002).

Contrastingly, unsupervised classifications, being solely based on the VI time series' information, are more operational without affecting the accuracy of the resulting maps, as shown by Cheema and Bastiaanssen (2010), Wardlow and Egbert (2008), Gumma et al. (2016). Hyperclustering techniques are commonly used, producing an initial large number of clusters that are then reduced by successive processing steps. This technique ensures that under-represented classes of interest are correctly isolated, and has the advantage of reducing the impact of the user-defined parameters (e.g. the number of clusters) on the results (Cihlar, 2000). However, the classification refinement techniques used in these studies, such as merging or splitting clusters, are highly expert knowledge-driven and often rely on external data (e.g. cropland masks, reference ground data), which limits their reproducibility.

In order to deal with the high intra-class spectral variability without depending on external data, cropping systems mapping over large areas would largely benefit from spectral variability reduction steps prior to classification. Geographic object-based image analysis (GEOBIA) techniques are widely used to reduce spectral variability by grouping individual pixels into thematically meaningful objects that summarize their spectral properties (Blaschke et al., 2014; Bunker et al., 2016; Hay and Castilla 2008). Another spectral variability reduction strategy, referred to as stratification, consists in partitioning the study regions into "biophysical-meaningful" sub-regions prior to classification (Cai et al., 2011; Cano et al., 2017; Franklin and Wulder, 2002; Vintrou et al., 2012; Wardlow and Egbert, 2008; Xiong et al., 2017). Traditionally, external maps are used (e.g. ecoregion, agro-ecological zoning maps), but some novel approaches have started to explore the potential of EO data-based landscape stratifications for large area land-cover classifications (e.g. Cano et al., 2017; Vintrou et al., 2012).

Based on this context, we hypothesize that the use of stratification prior to unsupervised classification improves the accuracy of the regional-scale cropping systems' classifications at the field-level by reducing the spectral variability of the EO data. To test our hypothesis, we propose a new semi-automatic GEOBIA clustering approach, and compare the classification accuracies obtained with two variants of the approach: an hyperclustering approach, which includes a field-level segmentation step prior to classification, and a landscape-clustering approach, which combines two spectral variability reduction steps prior to classification, a field-level segmentation and a landscape stratification.

This work contributes to the development of simple and generic tools to improve large-area cropping systems' mapping with new classification methods fully based on EO data. The approach benefits from the high spatial resolution of a Landsat mosaic and the high temporal resolution of a MODIS VI time series, and is tested on the Tocantins state for the 2014-2015 growing season. The classification results are evaluated with ground-truth data and the annual crop estimates of the Brazilian Institute of Geography and Statistics (IBGE) (IBGE, 2014).

## 2. Materials

### 2.1. Study region

The Tocantins state, located in the centre-north of Brazil (Fig. 1a), is an agricultural region of strong anthropic influence which has undergone fast land-use changes over the last decades. Since its creation in 1988, Tocantins experiences a growing population and a fast-evolving economy based on a developing agricultural sector (Pedroso da Silva and de Almeida, 2007; Fornaro, 2012). Apart from the extensive

development of the livestock systems which has been induced by the spreading of the transport network since the 1970s, this region undergoes, since the last decade, an expansion of the cropland mainly to the detriment of areas of degraded pasture, but also of some native vegetation areas (Barona et al., 2010; Silva, 2007; Fornaro, 2012). This trend is greatly driven by the global demand for Brazilian soybean, which is responsible for the expansion of the soybean cultivated areas in the latest agricultural frontier of Brazil, the MATOPIBA region, named after the states of Maranhão, Tocantins, Piauí and Bahia (Bolfe et al., 2016; Fornaro, 2012).

As a result, three soybean-based cropping systems, characterized by highly mechanized management practices and large field sizes, of around 100 ha, dominate the agricultural landscape of Tocantins. One of the cropping systems is a single cropping system, characterized by a single soybean crop cycle, and the other two cropping systems are characterized by a double-cropping pattern with two distinct crop cycles in the same growing season. This form of intensification is commonly adopted in this region since it allows increasing yields while reducing the need for expansion of the cropland areas (Spangler et al., 2017).

In one of the double cropping systems, henceforth referred to as the 'soybean-cereal double cropping system', a soybean crop is grown in the summer season, followed, in succession, by a cereal crop (mostly maize) during the end of summer to fall season. The other double-cropping system henceforth referred to as the 'rice-soybean double cropping system', is localized in the lowlands of the Formoso River basin and concerns a rainy season rice crop (grown under flood irrigation) succeeded by a late soybean crop.

### 2.2. Ground-truth data collection

Field surveys were conducted in the study area during October 2015, at the end of the 2014-2015 growing season, just before the rainy summer period and beginning of the next season. The sampling strategy aimed at collecting multiple observations for the three major cropping systems' classes and other LULC types over a wide variety of zones.

In particular, different agricultural zones, scattered throughout the study area were targeted and surveyed thoroughly, in order to capture the regional variability of the cropping systems. The selection of zones to survey was assisted by the NTSA (Thematic Unit on Agricultural Systems) of EMBRAPA (Brazilian Agricultural Research Corporation) Fisheries and Aquaculture, and was mainly focused on identifying zones where the crop fields are spatially "concentrated".

The survey equipment included a Trimble Yuma2 rugged tablet with an integrated GPS and Quantum GIS software, which allowed real-time navigation and direct recording of the waypoints in a spatial database. The tablet was equipped with the multispectral Landsat 8 OLI 2014-2015 mosaic (presented in Section 2.4.2), overlaid by the field-level segmentation (presented in Section 3.1.1). The sampling consisted in associating the observed class to its corresponding object in the segmentation layer. The commonly adopted no-tillage system allowed us to identify the cropping systems by the crop residues that cover the fields after the harvest of the crops.

A total of 900 GPS waypoints belonging to 35 different municipalities were registered (presented in Fig. 2), representing the three major cropping systems belonging to the annual cropland and other LULC types.

### 2.3. Area-based agricultural statistics

Area-based estimates of the soybean, second crop maize and rice crops were acquired from the Municipal Agricultural Production database (PAM) of the IBGE (2014). The 2015 harvested area estimates of the 139 municipalities of Tocantins were used to evaluate our classification results. The 'second crop maize' area estimates were used to evaluate the classification of the soybean-cereal double cropping
<figure><img src="image_84.png" /><figcaption>Fig. 1. (a) Location of Tocantins state, Brazil, (b) Administrative Micro-regions that compose the Tocantins state and limits of municipalities (grey lines).</figcaption></figure>

system since no statistics are available on the other cereals which are alternatively grown as cover or fodder crops after the soybean crop harvest (e.g. sorghum or millet).

## 2.4. Earth observation data acquisition and pre-processing

### 2.4.1. MODIS NDVI time series

An annual time series of the MODIS VIs product (MOD13Q1 collection 5) (Solano et al., 2010) was acquired from the NASA Land processes distributed active archive centre through the processing R package MODISsp (Busetto and Ranghetti, 2016). Two tiles (h13v09, h13v10) of the 16-day composited Normalized Difference Vegetation Index (NDVI) band were mosaicked pairwise for each of the composite periods to generate a total of 23 NDVI mosaic images, from the beginning of the growing season in October 2014 to October 2015, covering the study region at a 250 m spatial resolution.

Although the compositing process of the MOD13Q1 product already reduces noise effects by selecting the observation with the highest NDVI value and the smallest view angle among the 16 days' acquisitions (Solano et al., 2010), the dataset was further pre-processed to reduce remaining effects. The dataset was pre-processed in two steps (presented in detail in the '3.1. Pre-processing of the NDVI time series' section in Bellón et al. (2017)): First, a noise-reduction algorithm based on the Savitzky-Golay function-fitting smoothing method introduced by Chen et al. (2004) was applied to correct pixel values mainly affected by clouds and other atmospheric effects and improve the overall quality of the NDVI time series. Finally, the smoothed time series was linearly interpolated using the acquisition dates for each pixel retained during the compositing process, which are provided in the 'Composite Day of Year (CDOY)' band of the MOD13Q1 product, to reconstruct the NDVI

time series at regular 16-day intervals, and therefore avoid temporal discontinuities in the time series.

### 2.4.2. Landsat 8 OLI mosaic

Nineteen Landsat 8 Operational Land Imager (OLI) scenes were mosaicked to generate a seamless image covering the whole state of Tocantins at 30 m spatial resolution. The surface reflectance Landsat 8 OLI scenes were retrieved from the United States Geological Survey (USGS) Landsat archive via the ESPA Ordering interface (Jenkerson, 2013; Zanter, 2016). In order to obtain a cloud-free mosaic, we selected scenes which were cloud-free inside the Tocantins state's boundaries, acquired during July 2014 and July 2015. Table 1 presents the paths and rows numbers of the Landsat scenes and their corresponding acquisition dates.

## 3. Methods

The presented approach was conceived as a GEOBIA classification approach. The high resolution Landsat mosaic was segmented into thematically meaningful field level objects which were subsequently classified through clustering. Two clustering approaches were tested and compared: the hyperclustering approach and the landscape-clustering approach. Both approaches employed the k-means classifier and conducted clustering in two steps: a pre-clustering step partitioned the objects into a large number of initial sub-clusters, and a final-clustering step partitioned the resulting sub-clusters into the desired number of final clusters. The landscape-clustering approach, unlike the hyperclustering approach, uses landscape units, which are obtained though stratification of the study region with a landscape-level segmentation, for the pre-clustering step. The methodological framework
<figure><img src="image_85.png" /><figcaption>Fig. 2. Location and number of waypoints registered for each of the four classes in the spatial database during the ground-truth data collection field surveys in October 2015.</figcaption></figure>

used in this study is presented in Fig. 3 and described in the following sub-sections.

## 3.1. Segmentation at the field and landscape levels

### 3.1.1. Field-level segmentation

A field-level segmentation was performed on the Landsat mosaic to delineate the boundaries of the crop fields and other landscape features. The segmentation was carried out using the ‘Multiresolution Segmentation’ algorithm proposed by Baatz and Schäpe in the eCognition® Developer 9.0 software (Trimble Germany GmbH, 2014; Baatz and Schäpe, 2000).

The segmentation algorithm was applied to five spectral bands of the Landsat mosaic, which were equally weighted, and the segmentation parameters were empirically set after sequential testing of different parameter combinations and visual evaluation of the different segmentation results (Table 2 presents the final segmentation parameters). The objective of the segmentation being the delineation of the overall rectangular-shaped crop fields, we heavily weighted the ‘shape parameter’ compared to the ‘colour’ (radiometry), and the ‘compactness’ compared to the ‘smoothness’, resulting in compact, regular-shaped objects, which allowed an effective delineation of the agricultural fields over different parts of the study area (Fig. 4 presents subset images of the segmentation results). A total of 187,877 objects at a spatial resolution of 30 m were obtained over the entire study area.

### 3.1.2. Landscape-level segmentation

The study region was partitioned into landscape-level units following the stratification method proposed by Bellón et al. (2017) based

<caption>Table 1</caption>

Acquisition details of the Landsat 8 OLI scenes retrieved from the USGS archive to compose the cloud-free mosaic of Tocantins.

<table><thead><tr><th>Acquisition date</th><th>Path</th><th>Row</th><th>Sun elevation (degrees)</th><th>Sun azimuth (degrees)</th><th>Scene cloud cover (%)<sup>a</sup></th></tr></thead><tbody><tr><td rowspan="6">17/07/2015</td><td rowspan="6">223</td><td>64</td><td>50.52</td><td>46.03</td><td>0.06</td></tr><tr><td>65</td><td>49.35</td><td>45.00</td><td>0.05</td></tr><tr><td>66</td><td>48.18</td><td>44.04</td><td>0.00</td></tr><tr><td>67</td><td>46.98</td><td>43.15</td><td>0.00</td></tr><tr><td>68</td><td>45.78</td><td>42.33</td><td>0.00</td></tr><tr><td>69</td><td>44.56</td><td>41.56</td><td>0.00</td></tr><tr><td rowspan="6">07/07/2014</td><td rowspan="6">222</td><td>64</td><td>49.89</td><td>44.06</td><td>0.04</td></tr><tr><td>65</td><td>48.70</td><td>43.10</td><td>0.00</td></tr><tr><td>66</td><td>47.50</td><td>42.21</td><td>0.00</td></tr><tr><td>67</td><td>46.28</td><td>41.38</td><td>0.00</td></tr><tr><td>68</td><td>45.05</td><td>40.61</td><td>0.00</td></tr><tr><td>69</td><td>43.81</td><td>39.90</td><td>0.21</td></tr><tr><td rowspan="4">19/07/2015</td><td rowspan="4">221</td><td>66</td><td>48.39</td><td>44.46</td><td>0.00</td></tr><tr><td>67</td><td>47.20</td><td>43.56</td><td>0.00</td></tr><tr><td>68</td><td>46.00</td><td>42.72</td><td>0.00</td></tr><tr><td>69</td><td>44.79</td><td>41.94</td><td>0.44</td></tr><tr><td rowspan="3">25/07/2014</td><td rowspan="3">220</td><td>67</td><td>48.08</td><td>44.87</td><td>0.68</td></tr><tr><td>68</td><td>46.90</td><td>43.97</td><td>0.00</td></tr><tr><td>69</td><td>45.70</td><td>43.14</td><td>0.00</td></tr></tbody></table>

on the phenological information of a MODIS VI time series. This method consists in an automatic segmentation of the study area into units of land which present certain homogeneity in their physiognomy and land use that differs from adjacent units, while being

<footnote><sup>a</sup> The cloud covered areas in the scenes correspond to areas outside of the Tocantins state’s boundaries.</footnote>
<figure><img src="image_86.png" /><figcaption>Fig. 3. Methodological framework.</figcaption></figure>

heterogeneous in their structure, since they are extensive enough to contain a mosaic of patches with different land cover types, and thus are henceforth referred to as landscape units (see Bellón et al. (2017) for further details).

This landscape stratification method was originally conceived to produce an agricultural land-use systems' map at the landscape level by classifying the large-scale land units. In this study, the landscape stratification is otherwise used as a spectral variability reduction strategy, to assist regional-scale cropping systems' mapping at the field-level.

The landscape-level segmentation is performed using the 'Multiresolution Segmentation' algorithm previously presented (see Section 3.1.1), and the segmentation parameters were set according to the original method presented in Bellón et al. (2017) (shown in Table 2). A total of 120 landscape units were extracted following this method (Fig. 5).

## 3.2. Unsupervised classification

The two variants of the clustering approach were based on the temporal NDVI information contained in the MODIS time series. A common dataset was generated by extracting the median NDVI value per Landsat field-level object for each of the 23 dates in the time series. Due to the difference in spatial resolution between the MODIS pixels at 250 m and the Landsat objects at 30 m, only the values belonging to the pixels whose centroids were located inside of an object were used in the

calculation of the median NDVI value for that given object. The median value was chosen as the measure of central tendency since it is less sensitive to skewness in the distribution of values than the mean, and therefore is less affected by the influence of the mixed pixels (whose radiometry results from different cover types) situated on the object's boundary.

Both variants of the clustering approach were based on the k-means algorithm to classify the dataset of field-level objects containing the median NDVI value for each of the 23 dates in the MODIS time series. Originally formulated by MacQueen (1967), the k-means algorithm aims at partitioning the data into 'k' groups' such that the within-cluster sum of squares is minimized. This widely used partitional clustering algorithm has a number of variations. For this study, the clustering was performed by means of the algorithm proposed by Hartigan and Wong (1979), implemented in the 'kmeans' tool in the processing R package 'stats' (R Core Team, 2016).

### 3.2.1. Hyperclustering approach

In the hyperclustering approach, a pre-clustering was performed on the dataset of field-level objects. The number of sub-clusters *k* was set to 2400 to match the number of sub-clusters *k'* in the landscape-clustering approach presented in the next section (Section 3.2.2).

The mean NDVI temporal profile per sub-cluster was then calculated, resulting in a total of 2400 mean NDVI temporal profiles which were finally partitioned into ten final clusters. The number of final

<table><caption>**Table 2**
Field-level and landscape-level segmentation parameters (see Bellón et al. (2017) for further details on the landscape-level segmentation).</caption><thead><tr><th>Segmentation level</th><th>Input data</th><th>Spatial resolution</th><th>Bands (equally weighted)</th><th>Scale</th><th>Colour</th><th>Shape</th></tr></thead><tbody><tr><td>Field level</td><td>Landsat 8 OLI mosaic</td><td>30 m</td><td>Blue (B2)<br>Green (B3)<br>Red (B4)<br>NIR (B5)<br>SWIR 1 (B6)</td><td>110</td><td>0.2</td><td>0.8 (compactness = 1;<br>smoothness = 0)</td></tr><tr><td>Landscape level</td><td>Principal component (PC) images computed from the<br>MODIS NDVI time series</td><td>250 m</td><td>22 PC images</td><td>850</td><td>1</td><td>0</td></tr></tbody></table>
<figure><img src="image_87.png" /><figcaption>Fig. 4. Image subsets of the Landsat field-level segmentation results (field-level objects) over different landscapes. The background image corresponds to the false-colour Landsat mosaic.</figcaption></figure>

clusters was empirically set to 10, expecting 10 clusters to adequately discriminate the four targeted classes while still being easily interpretable.

Since random seeding was used for initial clustering in the k-means algorithm, ten replicates were conducted and the clustering result with the highest degree of difference between clusters and similarity within clusters was finally retained. This was justified by the highest ratios of the between-clusters sum of squares to the total sum of squares ($R^2$) and of the between-clusters sum of squares to the total within-clusters of squares ($R^2/(1 - R^2)$).

### 3.2.2. Landscape-clustering approach

In the landscape-clustering approach, the landscape units were first spatially overlaid by the field-level objects, and all of the objects whose centroid was located inside of the boundaries of a given landscape unit were grouped, resulting in one group of field-level objects per landscape unit. A pre-clustering was then performed on each group of objects, making as much partitions as number of landscape units (n = 120). The exact number of spectral classes per landscape unit being unknown, the number of sub-clusters per group of objects was empirically set to 20, which represents 5 times more clusters than final targeted classes, summing up to a total of 2400 *k'* sub-clusters (*k'* = 20 sub-clusters × 120 landscape units) for the whole study region.

The mean NDVI temporal profile was then calculated for each of the sub-clusters, resulting in a total of 2400 mean NDVI temporal profiles.

A final clustering was performed on the *k'* mean NDVI temporal profiles to obtain 10 final clusters.

Ten replicates of the final clustering on the mean NDVI temporal profiles were evaluated and the clustering result with the highest $R^2$ and $R^2/(1 - R^2)$ ratios was retained.

### 3.2.3. Labelling of final clusters

The mean NDVI temporal profiles of the ten final clusters of each approach were calculated and plotted. In general, the NDVI temporal profiles of cropping systems present fast NDVI increase and decrease trends related to the different stages of crop growth and harvest periods and thus of the crop cycles. This type of pattern allowed us to identify and label the final clusters through visual analysis of their NDVI variations over the growing season.

The clusters whose NDVI profile's shapes followed the cropping calendar of the main cropping systems of Tocantins were labelled with their corresponding cropping system class, while the rest were labelled as the 'other LULC types' class.

## 3.3. Accuracy assessment of the classification results

The accuracy of both the landscape clustering and the hyperclustering results was estimated and compared by using the ground-truth data and the agricultural statistics as a reference.

Six accuracy measures were calculated using the 900 reference
<figure><img src="image_88.png" /><figcaption>Fig. 5. The landscape-level segmentation result.</figcaption></figure>

objects in the ground-truth database: the overall accuracy with its associated confidence interval (Rosenfield and Fitzpatrick-Lins, 1986), the Cohen's Kappa coefficient (Cohen, 1960), the quantity

disagreement and allocation disagreement (Pontius and Millones, 2011), the McNemar's Chi-squared test (McNemar, 1947) and the F-measure (Labatut and Cherifi, 2012).

The overall accuracy corresponds to the total number of the correctly classified objects divided by the total number of objects being evaluated and is expressed as a percentage. The Kappa coefficient measures the difference between the observed agreement and the probability of agreement occurring by chance and ranges between 0 (no more agreement than expected by chance) and 1 (perfect agreement). The Kappa coefficient is often used as an overall measure of accuracy, but its use is being questioned by numerous articles (Strahler et al., 2006; Foody, 2002; Pontius and Millones, 2011). The quantity disagreement and allocation disagreement have therefore been included as complementary measures for the accuracy assessment. These two measures decompose the overall disagreement (which is the complement of the overall accuracy) into two types of disagreement and thus provide additional insight on the accuracy of the classification results.

In our case quantity disagreement expresses the amount of difference between the reference objects and the classified objects that is due to the less than perfect match in the proportions of the classes, and allocation disagreement expresses the amount of difference between the reference objects and the classified objects that is due to the less than optimal match in the spatial allocation of the classes given the proportions of the classes in the reference and classified objects. These two measures are expressed as a percentage.

The McNemar's Chi-squared test was used to evaluate if there is a statistically significant difference between the two classification results' accuracies. Finally, the class-specific accuracy was evaluated using the F-measure, which corresponds to the weighted harmonic mean of the producer's and user's accuracy and ranges from 0 to 1, being the lowest and the highest accuracy, respectively.

The crop area-based estimates in the agricultural statistics at the municipality level provided by the IBGE were finally used to evaluate the accuracy of the best classification result according to the ground-truth data accuracy assessment tests. A correlation coefficient was

<figure><img src="image_89.png" /><figcaption>Fig. 6. The temporal NDVI profiles' plots of the 10 final clusters obtained with (a) the hyperclustering approach, and (b) the landscape-clustering approach, and their corresponding labels. The continuous line curves represent the mean NDVI of all of the objects belonging to each cluster against the day of the year (DOY). The mean NDVI profiles of the clusters which seem to have captured the three main cropping systems of Tocantins are represented in separate plots (i), (ii) and (iii) with their standard deviation (discontinuous line curves), and the mean NDVI profiles of the clusters which seem to have captured other LULC types are represented in the same plot (iv).</figcaption></figure>
computed between the estimated crop area and the reference IBGE estimates at the municipality level.

# 4. Results

## 4.1. Clustering results

Since 120 landscape units were obtained, the number of pre-clusters for the whole study area was 2400 for both the landscape-clustering approach (20 pre-clusters per landscape unit) and the hyperclustering approach. The plots of the mean NDVI temporal profiles of the final ten clusters are presented in Fig. 6 for the hyperclustering (Fig. 6(a)), and the landscape-clustering approach (Fig. 6(b)).

Three distinct cropping patterns related to the main cropping systems of Tocantins can be identified among the NDVI temporal profiles of the final clustering results. These are particularly noticeable in three clusters' profiles obtained with the landscape-clustering approach (Fig. 6(b)(i–iii)), and in one of the clusters' profiles obtained with the hyperclustering approach (Fig. 6(a)(iii)), since they present fast NDVI increase and decrease trends. The rest of the cluster's profiles present remarkably smoother shapes.

Nevertheless, two clusters' profiles obtained with the hyperclustering approach which present slight NDVI oscillations seem to have captured the temporal variations related to the cropping systems and thus can be labelled accordingly, even if the general flatness of their profiles suggests mixture of classes in both these clusters (Fig. 6(a)(i) (ii)).

The first cropping pattern which is clearly recognizable in the temporal NDVI profile plot in Fig. 6(b)(i) exhibits a major single crop cycle that peaks around February, suggesting that the corresponding cluster has grouped the objects belonging to the soybean single cropping system. The cropping pattern in the temporal NDVI profile plot in Fig. 6(b)(ii) seems distinctive of the soybean-cereal double cropping system since the NDVI values are low around Mars, reflecting the soybean harvest, after which the NDVI increases as a result of the second cereal's growth until it peaks around May. Finally, the cropping pattern in the temporal NDVI profile plot in Fig. 6(b)(iii) seems distinctive of the rice-soybean double cropping system, with a long crop cycle during the rainy season corresponding to the rice crop cycle, followed by a shorter, late crop cycle from May to September. This late second soybean crop cycle is evidenced by the NDVI peak period around July.

## 4.2. Accuracy assessment results

Table 3 presents the overall accuracy, the Cohen's Kappa coefficient and the quantity and allocation disagreement of the classification results in four classes (the three main cropping systems and the 'other LULC types' class), derived from the hyperclustering and the landscape-clustering approaches.

The overall accuracy is higher for the landscape-clustering classification than for the hyperclustering classification i.e. the landscape-clustering approach correctly classified a larger number of objects than the hyperclustering approach, which suggests that the landscape-clustering approach performed better than the hyperclustering approach in the classification of the four classes. Accordingly, the overall disagreement (composed of the quantity disagreement and the allocation

disagreement) was lower for the landscape-clustering classification (6.8%) than for the hyperclustering classification (22.1%). Quantity disagreement was higher than allocation disagreement for both classification results, suggesting that the classification errors were related in a higher degree to a quantification error than to a location error.

The Kappa coefficient result for the landscape-clustering approach is also higher than that of the hyperclustering approach i.e. the difference between the observed agreement between the classified objects and the reference objects and the probability of agreement occurring by chance was higher for the landscape-clustering classification than for the hyperclustering classification, suggesting that the landscape-clustering classification's reliability is higher.

Furthermore, the results of the McNemar's test (presented in Table 4), showed a statistically significant improvement of the landscape-clustering classification results over the hyperclustering results ($p$-value < .0001).

Class-specific accuracy is presented in Table 5 through the F-measure results of the landscape-clustering and the hyperclustering approaches. The 'soybean single cropping system' class obtained the lowest F-measure value for both classifications (F-measure < 0.7), and contrastingly, the 'rice-soybean double cropping system' class was correctly classified by both approaches, obtaining the highest F-measure (F-measure = 1). Both the 'soybean-cereal double cropping system' and the 'other LULC types' classes obtained high F-measure values (F-measure > 0.8) for the landscape-clustering approach but only the 'other LULC types' class did for the hyperclustering approach, the 'soybean-cereal double cropping system' class presenting the greatest difference between both approaches. Except for the 'rice-soybean double cropping system' class which obtained the same F-measure for both approaches, the rest of the classes obtained lower scores with the hyperclustering approach, which indicates a better performance of the landscape-clustering approach.

The soybean crop area-based estimates in the municipality-level agricultural statistics provided by the IBGE are compared to the landscape-clustering approach's soybean area results in Fig. 7. The total area of soybean (considering all three cropping systems) estimated by the landscape-clustering approach accounts for 73% of the total area reported by the IBGE, obtaining a strong positive correlation ($r$) of 0.87 ($r^2$ = 0.75) at the municipality level.

Likewise, the second crop maize IBGE area estimates show a strong correlation with the landscape-clustering 'soybean-cereal double cropping system' classification results, with a correlation coefficient at the municipality level of 0.87 ($r^2$ = 0.76). Finally, the rice IBGE crop area estimates also show a strong correlation with the landscape clustering 'rice-soybean double cropping system', with a correlation coefficient of 0.99 ($r^2$ = 0.97).

## 4.3. Cropping systems' map obtained with the landscape-clustering approach

Fig. 8 presents the map derived from the classification results of the main cropping systems of Tocantins at the field level obtained with the landscape-clustering approach. The accuracy of the map corresponds to the accuracy of the landscape-clustering classification results presented in Section 4.2.

As shown in the map (Fig. 8), the rice-soybean double cropping system is exclusively localized in the southwest of Tocantins, while the

<caption>**Table 3**</caption>

The overall classification accuracy and its associated confidence interval (CI), the overall Kappa coefficient and the quantity and allocation disagreement results for 900 classified objects from each clustering approach.

<table><thead><tr><th>Clustering approach</th><th>Overall accuracy (%)</th><th>95% CI (%)</th><th>Overall Kappa</th><th>Quantity disagreement (%)</th><th>Allocation disagreement (%)</th></tr></thead><tbody><tr><td>Landscape-clustering</td><td>93.2</td><td>91.4–94.8</td><td>0.81</td><td>4.5</td><td>2.3</td></tr><tr><td>Hyperclustering</td><td>77.9</td><td>75.0–80.6</td><td>0.52</td><td>15.3</td><td>6.8</td></tr></tbody></table>
<caption>**Table 4**</caption>

Contingency table summarizing the outcomes of the two clustering approaches (number of objects correctly classified or misclassified) and McNemar's Chi-squared test ($\chi^2$) results.

<table><caption>**Table 5**</caption><thead><tr><th></th><th>Landscape-clustering approach<br>(correctly classified)</th><th>Landscape clustering approach<br>(misclassified)</th><th></th></tr></thead><tbody><tr><td>Hyperclustering approach (correctly<br>classified)</td><td>684</td><td>17</td><td rowspan="2">$\chi^2 = 110.72$, with degree of freedom = 1 (p-value < .0001)</td></tr><tr><td>Hyperclustering approach (misclassified)</td><td>155</td><td>44</td></tr></tbody></table>

Class F-measure results.

<table><thead><tr><th>Class</th><th>Landscape-clustering approach</th><th>Hyperclustering approach</th></tr></thead><tbody><tr><td>Soybean single cropping system</td><td>0.65</td><td>0.46</td></tr><tr><td>Soybean-cereal double cropping system</td><td>0.81</td><td>0.59</td></tr><tr><td>Rice-soybean double cropping system</td><td>1.00</td><td>1.00</td></tr><tr><td>Other LULC types</td><td>0.96</td><td>0.85</td></tr></tbody></table>

soybean single cropping system and the soybean-cereal double cropping system are locally concentrated in the south-east border of Tocantins and in the north-east border of the state, respectively, and present a scattered distribution over the rest of state.

## 5. Discussion

Globally, the results demonstrated that the new object-based clustering approach, based exclusively on EO data, can be used to produce field-level maps of cropping systems, described by the crop type and cropping pattern, at a regional scale.

The accuracy assessment results showed that both the hyperclustering and the landscape-clustering variants of the approach produced good classifications, scoring overall accuracies of more than 75%, even if the landscape-clustering approach clearly outperformed the hyperclustering approach. This is particularly evident in the results for the two cropping systems which are spatially dispersed, the single soybean single cropping system and the soybean cereal double cropping system,

the landscape-clustering approach obtaining an F-measure of around 0.2 higher than the hyperclustering approach for these two classes. The Kappa coefficient results further suggest a high reliability of the landscape-clustering classification, which is somewhat lower for the hyperclustering classification.

The significant improvement in the classification accuracy when using a stratification prior to clustering is in accordance with the results of Cano et al. (2017). This result suggests that partitioning the study region into sub-regions, in this case through landscape-level segmentation on the phenological information of a MODIS time series, is an effective spectral variability reduction strategy for improving LULC classification accuracy. This variability reduction strategy seems to particularly benefit the LULC types which are spatially dispersed and under-represented compared to more extensive LULC types since it maximizes the possibility of them being correctly isolated with initial clustering per landscape unit. On the contrary, the hyperclustering approach without stratification, even with a high number of initial clusters (in our case 2400), did not avoid the mixture of classes, thus making its final clusters hard to label and resulting in lower classification accuracies.

The soybean single cropping system obtained the lower accuracies with both approaches (F-measures of 0.65 and 0.46 for the landscape-clustering and the hyperclustering approaches, respectively). This was mostly due to its confusion with pasture areas with strong vegetative development during the rainy season, and thus presenting a high NDVI peak over the same period of the soybean crop's maximum growth stage. Contrastingly, all of the ground-truth objects for the rice-soybean double cropping system were correctly classified with both approaches. This excellent result may be related to the singular phenological pattern of this cropping system, with two distinct and temporary distant crop

<figure><img src="image_90.png" /><figcaption>**Fig. 7.** Comparison of Tocantins' municipal-level soybean area estimates in square kilometres between the remote sensing landscape-clustering results and the IBGE agricultural statistics for 2015.</figcaption></figure>
<figure><img src="image_91.png" /><figcaption>Fig. 8. Map of the main cropping systems of Tocantins at the field level in the 2014–2015 growing season obtained with the landscape-clustering approach.</figcaption></figure>

cycles, which does not resemble the phenological pattern of any other LULC class. Moreover, the landscape-clustering approach did not show an improvement in the classification of this particular cropping system, since all of the fields are localized in the same micro-region (the Rio Formoso micro-region) and so limited spatial variability is expected. The good location of this cropping system might partly explain the lower location disagreement in comparison to the quantity disagreement between the classifications of both approaches and the reference dataset.

The combination of spectral variability reduction steps (field-level and landscape-level segmentations) and a two-step unsupervised classification, appears to be an effective way of reducing the dataset's spectral variability with no need of complex expert knowledge-driven processing steps and external data. Covering a territory such as the Tocantins state (surface area of 277,621 km²) with a coarse spatial resolution dataset such as MODIS at 250 m represents an initial dataset of around 10 million pixels (per MOD13Q1 composite). Following the presented methodology, this high spatial dimensionality was successfully reduced to 187,877 field-level objects, then to 2400 pre-clusters, and finally to 10 clusters. Other than the user-defined segmentation parameters and the number of clusters, the processing is fully automatic, and the low number of final clusters allows labelling through visual analysis, yet requiring basic knowledge of the crop calendar of the cropping systems in order to relate the final clusters to the cropping systems of interest.

The comparison of the landscape-clustering results to the annual soybean, rice, and second crop maize area estimates of the IBGE showed strong positive correlations at the municipality level of around r = 0.9. The slight underestimation of the soybean area in the clustering results

is due to the misclassification of cropland fields with areas of grassland or pasture. These underestimations together with some overestimations are responsible for the dispersion obtained, which in the case of bordering municipalities is particularly evident since crop fields which are actually located inside the mapped state's boundary are sometimes censused on the neighbouring state's statistics. Nonetheless, although the number of clusters was empirically set to simplify the processing, the final classification may benefit from internal evaluation of different combinations of cluster numbers in order to obtain optimal clustering results from a statistical point of view (maximizing the variance between clusters and minimizing the variance within clusters) which may result in a better separability of the misclassified classes.

The inclusion of other processing steps may help improve the classification results at the expense of more computationally-intensive and time-consuming processing. In particular, radiometry of the field-level objects might still be affected by the border effect due to the disparity of spatial resolutions between the MODIS and the Landsat datasets. The presented methodology might thus benefit from spectral unmixing and spatial resampling of the MODIS pixels to the Landsat spatial resolution (e.g. Bisquert et al., 2015; Gevaert and García-Haro, 2015; Knauer et al., 2016). Potential errors also might be induced by the acquisition dates of the Landsat scenes used for the field-level segmentation. For example, some fields might have been divided into two different cropping systems during the growing season but had a homogeneous aspect at the Landsat summer acquisition date. The methodology might thus benefit from multi-temporal segmentation but at the expense of increased complexity and data requirements.

The reproducibility of the presented landscape-clustering approach is currently being tested in the neighbouring Brazilian state of Mato
Grosso for the 2015–2016 growing season and could be reproduced in other states of the MATOPIBA region to provide soybean-based cropping systems’ maps over the whole interstate region. This approach could also be regularly reproduced to provide timely insight on the evolution of the cropping systems in such dynamic agricultural areas, allowing the analysis of land-use trajectories like the ones presented by Spera et al. (2014) and Dubreuil et al. (2012).

Reproducibility of the presented approach in other Brazilian regions or other world regions with large crop field sizes may require adaptation of the scale parameter of the ‘Multiresolution segmentation’ algorithm (Baatz and Schäpe, 2000) to obtain equivalently scaled objects in adequacy to the agricultural landscapes of the targeted region: field-level objects adapted to the region’s fields and landscape units extensive enough to contain multiple field-level objects. The choice of the weight given to the other segmentation parameters may be guided by the indications provided (ex. heavily weighting the ‘shape parameter’ and the ‘compactness’ to obtain the field-level objects). The use of coarse spatial resolution data however prevents the approach from being reproducible in smallholder agriculture regions with fields smaller than a MODIS pixel. In such regions, higher spatial resolution datasets are needed to access the spectrally-pure information of the cropping systems.

## 6. Conclusions

The presented object-based clustering approach, being exclusively based on freely available MODIS and Landsat data, can be potentially reproduced to map cropping systems of other Brazilian regions or other world regions with large field sizes with minimal adaptation. The simplicity of the methods, which neither need external data nor computationally intensive processing demarcates from the data-intensive or highly-expert driven remote sensing approaches which are commonly used when addressing cropping systems’ mapping over large areas. This makes it particularly adapted for cropping systems’ classifications when field data is hard to obtain and external data is lacking.

The novel landscape-clustering approach successfully mapped the main cropping systems of the Brazilian state of Tocantins at the field level by combining two spectral variability reduction steps, a field-level segmentation and a landscape stratification, and an unsupervised classification. The inclusion of a landscape stratification step prior to unsupervised classification significantly improved the classification accuracy, from 78% overall accuracy to 93%, which confirms our initial hypothesis. The classification results of the landscape-clustering approach were strongly correlated to the official crop area-based estimates, with correlation coefficients of around 0.9. Landscape stratification seems to be a significant contribution towards large area cropping systems’ mapping, particularly in regions where cropping systems are spatially spread and under-represented compared to other LULC types.

The produced map can be used as fine-scale complementary information to the agricultural statistics estimates which are generally aggregated at administrative levels, and could eventually be used as an input to crop yield forecasting models. If reproduced yearly, this map could also contribute to the monitoring of cropping systems’ dynamics such as cropping frequency change and their implications for the environment.

## Conflicts of interest

None.

## Acknowledgements

This work was supported by the CIRAD (French Agricultural Research Centre for Development); the FP7 SIGMA (Stimulating Innovation for Global Monitoring of Agriculture) project under grant agreement no.603719; the French Space Agency (CNES) through the

funding of the TOSCA AGRIZONE project; and the CAPES/COFECUB program through the funding of the GeoABC Project (Methodologies and technological innovation for satellite monitoring of low carbon agriculture in support of Brazil’s ABC Plan), project No. 845/15.

## References

Arvor, Damien, Jonathan, Milton, Penello Meirelles, Margareth Simões, Dubreuil, Vincent, Durieux, Laurent, 2011. Classification of MODIS EVI time series for crop mapping in the state of mato grosso, Brazil. Int. J. Remote Sens. 32 (22), 7847-7871. http://dx.doi.org/10.1080/01431161.2010.531783.

Atzberger, Clement, 2013. Advances in remote sensing of agriculture: context description, existing operational monitoring systems and major information needs. Remote Sens. 5 (2), 949-981. http://dx.doi.org/10.3390/rs5020949.

Bégué, Agnès, Damien, Arvor, Camille, Lelong, Elodie, Vintrou, Margareth, Simoes, 2015. Agricultural systems studies using remote sensing. Land Resources Monitoring, Modeling, and Mapping with Remote Sensing. CRC Press, pp. 113-130.

Baatz, Martin, Schäpe, Arno, 2000. Multiresolution segmentation: an optimization approach for high quality multi-scale image segmentation. Angewandte Geographische Informationsverarbeitung XII 58. Wichmann-Verlag, Heidelberg, pp. 12-23.

Barona, Elizabeth, Ramankutty, Navin, Hyman, Glenn, Coomes, Oliver T., 2010. The role of pasture and soybean in deforestation of the Brazilian Amazon. Environ. Res. Lett. 5 (2). http://dx.doi.org/10.1088/1748-9326/5/2/024002.

Bellón, Beatriz, Agnès Bégué, Danny Lo Seen, Claudio Aparecido de Almeida, and Margareth Simões 2017. 'A Remote Sensing Approach for Regional-Scale Mapping of Agricultural Land-Use Systems Based on NDVI Time Series'. Remote Sensing 9 (6). Multidisciplinary Digital Publishing Institute: 600.

Bisquert, Mar, Bordogna, Gloria, Bégué, Agnès, Candiani, Gabriele, Teisseire, Maguelonne, Poncelet, Pascal, 2015. A simple fusion method for image time series based on the estimation of image temporal validity. Remote Sens. 7 (1), 704-724. http://dx.doi.org/10.3390/rs70100704.

Blaschke, Thomas, Geoffrey Hay, J. Hay, Maggi, Kelly, Lang, Stefan, Hofmann, Peter, Addink, Elisabeth, Feitosa, Queiroz, et al., 2014. Geographic object-based image analysis-towards a new paradigm. ISPRS J. Photogramm. Remote Sens. 87, 180-191. http://dx.doi.org/10.1016/j.isprsjprs.2013.09.014. International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS).

Bolfe, Édson L., De C. Vitória, Daniel, Contini, Elisio, Bayma-Silva, Gustavo, Spinelli-Araujo, Luciana, Gomes, Daniel, 2016. Matopiba Em Crescimento Agrícola Aspectos Territoriais E. Revista de Política Agrícola 25 (4), 38-62.

Brown, J., Christopher, Jude H., Camargo Coutinho, Alexandre, de Castro Victoria, Daniel, Bishop, Christopher R., 2013. Classifying multiyear agricultural land use data from mato grosso using time-series MODIS vegetation index data. Remote Sens. Environ. 130, 39-50. http://dx.doi.org/10.1016/j.rse.2012.11.009.

Bunker, Brian E., Tullis, Jason A., Cothren, Jackson D., Casana, Jesse, Aly, Mohamed H., 2016. Object-based dimensionality reduction in land surface phenology classification. AIMS Geosci. 2 (4), 302-328 Aims Press.

Busetto, Lorenzo, Ranghetti, Luigi, 2016. MODISTsp: an R package for automatic preprocessing of MODIS land products time series. Comput. Geosci. 97, 40-48 Elsevier.

Cai, Hongyan, Zhang, Shuwen, Bu, Kun, Yang, Jiuchun, Chang, Liping, 2011. Integrating geographical data and phenological characteristics derived from MODIS data for improving land cover mapping. J. Geog. Sci. 21 (4), 705-718 Springer.

Cano, Emmanuelle, Denux, Jean-philippe, Bisquert, Mar, Hubert-, Laurence, Chéret, Véronique, Cano, Emmanuelle, Denux, Jean-philippe, Bisquert, Mar, Hubert, Laurence, 2017. Improved forest-cover mapping based on MODIS time series and landscape stratification. Int. J. Remote Sens. 38 (7), 1865-1888. http://dx.doi.org/10.1080/01431161.2017.1280635. Taylor & Francis.

Cheema, M.J.M., Bastiaanssen, W.G.M., 2010. Land use and land cover classification in the irrigated indus basin using growth phenology information from satellite data to support water management analysis. Agric. Water Manage. 97 (10), 1541-1552. http://dx.doi.org/10.1016/j.agwat.2010.05.009. Elsevier B.V.

Chen, Jin, Jönsson, Per, Tamura, Masayuki, Gu, Zhihui, Matsushita, Bunkei, Eklundh, Lars, 2004. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter. Remote Sens. Environ. 91 (3), 332-344. http://dx.doi.org/10.1016/j.rse.2004.03.014.

Cihlar, J., 2000. Land cover mapping of large areas from satellites: status and research priorities. Int. J. Remote Sens. 21 (6), 1093-1114. http://dx.doi.org/10.1080/0143116002100902.

Cohen, Jacob, 1960. A coefficient of agreement for nominal scales. Educ. Psychol. Meas. 20 (1), 37-46 Sage Publications Sage CA: Thousand Oaks, CA.

Dubreuil, Vincent, Damien, Arvor, Nathan, Debortoli, 2012. Monitoring the pioneer frontier and agricultural intensification in mato grosso using SPOT vegetation images. Revue Française de Photogrammétrie et de Télédétection 200, 2-10.

FAO. 2009. 'How to Feed the World in 2050'. Insights from an Expert Meeting at FAO 2050 (1): 1-35. 10.1111/j.1728-4457.2009.00312.x.

Foody, Giles M., 2002. Status of land cover classification accuracy assessment. Remote Sens. Environ. 80 (1), 185-201. http://dx.doi.org/10.1016/S0034-4257(01)00295-4.

Fornaro, Alexandre Caselli, 2012. Logística E Agronegócio Globalizado No Estado Do Tocantins: Um Estudo Sobre a Expansão Das Fronteiras Agrícolas Modernas No Território Brasileiro.

Franklin, S.E., Wulder, M.A., 2002. Remote sensing methods in medium spatial resolution satellite data land cover classification of large areas. Prog. Phys. Geogr. 26 (2), 173-205 Sage Publications Sage CA: Thousand Oaks, CA.

Gevaert, Caroline M., García-Haro, F. Javier, 2015. A comparison of STARFM and an unmixing-based algorithm for landsat and MODIS data fusion. Remote Sens. Environ.
156, 34–44. http://dx.doi.org/10.1016/j.rse.2014.09.012.

Guan, Xudong, Huang, Chong, Liu, Gaohuan, Meng, Xuelian, Liu, Qingsheng, 2016. Mapping rice cropping systems in vietnam using an NDVI-based time-series similarity measurement based on DTW distance. Remote Sens. 8 (1). http://dx.doi.org/10.3390/rs8010019.

Gumma, Murali Krishna, Thenkabail, Prasad S., Teluguntla, Pardharsadhi, Rao, Mahesh N., Mohammed, Irshad A., Whitbread, Anthony M., 2016. Mapping rice-fallow cropland areas for short-season grain legumes intensification in south asia using MODIS 250 M time-series data. Int. J. Digital Earth 9 (10), 981–1003. http://dx.doi.org/10.1080/17538947.2016.1168489.

Hartigan, John A., Wong, Manchek A., 1979. Algorithm AS 136: a K-Means clustering algorithm. J. R. Stat. Soc. Ser. C (Appl. Stat.) 28 (1), 100–108 JSTOR.

Hay, Geoffrey J., Castilla, Guillermo, 2008. Geographic object-based image analysis (GEOBIA): a new name for a new discipline. Lecture Notes in Geoinformation and Cartography. pp. 75–89. http://dx.doi.org/10.1007/978-3-540-77058-9_4.

IBGE, 2014. Instituto Brasileiro de Geografia E Estatística, Produção Agrícola Municipal (PAM). http://www.ibge.gov.br/home/estatistica/economia/pam/2014/.

Calli, Jenkerson, 2013. User Guide: Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) on Demand Interface.

Knauer, Kim, Gessner, Ursula, Fensholt, Rasmus, Kuenzer, Claudia, 2016. An ESTARFM fusion framework for the generation of large-Scale time series in cloud-prone and heterogeneous landscapes. Remote Sens. 8 (5). http://dx.doi.org/10.3390/rs8050425.

Labatut, Vincent, Cherifi, Hocine, 2012. Accuracy measures for the comparison of classifiers. In: The 5th International Conference on Information Technology. Amman, Jordanie.

Leenhardt, Delphine, Angevin, Frédérique, Biarnès, Anne, Colbach, Nathalie, Mignolet, Catherine, 2010. Describing and locating cropping systems on a regional scale. A review. Agron. Sustain. Dev. 30 (1), 131–138. http://dx.doi.org/10.1051/agro/2009002.

MacQueen, James, 1967. Some methods for classification and analysis of multivariate observations. Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability 1, 281–297.

McNemar, Quinn., 1947. Note on the sampling error of the difference between correlated proportions or percentages. Psychometrika 12 (2), 153–157 Springer.

Mingwei, Zhang, Qingbo, Zhou, Zhongxin, Chen, Jia, Liu, Yong, Zhou, Chongfa, Cai, 2008. Crop discrimination in northern China with double cropping systems using fourier analysis of time-series MODIS data. Int. J. Appl. Earth Obs. Geoinf. 10 (4), 476–485. http://dx.doi.org/10.1016/j.jag.2007.11.002.

Pedroso da Silva, Antenor Roberto, de Almeida, Maria Geralda, 2007. Agronegócio E O Estado Do Tocantins: O Atual Estágio de Consolidação. Caminhos de Geografia 8 (21), s1.

Pontius, Robert Gilmore, Millones, Marco, 2011. Death to kappa: birth of quantity

disagreement and allocation disagreement for accuracy assessment. Int. J. Remote Sens. 32 (15), 4407–4429. http://dx.doi.org/10.1080/01431161.2011.552923.

R Core Team, 2016. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.

Rosenfield, George H., Fitzpatrick-Lins, Katherine, 1986. A coefficient of agreement as a measure of thematic classification accuracy. Photogramm. Eng. Remote Sens. 52 (2), 223–227.

Silva, Luís Antônio Guerra Conceição, 2007. Biomas presents No Estado Do tocantins. Consultoria Legislativa da Câmara dos Deputados.

Solano, Ramon, Kamel Didan, Andree Jacobson, and Alfredo Huete. 2010. 'MODIS Vegetation Index User's Guide (MOD13 Series)'.

Spangler, Keith R., Lynch, Amanda H., Spera, Stephanie A., 2017. Precipitation drivers of cropping frequency in the Brazilian Cerrado: evidence and implications for decision-making. Weather Clim. Soc. 9 (2), 201–213. http://dx.doi.org/10.1175/WCAS-D-16-0024.1.

Spera, Stephanie A., Avery, S. Cohn, VanWey, Leah K., Mustard, Jack F., Bernardo, Bernardo F., Risso, Joel, Adami, Marcos, 2014. Recent cropping frequency, expansion, and abandonment in Mato Grosso, Brazil had selective land characteristics. Environ. Res. Lett. 9 (6) IOP Publishing: 64010.

Strahler, Alan H., Luigi Boschetti, Giles M. Foody, Mark A. Friedl, Matthew C. Hansen, Martin Herold, Philippe Mayaux, Jeffrey T. Morisette, Stephen V. Stehman, and Curtis E. Woodcock. 2006. 'Global Land Cover Validation: Recommendations for Evaluation and Accuracy Assessment of Global Land Cover Maps'. Scientific and Technical Research Series. Vol. EUR 22156. 10.1080/01431160512331326521.

Tilman, David, Balzer, Christian, Hill, Jason, Befort, Belinda L., 2011. Global food demand and the sustainable intensification of agriculture. Proc. Natl. Acad. Sci. U. S. A. 108 (50), 20260–20264. http://dx.doi.org/10.1073/pnas.1116437108.

Trimble Germany GmbH, 2014. eCognition Developer© 9.0 User Guide. Munich, Germany.

Vintrou, Elodie, Desbrosse, Annie, Bégué, Agnès, Traoré, Sibiry, Baron, Christian, Lo Seen, Danny, 2012. Crop area mapping in West Africa using landscape stratification of MODIS time series and comparison with existing global land products. Int. J. Appl. Earth Obs. Geoinf. 14 (1), 83–93. http://dx.doi.org/10.1016/j.jag.2011.06.010. Elsevier B.V.

Wardlow, Brian D., Egbert, Stephen L., 2008. Large-area crop mapping using time-series MODIS 250 m NDVI data: an assessment for the U.S. central great plains. Remote Sens. Environ. 112 (3), 1096–1116. http://dx.doi.org/10.1016/j.rse.2007.07.019.

Xiong, Jun, Thenkabail, Prasad S., Gumma, Murali K., Teluguntla, Pardhasaradhi, Poehnelt, Justin, Congalton, Russell G., Yadav, Kamini, Thau, David, 2017. Automated cropland mapping of continental Africa using google earth engine cloud computing. ISPRS J. Photogramm. Remote Sens. 126, 225–244 Elsevier.

Zanter, K., 2016. Landsat 8 (L8) Data Users Handbook. LSDS-1574 Version 2.0.