JAXA High-Resolution Land Use/Land Cover Map for Central Vietnam in 2007 and 2017

Abstract

Robust remote monitoring of land cover changes is essential for a range of studies such as climate modeling, ecosystems, and environmental protection. However, since each satellite data has its own effective features, it is difficult to obtain high accuracy land cover products derived from a single satellite’s data, perhaps because of cloud cover, suboptimal acquisition schedules, and the restriction of data accessibility. In this study, we integrated Landsat 5, 7, and 8, Sentinel-2, Advanced Land Observing Satellite Advanced Visual, and Near Infrared Radiometer type 2 (ALOS/AVNIR-2), ALOS Phased Array L-band Synthetic Aperture Radar (PALSAR) Mosaic, ALOS-2/PALSAR-2 Mosaic, Shuttle Radar Topography Mission (SRTM), and ancillary data, using kernel density estimation to map and analyze land use/cover change (LUCC) over Central Vietnam from 2007 to 2017. The region was classified into nine categories, i.e., water, urban, rice paddy, upland crops, grassland, orchard, forest, mangrove, and bare land by an automatic model which was trained and tested by 98,000 reference data collected from field surveys and visual interpretations. Results were the 2007 and 2017 classified maps with the same spatial resolutions of 10 m and the overall accuracies of 90.5% and 90.6%, respectively. They indicated that Central Vietnam experienced an extensive change in land cover (33 ± 18% of the total area) during the study period. Gross gains in forests (2680 km²) and water bodies (570 km²) were primarily from conversion of orchards, paddy fields, and crops. Total losses in bare land (495 km²) and paddy (485 km²) were largely to due transformation to croplands and urban & other infrastructure lands. In addition, the results demonstrated that using global land cover products for specific applications is impaired because of uncertainties and inconsistencies. These findings are essential for the development of resource management strategy and environmental studies.

1. Introduction

Land use/cover change (LUCC) is increasingly impacting on the Earth’s surface biophysics, biogeochemistry, and biogeography at any rate or scale such as ecosystem services [1,2,3], water balance [4,5,6,7,8], climate [9,10,11,12,13,14], biodiversity conservation [15,16,17], and agriculture [18]. It means land use/cover information is important for natural resources planning and management [19,20]. In Central Vietnam, land cover has substantially altered as a result of rapid socioeconomic development activities over recent years [21,22,23]. Future changes are also anticipated to occur [24], since the region has an economic growth rate of approximately 10% a year, which is higher than the average of Vietnam [25]. The fast-growing economy has rapidly converted forest and agricultural lands into industrial or service zones [26]. The development also increases the region’s energy requirement, followed by the mass development of hydropower plants because of the Central Vietnam’s suitable geography, topography, and hydrological regime for the hydropower plants. These plants are changing the land cover around them [27]. On the other hand, natural disasters such as drought, floods, and typhoons are also causing land cover changes [28,29].

The LUCC has negatively affected a variety of resources such as biodiversity, carbon sequestration [21,26], and food security [30] over the region. Specifically, a decline of rice yields (by 30%), carbon storage (by 15%), and sequestration (by 12%) due to the expansion of infrastructure lands are predicted until 2100 [31]. The reduction of rice yields results in the concern of food security; Vietnam is the second largest exporter of rice [30]. The expansion of build-up land is projected to have influences on urban heat islands in several cities (e.g., Hanoi). Even though the LUCC insignificantly boosts the peak mean air temperature, the number of hot-spots is growing, particularly in the new infrastructure zones [32]. In addition, Vietnam has seen a gain in forest-cover that is estimated at 1696 million hectares [33], while Central Vietnam has seen a dramatic decrease of forest due to conversion into agricultural land, resulting in the increasing emission of carbon dioxide [34]. The conversions of forest into agriculture also have led to the increase of about 30% in surface runoff and approximately 55% in sediment yield from 2000 to 2008 in Dong Nai province [35,36]. Hence, monitoring of the LUCC is necessary for the sustainable management of natural resources and environment in the region and the achievement of Sustainable Development Goals (SDGs), especially the goal of Life On Land: “Protect, restore, and promote sustainable use of terrestrial ecosystems, sustainably manage forests, combat desertification, and halt and reverse land degradation and halt biodiversity loss [37].”

However, it is not easy to detect the land cover of a region with existing maps, i.e., global land cover maps. They seem to suffer from a low accuracy and coarse spatial resolution. Specifically, Global Land Cover 2000 (GLC2000) [38], GlobCover 2009 [39], the International Geosphere-Biosphere Program Data and Information System’s (IGBP) DISCover land cover [40], and the Moderate Resolution Imaging Spectroradiometer land cover (MCD12Q1) [41] have an overall accuracy of 68.6%, 67.5%, 66.9%, and 78.3% respectively. These maps have a spatial resolution of 300 m to 1000 m. Even with the Finer Resolution Observation and Monitoring—Global Land Cover (FROM-GLC) [42] has a spatial resolution of 30 m—it suffers from a lower overall accuracy (67.08%). This means that global land cover maps designed for particular global purposes are likely to be uncertain and inconsistent for specific local, subnational, or national applications while most counties lack national scale land cover mappings. Mapping land cover for a specific nation with the assistance of users in that nation is designed to meet particular user requirements. This is also the ongoing project of JAXA EORC Ecosystem Research Group, which are producing a country-by-country global land cover map.

The major problem of remotely sensed mappings can be because of insufficient satellite imageries available for land cover mappings, particularly due to suboptimal receiving agendas, limited data accessibility, and cloud cover. In fact, cloud covers 65% of the global surface and 70% of tropical surface in a year [43]. Landsat 5, 7, and 8 were designed to acquire data with a 16-day cycle but most areas have not been constantly imaged every 16 days due to the effects of seasonality, solar zenith angle, cloud cover, and because priority is given to the continental US [44]. The temporal resolution of Landsat 7, apart from the broken scan-line corrector (SLC-off; from 31/05/2003), is frequently longer than 16 days [45]. Data accessibility is also a serious problem, even though Landsat data are freely available via the Earth Resources Observation and Science (EROS) Center, users cannot acquire this data without a brief project description that must be approved [46]. To close these gaps, fusion of multiple remote sensing data with ancillary data is one of the best solutions.

The fusion of optical and radar satellite images (e.g., Landsat and L-band SAR) has recently been proven to be an advancement for monitoring land cover [47,48,49,50,51] and forests [52,53] in tropical areas. With the development of the recent European Space Agency Sentinel-1, -2, and -3 with high spatial resolution, the fusion of multiple sensors is more prevalent and effective. These new data have been effectively combined with Landsat for urban mapping [54,55] with the mosaics of Advanced Land Observing Satellite-2 Phased Arrayed L-band Synthetic Aperture Radar-2 (ALOS-2 PALSAR-2) for mangrove and forest monitoring [56,57]. Integration of multiple optical and radar sensors with different electromagnetic spectra can recognize various land cover features better than a single sensor [58,59]. However, most studies have integrated multiple sensors to map land cover in a certain time instead of temporal land cover changes; or to map a specific land cover type [60]. Fewer studies have estimated the advancement of data fusion among various land cover types, or used high-resolution data fusion, e.g., Advanced Land Observing Satellite Advanced Visible and Near Infrared Radiometer type 2 (ALOS AVNIR-2), for mapping temporal changes in land cover.

The aim of this study was to generate land cover maps over Central Vietnam during the period of 2007 to 2017, using a kernel density estimation and remotely sensed data from multiple sensors. This research shows the potential of combining multiple remotely sensed data and ancillary data for mapping large land cover dynamics. Results improve the understanding of land cover dynamics over Central Vietnam and can contribute to resource management and policy-maker decisions. The results also demonstrate the uncertainties of global land cover products.

2. Materials and Methods

2.1. Study Area

The research site is over Central Vietnam (13°00′–20°00′ N, 105°50′–109°12′ E; Figure 1a) surrounded by the ocean to the east, Laos and Cambodia to the west, Thanh Hoa province in the north, and Phu Yen and Dak Lak provinces in the south. Its total area is approximately 95,000 km² with three main areas: North Central Coast, South Central Coast, and Central Highlands (with highest elevation at 3142 m above sea level). They have a variety of landscapes from deltas, hill lands, mountainous regions or highlands, to coastal zones with a diverse climate from humid subtropical, monsoon to tropical savanna climates. The region experiences four seasons: spring (February to April), summer (May to July), fall (August to October), and winter (November to January). The mean annual rainfall is 700–5000 mm [61] and the mean annual temperature is 23.9–25.9 °C [28], which significantly controls the crop seasons over the region [30]. The diverse climate, complex topography, and various ethnicities lead to a complex geography and landscape with the dominant land cover types of rice paddy, crops, grassland, wetland, urban, forest, bare land, and mangrove.

2.2. Classification Scheme and Reference Data Design

A land cover and land use category system was established to identify the dominant land cover types for certain purposes in the study region [62,63,64]. Based on the local knowledge and the Land Cover Classification System [65], this research used nine land cover types including water, urban & built-up, paddy, other crops, grassland, barren, forest, and mangrove (

Table 1. Description of the land cover categories of Central Vietnam, Vietnam.

Code	Categories	Definition
W	Water	Oceans, seas, lakes, reservoirs, and rivers. They can be either fresh or saltwater bodies.
U	Urban	Land covered by buildings and other man-made structures.
P	Paddy	The cover type is rice paddy and is influenced by the presence of water.
C	Crop	Lands covered with temporary crops followed by harvest and a bare soil period (e.g., single and multiple cropping systems). Note that perennial woody crops will be classified as the appropriate forest or shrub land cover type.
G	Grass	Lands with herbaceous types of cover. Tree and shrub cover is less than 10%.
O	Orchard	An orchard is an intentional planting of trees or shrubs that is maintained for food production.
B	Bare land	Lands with exposed soil, sand, rocks, or snow and that have never had more than 10% vegetated cover during any time of the year.
F	Forest	Lands dominated by woody vegetation with a percent cover >60% and height exceeding 2 m. Almost all trees and shrubs remain green year round. Canopy is never without green foliage.
M	Mangrove	Mangroves are a group of trees and shrubs that live in the coastal intertidal zone.
Os	Others	The other land cover categories

); we also referred to a past paper [66] in order to produce a consistent land cover map for larger scales.

We found that imbalanced stratified random sampling obtains higher accuracy than balanced sampling, which also has an agreement with a past study [67]. This study, therefore, employed an imbalanced stratified random sampling to design reference data for training and testing the classifier. Approximately 3000 reference data were collected from the field survey and about 95,000 were extracted from Google Street View, Degree Confluence Project, Mapillary, and Google Earth by using a visual interpretation (Figure 1b,c). About 65% of the data was used to train the classification model, while the others were used for accuracy assessment. In order to achieve strict training sample standards, every extracted reference data must cover a homogeneous land cover type region with a diameter greater than 20 m. The data contains Geo-location, land cover category, observation time, and the GPS photo if available.

2.3. Data and Image Preprocessing

This study used a variety of multi-temporal satellite imagery from multiple sensors (

Table 2. Dataset organization, layer composition for each sensor type in each dataset, and the total number of images for each position.

Sensor Type	Year of Acquisition	Image/Position	Spatial Resolutions (m)	Temporal Resolution
Sentinel 2	2017	10	10 and 60	10 days
Landsat 8 OLI	2017	8	30	16 days
ALOS AVNIR-2	2007	5	10	46 days
Landsat 7 ETM+	2007	5	30	16 days
Landsat 5 TM	2007	5	30	16 days
ALOS PALSAR Mosaic	2007	1	25	1 year
ALOS-2 PALSAR-2 Mosaic	2017	1
SRTM 1 Arc-Second Global
2000	1	30	-
Open street map	-	1	-	-

). For the 2007 reference year, we used the mosaic images of ALOS AVNIR-2 and ALOS PALSAR of the Japan Aerospace Exploration Agency (JAXA); the calibrated top-of-atmosphere (TOA) reflectance products of Landsat 5 Thematic Mapper (TM) and Landsat 7 Enhanced Thematic Mapper Plus (ETM+) of the United States Geological Survey (USGS). For the year 2017, we used the Sentinel-2 Multi-spectral Instrument (MSI), Landsat 8 Operational Land Imager (OLI), and ALOS-2 PALSAR-2 mosaic. All the data were divided into 24 small square tiles with each of them having a magnitude of 1° # 1° longitude and latitude degrees or less. A bilinear interpolation resampling to World Geodetic System (WGS84) latitude-longitude coordinates with 18.94 s × 18.94 sec resolution (about 10 × 10 m²) was performed for individual data. For optical data, we preprocessed basic atmospheric corrections, cloud masking, and geometric corrections if the data had not been corrected by producers or errors were found, while the radar data showed masked slope effects and filtered speckle; both were used for classification.

For geometric correction, while most of the Landsat L1T and Sentinel 2 products had geometric accuracy of within ± 1 pixel and ± 0.3 pixels, respectively, few images were less accurate. The less accurate images were checked and removed based on control points derived from the Global Land Survey (GLS) 2000 data [68]. The ALOS AVNIR-2 level 1B2 product has geometrically been corrected by JAXA.

For clouds and cloud shadow masking, we used the function of mask (Fmask) for Landsat images. Fmask employed temperature bands to identify clouds at various altitude [69,70]. The elevation of the cloud was later applied in shadow detection. The shadow was identified thanks to the cloud projection and additional determination of shadow pixels. The Sentinel-2 and ALOS AVNIR-2 images, however, did not have such thermal bands that most of the cloud masking methods rely massively on for clouds and cloud shadows masking. They faced difficulty for clouds and cloud shadows masking [70]. Hence, RGB and NIR (and SWIR for Sentinel 2) bands were used to detect cloud; the ratio of blue and green reflectance was employed to identify shadow [71]. Although the method did not use thermal bands, its performance reached similar accuracies to the VIIRS Cloud Mask [72,73] and VIIRS I-Band Cloud Mask [74] methods which used thermal bands [71].

To enhance the accuracy of the optical images, we calculated optical indices that were primarily severed to identify different land types in this research. While the Enhanced Built-Up and Bareness Index (EBBI; Equation (1) [75]), the Normalized Difference Built-Up Index (NDBI; Equation (2) [76]), the Urban Index (UI; Equation (3) [77]), and the Normalized Different Bareness Index (NDBaI; Equation (4) [78]) were employed for distinguishing built-up and bare land. The Normalized Difference Vegetation Index (NDVI; Equation (5) [79]), Enhanced Vegetation Index (EVI; Equation (6) [80]), Soil Adjusted Vegetation Index (SAVI; Equation (7) [81]), and Normalized Difference Water Index (NDWI; Equation (8) [82]) can be used for the sake of promptly monitoring vegetated regions in a complex heterogeneous landscapes to discriminate water, croplands, plantations, and forests [48,52,83,84,85,86].

$$EBBI=\frac{SWI{R}_1-NIR}{10\sqrt{SWI{R}_1+TIR}}$$

(1)

$$EDBI=\frac{SWI{R}_1-NIR}{SWI{R}_1+NIR}$$

(2)

$$UI=\frac{SWI{R}_2-NIR}{SWI{R}_2+NIR}$$

(3)

$$NDBaI=\frac{SWI{R}_1-TIR}{SWI{R}_1+TIR}$$

(4)

$$NDVI=\frac{NIR-VRed}{NIR+VRed}$$

(5)

$$EVI=2.5\frac{NIR-VRed}{NIR+6VRed-7.5VBlue+1}$$

(6)

$$SAVI=1.5\frac{NIR-VRed}{NIR+VRed+0.5}$$

(7)

$$NDWI=\frac{VGreen-NIR}{VGreen+NIR}$$

(8)

where the pixel value used in each channel is the digital number corresponding with each optical satellite imagery; the bands refer to

Table 3. The reference of spectral band information Landsat 5, 7, 8 [87]⁠, Sentinel-2 [88]⁠, and ALOS AVNIR-2 [89]⁠ used for the calculation of optical indices.

Data	Channel	Spectral Range (μm)	Electromagnetic Region
Landsat 8	Band 1	0.435–0.451	Coastal Aerosol
	Band 2	0.452–0.512	Visual Blue (VBlue)
	Band 3	0.533–0.590	Visible Green (VGreen)
	Band 4	0.636–0.673	Visible Red (VRed)
	Band 5	0.851–0.879	Near Infrared (NIR)
	Band 6	1.566–1.651	Short Wave Infrared (SWIR1)
	Band 7	2.107–2.294	Short Wave Infrared (SWIR2)
	Band 10	10.60–11.19	Thermal Infrared (TIR)
Landsat 5 and 7	Band 1	0.45–0.52	Visual Blue (VBlue)
	Band 2	0.52–0.60	Visible Green (VGreen)
	Band 3	0.63–0.69	Visible Red (VRed)
	Band 4	0.77–0.90	Near Infrared (NIR)
	Band 5	1.55–1.75	Short Wave Infrared (SWIR1)
	Band 6	10.40–12.50	Thermal Infrared (TIR)
	Band 7	2.09–2.35	Short Wave Infrared (SWIR2)
Sentinel 2	Band 1	0.433–0.453	Coastal Aerosol
	Band 2	0.458–0.522	Visual Blue (VBlue)
	Band 3	0.543–0.578	Visible Green (VGreen)
	Band 4	0.650–0.680	Visible Red (VRed)
	Band 8	0.785–0.899	Near Infrared (NIR)
ALOS AVNIR-2	Band 1	0.42–0.50	Visual Blue (VBlue)
	Band 2	0.52–0.60	Visible Green (VGreen)
	Band 3	0.61–0.69	Visible Red (VRed)
	Band 4	0.76–0.89	Near Infrared (NIR)

.

The ALOS PALSAR and ALOS-2 PALSAR-2 Mosaic included dual-polarized (HH and HV) channels. In order to increase the effectiveness of land cover classification [90], each channel was filtered speckle using a Refined Lee filter with the European Space Agency Sentinel Application Platform Toolbox v.6.0.0 (ESA SNAP; available online: http://step.esa.int/main/toolboxes/snap/). The digital number (DN) of HH and HV channels were converted into sigma naught values in decibel units (dB) using the following Equation (9):

$${\sigma }^0=10.{\log }_{10}\left(D{N}^2\right)+CF$$

(9)

where σ⁰ is the radar backscatter per unit area, and CF is the calibration factor (CF = 83.0 dB [91]).

For the purpose of anthropogenic activity detection, we took advantage of SRTM 1 Arc-Second Global and OpenStreetMap (OSM) data collected from National Aeronautics and Space Administration (NASA) and OpenStreetMap Foundation (OSMF), respectively. The SRTM 1 Arc-Second Global data were used to extract topography information, such as slope, whereas the level of human activity was identified by extracting a rasterized distance map to road network using the OSM.

All the preprocessing steps were automatically completed by employing C++ and python with the support of the Geospatial Data Abstraction Library (GDAL; available online: https://www.gdal.org/), Geographic Resources Analysis Support System, Geographic Information System v.7.2 (GRASS—GIS; available online: https://grass.osgeo.org/), and Quantum Geographic Information System v.2.18 (QGIS; https://www.qgis.org/en/site/).

2.4. Classification Method

Our overall flowchart included three main stages: image preprocessing, image classification and accuracy assessment, and change analysis (Figure 2). Following a previously described method [92], this study employed Bayesian logic together with a kernel density estimation (KDE [93]). More specifically, we estimated the probability density function of features (such as vegetation indices, band reflectance, etc.) for each category by generating small Gaussian functions around training data (which is KDE) in the feature space, and combined them to generate posterior probability using Bayes’ theorem. After the computation of the posterior probability of every land cover category for each image, we generated the joint posterior probability from all the overlapped images acquired from different periods and sensors. The final choice of category was based on the highest joint probability among all land cover categories. This method is suitable for detecting the seasonal change of land surface (phenology) as a classification key in a large area due to its fully-automatic robustness. More importantly, this KDE approach is more accurate than support vector machines and the maximum likelihood classification [92]. We offer a detailed description of the processes as follows:

The posterior probability of a category was computed based on a D-dimensional vector of input data x, i.e., spectral bands of reference and two-dimensional values representing observation date [t₁, t₂] (Equation (10)).

$$[t_1,t_2]=[\cos (2\pi \frac{DOY}{DO{Y}_{max}}),\sin (2\pi \frac{DOY}{DO{Y}_{max}})$$

(10)

where DOY is the date of the observation year (Julian day), and DOY_max (=365.25) is the maximum day in the year.

For every image, the posterior probability of a category C_k (k = 1, 2, …, M; M is the number of land cover categories; M = 9) was determined by using the Bayesian rule based on the input data x (Equation (11)).

$$p\left(C_k\mid x\right)=\frac{p\left(C_k\right)p\left(x\mid C_k\right)}{p\left(x\right)}=\frac{p\left(C_k\right)p\left(x\mid C_k\right)}{{{\displaystyle \sum }}_{k=1}^{M}p\left(C_k\right)p\left(x\mid C_k\right)}$$

(11)

where p(C_k) is the prior probability of C_k (which is assumed to be a uniform distribution), and p(x|C_k) is a category-conditional probability of x; p(x|C_k) was estimated based on the training data using kernel density estimation (KDE). KDE is used to compute the probability distribution of data as the sum of kernel functions that are of the same form and centered on each training data, by using the Gaussian kernel (Equation (12)) and Scott’s rule of thumb (Equation (13)) as follows.

$$p\left(x\mid C_k\right)=\frac{1}{N_k}{\displaystyle \sum }_{n=1}^{N_k}\{{\displaystyle \prod }_{d=1}^D\frac{1}{h_d}K(\frac{x_d-x_{n,d}}{h_d})\}$$

(12)

$$K\left(u\right)=\frac{1}{\sqrt{2\mathsf{\pi }}}\exp (\frac{-u^2}{2})$$

(13)

$$h_d=N^{-1/\left(D+4\right)}.{\mathsf{\sigma }}_d$$

(14)

where, N_k is the number of training data of a category C_k, h_d is a bandwidth parameter estimated by Equation (14), N is the total number of training data (N = N₁ + N₂ + ... + N_M), and σ_d denotes the standard deviation of d-th dimension of training data {x_n,d|1 ≤ n ≤ N}.

In the next step, at each pixel for each category, we integrated the posterior probability of all overlapped images by multiplying the posterior probability of all images. In general, the higher the joint posterior probability of a category, C_k, the more possible the land cover category C_k is. However, in reality, even if the true land cover is category C_k, the p(C_k|x) of one image sometimes might be close to or equal to zero due to noise, cloud, or insufficient training data. If it occurs, it would make the joint of the posterior probability of category C_k also close to or equal to zero, because multiplying by zero always gives zero. This means that even if the p(C_k|x) of the most overlapped images are as high as 1, the final prediction of the land cover cannot be in the category C_k. To overcome this issue, the posterior probability of each image should not be too close to zero. To realize it, we used Equation (15) [94]. The final probability of a category C_k (p’(C_k)) was estimated by Equation (16).

$$p^{\prime }\left(C_k\mid x\right)=ap\left(C_k\mid x\right)+\frac{1-a}{M}$$

(15)

$$p^{\prime }\left(C_k\right)={\displaystyle \prod }_{i=1}^S{p}_i^{\prime }\left(C_k\mid x_i\right)$$

(16)

where a is a constant value (a = 0.7), and S is the number of images.

The final choice of a category is the category with the highest joint probability among all land cover categories. Supposedly, at a pixel r of a classified land cover map with two categories, i.e., water and urban, has the joint probability of water: p’(C_water) = 0.6 and the joint probability of urban: p’(C_urban) = 0.4. The highest joint probability of pixel r is 0.6 and the land cover of pixel r is water.

To shorten data processing and complex landscapes, we extracted imagery data into 24 small square tiles, each of them having a magnitude of 1° # 1° with longitude and latitude degrees being lesser (Figure 2). The classification process was performed independently for each tile using Saclass Software version 1.7 developed by the University of Tsukuba and JAXA [94,95,96].

2.5. Accuracy Assessment

Assessment of classification accuracy of 2007 and 2017 maps was the most important part in order to determine the quality of these maps. While there was a large number of accuracy measures used to estimate algorithm performance, it is crucial to carry out accuracy estimation for any category if the classification results are valuable for changing detection [97]. We, therefore, chose a confusion matrix [98] because it is one of the most general methods that is easy to understanding and has useful values [99]. An imbalanced stratified random approach was used to extract approximately 27,000 reference data from ground truth data and visual interpretations. These reference data represented all the designed land cover categories of the region. In addition, a Kappa coefficient of agreement was used to measure the range of single classification accuracy [100]. It is presented by the following formula (Equation (17)).

K = [P_o − P_e]/[1 − P_e]

(17)

where P_o is the observed proportional agreement between actual and predicted categories defined as $P_o=\frac{1}{n}{\displaystyle \sum }_{i=1}^g{f}_{ii}$ and P_e is the expected agreement by chance defined as $P_o=\frac{1}{n^2}{\displaystyle \sum }_{i=1}^g{f}_{i+}{f}_{+i}$, where fi+ is the total for the i th row and f+i is the total for the i th column in the confusion matrix.

3. Results

The classified maps and areas of land cover change within the 10-year period are shown in Figure 3 and their accuracy assessments (confusion matrix) are shown in

Table 4. Accuracy assessment of the 2007 land cover maps over Central Vietnam, using a confusion matrix.

Predicted Category
Actual Category		W	U	P	C	G	O	B	F	M	Total	PA (%)
	W	644	0	4	0	0	0	0	0	11	659	97.8
	U	0	1005	3	64	31	76	2	0	0	1181	85.1
	P	0	6	1138	12	3	2	15	2	58	1236	92.1
	C	0	12	8	1123	73	105	6	6	2	1335	84.2
	G	0	6	0	4	377	5	1	0	0	393	96
	O	0	8	2	34	11	637	0	45	0	737	86.5
	B	0	7	0	1	1	1	496	0	0	506	98.1
	F	1	0	1	0	1	60	0	986	0	1049	94
	M	23	0	16	0	2	0	0	0	495	536	92.4
	Total	668	1044	1172	1238	499	886	520	1039	566	7632	91.8
	UA (%)	96.5	96.3	97.1	90.8	75.6	71.9	95.4	94.9	87.5	89.6	90.5
	Ka	0.02	0.05	0.02	0.02	0	0	0.01	0.03	0	0.15	0.9

UA: Accuracy for users; PA: Accuracy for producers; K_a: Kappa coefficient; W: Water; U: Urban; P: Paddy; C: Crop; G: Grass; O: Orchard; B: Bare land; F: Forest; M: Mangrove.

and

Table 5. Accuracy assessment of the 2017 land cover maps over Central Vietnam, using a confusion matrix.

Predicted Category
Actual Category		W	U	P	C	G	O	B	F	M	Total	PA (%)
	W	2636	2	1	0	4	0	0	0	2	2645	99.7
	U	0	1799	1	70	58	75	17	0	0	2020	89.1
	P	5	0	2971	17	0	11	0	28	19	3051	97.4
	C	2	33	157	2454	181	171	0	28	48	3074	79.9
	G	84	74	7	107	1821	135	0	46	7	2281	79.9
	O	1	12	4	153	72	765	0	54	8	1069	71.6
	B	12	60	5	16	7	3	1264	0	0	1367	92.5
	F	0	0	8	3	1	10	1	3291	0	3314	99.4
	M	7	3	5	2	2	1	1	22	781	824	94.8
	Total	2747	1983	3159	2822	2146	1171	1283	3469	865	19645	89.4
	UA (%)	96	90.8	94.1	87	84.9	65.4	98.6	94.9	90.3	89.1	90.6
	Ka	0.02	0.01	0.02	0.02	0.01	0	0	0.03	0	0.13	0.9

UA: Accuracy for users; PA: Accuracy for producers; K_a: Kappa coefficient; W: Water; U: Urban; P: Paddy; C: Crop; G: Grass; O: Orchard; B: Bare land; F: Forest; M: Mangrove.

. The overall accuracies of the maps for 2007 and 2017 are 90.5% (kappa coefficient of 90%) and 90.6% (kappa coefficient of 90%), respectively. Most categories have accuracy for users and producers greater than or close to 90%, except for grass and orchards. Water, bare land, paddy, and forest have the highest accuracies that are over or close to 95%, followed by urban and crops that account for approximately 91% and 90%, respectively. Orchards and grassland have the lowest accuracies (<85%) in the two maps. The reason for misinterpreted classification of orchard and grassland may be the correspondent spectral characteristics between orchards, grass, and the other categories.

Central Vietnam experienced an extensive change in land cover from 2007 to 2017 (Figure 3c). A total of 31,380 ± 16,920 km² (33 ± 18% of total area) underwent changes with the major changes occurring in orchards, forests, and croplands in the coastal areas and central highlands. In 2007, central Vietnam covered an estimated area of about 94,000 km² with 93% of the total area being vegetated areas and the other being water bodies, bare land, and urban & built-up. The vegetated area included dense tree cover, i.e., forest (about 51,000 km²) and mangrove (1000 km²), and dynamic land cover, i.e., paddy fields (6500 km²), grassland (2800 km²), cropland (9900 km²), and orchards (15,500 km²). Over the recent decade, water, urban & built-up, cropland, grassland, forest, and mangrove areas increased by approximately 560, 40, 1100, 2680, and 930 km², respectively, while bare land and paddy fields decreased the similar amount of about 500 km². Surprisingly, the period witnessed a sharp decline in orchards (4600 km²).

4. Discussion

It is not easy to create an accurate map for a cloudy area such as Central Vietnam, a tropical region, where clouds cover 70% of the area over the year [43]. Our methods were customized to Vietnam’s unique conditions by integrating a variety of satellite images from multiple sensors with ancillary data. This is an effective approach for land cover mapping for the cloudy and large region as the result of the availability of free satellite images and the development of image processing and computational power. For instance, by using an automated image preprocessing approach, with an automated classification method, and support of servers, we could produce a 10 m spatial resolution land cover map over Central Vietnam (94,000 km²) within approximately 3 days. Much more effort was frequently required for the careful work of data organization, field survey, and land cover change estimation. Further effectiveness and efficiency of this approach can be explained as follows.

Our research showed the potential of combining multisensor remote sensing data for land cover classification and change detection in the tropical area. The outperforming of data fusion (e.g., Landsat and L-band SAR) over individual sensors for improving overall accuracies has been used for land cover monitoring in tropical regions such as in Indonesia [49,51] and West Africa [50], which obtained similar overall classification accuracies. Nevertheless, these researches mapped land cover for smaller areas (<900 km²) compared with the current study (94,000 km²). Note that this study corroborated the research by Hoang et al. [66], which used the combination of multisensor data (i.e., Landsat 5 and 8 and ASTER-VA version) for analyzing land cover in Northern Vietnam that demonstrated the potential of the fusion of multiple sensors for mapping heterogeneous landscapes. However, our maps are better than their maps in term of the overall accuracy and spatial resolution (90.5% and 10 m vs. 81.0% and 15 m). A possible explanation for this might be that we used finer spatial resolution images (10 m against 15 m), more reference data (98,000 against 65,000), and the improvement of posterior probability integration (Equation (15)). Another explanation is that, instead of using 4 bands (Blue, Green, Red, and Infrared), we used the best combination of bands and a set of spectral indices.

Surprisingly, while most categories have users’ accuracy (UA) and producers’ accuracy (PA) higher than 90%, few limitations can be found in the accuracy of grass, crops, and orchard categories (Table 4 and 5). Specifically, grass and orchard have the lowest UA while misclassification between grass, crop, and orchard has occurred for both years. A possible explanation for this might be that orchard and crop categories include a large variety of orchard and crop types leading to a significant variance of spectral reflectance patterns. Another probable reason may be that several grasslands for raising cattle are cultivated as croplands (temporary crops followed by harvest and a bare soil period), resulting in confusion over whether they are grassland or cropland. The misclassification between crop, orchard, and urban may be due to the specialty of the traditional Vietnamese farm, a form of domestic agriculture in which food gardening, fish rearing, and animal husbandry are wholly combined [101,102]. However, these systems are frequently used in small complex areas of mixed land cover/land use, causing difficulty in distinguishing land-cover. In order to overcome the misclassification, it may be worth considering the use of very high-resolution remotely sensed data (e.g., Ikonos, QuickBird, and Kompsat-2 [103,104]) or classifying these categories into different subcategories before integrating them into a single one. The issue of the misclassification between the mixed and complex land covers is an intriguing one which could be explored in further research.

4.1. Uncertainties of Global Land Cover Maps over Central Vietnam

With a spatial resolution of 10 m and accuracy of 90.5% (Kappa coefficient: 0.9), our maps are better than the existing maps which are considered to be of coarse spatial resolution (30–1000 m) and low accuracy (<80%; Figure 4). Figure 5 compares our maps with (a) Climate change initiative (CCI) 300-m land cover V2 for the year 2015 released by the European Space Agency (ESA), (b) the GlobeLand30 map for the year 2015 published by the National Geomatics Center of China, (c) the MCD1Q1 0.5 km MODIS-based global land cover climatology for the year 2001–2010 published by the USGS, (d) the Global PALSAR-2 25-m Forest/Non-forest map for the year 2017, and (e) the Global PALSAR 25-m Forest/Non-forest map for the year 2007 from JAXA based on visual interpretation.

Results show that global land cover maps seem to be uncertain and inconsistent. To be specific, the CCI map tended to misclassify most inland water, whereas underestimated urban & built-up areas and overestimated croplands (Figure 5a). Although the GlobeLand30 could detect well inland water, it was likely to underestimate urban & other infrastructure lands and overestimate grassland (Figure 5b). These issues may be the result of using coarse spatial resolution satellite images, cloud cover, or reference data shortcoming. The MCD1Q1 0.5 km MODIS-based global land cover climatology had a tendency to misinterpret most inland water as wetland regions and overestimate cropland (Figure 5c), perhaps because of the difference of land cover type definition or the use of very coarse spatial satellite imagery (500 m). For forest estimation, based on Google Earth view and the 2017 forest map in this research, we found that FNF maps probably misclassified some forest areas and could not accurately detect inland water, since a large number of reservoirs disappeared in the map (Figure 5d,e). These problems may be the result of using only SAR images (ALOS PALSAR or ALOS-2 PALSAR-2). Although the SAR images are not blocked by clouds or cloud shadows, it often suffers from speckle which can be reduced by using noise reduction filters, however still constraining classification accuracies [105,106,107]. In summary, global land cover maps contain large uncertainties for environmental studies.

4.2. Ten-Year Land Cover Change over Central Vietnam

Central Vietnam has a heterogeneous landscape that experienced rapid and extensive changes during the period of 2007 to 2017. To observe the mass conversion of land use in Central Vietnam due to recent socioeconomic transformation [108,109], three testing locations were chosen to detect qualitative changes between 2007 and 2017. The testing locations are sites A, B, and C in Thua Thien Hue, Quang Nam, and Thanh Hoa provinces, respectively (Figure 6). In site A, a number of reservoirs have been constructed, which is a common phenomenon over Central Vietnam. These reservoirs converted orchard to water surface while a number of neighboring forests changed to crops. This finding is in agreement with other satellite analysis [110,111] and can explain why orchards decreased while croplands increased over the recent decade. Site B illustrated a shift from croplands to forests, probably the result of recent government policy to reforest some parts of Vietnam by providing financial and technical resources [112,113]. This forest gain also agrees with other satellite analysis [23] and demographic statistics [114] showing the forest area increased by 1.696 million hectares on the national scale from 2005 to 2015. Site C presents the change of paddy to crops or urban & built-up areas. This could explain the decrease of paddy fields and the increase of croplands in the region. Another reason for the decrease of paddy fields can be from conversion to aquaculture because of decreasing rice productivity as the result of the intrusion of saltwater [27].

The region has been experiencing extensive changes, particularly in the decrease of paddy fields and the increase of inland water surface that can be detected easily based on this study. These changes may have generated unprecedented new ecosystems that impact environmental sustainability and food security. Results show that more than 21 huge dams have been constructed at upstream rivers (e.g., Huong, Vu Gia—Thu Bon, Dong Nai, and Sre Pok) and many more are now planned. These dams can block suspended sediment from upstream areas, which may cause large-scale shoreline erosion and land loss. Also, the construction of upstream dams restricts downstream river flow leading to a decrease in water level at estuaries, while the sea water level is expected to rise [115,116,117]. This can also result in severe erosion and intense saltwater intrusion in lowland areas, followed by the expansion of salinity effects on plant growth and yield such as rice [118,119], and the conversion of rice to aquaculture or other lands leading to the decrease of rice productivity. Because Vietnam is the second largest exporter of rice, domestic food production and international rice trade are likely to be at risk unless a sustainable development strategy is considered in the near future.

4.3. Potential Application and Future Work

Even though the fusion of remote sensing data has been applied in recent LULC mapping studies [49,50,51,66], our study has made significant contributions to the field. As in the above discussions, Central Vietnam, a cloudy area, has been a challenging area for mapping with LULC, especially with high spatial resolution, due to the scarcity of data availability, which can be customized by our approach. Such high spatial resolution maps of the heterogeneous and large-scale areas have rarely been carried out, particularly in the year before 2015 (when Sentinel-2 was launched). This study made use of ALOS AVNIR-2 to generate a past-time 10-m LULC map, which may serve as a baseline map to compare observed changes or critical data for local or national long-term land use planning.

The findings and maps presented in this research can be used for multipurpose applications. First, the construction of hydroelectric reservoirs can result in problems such as shoreline erosion, salinity intrusion, extreme water level variations, and sediment delivery issues, which have been occurring in the region. Also, the region is sensitive and vulnerable to the influences of climate change and consequent sea level rise. These problems link to land cover/land use changes. While several mapping projects have been using coarse resolution satellite imagery for land cover change detection, they may not be effective for a complex and fragmented landscapes such as Central Vietnam [120,121], and our maps seem to be more suitable for land cover change analysis over Central Vietnam, and can provide the policy-makers and scientific associations with input data for the further discussion of environmental management, in particular water balance, sediment estimation, and food security. In addition, the maps can be a critical data for managing ecosystems and biodiversity such as for Global Forest Resources Assessment (FRA) of FAO. They may also serve as baseline maps for other land cover/land use projects such as the “Land Use Status, Change, and Impacts in Vietnam, Cambodia, and Laos” of NASA. Although our maps achieved a certain level of overall accuracy, they may be insufficient for a quantitative analysis of changes with a certain level of statistical significance. In order to improve the accuracy, new satellite data (e.g., Sentinel-1), or full polarimetric SAR data should be considered for next steps. On the other hand, the combination of this approach with others may be necessary. Finally, due to the uncertainties of global land cover products, it may be better to create individual national-scale maps and combine them into a global map instead of generating a whole global map.

5. Conclusions

Based on the kernel density estimation, we produced land cover maps of the over Central Vietnam between 2007 and 2017 using high-resolution remotely sensed data from multiple sensors. These maps have a spatial resolution of 10 m and an overall accuracy of 90.6% (kappa coefficient 0.9). This accuracy and spatial resolution are higher than that of existing land cover maps which tend to have a coarse resolution (30 m to 1000 m) and low accuracy (<80%), causing uncertainties for users. This study indicates the potential of multisensor fusion for monitoring land cover dynamics in a cloud and large area.

Our results show that although global land cover products are fundamental variable for global specific applications, there remains a considerable amount of uncertainties and inconsistencies for particular applications at local and national scales which could be solved by using products in this study.

Anthropogenic pressures on the land cover system over Central Vietnam are growing because of the rapid socioeconomic development process. Over the recent decade, forest areas have significantly expanded due to government efforts to reforest by changing policy and providing technical resources. However, the quality of forest in central Vietnam still remains a mystery, and a major concern in ensuring forest management in Vietnam or the Sustainable Development Goal (SDG) 15.1.1. Urban & other infrastructure areas have expanded around crowded cities such as Thanh Hoa, Vinh, Hue, and Da Nang due to population growth and the movement of citizens from rural areas to urban regions. Population growth is also accompanied by an increasing demand for water, domestic and industrial irrigation, and hydropower, resulting in the expansion of inland water surface. These changes may damage environmental sustainability, particularly by shoreline erosion, land loss, and salinity intrusion. The findings of land cover dynamics together with an interpretation of driving factors can provide policy-makers and scientific associations with appropriate input data for the further discussion of land environment management.

For further applications or other interests, readers can refer to the supplementary materials or/and download the land cover map results in this study on the JAXA/EORC website: http://www.eorc.jaxa.jp/ALOS/en/lulc/lulc_vnm_v1807.htm.