Revealing Cropping Intensity Dynamics Using High-Resolution Imagery: A Case Study in Shaanxi Province, China

Abstract

Reliable and continuous information on cropping intensity is crucial for assessing cropland utilization and formulating policies regarding cropland protection and management. However, there is still a lack of high-resolution cropping intensity maps for recent years, particularly in fragmented agricultural regions. In this study, we combined Landsat-8 and Sentinel-2 imagery to generate cropping intensity maps from 2019 to 2023 at a 10 m resolution for Shaanxi Province, China. First, the satellite imagery was harmonized to construct 10-day composite enhanced vegetation index (EVI) time series. Then, the cropping intensity was determined by counting the number of valid EVI peaks within a year. Assessment based on 578 sample points showed a high level of accuracy, with overall accuracy and Kappa coefficient values exceeding 0.96 and 0.93, respectively. We further analyzed the spatiotemporal patterns of cropping intensity and generated a map of abandoned cropland in Shaanxi. The results indicated that cropland in Shaanxi Province was mainly utilized for single-cropping (52.9% of area), followed by double-cropping (35.2%), with non-cropping accounting for 11.9%. Cropping intensity tended to be lower in the north and higher in the south. Temporally, the average cropping intensity of Shaanxi increased from 1.1 to over 1.3 from 2019 to 2023. Despite this upward trend, large areas of cropland were abandoned in northern Shaanxi. These results demonstrate the potential of utilizing Landsat-8 and Sentinel-2 imagery to identify cropping intensity dynamics in fragmented agricultural regions and to guide more efficient cropland management.

1. Introduction

Understanding the spatiotemporal dynamics of cropping intensity is pivotal for effective cropland management and sustainable agricultural development. Cropping intensity, defined as the cultivation frequency on a cropland plot within a year, serves as a crucial indicator of cropland utilization efficiency and productivity [1]. In China, the multiple cropping system, i.e., harvesting two or more crops within a single year, has a long history dating back to the first century [2]. This strategic approach to land use aims to increase crop production within the constraints of limited cropland area, thereby maximizing the utilization of agroclimatic resources [3].

Shaanxi Province, located in northwestern China, is a typical agricultural region characterized by diverse topography and climatic conditions. Shaanxi was one of three pilot provinces at the beginning of the Grain for Green program in 1999 in which steeply sloping cropland was converted to forest or grassland [4]. Consequently, there has been a substantial reduction in cropland area in Shaanxi [5,6]. Yet, large portions of abandoned cropland have also been observed in gently sloping areas [7]. Since 2020, Shaanxi has implemented several measures aimed at the recultivation of abandoned cropland to increase food supply. Thus, the period from 2019 to 2023 marked a critical juncture in the trajectory of cropland management in Shaanxi.

Remote sensing technology has played a crucial role in agricultural land use mapping since the 1970s [8]. Over the past few decades, research in crop mapping using remote sensing has grown significantly, with an increasing reliance on time series imagery to monitor crop phenology and seasonal variability [9,10]. Medium-resolution imagery, such as from the Moderate Resolution Imaging Spectroradiometer (MODIS), has been extensively used for generating cropping intensity assessments [11,12]. For example, [12] produced a global cropping intensity dataset from 2001 to 2019 using MODIS enhanced vegetation index (EVI) time series with a 16-day temporal resolution.

While earlier studies have made significant progress in crop intensity mapping using remote sensing, many still rely on data from a single sensor [8], which involves trade-offs between spatial and temporal resolution due to the physical limitations of individual sensors [13]. For example, Landsat imagery provides a high spatial resolution (30 m), but its temporal resolution is often insufficient to capture cropping patterns [14]. To address these limitations, spatiotemporal fusion techniques have been developed, combining fine- and coarse-resolution imagery to generate datasets with high spatiotemporal resolutions [15,16]. By fusing Landsat and MODIS data, it is possible to create 8-day-interval time series with a 30 m spatial resolution [17]. Additionally, the fusion of Landsat-8 and Sentinel-2 imagery offers even richer observations, achieving a temporal resolution of less than 5 days [18], enabling more detailed analysis of seasonal cropland dynamics. Given the substantial data volume and computational complexity involved, cloud computing platforms like Google Earth Engine (GEE) have become indispensable, offering planetary-scale processing capabilities for handling large datasets efficiently [19]. By leveraging its multi-petabyte data catalog and freely accessible computational service, GEE has been widely used for large-scale crop mapping and agricultural monitoring [20,21,22].

Despite these advancements, the potential of high-spatiotemporal-resolution imagery to accurately reveal cropping intensity remains underexplored, particularly in fragmented agricultural regions with complex cropping patterns like Shaanxi Province [23]. Although some studies have employed fused data to develop cropping intensity maps for individual years [14,24,25], these is still a notable gap in the availability of continuous, high-resolution cropping intensity maps for such challenging regions.

The primary objective of this study was to create a continuous, high-resolution cropping intensity dataset for Shaanxi Province from 2019 to 2023. Using a combination of Landsat-8 and Sentinel-2 imagery processed through GEE, we aimed to capture the spatiotemporal dynamics of cropping intensity and identify abandoned cropland in this region. Our findings will provide valuable insights to inform more effective cropland management strategies, and the methodology developed in this research offers the potential for broader applicability, facilitating continuous monitoring of the cropping intensity in other regions by leveraging multi-source remote sensing data.

2. Materials and Methods

2.1. Study Area

This study was conducted for Shaanxi Province (105°29′–111°15′E, 31°42′–39°35′N) in China, which has a total land area of 20.56 × 10⁴ km² (Figure 1a). Shaanxi exhibits a varied range of altitudes and climates. Geographically, it is characterized by three natural zones from north to south: the Loess Plateau, the Guanzhong Plain, and the Qinba mountainous region. Shaanxi spans three climate zones: temperate, warm temperate, and northern subtropical regions. In recent years, the annual accumulated temperature above 0 °C ranged between 4000 and 6000 °C (Figure 1b), while annual precipitation increased from 400 mm in the north to 1400 mm in the south (Figure 1c). Cropland is dispersed throughout the entire province, with the majority concentrated in central Shaanxi. The cropping system in Shaanxi transitions from single-cropping in the north to double-cropping in the south.

2.2. Workflow

The data and workflow for generating cropping intensity are presented in Figure 2. First, Sentinel-2 and Landsat-8 imagery were processed and harmonized to create a high-spatiotemporal-resolution dataset. Next, we constructed a consistent EVI time series through compositing, gap-filling, and smoothing techniques. Subsequently, cropping intensity was derived by counting the number of valid EVI peaks within a year. The accuracy of the generated product was then assessed using a large number of sample points.

2.3. Satellite Imagery

To gather cropland planting information, we utilized Landsat and Sentinel time series data spanning from 2019 to 2023. On a globe scale, the dual Sentinel satellite constellations (S2A and S2B) revisit every 3.8 days, and incorporating Landsat-8 data will further reduce the interval to 2.9 days on average [18]. The surface reflectance products of Landsat-8 Operational Land Imager (OLI) and Sentinel-2 MultiSpectral Instrument (MSI) were processed through GEE [19]. Low-quality observations were identified and excluded from further processing. We used the quality assurance layers QA60 and QA_PIXEL for Sentinel-2 and Landsat-8, respectively, to remove clouds and other low-quality observations.

Imagery harmonization was conducted to merge Landsat-8 and Sentinel-2 time series. We used the red, blue, and near-infrared bands from two sources: OLI bands 2, 4, and 5, and MSI bands 2, 4, and 8A [26]. First, given the spectral response similarities of OLI and MSI sensors, as shown in

Table 1. The spectral band characteristics of Landsat-8 OLI and Sentinel-2 MSI sensors and the transformation equations.

Spectral Band	OLI		MSI-2A		MSI-2B		Transformation Equation
	Central Wavelength (nm)	Bandwidth (nm)	Central Wavelength (nm)	Bandwidth (nm)	Central Wavelength (nm)	Bandwidth (nm)
Blue	482.0	60	492.7	65	492.3	65	OLI = 0.0003 + 0.9570 MSI
Red	654.6	37	664.6	30	664.9	31	OLI = 0.0041 + 0.9533 MSI
NIR	864.7	28	864.7	21	864.0	21	OLI = 0.0077 + 0.9644 MSI

and Figure 3, we employed linear models described by [27] to transform MSI data into OLI-like data. Then, to reconcile the difference in spatial resolutions, bands with 20 m or 30 m resolution were downscaled to 10 m. Finally, two normalized datasets with comparable spectral ranges and spatial resolutions were merged into a single dense dataset.

We counted the total number of observations and the number of good-quality observations per pixel in Shaanxi throughout 2021 as an example (Figure 4). On average, the Sentinel-2 constellations and Landsat-8 exhibited revisit cycles of approximately 2.9 days and 13 days, respectively. Their combination resulted in a more frequent observation interval, reduced to 2.4 days. After the removal of cloud-affected data, the frequency of good-quality observations remained reasonably high. To provide a comparative perspective, the revisit cycles for Landsat-8, Sentinel-2, and the combined Landsat-8 and Sentinel-2 datasets were approximately 21.2 days, 4.8 days, and 3.9 days, respectively.

2.4. WorldCover Product

We used the European Space Agency (ESA) WorldCover land cover product to prepare the cropland basemap (Figure 5a). This product was the first global land cover product at 10 m resolution based on both Sentinel-1 and Sentinel-2 data. The WorldCover 2021 product [30] represents the global land cover for 2021 and is the successive version of the initial product obtained for 2020. It was generated with an improved algorithm (v200 vs. v100), leading to an overall improved retrieval and an increased accuracy (76.6% vs. 74.7%). As shown in Figure 5a, extensive croplands in Shaanxi are predominantly situated within the Guanzhong Plain, characterized by the prevalence of moderate- and large-sized crop fields. Conversely, the Loess Plateau in northern Shaanxi, as well as the Qinba mountainous area of southern Shaanxi, exhibit small-sized crop fields.

The Third National Land Survey (TNLS) data of Shaanxi Province, released in 2022, served as the benchmark for evaluating the derived cropland basemap. This dataset resulted from tremendous efforts integrating field surveying, drones, and satellite imagery with a resolution of less than or equal to one meter [31]. By collecting the cropland area figures from the TNLS result for each city, ranging from 85.6 Kha to 926 Kha, we assessed the accuracy of the WorldCover cropland product. To achieve this, we calculated the cropland area for each city using the cropland basemap and compared it with the TNLS data, employing linear regression for analysis.

For cropland area at the city level, the derived cropland map from the WorldCover product exhibited reasonably good agreement with the TNLS dataset. The correlation coefficient between these two datasets was 0.85, and the regression line closely approximated the 1:1 line, as depicted in Figure 5b. For most cities, the data points were near to the 1:1 line. As a whole, the cropland area estimated from the WorldCover product was higher than the area from the TNLS dataset by about 354 Kha. Consequently, the WorldCover product shows promise as a reliable source for generating a land use map for Shaanxi, particularly for extracting the distribution of cropland.

2.5. EVI Calculation

EVI has shown a strong correlation with canopy structural attributes, such as leaf area index and chlorophyll content, making it a prevalent indicator of green vegetation [32]. EVI time series have been employed extensively to effectively identify crop phenological traits, reflecting various growth stages, such as time of planting and harvesting. Notably, different cropping intensities exhibit distinct EVI patterns within one year [11]. Therefore, EVI has been used in several studies for crop intensity mapping [11,12]. The EVI calculation formula is as follows:

$$\begin{array}{c}EVI=2.5\times {\displaystyle \frac{NIR-Red}{NIR+6\times Red-7.5\times Blue+1}}\end{array}$$

(1)

where NIR, Red, and Blue represent the surface reflectances in the near-infrared, red, and blue bands, respectively.

2.6. Time Series Construction

The minimum blue compositing method [11] was employed to create EVI products with uniform intervals. Within each compositing period, the EVI value with the lowest reflectance in the blue band was extracted. This approach is particularly effective in identifying the clearest conditions owing to the blue band’s sensitivity to aerosol contamination [33]. Considering the good-quality satellite observations and typical crop phenology in Shaanxi, a 10-day compositing period was chosen to generate 10-day composite EVI time series. Still, the composite time series might contain gaps due to a continuous absence of good-quality observations during the compositing period. Linear interpolation was utilized to fill these gaps by referencing EVI values before and after the respective time of the data gap.

We utilized the Savitzky–Golay filter [34] to smooth the composite EVI time series. Disturbances in the observations caused by snow and other factors often result in dramatic fluctuations within EVI curves. These fluctuations, theoretically, should not normally manifest themselves due to the limited duration of crop growth. Therefore, it was imperative to apply smoothing techniques to the composite time series before further analysis. The Savitzky–Golay filter is a widely used approach for reconstructing vegetation index time series [35,36]. This method performs polynomial least squares fitting over successive sets of adjacent values within a temporal window. The window size dictates the number of observations considered, while the polynomial degree regulates the order of the smoothing function. We implemented a moving window comprising nine composite dates, considering that the crop growth period typically exceeds 90 days in Shaanxi. A second-order polynomial was chosen to preserve subtle information related to crop phenology. Consequently, we generated a gap-filled and smoothly refined 10-day composite EVI time series for subsequent analysis (Figure 6).

2.7. Cropping Intensity Identification

The determination of cropping intensity relied on the extraction of phenological characteristics from smoothed EVI time series. Throughout a crop growing season, EVI values take a distinctive pattern: a gradual increase during the vegetative growth stage, reaching a peak when the canopy closes, and declining through the reproductive growth stage until harvest. Within one year, single-cropping systems will exhibit one distinctive cycle, whereas double-cropping systems will display two complete cycles. In contrast, non-cropping areas or fallow lands tend to have consistently low EVI values or significantly lower EVI peaks during the growing season. Therefore, the count of EVI peaks serves as a robust indicator of cropping intensity (Figure 6b). In this study, potential peaks in the EVI time series were identified as local maxima using a moving window spanning nine composite dates. Peaks with an EVI value lower than 0.35 were disregarded as spurious peaks [17]. Furthermore, peaks occurring after late October were excluded because winter wheat could exhibit a peak before entering dormancy. Cropping intensity was quantified by tallying the total count of peaks observed within the calendar year. Finally, we obtained the annual cropping intensity map for Shaanxi during 2019–2023.

2.8. Accuracy Assessment

The accuracy of cropping intensity quantification was assessed at the pixel scale using a stratified random sampling dataset. Sample points were randomly collected from each cropping intensity region estimated for the year 2019. The actual cropping intensity was independently determined through visual interpretation of EVI time series, very high-resolution imagery, and true color Sentinel-2 images. This process utilized tools such as Collect Earth v1.20.6 [37] and the Copernicus Browser [38]. Notably, during the interpretation, the points whose buffer zones contained noisy features such as ridges and trees were deliberately excluded. As a result, a total of 578 individual sample points were retained for the accuracy assessment (Figure 7). To enhance the independence of the validation, only sample points for the years from 2020 to 2023 were used for assessment, with each year containing more than 500 sample points. Subsequently, the confusion matrix and various metrics were calculated, such as overall accuracy and the Kappa coefficient. The overall accuracy serves as a fundamental metric for assessing the performance in terms of correct classifications across the entire dataset. It is determined as follows:

$$\begin{array}{c}{p}_o={\displaystyle \frac{t}{n}}\end{array}$$

(2)

where p_o represents the overall accuracy, t is the number of correctly classified samples, and n is the total number of samples. The Kappa coefficient is a statistic that measures the performance considering class imbalances and the possibility of agreement occurring by chance. It is calculated as follows:

$$\begin{array}{c}k={\displaystyle \frac{p_o-p_e}{1-p_e}}\end{array}$$

(3)

where p_e is the expected proportion of samples correctly classified by chance. The Kappa coefficient (k) ranges from −1 and 1, with a value of 1 indicating perfect agreement between the truth and predictions.

3. Results

3.1. Spatial Distribution of Cropping Intensity

The spatial distribution of cropping intensity and the corresponding area in Shaanxi Province during 2019–2023 are shown in Figure 8. In agricultural practice, single-cropping dominates in Shaanxi, followed by double-cropping, with non-cropping being the least common. In terms of the percentage of cropland area occupied by each cropping intensity, single-cropping consistently maintained the highest proportion of the total cropland area from 2019 to 2023, averaging approximately 52.9%. Double-cropping accounted for 35.2% of the total cropland area, while non-cropping occupied 11.9% of the total cropland area.

A considerable amount of cropland in northern Shaanxi has remained idle (Figure 9). For instance, in 2019, nearly 1/3 of the cropland in northern Shaanxi (366.4 Kha) had a cropping intensity of 0, constituting 72.3% of the non-cropping cropland in the entire province. This proportion remained consistently high, reaching or exceeding 90% from 2020 to 2023, peaking at 96.1% in 2021. This indicates that the majority of idle cropland was concentrated in northern Shaanxi, which is also evidenced by Figure 8. In contrast, there was very little idle cropland in central and southern Shaanxi, both in terms of absolute area and proportion, which were far lower than observed for the northern region. For example, in central Shaanxi, the proportion of idle cropland was 7.4% in 2019, the highest in the past five years, and subsequently remained below 2%.

3.2. Temporal Dynamics of Cropping Intensity

Temporally, the cropping intensity of most cropland in Shaanxi has undergone transitions over the past five years (Figure 10a). Only about 1/3 of cropland maintained a consistent cropping intensity during this period, including 100.4 Kha of non-cropping, 614.1 Kha of single-cropping, and 390.7 Kha of double-cropping cropland. In adjacent years, transitions primarily occurred between 1 and 2 or between 0 and 1. For example, from 2019 to 2020, the cropping intensity of 368.5 Kha of cropland increased from 1 to 2, while 267.0 Kha declined from 2 to 1. Additionally, 188.8 Kha and 160.9 Kha transitioned from 0 to 1 and from 1 to 0, respectively. When comparing 2019 and 2023, approximately 62% of croplands maintained their cropping intensity, with the remaining 656.8 Kha transitioning from 1 to 2, 233.1 Kha transitioning from 0 to 1, and 183 Kha and 104 Kha transitioning from 2 to 1 and from 1 to 0, respectively (Figure 10b).

The transition of cropping intensity has led to an increasing trend in the utilization of cropland in Shaanxi Province (Figure 10c). From 2019 to 2022, the area of double-cropping cropland steadily increased from 844.9 Kha to 1443.7 Kha. Although there was a slight decrease in 2023, the decline was relatively small. In contrast, the area of single-cropping cropland decreased annually from 2019 to 2023. The area of idle cropland exhibited a significant decrease in 2022, reduced by 45% compared to the preceding year. Consequently, the average cropping intensity of Shaanxi, which was 1.10 in 2019, increased to 1.37 in 2022 and 1.32 in 2023.

3.3. Distribution of Abandoned Cropland

The abandoned cropland map of Shaanxi Province was derived from cropping intensity maps, as shown in Figure 11. The definition of abandoned cropland varies across previous studies, particularly in terms of the duration of time for which the land remains uncultivated [39,40,41]. In this study, we followed the FAO’s definition, considering cropland that has been idle for at least five consecutive years as abandoned [42]. Cropland with a cropping intensity of zero for the period 2019–2023 was recognized as abandoned. In total, Shaanxi Province has 100.4 Kha (about 3%) of abandoned cropland, primarily concentrated in northern Shaanxi. At the city level, Yulin and Yan’an account for 90.0 Kha and 8.4 Kha of the abandoned cropland, respectively. Of the other regions, Weinan has the largest area of abandoned cropland, at approximately 0.58 Kha.

3.4. Accuracy of Cropping Intensity Based on Sample Points

The accuracy assessment, conducted on 578 sample points, exhibited a high level of accuracy in extracting cropping intensity. Detailed confusion matrix and accuracy metrics for estimating the cropping intensity from 2020 to 2023 are presented in

Table 2. The confusion matrix and accuracy metrics for estimating cropping intensity for Shaanxi Province in 2020–2023.

2020					2021					2022					2023
		Truth					Truth					Truth					Truth
		0	1	2			0	1	2			0	1	2			0	1	2
Prediction	0	76	4	0	Prediction	0	73	4	0	Prediction	0	30	3	0	Prediction	0	43	3	0
	1	1	255	4		1	0	247	2		1	4	265	1		1	1	252	1
	2	0	10	175		2	0	5	200		2	0	11	199		2	0	10	201
Overall Accuracy: 0.964					Overall Accuracy: 0.979					Overall Accuracy: 0.963					Overall Accuracy: 0.971
Kappa Coefficient: 0.940					Kappa Coefficient: 0.966					Kappa Coefficient: 0.933					Kappa Coefficient: 0.949

. On average, the overall accuracy and Kappa coefficient values reached 0.969 and 0.947, respectively. Moreover, for each year, the overall accuracy and Kappa coefficient values exceeded 0.96 and 0.93, respectively. Thus, the proposed approach could accurately detect the cropping intensity across multiple years in Shaanxi.

4. Discussion

4.1. Distribution of Cropping Intensity and Influencing Factors

Cropping intensity in Shaanxi generally exhibited a distribution pattern of lower in the north and higher in the south, aligning with the spatial pattern of thermal conditions. We examined the annual accumulated temperature above 0 °C from 2019 to 2023, a reliable indicator for assessing the suitability of multiple cropping systems [1,43,44]. This parameter was calculated as the sum of daily average air temperature consistently above 0 °C, with the beginning and end dates determined using a moving average method over five days. Based on the average values during 2019–2023 (Figure 1b), most parts of Shaanxi had accumulated temperatures $\ge 0$ °C ranging between 4000 and 6000 °C, with less in the north and more in the south. Recent research by [45] suggested that the minimum accumulated temperature requirement was about 4400 °C for double-cropping systems, and 6000 °C for triple-cropping systems. In northern Shaanxi, accumulated temperatures mostly ranged between 4000 and 4500 °C, making the area suitable for single-cropping, while central and southern Shaanxi had accumulated temperatures ranging between 4500 and 6000 °C, sufficient to support double-cropping, but inadequate for triple-cropping.

Topography likely plays the dominant role in shaping the distribution of idle cropland in Shaanxi. Topographical attributes have been demonstrated to have a significant influence on agricultural production, affecting the spatial arrangement of soil nutrients, water resources, and cultivation practices [46,47]. The abandoned cropland map (Figure 11) illustrates that idle cropland in Shaanxi is predominantly concentrated in high-altitude regions of northern Shaanxi. To delve deeper, we conducted a statistical analysis of the area of idle cropland under different relief degrees of land surface (RDLSs) (Figure 12). The RDLS, originally proposed by [48], is a synthetic representation of the regional altitude and surface cutting degree. Leveraging the RDLS dataset for China [49], we projected it to align with the cropping intensity map prior to analysis. The findings suggested that during 2020–2023, more than 80% of idle cropland was concentrated in regions where the RDLS ranged from 0.9 to 1.6, with a sparse distribution observed in flat areas.

4.2. Implication for Cropland Management

The rise in cropping intensity in Shaanxi was closely linked to the strong cropland protection policies implemented in recent years. China has devoted great efforts to cropland protection, and the 14th Five-Year Development Plan has reinforced the strictest cropland protection system to ensure national food security [41]. Guidance from the Ministry of Agriculture and Rural Affairs emphasized the urgent need to rehabilitate abandoned cropland in plain areas for cultivation. Our findings revealed a notable decline in idle cropland area in central Shaanxi in 2020 (Figure 8 and Figure 9). Additionally, the increasing trend in the utilization of cropland over the past five years reflects the efficacy of the “storing grain in land” strategy.

We found that there were still some idle croplands in the hilly and mountainous areas of northern Shaanxi, which should be the areas of focus for the recultivation of abandoned cropland. This region has complex terrain, with alternating ridges and hills, numerous gullies, and serious fragmentation of cropland. Additionally, there are some areas affected by salinization and wind erosion. Therefore, we recommend the development of diverse agricultural production based on the local conditions, and thereby growing suitable crops for different areas, such as staple crops, miscellaneous grains, and forage grasses, to increase the supply of diversified products.

The northern part of central Shaanxi Province also exhibited cropping intensity gaps. This region is primarily dominated by single-cropping cultivation, and there are also scattered idle croplands. In terms of the thermal conditions, the accumulated temperature ≥0 °C in this region has exceeded 4500 °C in recent years, which is sufficient to meet the requirements for double-cropping (Figure 1b). Additionally, the annual precipitation in this region has also surpassed 600 mm in recent years (Figure 1c). Furthermore, climate change may increase the production potential of this area in the future [50]. However, the lack of water infrastructure makes this region a typical rainfed agriculture area. On the one hand, rainfed farming techniques and water-saving irrigation can be promoted to elevate non-cropping to single-cropping. On the other hand, improving irrigation conditions through construction of water infrastructure could help further increase cropping intensity, thereby fully exploiting the production potential of thermal and water resources.

4.3. Potential of Landsat-8 and Sentinel-2 Imagery to Identify Cropping Intensity

The combination of Landsat-8 and Sentinel-2 has demonstrated significant potential for mapping cropping intensity. Globally, the median revisit interval of their combination is 2.9 days [18]. For Shaanxi Province, our results showed an even shorter interval of 2.4 days. After eliminating cloud-affected, low-quality data, the fusion still offers frequent observations at a 3.9-day interval. This short revisit cycle allows for the creation of dense time series, enabling the capture of seasonal variability and the assessment of cropping intensity. Additionally, the high spatial resolution enhances the ability to monitor fragmented cropland, which is common in Shaanxi and similar agricultural regions, as compared to previous studies [11,12]. The long historical archive of Landsat-8 and Sentinel-2 also supports continuous temporal analysis, facilitating the identification of cropland abandonment. The combination delivers the data needed for high-accuracy and consistent cropping intensity mapping, ultimately supporting more effective cropland management strategies.

4.4. Potential Sources of Uncertainty

Several factors may affect the accuracy of cropping intensity maps generated from the proposed algorithms. Firstly, we utilized the WorldCover product to derive the cropland basemap of Shaanxi Province. While the extracted cropland areas were largely consistent with the TNLS results, there were still some disparities. For example, the cropland areas appeared slightly smaller in southern Shaanxi (Figure 5b). This could be attributed to the high degree of fragmentation in these regions [51], where individual cropland plots may be smaller than the size of a single pixel, leading to misclassification as other land cover types. Secondly, in some areas of Shaanxi, the two-year–three-crop rotation system is practiced [52]. However, this study did not specifically distinguish this category, but rather broke it down into a single-cropping and a double-cropping classification. While this approach may slightly impact the calculated cropping intensity, its effect was limited due to the relatively small area of these plots.

4.5. Limitations and Future Directions

At present, there are numerous spectral indices available in addition to the EVI used in this study. Although EVI alone showed a high accuracy, incorporating additional indices could potentially enhance the robustness of the results. For example, the land surface water index may improve the differentiation of cropping cycles [23]. Furthermore, the cropping intensity generation method in this study relies on the shape characteristics of the EVI time series, which can be effectively extracted using deep learning techniques [8]. Applying deep learning to automate the feature extraction process could significantly improve the efficiency of cropping intensity assessments. In summary, while this study demonstrates the effectiveness of EVI in mapping cropping intensity, future research should consider integrating multiple spectral indices and deep learning techniques to provide more robust and reliable assessments.

5. Conclusions

We explored the potential of fused high spatiotemporal remote sensing imagery for revealing cropping intensity dynamics. Using Shaanxi Province, China, as a case study, we combined Landsat-8 and Sentinel-2 imagery to generate cropping intensity maps at a 10 m resolution from 2019 to 2023. Subsequently, we investigated the spatiotemporal characteristics of cropping intensity and the distribution of abandoned cropland. Our findings revealed that single-cropping dominated crop production in Shaanxi (52.9% of cropland area), followed by double-cropping (35.2%), with non-cropping accounting for 11.9%. Cropping intensity tended to be lower in the north and higher in the south. While the transition in cropping intensity led to an upward trend in cropland utilization, large areas of cropland were abandoned, particularly in northern Shaanxi. These findings demonstrate the capability of high-resolution imagery to identify cropland utilization dynamics in fragmented agricultural regions and offer valuable guidance for effective cropland management practices.