Land use/land cover change detection combining automatic processing and visual interpretation

ABSTRACT This article presents a hybrid classification method combining image segmentation, GIS analysis, and visual interpretation, and its application to elaborate a multi-date cartographic database with 23 land use/cover (LUC) classes using SPOT 5 imagery for the Mexican state of Michoacan (~60,000 km2). First, the resolution of an existing 1:100,000 LUC map produced through visual interpretation of 2007 SPOT images was improved. 2007 SPOT images were segmented, and each segment received the “majority” LUC category from the 1:100,000 map. Segments were characterized from the images (spectral indices) and the map (LUC class). A multivariate trimming was applied to detect “uncertain” segments presenting discrepancy between their spectral response and the LUC class assigned from the map. For these uncertain segments, a category was determined by digital classification, but a definitive category was assigned through visual interpretation. Finally, accuracy of the resulting LUC map was assessed. The same procedure was applied to downgrade (2004) and to update (2014) the map. The implemented method enabled us to improve the scale of an existing 2007 LUC map and to detect land use/cover changes in previous (downgrading) and later (updating) dates with an overall accuracy of 83.3% ± 3.1%.


Introduction
The importance of updated and consistent land use/ cover (LUC) cartography is widely recognized. Land use/cover changes (LUCC) are one of the most evident and fastest processes changing both the characteristics of land surface and the potential of natural ecosystems for provision of environmental services (Herold, Achard, deFries, & Mollicone, 2009). LUCC are one of the main causes of climate change through the modification of carbon, water, and energy cycles. Quality data about LUC and LUCC are crucial to understand current and future dynamics of deforestation and its causessuch as wildfires, urbanization, cropping, loss of biodiversity, and the influence of LUCC on climate change (Herold et al., 2009;Radoux et al., 2014). LUC thematic maps are of great value, because they represent an important tool for designing public policies about urban planning, forest conservancy, risk assessment, and other major concern issues. Remote Sensing data have been widely used in the last decades to elaborate LUC maps (Manakos & Braun, 2014;Millington & Alexander, 2000;Thenkabail, 2015). The increasing availability of remote sensing data and the constant improvement in change detection techniques makes possible to assess dynamics such as forest extent and deforestation (Lu, Mausel, Brondizio, & Moran, 2004). Detailed assessment of deforestation requires high spatial and temporal resolution images, such as Landsat, ASTER, SPOT, and Sentinel (Desclée, Bogaert, & Defourny, 2006;Mas, Lemoine-Rodríguez, & Taud, 2016b). In spite of the improvement of remotely sensed images in the last decades, updating LUC databases is not an easy task (Radoux & Defourny, 2010).
One of the main factors making deforestation monitoring difficult is spectral confusion: different types of LUC can present similar spectral responses and the same type of LUC can present different spectral responses depending on phenology, conservation state, and density. Due to rapid LUCC, frequent updating of maps is needed. Moreover, the elaboration of multi-date cartographic databases is required to assess change (Radoux & Defourny, 2010). Different approaches can be used to elaborate and update existing LUC maps. Visual interpretationoften computer aidedhas been widely used to elaborate LUC cartography (Afrasinei et al., 2017;Asner, Keller, Pereira, & Zweede, 2002;Gurgel, Farias, & Oliveira, 2017;Sunar, 1998), including cartography of large areas, such as Europe (Feranec, Hazeu, Christensen, & Jaffrain, 2007), Africa (Disperati & Virdis, 2015), China (Zhang et al., 2014), and India (Roy et al., 2015). Visual interpretation enables map producers to include many classification criteria, such as texture, shape, pattern, size, and proximity between objects, and the interpreter's knowledge about the study area (Lu et al., 2004;NRSC/ISRO, 2010). Although this technique is time-consuming and requires skilled analysts, visual interpretation is still widely used because it often achieves more accurate results than spectral-based digital approaches (Mas & Ramirez, 1996;Palacio Prieto & González, 1994;Ruelland, Tribotte, Puech, & Dieulin, 2011;Sader, Stone, & Joyce, 1990;Van Den Broek, Smith, & Toet, 2004). When visual interpretation is used to update existing cartography, LUCC are accurately extracted (Disperati & Virdis, 2015;Zhang et al., 2014). On the other hand, fully automatic processing approaches based on digital classification of spectral data (eventually combined with ancillary information) allow for a faster analysis, but they often present more classification errors than visual interpretation. Due to the increasing amount of remote sensing data to analyze, operational monitoring systems require incorporating automated methods. By using an automated procedure, the timeconsuming human interpretation by visual interpretation could be limited to class labeling of identified changed areas.
This paper presents a hybrid semi-automatic classification method combining image segmentation, GIS analysis, and visual interpretation. The method was applied to elaborate a multi-date cartographic database from remote sensing data for the state of Michoacan, Mexico (with an area of about 60,000 km 2 ) based on a classificatory system of 23 LUC classes, in order to produce detailed cartography (1:50,000 scale) and assess LUCC. Finally, two accuracy assessments of the cartography were carried out.

Study area
The state of Michoacan is located in central western Mexico between latitudes 17°55ʹ and 20°24ʹ north, and longitudes 103°44ʹ and 100°04ʹ west (Figure 1), having an area of about 60,000 km 2, representing 3% of the total area of Mexico. The climates are tropical rainy, dry, and temperate (García, 1981). The elevation ranges from sea level to 3,840 m.a.s.l., because of which the vegetation in the state is one of the most diverse in Mexico (Bocco, Méndoza, & Masera, 2001). The temperate region of the state presents oak, pine, fir, and mountain cloud forests, and deciduous and semi-deciduous tropical forests cover the tropical region. These forests have suffered important processes of LUCC, most of them as consequence of the establishment of cropped grasslands and perennial crops, resulting in a complex mosaic of LUC patches (Bocco et al., 2001).

Datasets and software
For this study 32 multispectral SPOT 5 images for 2004, 2007, and 2014 with a 1B processing level and a spatial resolution of 10 m were employed. A LUC map scale 1:100,000 with 23 LUC classes (Table 1) obtained by visual interpretation of the 2007 SPOT 5 images was also used. Field plots from the National Forest and Soil Inventory (INFyS in its Spanish acronym) collected by the National Forestry Commission (CONAFOR in its Spanish acronym) were used to carry out a qualitative assessment of the cartography. These plots, with an area of 1 ha, were selected following a systematic stratified sampling and include information about land cover, tree species, tree height, and tree diameter. Image processing was carried out using the spatial modeling platform DINAMICA EGO (Soares-Filho, Cerqueira, & Pennachin, 2002), BIS Cloud (Berkeley Image Segmentation on the cloud; Berkeley Image Segmentation, 2015), the Geographic Information System QGIS (QGIS Development Team, 2015), the QGIS Plugin AccurAssess , the environment for statistical analysis R (R Core Team, 2013), the R package for decision trees, and rule-based models for pattern recognition "C50" (Kuhn, Weston, & Coulter, 2015), and the R package "np" for nonparametric kernel smoothing methods for mixed data types (Hayfield & Racine, 2008).

The hybrid method
The hybrid classification method is an adaptation from the methods for updating cartography and detecting LUCC proposed by Radoux and Defourny (2010), Radoux et al. (2014), and the interdependent classification method of FAO (1996). This approach is based on the comparison between a map and a satellite image in order to detect discrepancies between map category and spectral response. When the map and image belong to the same date, the method can be applied to improve the resolution of the map. When the image is more recent than the map, the method enables to detect LUCC and to update the map. First, the resolution of a 1:100,000 scale LUC map (Figure 2(a)) obtained by visual interpretation of SPOT 5 2007 images with 23 LUC categories was improved (Table 1). Then, the resulting map was used as base map to produce 2004 and 2014 maps and to assess LUCC. And finally, the accuracy of the resulting LUC maps was assessed. The method is described in detail in the following sections.
Image pre-processing The method requires only few pre-processing steps to prepare the data for the subsequent analysis. It is necessary that all the images that will be used in the process and the LUC base map match geographically. For this study, 2004 and 2014 SPOT 5 images were corrected geometrically using the 2007 SPOT 5 images as reference. 2007 SPOT 5 images were used as reference because the photo-interpreted map was based on those images.
The 2007 SPOT 5 images were segmented by the region growing algorithm, which segments the image by examining the neighboring pixels of a set of points known as seed pointsand determining whether the pixels could be classified to the seed point cluster or not (Adams & Bischof, 1994). This process allowed for obtaining groups of spatially continuous pixels (objects) that are spectrally homogeneous, with a minimum mapping area of one hectare (Figure 2(b)).
Map to segments labeling Afterwards, the 1:100,000 scale LUC map was overlaid with the derived segments of the 2007 SPOT images in order to assign a label (LUC category) from the map to each segment by a majority operation ( Figure 2(c)), consisting in assigning to segments the label of the map category that overlaps the largest proportion of each segment.
Image to map discrepancy detection A density function was calculated from the spectral responses of the segments of each category of LUC. This function indicates the probability for each segment to belong to a certain LUC class as a function of its spectral response (Figure 2(d)). The density of this probability was estimated using the Kernel estimation method proposed by Aitchison and Aitken (1976). By that procedure, the segments presenting a spectral response not matching the LUC category received from the map (outlier) were labeled as uncertain (Figure 2(e)). Additionally, an automatic classification using a C5 tree algorithm (Quinlan, 1993) was applied, which assigned a LUC category to each segment based only on its spectral response. The LUC classes assigned by the automatic classification were used as reference of the spectral confusion between classes during visual interpretation.

Visual interpretation
The segments labeled as uncertain were inspected in pan-sharpened 2007 SPOT images with a spatial resolution of 2.5 m and in very high resolution images from GoogleEarth, using the visual interpretation criteria: texture, color, pattern, size, form, and tone. The discrepancies between the 1:100,000 scale map and the 2007 SPOT image were polygons that were not represented in the map due to scale. The labeling of these segments enabled the interpreter to improve the scale of the map from 1:100,000 to 1:50,000 (Figure 2(f)). Finally, the image boundary sections where SPOT images overlap were manually edited Oak forest (primary) 7 Oak forest (secondary) 8 Fir forest (primary) 9 Fir forest (secondary) 10 Pine forest (primary) 11 Pine forest (secondary) 12 Mountain cloud forest (primary) 13 Mountain cloud forest (secondary) 14 Pine-oak forest (primary) 15 Pine-oak forest (secondary) 16 Tropical deciduous forest (primary) 17 Tropical deciduous forest (secondary) 18 Tropical semi deciduous forest (primary) 19 Tropical semi deciduous forest (secondary) 20 Water bodies 21 Mangrove 22 Herbaceous wetland vegetation 23 Areas without vegetation in order to avoid unreal straight lines on the segments caused by image borders. The same procedure was carried out to generate the 2004 (downgrade) and 2014 (update) LUC maps, using as base map the 2007 1:50,000 scale map, and the SPOT images of each corresponding date were used to obtain segments and spectral values, the discrepancies between each pair of maps corresponds to LUCC areas.

Accuracy assessment
Accuracy assessment of the cartography was applied to the 2007 1:50,000 LUC map, because this was the base map used to generate the 2004 and 2014 maps. 946 reference sites with areas of 1 ha were selected based on stratified random sampling using the map categories as strata. Stratified sampling based on map categories enabled us to select the number of reference sites for each map category. Reference sites were photo-interpreted in pan-sharpened 2007 SPOT images with a spatial resolution of 2.5 m and in very high resolution images from GoogleEarth. In cases where more than one category was present, the category covering the majority of the plot was selected. Fieldwork was carried out to identify sites difficult to interpret in the images.
In the stratified sampling, the number of samples for each map category was not proportional to the area covered by each category on the map. This lack of proportion should be taken into account when computing accuracy indices. According to Card (1982) and Olofsson et al. (2014), the confusion matrix was adjusted by weighing the number of sites using the area of each category on the map. This method enables assessment of various accuracy indices and their certainty (confidence interval). The calculated indices were: (1) overall accuracy (proportion of map correctly classified), (2) user accuracy (one minus the commission error of the category), and (3) producer accuracy (one minus the omission error of the category). Each category assigned by photo-interpretation was compared with the corresponding map category in order to obtain a confusion matrix. In order to carry out a second assessment of the map, the 2007 scale 1:50,000 map was compared with field data from the INFyS. This was a qualitative assessment consisting in comparison of field data with the map, and assessment of the matches and discrepancies. Accuracy indices were not computed from this comparison, because INFyS plots were sampled using as strata a LUC map from INEGI and, therefore, sampling for the map under assessment was not probabilistic (Stehman, 1999).

Results
The 2007 Figure 4). This multi-date data enabled us to quantify LUCC and to detect deforestation hot spots in the state of Michoacan.
According to our cartography, during the last decade, annual deforestation rates in Michoacan have decreased from 0.16% for temperate forests and 0.17% for tropical forests in the 2004-07 period to 0.09% and 0.06%, respectively in 2007-14. Most of the clearing of temperate forests has been due to implementation of crops, mainly perennial cropland, while tropical forests have been replaced by grassland and rainfed crops.
Accuracy assessment was based on 40-50 verification sites by category, except for two categories with a small area. During accuracy assessment, category determination by photo-interpretation of the images was difficult in 110 reference sites that were labeled by ground truth. Among these sites, 23 were incorrectly labeled in the map, 12 were confusions between primary and secondary class of the same category, and 66 were correctly classified. In the raw confusion matrix (Table 2), rows correspond to LUC categories in the 2007 1:50,000 map and columns correspond to true LUC categories (reference plots). The accuracy assessment of the 2007 LUC map indicated an overall accuracy of 83.3% with a confidence interval of 3.1%. Table 3 shows per class accuracy indices.

Discussion
A large number of studies intended to capitalize on the strengths of both visual interpretation and digital classification to map land cover and assess land cover change. For instance, Haag and Haglund (2002) elaborated a land systems map by visual interpretation and carried out independent supervised digital classifications for each of these land systems in order to reduce spectral confusion. Zhou, Cadenasso, Schwarz, and Pickett (2014) went further in this nested approach in order to quantify the fine-scale heterogeneity in urban landscapes. These authors first delineated, by visual interpretation, patches that contain a mix of land cover which served as pre-defined boundaries for finer-scale segmentation and classification. Patches were then classified using the within-   patch proportion land cover categories. Liu et al. (2017) segmented images using different scale parameters and assessed the segmentation results by visual inspection to obtain segments that coincide with landscape features. When updating the database, they used the objects based on previous images as pre-defined boundaries to segment the earlier image. Zhou, Troy, and Grove (2008) carried out an objectbased post-classification change detection using "ifthen" rules which take into account shape, adjacency and from-to category in order to reduce the propagation of attribute and position errors in the two maps used in the post-classification comparison.
The change detection method that we implemented in this study is not based on the comparison of the spectral response of the same object observed in two images acquired in different dates. The methods based on such comparisons need radiometrically corrected images and data of the same (or close) date of the year to avoid differences in the radiometric response that are not due to land change, such as atmospheric effects, differences in illumination and observation angles, and differences in sensor calibration or phenology (Xian & Homer, 2010) . Instead, the method we applied was based on the comparison between a map (or a classified image) and multispectral images, and the detection of change was based on the discrepancy between the spectral responses of different categories of objects in the same image. In the following paragraphs we discuss the potential advantages and limitations of this approach.
The first advantage of comparing a map and a spectral image is that radiometric corrections are not a critical issue. Corrections are necessary only if spectral indices are used, such as the normalized difference vegetation index. For the same reason, the date can vary between images, because the method is more robust in case of change in the phenology of vegetation cover. However, images acquired during certain periods of the year are likely to enable distinguishing particular categories. For instance, images captured at the beginning of the dry season will likely be more useful to discriminate dry tropical forest and pasture land than images taken at the end of the dry season, because at the end of the dry season both covers will not present any vegetative response.
The second advantage is that the method is based on an updating of a previous map, which guarantees temporal congruency. Change detection methods based on post-classification comparison tend to overestimate changes, because due to spectral confusion small spectral variations might lead to classifying the same object into different categories (Mas et al., 2016a). In the case of the present study, the decision of mapping a change is taken through a visual interpretation of the images, which enables taking into account visual interpretation criteriasuch as texture, size, shape, and distribution of the objecttogether with the knowledge of the interpreter. This method is then closer to the interdependent interpretation of multi-date images proposed by FAO (1996), which consists in interpreting the new remote sensing data with reference to the former, therefore ensuring compatibility between historical and recent image classifications, and reducing the error associated to the estimation of changes. The difference between the method we herein propose and the interdependent interpretation is that areas with plausible changes are first detected by image segmentation and trimming. This automated detection step enables the photointerpreter to focus in the areas likely to have changed. It is worth noting that the automated detection step can be fine-tuned to focus in the detection of a particular class and to increase or decrease the number of polygons that will be labeled as uncertain. The number of uncertain polygons can be adapted to the expected changes or to the time available for carrying out the assessment. Increasing the sensibility of the automated step (by increasing the number of polygons labeled as uncertain) will increase the commission error during the first step (polygons labeled as likely to have changed that after visual inspection are finally labeled as not changed) but will decrease the omission error (true changed polygons not detected as likely to have changed and, therefore, notor superficiallyvisually inspected by the interpreter are finally not depicted as change). The hybrid method can therefore be seen as a numerically improved version of the interdependent interpretation method.
In some particular cases, the advantages of the hybrid method can appear as limitations. For example, in the case of categories having similar spectral responses, segments that changed from one category to another will not be detected because they will not be identified as outliers. It is worth noting that such cases of spectral confusion will affect any method based on spectral responses of objects. Using the method we propose here, such changes can be detected through visual interpretation, but the automated step of detection will not be useful. The method will also not be useful in regions that experience a very large amount of change, because segments of change are detected as outliers because they are compared to a vast majority of segments that remain in the same categories over time.
The development of wall-to-wall land cover monitoring system is not a simple task in the cases of large, diverse, and complex countries such as Mexico. An operational system should be able to produce accurate thematic maps and accurate estimations of changes. In the case of the monitoring, reporting, and verification system under the international Reduced Emissions from Deforestation and Degradation plus (REDD+) program, the system should be able to detect subtle changes as forest degradation. The method we propose in this paper enabled us to produce, in a reasonably short period of time, time series of maps exhibiting a high consistency over time and providing accurate figures of deforestation area. However, at the current stage of development, such maps do not give an accurate assessment of forest degradation, and the method must be strengthened with visual interpretation by local expert knowledge and complemented with additional data (field data, transects with LIDAR, or very high resolution imagery).

Conclusions
The classification method we implemented allowed us to detail and generate the 2004, 2007, and 2014 scale 1:50,000 maps in a relatively short time (one month per date by three GIS/RS technicians) and with a high level of accuracy. This was possible because the method benefits from automatic processing (analyzing large amounts of data in short time) and takes advantage of traditional cartographic methods (visual interpretation) in order to supervise the process.
It is important to say that, except for image segmentation, all processing was carried out with open source (R, QGIS) and free software (DINAMICA). However, image segmentation can also be executed in open source software such as, for example, SPRING (Câmara, Souza, Freitas, & Garrido, 1996). Although we obtained the limits of the objects by automatized processing, all cases of LUCC presented in the maps were visually inspected and approved using expert knowledge. Moreover, this interdependent way to update maps allowed us to obtain consistent maps without the false changes usually generated by independent classification. The scale (1:50,000) of the cartography enables detecting deforestation at a detailed level. This is particularly important in countries where, as in the case of Mexico, deforestation mainly occurs in small patches.

Acknowledgments
This study was supported by the projects Monitoreo de la cubierta del suelo y la deforestación en el Estado de Michoacán: un análisis de cambios mediante sensores remotos a escala regional (Fondo Mixto Conacyt-Gobierno del Estado de Michoacán, grant number 192429) and ¿Puede la modelación espacial ayudarnos a entender los procesos de cambio de cobertura/uso del suelo y de degradación ambiental? (SEP-CONACYT, grant number 178816). SPOT images were obtained through the ERMEXS-UNAM agreement. We also thank the experts of the National Institute of Statistics and Geography (INEGI) for their assessment of our cartography. Final corrections of the paper were done by the first author during a sabbatical stay at Universidade Federal da Bahia and Universidade Estadual de Feira de Santana (Brazil) within the project INCT em Estudos Interdisciplinares e Transdisciplinares em Ecologia e Evolução (IN-TREE) with the support of the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (Brazil) and the Programa de Apoyos para La Superación del Personal Académico de La UNAM (PASPA, UNAM, Mexico).

Disclosure statement
No potential conflict of interest was reported by the authors.