Mapping of Vegetation Using Multi-Temporal Downscaled Satellite Images of a Reclaimed Area in Saemangeum, Republic of Korea

Abstract

The aim of this study is to adapt and evaluate the effectiveness of a multi-temporal downscaled images technique for classifying the typical vegetation types of a reclaimed area. The areas reclaimed from estuarine tidal flats show high spatial heterogeneity in soil salinity conditions. There are three typical vegetation types for which the distribution is restricted by the soil conditions. A halophyte-dominated vegetation is located in a high saline area, grass vegetation is found in a mid- or low saline area, and reed/small-reed vegetation is situated in a low saline area. Multi-temporal satellite images were used to classify the vegetation types. Landsat images were downscaled to take into account spatial heterogeneity using cokriging. A random forest classifier was used for the classification, with downscaled Landsat and RapidEye images. Classification with RapidEye images alone demonstrated a lower level of accuracy than when combined with multi-temporal downscaled images. The results demonstrate the usefulness of a downscaling technique for mapping. This approach can provide a framework which is able to maintain low costs whilst producing richer images for the monitoring of a large and heterogeneous ecosystem.

1. Introduction

The accomplishment of many reclamation projects has changed tidal flats into agricultural areas on the West Sea coast of the Korean Peninsula. The Saemangeum is one of the largest reclamation projects in the Republic of Korea. The reclaimed areas have typical ecological characteristics. The soil salinity of the land is usually high, which leads to harsh conditions for plants, restricting the colonization of new plants and the development of vegetation. Until the soil salinity reaches a normal level, it is a strong driving force in regulating the distribution of vegetation communities within the area. The soil salinity in the area is continuously changing and heterogeneous [1,2], and as a result, it shapes the unique and dynamic characteristics of the ecosystems.

To understand the development and change of these vegetation communities, it is critical to determine the distribution of vegetation types that have different resistances to salinity. Studies of vegetation in reclaimed areas are relatively rare [1,2,3,4,5] in Korea compared to studies in salt marshes [6,7,8,9,10]. The studies have primarily focused on the relationship between the soil conditions and physiological traits of plants in reclaimed areas. Lee et al. [4] surveyed four tidal reclamation project areas and concluded that the distribution of plants was strongly restricted by the soil salinity. The clear relationship between the soil salinity and the distribution of vegetation was also found in the Saemangeum area [1,3]. Kim et al. [3] classified the vegetation types of the Saemangeum area as: (1) halophyte vegetation located in a high saline area (mean electrical conductivity (EC): 14 dSm-1), dominated by Salicornia europaea, Suaeda asparagoides, and S. japonica; (2) mixed halophyte vegetation located in mid saline areas (mean EC: 6.7 dSm-1), dominated by Phragmites communis, Puccinellia nipponica, and Carex pumila; and (3) low-saline vegetation located in low saline areas (mean EC: 3.0 dSm-1), dominated by diverse plants, including Aster subulatus, A. tripolium, and Echinochloa spp. Despite the physiological relationships, there is a limited understanding of the ecosystem conditions in a spatial context because these studies were usually conducted at experimental sites [1,2,3,5].

Remote sensing is an effective method for estimating or mapping target quantities in a large area. It is cost-effective, rapid, and able to cover a large area. Remote sensing has also been extensively used to estimate vegetation in salt marshes [6,7,11,12,13,14,15,16,17], but there is little research on using remotely sensed data for reclaimed areas [18]. In salt marsh studies, many researchers have tackled classifying halophyte plant and/or vegetation types that share many common characteristics of vegetation in an area reclaimed from estuarine tidal flats. Some studies have focused on (multi-temporal) spectral differences [6,7,11,16,17], whilst others have concentrated on classification methods such as Neural Networks Classification [12] or the random forest method [19].

For the application of remote sensing techniques to be effective, remotely sensed data should have a strong connection to the target quantity or characteristics, in addition to the proper resolution, in order to take into account spatial variability. Vegetation in the Saemangeum is dominated by grass/herbaceous plants, and the biomass is low. The vegetation types show similar spectral responses, which makes it difficult to classify the types using spectral characteristics alone. Additionally, the vegetation patches are diverse in terms of size and shape. The patch sizes vary from several m² to thousands of m², and the shapes are irregular. To take into account spatial heterogeneity, it is important to select the proper resolution of the image.

It is well known that the phenological traits of vegetation can be used to improve the accuracy of vegetation or plant species mapping [11,20,21,22]. The plants in the study area primarily appear after the rainy season (mid-June to mid-July), due to the dry soil conditions in the spring. The biomass reaches a maximum in late-August to mid-September. Gilmore et al. [11] found that information on phenological variability in the growing season was useful for distinguishing dominant marsh plant species. However, it is difficult to obtain multi-temporal images with a high resolution during the short growing season. If images with short revisit frequencies and coarse resolutions can be downscaled (i.e., increased in spatial resolution), the utility of the images can be effectively increased.

The downscaling of imagery has gained attention in recent years, and many approaches have been developed (see Atkinson [23]). Among them, the use of cokriging for image sharpening provides several advantages. Cokriging has a well-established theoretical model and uses semivariograms and cross-variograms of two or more images of different resolutions [24,25,26]. Because cokriging is a generic tool, it is possible to include diverse types of data, such as topographic maps, thematic maps, or experimental data. Additionally, the downscaling cokriging method has demonstrated a lower error (mean error and mean squared error) than the method of using a high pass filter [24].

The objective of this study was to propose a modified approach to map vegetation in a reclaimed area using multi-temporal downscaled images. This study used cokriging methods to downscale Landsat imagery to the resolution of a RapidEye image. This study provides an effective method for creating an accurate vegetation map, which is essential for monitoring and managing the ecosystem of the reclaimed Saemangeum area.

2. Materials and Methods

2.1. Study Area

The Saemangeum is a reclaimed area that was created by the construction of a sea dike (33 km) across an estuarine tidal flat on the West Sea coast of South Korea (completed in April 2010). The dike encloses a total area of 401 ${\mathrm{km}}^2$ (281 ${\mathrm{km}}^2$ of land and 120 ${\mathrm{km}}^2$ of fresh water reservoir). Although there are many reclaimed areas in Korea, the Saemangeum is unique because of its large scale and relatively long development period. These conditions allow large areas to remain under natural conditions, with very limited anthropogenic disturbance to the vegetation. It is worth noting that the soil is not imported from outside of the reclamation area, due to its large scale. The reclaimed land was created by lowering the water elevation level (EL). Since the end of 2010, the EL has been maintained at about −1.5 m (below sea level, Figure 1). At the EL, a large area of estuarine tidal flat (ca. 180 ${\mathrm{km}}^2$) became land [27], while some portions of the area remained natural. The study area was selected, based on the criteria of: (1) minimal anthropogenic disturbance; and (2) heterogeneity and typical vegetation type of the reclamation area (Figure 1).

2.2. Collecting Field Data

The study area was approximately 3.5 km². A total of 127 sample plots (each plot: 5 $\mathrm{m}\times 5\mathrm{m}$) were selected by simple random sampling (Figure 2) within the study area. A field survey was conducted in September 2014. We recorded the coverage of each species and dominant vegetation type within 25 m² of the sample points and classified the typical vegetation types. There were three typical vegetation types: halophyte vegetation (HV), mid- to low saline vegetation (MLV), and reed/small-reed vegetation (Reed, Figure 3). A total of three randomly selected soil core samples (0–5 cm depth) were collected from each sample point. The soil samples, which were air dried, crushed, and uniformly mixed, were sifted through a 2 mm sieve for EC and pH analysis and a 0.5 mm sieve for organic matter content analysis. The soil EC and pH were measured in saturated paste extract and saturated paste, respectively, according to the methods described by the NAAS [28]. The organic matter content of the soils was measured by Tyurin’s method.

2.3. Satellite Data Acquisition and Processing

RapidEye and Landsat OLI (Operational Land Imager) satellite images were used. A RapidEye image was acquired on 1 September 2014. The RapidEye imagery is multispectral, with a red edge (690–730 nm) and blue (440–510 nm), green (520–590 nm), red (630–685 nm), and near-infrared (NIR, 760–785 nm) bands. A level 2B RapidEye image was used, and the image was georeferenced to the Universal Transverse Mercator (UTM) coordinate system using the WGS 1984 datum with an accuracy of less than the root mean square of 0.01 pixels. The image was resampled to a 6 $\mathrm{m}\times 6\mathrm{m}$ spatial resolution. Image processing was conducted with ERDAS Image software (version 9.1).

Landsat 8 OLI images were acquired on 1 July and 5 October 2014. The Landsat 8 OLI blue (450–515 nm), green (525–600 nm), red (630–680 nm), and NIR (845–885 nm) bands were used for downscaling. The digital number (DN) values of the image were converted to reflectance with a multiplicative rescaling factor and an additive rescaling factor. Topographic correction was not necessary because the study area was very flat. The atmospheric correction was carried out with the Second Simulation of Satellite Signal in the Solar Spectrum (6S) model [29], using GRASS software [30].

2.4. Downscaling Cokriging

The following description of downscaling cokriging is mainly based on the work of Pardo-Iguzquiz et al. [24]. All of the formulas and notations followed Pardo-Iguzquiz et al. [24]. The cokriged finer-spatial-resolution image (downscaled Landsat) of band k, calculated from a Landsat (band k) image and RapidEye (band l), is given by:

$${\widehat{Z}}_u^k\left({\mathit{x}}_0\right)={\displaystyle \sum }_{i=1}^N{\lambda }_i^0{Z}_V^k\left({\mathit{x}}_i\right)+{\displaystyle \sum }_{j=1}^M{\beta }_j^0{Z}_u^l\left({\mathit{x}}_j\right)$$

(1)

where:

${\widehat{Z}}_u^k\left({\mathit{x}}_0\right)$ is a random variable (RV) of a pixel of areal size u (RapidEye), with the spatial location ${\mathit{x}}_0$ and spectral band k estimated by cokriging.

$Z_V^k\left({\mathit{x}}_i\right)$ is a RV of the pixel of the coarse spatial resolution image with areal size V (Landsat OLI) and spectral band k. N of these pixels are used. The weight assigned to the random variable of the i-th pixel is ${\lambda }_i^0$. The number of window pixels for Landsat (N) was nine (=3 $\times$ 3).

$Z_u^l\left({\mathit{x}}_j\right)$ is a RV of the pixel of the fine spatial resolution image with areal size u (RapidEye) and spectral band l. M of these pixels are used. The weight assigned to the random variable of the j-th pixel is ${\beta }_j^0$. The number of window pixels for RapidEye (M) was 25 (=5 $\times$ 5).

The two sets of weights {${\lambda }_i^0$; i = 1,…,N} and {${\beta }_j^0$; j = 1,…,M} are obtained through cokriging. To be optimal in the sense of giving a minimum variance unbiased estimator, the sum weights of the variable $Z_V^k$ must be 1 (${\sum }_{i=1}^N{\lambda }_i^0=1$), and the sum of the weights of the variable $Z_u^l$ must be 0 (${\sum }_{j=1}^M{\beta }_j^0=0$). The weights are calculated by formulizing the spatial structure of the image data and minimizing the unbiased estimation variance.

There are two ways to address the local mean in the downscaling cokriging process. One is to consider non-stationarity (taking into account the variation of local mean) [26], and the other is to use a global model [24]. This study used a global model due to its ease of application and statistical reliability [24]. The point spread function (PSF) is important to consider when factoring in area-to-point kriging. A uniform PSF was applied in this study because there was little information available to formulate a specific PSF. The Fortran source code for downscaling cokriging (DISKORI) was obtained via the Computers and Geosciences website [31].

2.5. Classification of Vegetation Type

The random forest classifier (RFC) proposed by Brieman [32] was applied to image classification because it is widely used in a range of fields and often yields good results [19,33,34,35]. The random forest algorithm (RF) is a type of decision tree classification. The RFC only has two parameters: the number of trees in the forest (number of bootstrap samples from the original data) and the number of random variables at each node (the maximum number of this parameter is the number of predictors). The RFC can use diverse types of data with little restriction because it does not make specific assumptions about the data distribution. Additionally, it can easily handle a large number of input variables.

R software [36] was used for image classification. The RFC was conducted with the ‘randomForest’ package (version 4.6–10), which was implemented by Liaw and Wiener [37] based on the algorithm proposed by Brieman [32].

Classification was performed for different combinations of images (

Table 1. Image bands used for classification.

Input Images		RF Model
		1	2	3	4	5	6
RapidEye 1 Sepetember	Blue	O	O	O	O	O	O
	Green	O	O	O	O	O	O
	Red	O	O	O	O	O	O
	Red edge	O	O	O	O	O
	NIR	O	O	O	O	O	O
	NDVI	O	O	O		O	O
Downscaled Landsat 1 July	Blue		O		O	O	O
	Green		O		O	O	O
	Red		O		O	O	O
	NIR		O		O	O	O
	NDVI		O			O	O
Downscaled Landsat 5 October	Blue			O	O	O	O
	Green			O	O	O	O
	Red			O	O	O	O
	NIR			O	O	O	O
	NDVI			O		O	O

). In addition to the reflectance of image bands, the normalized difference vegetation index (NDVI) was used. The NDVI was computed with the following equation:

$$\mathrm{NDVI}=\frac{{\rho }_{nir}-{\rho }_{red}}{{\rho }_{nir}+{\rho }_{red}}$$

(2)

where ${\rho }_{red}$ and ${\rho }_{nir}$ are the surface reflectance at the red and near infrared band, respectively.

Training samples were selected, based on the 127 randomly selected field survey points. Pixels with centers within a 9 m radius of the training samples were assigned the same vegetation type as the training samples (Figure 4). The radius criterion was determined by only selecting adjacent pixels (considering the 6 m resolution of the RapidEye and downscaled Landsat imagery) to the sampling points. Diagonal pixels were selectively included, according to the distance from the sampling point. The total sample size was 915. The vegetation classification was trained using randomly selected 200 pixels from each of the three vegetation types (600 training data = 200 $\times$ three types). After building RF models, the models were used to predict the vegetation type of the remaining 315 pixels. An accuracy assessment of classifications was carried out by comparing the predicted vegetation type and field data. The confusion matrices for each model were created after classification. The overall accuracy, producer’s accuracy, user’s accuracy, and Cohens’ kappa [38] were calculated, as suggested by Congalton [39].

2.6. Characterisics of Vegetation Distribution

After creating the vegetation map, the spatial characteristics of the vegetation were analyzed with R software [36]. The mean NDVI of each vegetation type was calculated with a ‘raster’ package (version 2.5–8). The patch statistics of the vegetation types were calculated with an ‘SDMTools’ package (version 1.1–221).

3. Results

3.1. Characteristics of Vegetation Types

HV was located in the high saline area and was dominated by Suaeda spp. and/or Salicornia europaea. The plant density was generally low and very low in some locations. Although there were some areas that were rarely covered by plants, a bare soil class was not included, because it is difficult to clearly distinguish between a bare soil class and a halophyte-dominated vegetation class, due to a low plant density and biomass. Additionally, it is worth noting that the area is not degraded, but is in the process of being colonized by halophyte species. This means that the areas of low plant density will be colonized by HV in the next few years, until the condition has been eased for non-halophyte species colonization. The height of HV was less than 0.5 m. The MLV was located in the mid- or low saline areas and was dominated by diverse grass species, with an admixture of Conyza Canadensis and/or Carex scabrifolia. The height of MLV was 0.5~1.5 m. The plant density of MLV was higher than HV, but generally lower than Reed vegetation. Reed vegetation was located in the low saline area and was dominated by Phragmites communis and/or Calamagrostis epigeios. The plant density of this type was usually higher than the other types. The height of Reed vegetation was 1.5~2.5 m.

3.2. Soil Environment of Vegetation Types

The texture of all soil samples was clayey. The soil EC, pH, and organic matter (OM) of the vegetation are summarized in

Table 2. Environmental conditions of sample points according to the vegetation types. Means with different superscript letters indicate significant difference between vegetation types (p < 0.05). (Mean: Arithmetic mean, SD: Standard deviation)

	HV		MLV		Reed
	Mean	SD	Mean	SD	Mean	SD
EC (dS/m)	9.7 a	5.5	1.3 b	2.1	0.5 b	0.4
pH	7.3 b	0.5	7.8 a	0.6	7.6 a	0.4
OM (%)	1.7 a	0.4	1.6 a	0.4	2.0 a	1.1

. The Wilcoxon rank sum test was used for the multiple comparison of mean, because the distributions were non-normal. The HV type showed the highest mean and standard deviation EC. The mean EC of the MLV and Reed vegetation types was significantly lower than HV. The high soil EC restricts the colonization of low- or non-saline plants. The pH was close to neutral in all vegetation types. However, the mean pH of HV was significantly lower than the other types. The mean OM was less than 2%. The low OM represents the poor soil conditions of the area.

3.3. Spectral Characteristics of Classification Classes

The spectral signatures of the RapidEye image are shown in Figure 5. The reflectance of the MLV and Reed classes was almost identical in the blue, green, and red bands. There was a large overlap between the HV and MLV classes in the NIR band. Although the red edge band was expected to be useful for the detection of fine differences in the vegetation condition, the three vegetation types showed a similar reflectance value in this band. The similarities in spectral reflectance among classes can reduce classification separability. However, the water class very clearly showed a distinguishable signature in the mean and standard deviation, especially in the red, red edge, and NIR bands.

3.4. Downscaling Cokriging

In this study, Landsat images with a 30 m resolution were downscaled to a 6 m resolution, to account for the spatial heterogeneity of the vegetation distribution. To conduct downscaling cokriging, the variogram and cross-variogram were fitted. An isotropic exponential model was applied as the fitting model for all variograms, because there was little difference in the direction of the experimental variogram. A nested model with two exponential structures was used, because one experimental structure was not sufficient for fitting the experimental variogram. The variogram and cross-variogram of the NIR band are shown in Figure 6. The fitted (cross) variogram models are summarized in

Table 3. Variogram models used in downscaling cokriging.

		RapidEye (1 Sepetember)			Landsat (1 July)			Cross-Variogram
		Nugget	First Structure	Second Structure	Nugget	First Structure	Second Structure	Nugget	First Structure	Second Structure
Common Range			40	300		40	300		40	300
Partial Sill	Blue	0	195	270	0	10	270	0	5	260
	Green	0	340	300	0	30	490	0	30	350
	Red	0	470	390	0	60	800	0	30	510
	NIR	0	1540	3420	0	180	810	0	75	1660
		RapidEye (1 Sepetember)			Landsat (5 October)			Cross-Variogram
		Nugget	First Structure	Second Structure	Nugget	First Structure	Second Structure	Nugget	First Structure	Second Structure
Common Range			40	300		40	300		40	300
Partial Sill	Blue	0	195	270	0	10	110	0	20	185
	Green	0	340	300	0	25	230	0	45	320
	Red	0	470	390	0	40	285	0	85	370
	NIR	0	1540	3420	0	155	1350	0	190	2325

. As the nested structure was applied, two ranges were identified. The range of the first structure was 40 m and that of the second was 300 m. In the case of the RapidEye image, the partial sills of the two structures were relatively similar. In contrast, the partial sills of the second structure were much larger than those of the first one in the Landsat image. This difference is related to the resolution of the images. For example, the small variation of the first structure can be explained by its short range. A pair of two points, separated by a smaller distance than the first range (40 m), was usually located in adjacent pixels or the same pixel. This leads to little difference in the values of the two points.

The plots of the reflectance value of the RapidEye image versus the downscaled Landsat images are shown in Figure 7. The reflectance of the July image was higher than that of the one produced in October, in all bands. However, it was lower than that of the RapidEye image in the blue and NIR bands (below of 1:1 line). The reflectance of the October image was more highly correlated with the RapidEye image in all bands.

Figure 8 and Figure 9 show the original Landsat images and downscaled images. Downscaled images provide spatially detailed information, especially in terms of the linear shape. The difference in the phenological traits of vegetation is easily found in these multi-temporal images. The color change of the images was variable, according to the location. Some areas remained unchanged or underwent little change, while other areas changed in color.

3.5. Result of Image Classification

The RFC was used to test whether the number of trees (the size of the forest) or the number of variables investigated at each node affected the classification accuracy. The former did not create a significant change in the results when the number of trees was larger than 1000. The accuracy was also not substantially affected by the number of variables investigated at each node when the number of trees was larger than 1000. The RFC was carried out with 1000 trees and two variables were investigated at each node.

A total of six RF models (Table 1) were tested for vegetation classification accuracy. The results are shown in

Table 4. Confusion matrix of the RF model using RapidEye (RF model 1).

$\widehat{\mathit{K}}$ = 0.72		Reference
		HV	MLV	Reed	User’s Accuracy
Prediction	HV	112	14	1	88.2%
	MLV	7	108	8	87.8%
	Reed	0	25	40	61.5%
	Producer’s accuracy	94.1%	73.5%	81.6%	82.5%

,

Table 5. Confusion matrix of the RF model using RapidEye and Landsat (1 July) (RF model 2).

$\widehat{\mathit{K}}$ = 0.81		Reference
		HV	MLV	Reed	User’s Accuracy
Prediction	HV	115	15	0	88.5%
	MLV	3	120	6	93.0%
	Reed	1	12	43	76.8%
	Producer’s accuracy	96.6%	81.6%	87.8%	88.3%

,

Table 6. Confusion matrix of the RF model using RapidEye and Landsat (5 October) (RF model 3).

$\widehat{\mathit{K}}$ = 0.81		Reference
		HV	MLV	Reed	User’s Accuracy
Prediction	HV	113	12	1	89.7%
	MLV	6	120	4	92.3%
	Reed	0	15	44	74.6%
	Producer’s accuracy	95.0%	81.6%	89.8%	87.9%

,

Table 7. Confusion matrix of the RF model using all of the bands, excluding NDVI (RF model 4).

$\widehat{\mathit{K}}$ = 0.87		Reference
		HV	MLV	Reed	User’s Accuracy
Prediction	HV	115	13	0	89.8%
	MLV	4	129	4	94.2%
	Reed	0	5	45	90.0%
	Producer’s accuracy	96.6%	87.8%	91.8%	91.7%

,

Table 8. Confusion matrix of the RF model using all of the bands (RF model 5).

$\widehat{\mathit{K}}$ = 0.87		Reference
		HV	MLV	Reed	User’s Accuracy
Prediction	HV	114	12	0	90.5%
	MLV	5	131	4	93.6%
	Reed	0	4	45	91.8%
	Producer’s accuracy	95.8%	89.1%	91.8%	92.1%

and

Table 9. Confusion matrix of the RF model using all of the bands, excluding Red Edge band (RF model 6).

$\widehat{\mathit{K}}$ = 0.88		Reference
		HV	MLV	Reed	User’s Accuracy
Prediction	HV	114	10	0	91.9%
	MLV	5	133	5	93.0%
	Reed	0	4	44	91.7%
	Producer’s accuracy	95.8%	90.5%	89.8%	92.4%

. The RF model (Model 1) that used a single RapidEye image yielded a poorer result (Overall Accuracy (OA) = 82.5%) than a multi-temporal image, which improved accuracy. Two images (RF model 2 (OA = 88.3%) and RF model 3 (OA = 87.9%)) were more efficient at producing accurate results. However, the best result was shown in the later three RF models (RF model 4–6) that used three image data, especially the RapidEye and Landsat images acquired in July and October, without the Rededge band (RF model 6 (OA = 92.4%) Table 9). There was no substantial difference in OA and $\widehat{K}$ among the three RF models.

The classification map of the RF model 6 is shown in Figure 10. The Reed patches were surrounded by MLV, and HV shared most of its boundaries with MLV. This pattern of vegetation distribution implies that there is a strong environmental gradient that has been spatially structured and regulates the vegetation distribution. Additionally, there were many linear patches located along narrow channels that were naturally formed when the area was an estuarine tidal flat.

The percentage of each vegetation type was 39.8% (HV), 38.9% (MLV), and 21.3% (Reed) (

Table 10. Zonal statistics of each vegetation type.

		HV	MLV	Reed
Area (%)		39.8	38.9	21.3
Mean NDVI	1 July (Downscaled Landsat 8)	0.188	0.281	0.463
	1 Sepetember (Rapid Eye)	0.262	0.441	0.527
	5 October (Downscaled Landsat 8)	0.216	0.389	0.491

). Almost 40% of the total area was covered by halophyte-dominated vegetation. Another 40% was colonized by diverse grass/herbaceous plants. The rest of the area was occupied by relatively dense vegetation, such as reed and small-reed plants. However, the Reed class is prone to expand due to its rhizome roots system, which is competitive in grass/herbaceous communities.

The pattern of NDVI reflects the plant phenology. Early July marks the beginning of the growing season. The leaf biomass increases until late summer or early autumn. After the peak, the NDVI decreases, due to the loss of leaf biomass. The pattern was the same for all vegetation types (Table 10). The Reed type had the highest NDVI, meaning that the Reed type may substantially contribute to the above ground biomass of the area, regardless of its relatively small coverage.

The area of the vegetation patch was calculated and the frequency is shown in Figure 11. There were many small patches. The number of patches smaller than 450 m² (half of Landsat support area: 30 m $\times$ 30 m) was 569. These small patches are difficult to identify using the Landsat image for the vegetation classification. Each pixel of Landsat was composed of 25 subpixels that contained information on the vegetation types. Among the Landsat pixels, only 47.7% were composed of a single vegetation type, 44.6% were composed of two vegetation types, and 7.7% were composed of three vegetation types. The highly diverse vegetation distribution needs finer image data to figure out the state of vegetation distribution.

4. Discussion

The study area became permanent land after the water level was lowered in 2010, and the soil condition has been changing. Kim et al. [3] reported the mean soil EC of halophyte vegetation as 14 dSm-1, mixed vegetation as 6.7 dSm-1, and reed dominant vegetation as 3.0 dSm-1 in 2010. Compared to an earlier study, the soil EC showed lower levels in 2014. The mean soil EC of HV was 9.5 dSm-1, MLV was 1.0 dSm-1, and Reed was 0.5 dSm-1. The rapid changes in the soil conditions during this relatively short period also influenced the distribution of vegetation.

Although the characteristics and distribution of vegetation are crucial for determining the condition of the ecosystem, it is hard to create accurate vegetation maps of grass/herbaceous plant-dominated areas. There have been a number of studies that have tried to improve the accuracy of vegetation mapping using the phenological characteristics of plants [11,20,21,22]. Gilmore [11] suggested that it is possible to discriminate Phragmites spp. using single-date images acquired in the autumn season. However, it was difficult to classify Phragmites-dominated vegetation using timely acquired RapidEye imagery. More than 20% of the Reed class was misclassified (Table 4). The accuracy increased when the multi-temporal images were used in classification (Table 5, Table 6, Table 7 and Table 8).

The goal is to acquire multi-temporal images, especially in the growing season. This research applied cokriging to downscale Landsat imagery, which is easy to obtain, but has limited utility due to its relatively low spatial resolution. Sharpened Landsat imagery taken at different time periods improved the accuracy of the vegetation map. Cokriging is a useful approach because it is based on a well-established geostatistical theory and provides a framework to combine different types of data [23,24]. However, the application of the approach has been limited, until recently [40]. The results show that image downscaling using cokriging can be applied in vegetation studies and is useful for the analysis of multi-temporal data.

In this study, a nested variogram model with two structures was used in cokriging. A nested model is useful for modelling the spatial structure with different scales [41,42]. Pardo-Iguzquiz et al. [24] also used a nested model for downscaling two images with different resolutions. The nested structure model fitted the data well, with a 0 nugget value. The partial sill of the first structure (short range) was much smaller than that of the second structure in a variogram model of a low resolution image. It is consistent with the results of this study. However, in both cases, further research is needed to quantify the uncertainty of using the nested model for downscaling images.

Multi-temporal and multispectral data analysis raises the problem that the number of input variables can rapidly increase. The RFC can easily and effectively handle a large number of input variables [19,33,34,43]. For that reason, it was used in this study and yielded good results.

5. Conclusions

An area reclaimed from estuarine tidal flats forms a very unique and dynamic ecosystem. The creation of a large reclamation area in itself is a difficult challenge for ecosystem monitoring and management.

This paper presents a novel approach for the classification of vegetation in reclaimed areas, using multi-temporal downscaled Landsat imagery along with RapidEye imagery. The study area has a typical environmental gradient of soil salinity and vegetation types, of which the distribution is mainly controlled by the environmental conditions. This study provides a framework to reduce spatial and temporal resolution by downscaling images with a low spatial and high temporal resolution, with an image that has a high spatial and low temporal resolution. Multi-temporal downscaled images were compared to single high resolution images, to classify the vegetation type. The accuracy of the vegetation map using the multi-temporal downscaled image was approximately 10% higher than the accuracy of the map using a single high resolution image.

These results demonstrate that the downscaling technique is useful when a high resolution image is not available. This technique is also capable of increasing the utility of relatively low spatial resolution images, such as Landsat imagery. The results also confirmed the performance of the RFC. The analysis of multi-temporal multi-band data causes difficulties in handling high dimensional data. The RFC was very versatile and demonstrated a high performance when classifying the vegetation type with a large number of variables.

The results of this study can be utilized in ecosystem monitoring for target objects that are difficult to distinguish with single, high spatial resolution imagery. The method can not only expand the time span of the data, but also increase the accuracy of mapping. It is also possible to apply this approach to the monitoring of ecosystems that witness rapid changes.