Geo-Parcel Based Crop Identification by Integrating High Spatial-Temporal Resolution Imagery from Multi-Source Satellite Data

Abstract

Geo-parcel based crop identification plays an important role in precision agriculture. It meets the needs of refined farmland management. This study presents an improved identification procedure for geo-parcel based crop identification by combining fine-resolution images and multi-source medium-resolution images. GF-2 images with fine spatial resolution of 0.8 m provided agricultural farming plot boundaries, and GF-1 (16 m) and Landsat 8 OLI data were used to transform the geo-parcel based enhanced vegetation index (EVI) time-series. In this study, we propose a piecewise EVI time-series smoothing method to fit irregular time profiles, especially for crop rotation situations. Global EVI time-series were divided into several temporal segments, from which phenological metrics could be derived. This method was applied to Lixian, where crop rotation was the common practice of growing different types of crops, in the same plot, in sequenced seasons. After collection of phenological features and multi-temporal spectral information, Random Forest (RF) was performed to classify crop types, and the overall accuracy was 93.27%. Moreover, an analysis of feature significance showed that phenological features were of greater importance for distinguishing agricultural land cover compared to temporal spectral information. The identification results indicated that the integration of high spatial-temporal resolution imagery is promising for geo-parcel based crop identification and that the newly proposed smoothing method is effective.

1. Introduction

With the development of agricultural management, agro-ecology studies and agricultural policy making, especially in precision agriculture, there have been increasing demands for crop distribution information on a land plot scale. Field survey methods to acquire crop types are time-consuming and laborious. In contrast, remote sensing has been widely used in agricultural management because it can be effectively applied on a large scale. However, due to being limited by data acquisition and satellite reception issues, remote sensing imagery has difficulty supplying high spatial and high temporal resolution at the same time.

A number of studies have identified different crops, based on their phenology, using high temporal resolution satellite data, particularly from coarse-resolution images, such as those from Advanced Very High Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS). With a short revisit interval, they help to continuously observe land cover changes and detect vegetation phenological development. Many studies focused on constructing a vegetation index (VI) time-series using these kinds of images. Reed et al. [1] derived a suite of 12 metrics from 4 years of AVHRR normalized difference vegetation index (NDVI) time-series, and they found that the metrics were strongly consistent with various land cover types and their predicted phenological characteristics. Hill et al. [2] disaggregated quantitative metrics from AVHRR NDVI time-series for 8 years, analyzed the temporal and spatial patterns of key NDVI metrics, and made classifications. Sibanda et al. [3] derived phenological metrics from a 16-day MODIS NDVI time-series, compared the greenness change rate of cotton, maize and sorghum over the onset of green-up, green-peak and the senescence growth period, and significantly (p < 0.05) distinguished cotton fields from maize and sorghum fields. As shown above, the vegetative index (VI) time-series can apparently respond to vegetative phenology. However, due to the sensor observation angle, response accuracy, sunlight conditions, cloudy weather and atmospheric factors, the VI time-series showed an irregular status, which made it difficult to extract phenological metrics. Many studies have proposed various VI smoothing methods, and the commonly used VI curve filtering methods are the Best Index Slope Extraction Method (BISE) [4], Mean Value Iteration (MVI) [5], Savitzky–Golay Filter (S–G) [6], Asymmetric Gaussian Fitting (AG) [7], Whittaker smoother (WT) [8] and Double Logistic function fitting algorithm (DL) [9]. Hird et al. [10] selected six kinds of NDVI filtering methods for a quantitative comparative analysis and explained that the purpose of filtering was to maximize the elimination of abnormal impacts while maintaining true values. These existing methods have been applied to smooth VI time-series and eliminate unnecessary fluctuations, especially for crop single growth seasons. However, crop rotation practices demand a smoothing method for fitting a time profile with multiple peaks. Our study presents an adaptive smoothing method for identifying crop rotations.

Although it is advantageous to form continuous observations of land cover using coarse-resolution images, because of low spatial resolution and resulting mixed pixels, pixel-based analyses often lead to inaccurate results. Consequently, coarse-spatial resolution is unsuitable for identifying crops in small agricultural field plots that have high spatial heterogeneity. With the development of high-resolution technology at the beginning of 21st century, high spatial resolution images, such as images from SPOT-5, IKONOS, and Worldview satellite, are capable of depicting detailed geometric and textural information, which is often negligible in low- and medium-resolution images. Turker et al. [11] identified geo-parcel based crop types separately, using SPOT-4, SPOT-5, IKONOS and Quickbird images for similar climatic conditions, and they compared the classification accuracies between filtered and unfiltered images. As a whole, the resulting accuracies were above 83%. Immitzer et al. [12] used the Random Forest algorithm to classify 10 tree species in Austria, using high spatial resolution imagery, Worldview-2 with eight bands, and an overall accuracy of approximately 82%. To effectively utilize the geometric/textual information of high-resolution imagery, the object-based analysis method was developed. The performances of pixel-based versus object-based approaches were compared by many studies, and scholars found that the latter approach outperformed the former approach [13]. However, the revisit cycle of high spatial resolution imagery is too long to capture the optimal phase for crop identification, not mentioning the frequent interference of cloud cover during the acquisition time.

The contradiction between the spatial and temporal resolutions of single-source satellite data may limit its potential for many applications that demand observations with high spatial and temporal resolution. To overcome the limits of single-source data, a trend of synchronously utilizing data from multiple satellite sensors has recently been observed. One of these types of techniques is Spatial and Temporal Data Fusion (STDF), which aims to construct relatively intensive time-series imagery with fine spatial resolution by blending coarse and fine spatial resolution images [14,15,16]. Hilker et al. [15] used the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM) to blend Landsat images with MODIS data for crop monitoring in central British Columbia, Canada. In their study, 18 synthetic time-series Landsat images were obtained using the STARFM algorithm, and the results showed that the crop growing patterns can be captured well by the synthetic time-series data. However, resolutions of these Landsat-like synthetic images are still relatively coarse (mostly 30 m resolution), and they are mainly used for pixel-based applications. Recently, some object-based STDF methods were proposed, which integrate objects segmented from high-resolution imagery with spectral/temporal information from coarse imagery. In a study of crop rotation classification, Conrad et al. [17] extracted the land plot objects through segmenting the SPOT-5 data and then assigned tasseled cap indices, derived from the bi-temporal ASTER images, to these objects. Their results demonstrated the effectiveness of this method. Singha et al. [18] obtained synthetic time-series dataset with a 16-day revisiting cycle and 30 m spatial resolution from MODIS and HJ-1 A/B CCD by using the STARFM algorithm and used a multi-resolution segmentation algorithm on an HJ-1 dataset to acquire objects. After the selection of key temporal features, a decision tree classifier produced a satisfactory performance (overall accuracy of >84%) in paddy rice mapping. Although the object-based STDF methods could achieve relatively high accuracy in crop identification, objects that were segmented from imagery, in the abovementioned resolution, may be inconsistent with the real boundaries of farm land plots, and the only single sensor for medium-resolution data may lead to insufficient temporal coverage. Therefore, higher-resolution spatial and temporal data from multi-source data should be integrated, to improve the geometric accuracy of land plot objects and expand the temporal observation cover range.

This paper aimed to solve the limitations of single-source satellite data in geo-parcel based crop identification. Fine spatial resolution images with 0.8 m spatial resolution from the GF-2 satellite were used to obtain precise agricultural field plots, which were considered to be the smallest analytical elements. Multi-source medium-resolution data from GF-1 and Landsat 8 OLI were utilized together, in order to acquire intensive observations. To improve the EVI time profile fitting effect, we developed a piecewise EVI time-series smoothing method that worked well, especially for rotated crop cases. Then, we constructed a feature set that contained phenological metrics that were derived from the smoothed EVI time profile, and multi-temporal spectral features, and made classifications using Random Forest. The results showed that the proposed piecewise fitting method performed better than the S–G and the Harmonic Analysis of Time Series (HANTS) methods for crop classification, and the integration of the two kinds of features could significantly improve classification accuracy.

2. Materials and Methods

2.1. Study Area

Lixian, with an area of 2107.3 km², is located in the northwest part of the Hunan Province (Figure 1). It is named after Lishui River, which flows through the entire territory. Its middle part is situated in the largest plain in Hunan, the famous Liyang Plain, and the northwest part of the province is mountainous. It has a subtropical humid monsoon climate, with an average temperature and rich water resources. The average annual temperature is 16.7 °C, and the annual precipitation is 1200 to 1900 mm. The climatic conditions are suitable for crop growth, including rice, cotton and cole crops.

2.2. Data

2.2.1. Remotely Sensed Imagery

Two disparate spatial resolution images were used in our study. One of them was fine spatial resolution imagery, obtained from GF-2 satellite data. Due to its narrow width of 45 km, 12 images in total covered the study area completely. The Gram–Schmidt Pan Sharpening method was used to sharpen multispectral (4 m) bands using high-resolution panchromatic (1 m) bands. This created a high-resolution multispectral GF-2 dataset, and the resolution of sub-satellite points reached 0.8 m. Then, Landsat 8/OLI imagery, provided by the United States Geological Survey (USGS), was used as reference data during orthorectification and geometric correction. The second-order polynomial interpolation method was used, and approximately 50 ground control points were selected for precise geometric registration, resulting in a controlled correction error within 0.5 pixels.

The medium-resolution imagery was mainly from the GF-1 WFV (Wide Field View) sensors. Four WFV sensors had the ability of acquiring imagery with 16 m spatial resolution and a four-day revisit cycle. The spectral bands ranged from visible to near-infrared wavelengths. Technical specification of the GF1 dataset is provided in

Table 1. Technical specifications of the GF-1 WFV sensors.

Satellite	Payloads	Bands No.	Spectral Range (μm)	Spatial Resolution (m)	Swath Width (km)	Repetition Cycle (Day)
GF-1	WFV	1	0.45–0.52	16	800 (four cameras combined)	4
		2	0.52–0.59
		3	0.63–0.69
		4	0.77–0.89

. Due to cloudy and rainy weather, valid GF-1 images were always insufficient in continuous observation, as shown in

Table 2. Information for the medium-resolution images that were used in the study (WFV1 for NO.1 WFV sensor in GF-1, and so on).

Satellite	Sensor	Acquisition Time	Day of Year
GF-1	WFV2	8-02-2016	38
GF-1	WFV1	28-03-2016	88
GF-1	WFV3	18-04-2016	109
GF-1	WFV4	05-05-2016	126
GF-1	WFV1	12-05-2016	133
GF-1	WFV1	16-05-2016	137
GF-1	WFV2	14-06-2016	166
GF-1	WFV1	14-06-2016	166
GF-1	WFV4	23-06-2016	175
GF-1	WFV2	09-07-2016	191
GF-1	WFV3	09-07-2016	191
GF-1	WFV1	25-07-2016	207
GF-1	WFV1	29-07-2016	211
GF-1	WFV3	15-08-2016	228
GF-1	WFV3	15-08-2016	228
GF-1	WFV3	03-10-2016	277
GF-1	WFV1	04-11-2016	308
GF-1	WFV4	26-11-2016	330
GF-1	WFV4	04-12-2016	339
GF-1	WFV1	15-12-2016	350
Landsat 8	OLI	30-07-2016	212
Landsat 8	OLI	16-09-2016	260

. To construct intensive time-series, Landsat 8 OLI multispectral bands were used as supplementary data when a critical phase for crop identification was missing or key time point should be caught. Table 1 lists the medium-resolution data that were used in the experiment.

2.2.2. Ancillary Data

Digital elevation model (DEM) data were used for the orthorectification of GF-1 data. To assist with crop classification and accuracy assessment of the classification results, an extensive field survey was conducted in March, August and November 2016. Samples from 500 farming plots were collected using GPS to record their geometric positions. Other ancillary data, including administrative boundaries, road networks, and phenological information about crops, were also collected.

2.3. Methods

To identify farmland crop types, much pre-processing work on multi-source satellite data was done. We proposed a piecewise EVI time profile smoothing method for fitting the EVI time series, and we made classifications using a Random Forest classifier.

The workflow of the farmland-based crop mapping method is shown in Figure 2.

2.3.1. Geo-Parcels Extraction

The first step of the object-based STDF analysis involves object extraction to produce separate and homogenous geo-parcels. Geo-parcel extraction was based on GF-2 0.8 m fusion data, which accurately depicted the surface information in detail (Figure 3). It was performed manually along the real boundaries of farm land plots, as shown in the image. Each plot was considered to be the smallest analytical element. The total number of geo-parcels was approximately 590,000.

2.3.2. Pre-Processing of Medium-Resolution Images

GF-1 WFV data were radiometrically calibrated, atmospherically corrected, orthorectified and accurately geographically registered. The geographic correction used Landsat 8 OLI imagery as a reference, making the correction error controlled within 0.5 pixels. Landsat 8 OLI images were also radiometrically calibrated and atmospherically corrected. We used the FLAASH module in ENVI software for the atmospheric correction of the Landsat 8 OLI images. After that, image sharpening was applied to OLI images, using panchromatic and multispectral bands, which produced images with 15 m resolution. Due to the sensor performance and satellite imaging pose, it is essential to normalize different sensors in radiometric features. Radiometric inter-calibration was performed, using the method presented by Huang et al. [19]. The surface reflectance of other GF-1 WFV sensors and OLI was calibrated, using the reflectance of the GF-1 WFV1 sensor. Satellite images had many invalid pixels caused by cloudy and hazy weather conditions. Discarding the invalid pixels might lead to an inadequate observation frequency when constructing VI time-series. The strategy adopted here was to remove cloud and haze parts, fragment images and preserve the valid parts.

2.3.3. Construction of the Vegetation Index Time-Series

The VI time-series fluctuates according to physiological states during growing stages. Land cover types can be identified by analyzing the VI time-series and excavating its changing pattern [20,21,22]. Spectral information from the multi-temporal images was assigned to geo-parcel objects. Pixels that intersected the geo-parcel boundary were joined with calculations of temporal-spectral values. The average of several VI values that were derived from the adopted pixels represented the VI value of the field plot. The spectral value of each band was also calculated using this method.

In addition to the commonly used NDVI, EVI can also characterize the coverage of vegetation. Hongjun et al. [23] used mathematical statistics and geoscience methods to compare MODIS NDVI and EVI, and found that EVI was more sensitive during the vegetation growth period, whereas NDVI more easily reached saturation. Moreover, the standard deviation and value range of EVI was greater than NDVI at the same spatial resolution. We compared EVI to NDVI during crop growth seasons. Figure 4 depicts the NDVI and EVI time profile of typical land cover types through 2016. The comparison results showed that EVI provided a higher capability for extracting phenological metrics compared to NDVI. In the current study, EVI was chosen to characterize vegetation growth status.

2.3.4. EVI Time Series Smoothing Method

The planting patterns in the study area were very complex. It is common to rotate different crops multiple times on the same farm plot. The EVI time-series of rotated plot fields showed multiple peaks. The complex shapes of time profiles made it more difficult to smooth curves and derive phenological metrics. According to Jönsson et al. [24], the global function can be reconstructed by integrating several local functions [25]. In this study, we developed a piecewise fitting method, in which the global time profile was disaggregated into local profiles. We then performed fitting methods on each part, and finally connected them, so that they were continuous.

The implementation method was as follows:

(1)

First, the linear interpolation method was applied, to interpolate the missing points along the time profile. The time interval was set as five days, which produces a dense EVI curve.

(2)

Then, the S–G filtering method was used to eliminate small unnecessary fluctuations that are caused by system factors. The S–G filter method was applied twice, to enhance the smoothing effects. This procedure aims to highlight key points, namely summit points and bottom points. Here, the S–G parameter window size was set as five, and polynomial degree was set as three.

(3)

Based on the slight S–G smoothing process, it became easier to distinguish the summit and bottom points. Several groups of key points were selected. Each group was used to represent a vegetative greenness period. Every group consisted of a summit point (t_C) and two bottom points (t_L and t_R) at t_C side, the difference of which was required to be greater than 0.2. A difference that was less than 0.2 was considered to be an abnormal fluctuation, so the corresponding group was discarded. After that, the global time profile was divided into several temporal parts, and each part was from t_L to t_R.

(4)

Furthermore, based on temporal segments, the Gaussian and Polynomial fitting methods were used to smooth the local EVI time profile (Equation (1)), from t_L to t_R. Figure 5a shows the original data and its fitting result (f_center) between two bottom points (t_L and t_R). The Gaussian and Polynomial fitting methods were also conducted between every two neighboring summit points. The original data and their fitted functions are shown in Figure 5b,c. To connect fitted local functions, the function fitting method (Equation (2)) was developed, as shown in Figure 5d.

$$\mathrm{f}\left(t\right)=a{e}^{-\frac{z^2}{2}}+d{t}^2+e+f,\mathrm{where}z=\frac{t-b}{c}$$

(1)

$$\{\begin{array}{c}{f}_{\mathrm{left}}{t}_L<t<\frac{t_L+t_C}{2}\\ f_{\mathrm{center}}\frac{t_L+t_C}{2}<t<\frac{t_R+t_c}{2}\\ f_{\mathrm{right}}\frac{t_R+t_c}{2}<t<t_R\end{array}$$

(2)

2.3.5. Calculation of EVI Time-Series Metrics

Before the calculation of phenological metrics, the smoothed-and-connected global EVI time profile was disintegrated into temporal parts to recognize crop rotation. Crop rotation reflects the utilization of land use [26]. In our present study, the temporal segments with a fluctuation amplitude greater than 0.2 were considered to be the growth seasons. Then, 11 phenological metrics were derived from the segmented EVI time profile, including the time of onset of greenness (OnT), value at OnT (OnV), time of maximum EVI (maxT), value at maxT (maxV), time of end of greenness (EndT), value at EndT (EndV), growth rate (GR), senescence rate (SR), duration of greenness (DT), integrated area under segmented EVI time profile (IntegratedEVI), and the amplitude of EVI through greenness (Growth Amplitude, or GA). Many of the metrics were adapted from previous studies ([1,27]). Figure 6 shows the derivation of these metrics. The principle task of metric extraction is the identification of OnT and EndT. Their calculation used a threshold strategy, in which the threshold was assigned as an EVI value greater than the bottom value by 20% of the amplitude. Such a threshold value signified the soil background level. The times of the beginning and end of vegetative activity were then acquired, according to the mapping relationship between the threshold and the time axis. GR was estimated from the ratio of amplitude and duration time from beginning to maximum, and similarly, SR was defined. The integrated EVI indicated overall conditions through the growth season, and some studies demonstrated that cumulated indices had a good correlation with potential biomass [28].

2.3.6. Crop Identification Using Random Forest

According to empirical knowledge, typical staple crops in the study area could be planted no more than three times in four seasons of the year. We assumed that crop planting activities occurred three times per field plot. After the global EVI time profile was divided into segmented parts, the 11 phenological metrics of each part were derived. The global time profile that represented the phenological changes throughout the year contained 33 metrics at most. The phenological index of the field, which was rotated less than three times, was set as −1.

In addition to phenological characteristics, spectral information was also of great importance for crop identification. The monthly average reflectance of visible and near infrared bands was utilized to expand the feature set. Moreover, the monthly average EVI and monthly average Normalized Difference Water Index (NDWI) [29] were also combined into the feature set.

With large numbers of accessible features, RF was used to identify crops in the study. RF uses bagging integration, based on a Decision Tree (DT), and introduces a random feature selection process [30]. RF draws bootstrap samples from the original dataset to construct several individual DTs, and the bootstrap process makes it relatively difficult to over-fit. DT randomly selects a subset of variables and chooses the best-splitting variable as the splitting node. This step repeats until the individual DT ends. RF integrates a great many of the abovementioned DTs and finally generates classification results from the majority of the votes. In this study, the feature set was expanded to 105 dimensions, and the RF classifier method was considered to be an ideal choice, because it is easy to implement and has a low computational cost.

3. Results

This analysis was performed on agricultural fields in Lixian, and more than 590,000 field plots were classified. Multi-source medium-resolution images were organized in an orderly way and were transformed into EVI time-series for each field plot. The phenological metrics were derived, and the crop types were recognized. There were three dominant crops: paddy rice (Oryza sativa), cotton (Gossypium hirsutum) and cole (rapeseed). Figure 7 shows the EVI time profiles for those crops in the study area, during 2016. However, the crop planting patterns were varied, due to farming habits or production demands. This part of the paper shows how this method worked in Lixian.

3.1. Time Profile Smoothing

The common agricultural farming patterns are displayed in Figure 8. The following figures depict smoothed EVI time profiles and original data. The smoothed global functions are tightly aligned with the original EVI time-series and preserve vegetative activities that should be retained. The piecewise fitting method showed satisfactory performance for smoothing these kinds of EVI time profiles of complex planting patterns.

3.2. Rotation Frequency Judgement

By means of the rotation judgement technique, the rotation frequency of all field plots could be recognized. Figure 9 shows the spatial distribution and indicates that most farming fields were planted only once throughout the year, the planting structure of which was simple. They were spatially distributed across the entire area. Approximately 16.39% of the 59,000 plots were rotated, most of which were distributed in the east part of Lixian. This phenomenon can be explained by the natural terrain and climatic conditions. Eastern Lixian is located in the Dongting Lake Plain, which has an average elevation of less than 50 m. Rivers run through the plain, making it convenient for agricultural irrigation. However, the north part of Lixian is mountainous. Agricultural fields are scattered, and the majority of them are planted without rotation.

3.3. Phenological Features Extraction

We studied the phenological characteristics of typical crops. According to regular planting patterns, paddy rice that is planted once in a year is called single season rice, and rice that is planted twice is called double season rice. The greenness of single season rice usually lasts from May to October. The cultivation of double season rice is stricter. Transplanted in mid-April, the first season rice (early rice) is harvested in mid-July, and the second season rice (late rice) is planted soon after, to be harvested in late October. The cotton greenness lasts from May to October, and it has five growth stages: sowing, seedling, budding, flower blooming and boll opening periods. In the autumn sowing season, cole is transplanted in late October, then it hibernates through the cold winter weather, and turns green in February of the next year to be harvested in mid-May.

In our study, phenological metrics were derived after the smoothing procedure for the EVI time profiles. Figure 10 shows the statistics of the 11 phenological metrics in the form of bar charts with error bars. In the following charts, the averages and standard deviations describe the given variables. As shown by these charts, phenological characteristics that were related to time (e.g., OnT, EndT, MaxT, DT and integrated EVI) had higher distinguishability than others, because they had more distinct ranges with shorter error bars (especially OnT, EndT and MaxT). Phenological metrics, such as OnV, EndV, MaxV, GR and SR, did not make an obvious distinction for specific crops, as they had similar values between the different crops, and the ranges were broad.

3.4. Crop Classification and Spatial Distribution

By using the RF classifier, more than 590,000 agricultural plots were classified into seven farming types, which were single season paddy rice, paddy rice-cole, cole-paddy rice, double season rice, paddy rice-cotton, cole-paddy rice-cole, and other crops. The ‘-’ indicates crop rotation in sequenced seasons.

The integration of EVI time profile features and multi-temporal spectral features led to a satisfactory performance, the overall accuracy of which was 93.27% (

Table 3. Confusion matrix of classification using phenological features and spectral features.

All Features	Paddy Rice	Paddy Rice-Cole	Cole-Paddy Rice	Double Season Paddy Rice	Cole-Cotton	Cole-Paddy Rice-Cole	Other Crops	User’s Accuracy
Paddy Rice	79	0	0	0	0	0	29	0.7315
Paddy Rice-Cole	0	170	0	0	0	1	0	0.9942
Cole-Paddy Rice	0	0	218	0	21	1	0	0.9083
Double Season Paddy Rice	0	0	1	115	2	0	0	0.9746
Cole-Cotton	0	0	12	0	121	0	0	0.9098
Cole-Paddy Rice-Cole	0	0	0	0	0	245	3	0.9880
Other Crops	3	1	7	0	0	2	202	0.9395
							Overall Accuracy	0.9327
							Kappa	0.9201

). Additionally, it indicated that double season rice, cole-paddy rice-cole, and paddy rice-cole could be easily recognized, with accuracies higher than 95%. Single season paddy rice was most likely to be misidentified, due to its simple greenness without regular rotation. The results implied that plots that were rotated multiple times had higher accuracies because agricultural activities usually abided by natural law, and there were more limits on successive planting activities.

The crop spatial distribution was acquired, as shown in Figure 11. The distribution indicated that, in the plains region, farming plots were cultivated adequately for agricultural production. According to classification results, rice occupied the majority of field plots, accounting for approximately 44.32%. There were approximately 46,000 plots that were rotated with single season rice and cole, 11,919 plots were planted with double season rice in 2016, accounting for approximately 2%, and most of them were planted in central Lixian in a scattered distribution. As an autumn-sown crop, cole was likely to be rotated with rice or cotton.

From another perspective, the following figures (Figure 12 and Figure 13) show thematic mappings for specific crop types. Figure 12 displays a map of paddy rice and cotton. Paddy rice includes double season rice, single season rice without rotation, and rice rotated with cole. Cotton was mostly distributed in eastern Lixian, and some was planted in the south. Figure 13 presents a map of autumn-sown cole. Cole was mostly planted in eastern Lixian, with some distributed in the south, and some spread across the study area, even in mountainous terrain.

4. Discussion

4.1. Comparison of the S–G and HANTS Smoothing Methods

Comparing the S–G and HANTS smoothing method, the results show that the piecewise fitting analysis outperformed the others. Figure 14a shows the EVI time profile of farmland planted with single season rice. The fitting results of the three methods revealed that the HANTS method and our proposed piecewise fitting method were able to fit the temporal curve to a smooth one, whereas S–G conformed tightly to its original shape. To summarize, these methods all work well for fitting an EVI time profile of farmland that is not rotated. Figure 14b shows the EVI time profile with three summits of a rotated farmland. In this figure, the first summit indicates vegetative activity from the beginning of the year, and the third summit shows newly planted crops in autumn. The middle sub-period with the highest peak value indicates a complete growth season. It shows that the effect of the S–G filter was not satisfactory, because it was unable to eliminate unnecessary fluctuations, which causes difficulties in rotation judgement. Another method, the HANTS filter, which fits the time profile using various harmonic functions with different frequencies, can smooth curves, while preserving multiple peaks. However, in regard to time profiles with sub-periods that represent incomplete growth seasons, such as in Figure 14b, it performs poorly. As indicated, it maintains the basic shape of the curve but neglects important partial periods, and the fitted data seriously deviates from the original data. In contrast, our piecewise fitting method can effectively eliminate abnormal fluctuations, which helps to judge crop rotation, and preserves complete or incomplete greenness periods as well.

4.2. Classification Accuracy Comparison

Apart from visible distinctions between S–G, HANTS and the piecewise fitting methods in fitting the EVI time profile, the classification accuracy is also improved (Figure 15). We applied three groups of classification experiments using phenological features that were independently derived from the S–G, HANTS and piecewise fitting methods (PF). The overall accuracies were 82.78%, 77.42% and 86.72%, respectively. Phenological metrics of the HANTS fitting method led to the lowest accuracy, whereas the piecewise fitting method had the highest accuracy, but was still below 90%. We performed another three groups of experiments that utilized multi-temporal spectral features (SF) and phenological features, and the accuracies increased by different degrees. The accuracy of S–G + SF (phenological metrics derived from S–G and spectral features) achieved 85.62%, HANTS + SF (phenological metrics derived from HANTS and spectral features) achieved 82.44%, and PF + SF (phenological metrics derived from the piecewise fitting method and spectral features) achieved 93.27%. The results showed that the S–G smoothing method outperformed the HANTS in crop classification, and the piecewise fitting method worked best. Moreover, the integration of SF and phenological metrics significantly improved the classification accuracies.

To analyze the effects of feature set integration, we selected the well-performing piecewise fitting method to compare classifications while using and not using SF. In contrast, the performance of classification using 72 multi-temporal spectral features was unsatisfactory, with an accuracy of 81.21% (

Table 4. Confusion matrix of classification using spectral features.

Spectral Features	Paddy Rice	Paddy Rice-Cole	Cole-Paddy Rice	Double Season Paddy Rice	Cole-Cotton	Cole-Paddy Rice-Cole	Other Crops	User’s Accuracy
Paddy Rice	80	2	5	0	0	1	20	0.7407
Paddy Rice-Cole	60	78	2	0	0	30	1	0.4561
Cole-Paddy Rice	10	0	201	0	21	8	0	0.8375
Double Season Paddy Rice	1	0	0	114	0	0	3	0.9661
Cole-Cotton	1	0	13	0	118	0	1	0.8872
Cole-Paddy Rice-Cole	1	1	19	0	4	224	1	0.8960
Other Crops	13	1	8	0	5	0	188	0.8744
							Overall Accuracy	0.8121
							Kappa	0.7777

). Without support of phenological features, rotated crops were easily misclassified. For instance, paddy rice-cole was commonly identified as cole-paddy rice-cole, and its accuracy was only 45.61%. In comparison, the accuracy of crop identification using phenological features was better, with an accuracy of 88.73% and kappa coefficient of 0.778 (

Table 5. Confusion matrix of classification using phenological features.

Phenological Features	Paddy Rice	Paddy Rice-Cole	Cole-Paddy Rice	Double Season Paddy Rice	Cole-Cotton	Cole-Paddy Rice-Cole	Other Crops	User’s Accuracy
Paddy Rice	102	0	0	0	0	0	6	0.9444
Paddy Rice-Cole	0	170	0	0	0	1	0	0.9942
Cole-Paddy Rice	0	0	168	0	71	1	0	0.7000
Double Season paddy Rice	0	0	0	113	3	2	0	0.9576
Cole-Cotton	0	0	17	0	116	0	0	0.8722
Cole-Paddy Rice-Cole	0	0	0	0	5	243	0	0.9798
Other Crops	23	0	5	3	0	2	182	0.8465
							Overall Accuracy	0.8873
							Kappa	0.8672

). However, without the support of spectral information, the accuracy of cole-paddy rice recognition declined from 90.83% (in Table 3) to 70%, and it was most commonly misclassified as cole-cotton, due to close cultivation times between the rice and cotton.

To summarize, compared with the above results, our study revealed that the integration of methods could significantly improve identification performance (Table 3). Phenological metrics extracted from the EVI time profile reflect vegetative changes during greenness, and multi-temporal spectral information plays an important role when distinguishing specific crops.

4.3. Features Significance Ranking

The RF classifier generated a significance ranking of the 105 features that were integrated from phenological features and multi-temporal spectral features (Figure 16). In Figure 16, features derived from the EVI time profile were named as ‘OnT_1’, where ‘1’ represented the first greenness. Features extracted from the multi-temporal medium-resolution images were named as ‘B1_3’, where ‘3’ represented blue band reflectance (B1 for blue, B2 for green, B3 for red, and B4 for the NIR band) in March. The results showed that among the feature set, phenological features had great significance for recognizing crops. In the top 20 important features, 15 features were phenological metrics, accounting for a large proportion. Additionally, time-related phenological metrics of field second time greenness made the greatest difference in classification (for instance, OnT_2, EndT_2, MaxT_2 and DT_2).

5. Conclusions

Crop mapping is essential for management of precision agriculture. There is much previous research regarding crop mapping that uses satellite imagery from a single source. Because of the spatial and temporal resolution constraints that are limited by single source satellite sensor performance, problems exist in mixed pixels or salt-and-pepper noise. Based on this issue, our study focused on identifying geo-parcel based crop types by collectively using high spatial and temporal resolution images.

We acquired field plot boundaries from fine-resolution images from GF-2 and constructed EVI time profiles, using medium-resolution images from Landsat 8 OLI and GF-1 WFV sensors. Considering that an analysis of time-series requires high consistency in geometric and radiometric features, we did much pre-processing work of data. Then, multi-temporal spectral information of medium-resolution images was assigned to farming plots. To smooth irregular EVI time-series, we presented a piecewise EVI time profile smoothing method to meet the requirements of rotated farming plots. This could effectively eliminate abnormal fluctuations, which helped with judgement of crop rotation and preserved greenness periods. Following that, the Random Forest classifier was used for crop type identification by integrating phenological metrics that were derived from the smoothed EVI time profile with multi-temporal spectral features, which produced satisfactory accuracies. We demonstrated the significance of phenological features in crop recognition.

In comparison with the S–G and HANTS smoothing methods, the three methods performed well in smoothing the EVI time profile of cultivated land without rotation. When smoothing rotated farmland, the proposed piecewise smoothing method maintained the basic shape of the curve while retaining greenness periods. However, S–G and HANTS behaved unsatisfactorily. The comparison of classification accuracies demonstrated that the S–G smoothing method outperformed the HANTS method in crop classification, and the piecewise fitting method worked best. Additionally, the integration of SF and phenological metrics improved the classification accuracies significantly. So, PF + SF led to the highest classification accuracy.

The limitations of this method mainly exist in the following two aspects. On one hand, for tiny and fragmented plots, the assignment of pixels from medium-resolution images to farmland objects could affect temporal spectral values. This had effects on the EVI time-series and spectral indices, which would eventually impact the classification accuracy to some extent. On the other hand, the threshold of segmenting the global time profile into local parts and judging rotation activity might be optimized in specific situations. These issues will be studied in our future work.