A 30-m annual maize phenology dataset from 1985 to 2020 in China

. Crop phenology indicators provide essential information on crop growth phases, which are highly required for agroecosystem management and yield estimation. Previous crop phenology studies were mainly conducted using coarse- 10 resolution (e.g., 500 m) satellite data, such as the moderate resolution imaging spectroradiometer (MODIS) data. However, precision agriculture requires higher resolution phenology information of crops for better agroecosystem management, and this requirement can be met by long-term and fine-resolution Landsat observations. In this study, we generated the first national maize phenology product with a fine spatial resolution (30m) and a long temporal span (1985-2020) in China, using all available Landsat images on the Google Earth Engine (GEE) platform. First, we extracted long-term mean phenological 15 indicators using the harmonic model, including the v3 (i.e., the date when the third leaf is fully expanded) and the maturity phases (i.e., when the dry weight of maize grains first reaches the maximum). Second, we identified the annual dynamics of phenological indicators by measuring the difference in dates when the vegetation index in a specific year reaches the same magnitude as its long-term mean. The derived maize phenology datasets are consistent with in-situ observations from the agricultural meteorological stations and the PhenoCam network. Besides, the derived fine-resolution phenology dataset agrees 20 well with the MODIS phenology product regarding their spatial patterns and temporal dynamics. Furthermore, we observed a noticeable difference in maize phenology temporal trends before and after 2000, which is


Introduction
Accurate and timely crop phenology information, which contains multi-phase growth information from sowing to harvest, is highly required for precision agriculture management (Gao and Zhang, 2021;Zeng et al., 2020), such as irrigation schedules 5 and pest control. The agriculture management schemes should be precisely scheduled according to different growth phases, during which period the water requirements and the possibilities of pest and disease events are different (Yang et al., 2021;Xiao et al., 2020). Besides, the effect of climate change on crop phenology has been widely reported (Abbas et al., 2017;Zhang and Tao, 2013;Tao et al., 2012). Given that the altered growth phases of crops will influence crop production, thus, further research into the response of crop phenology to global warming is necessary, which requires long-term records of phenology 10 change. In addition, information on crop phenology is also helpful for crop mapping because different crops vary in their growth phases (Sakamoto et al., 2014;Zhong et al., 2014;Zhang et al., 2014;Huang et al., 2019b).
Remote sensing has become a profound tool for deriving crop phenology at a large scale (Pan et al., 2015;Liu et al., 2018).
The annual variations of crop phenology are affected by many factors, including climate conditions, soil properties, and anthropogenic activities (e.g., sowing dates) (He et al., 2020). The traditional in-situ based crop phenology recording is time -15 consuming and focused on limited sites (Gao and Zhang, 2021). These limitations have been considerably mitigated by satellite images, which provide revisit observations of crop growth at regional and global scales (Shanmugapriya et al., 2019;Zhang et al., 2003;Cao et al., 2015). Different phenological indicators (such as the start of season and the end of season) are retrieved for crop growth monitoring using satellite observations, including the moderate resolution imaging spectroradiometer (MODIS) data (Sakamoto et al., 2010), the advanced very high resolution radiometer (AVHRR) data (Zhang et al., 2014;Gim et al., 20 2020). The retrieved multiple phenological indicators can delineate the development stages of crops from sowing to harvest at a regional and global scale.
Fine resolution Landsat satellite data show great potential in providing crop phenological indicators with a fine resolution and a long-term span. Despite that those coarse satellite data (such as MODIS and AVHRR) have a fine temporal resolution, which 3 is helpful to depict the crop growth phases, they are limited in their spatial resolution. Recently, several attempts have been made at deriving phenology datasets using fine resolution satellite data, such as Landsat Senf et al., 2017), Sentinel-2 (Bolton et al., 2020), and the harmonized Landsat8 and Sentinel-2 (HLS) data (Claverie et al., 2018;Bolton et al., 2020). Compared with medium-resolution satellite data such as MODIS, the Landsat satellite data can provide numerous land surface records from 1985 to the present, which helps derive the long-term crop phenology dynamics. Unfortunately, limited 5 attempts have been made using Landsat data to map the crop phenology with a fine resolution and a long-term span in China due to the complex planting patterns (Luo et al., 2020;Wu et al., 2010). Also, the computing resources required for such a mapping project are a huge challenge (Dong et al., 2016).
The advent of the Google Earth Engine (GEE) platform relieves the huge stress of data storage and computing at regional and global scales. The GEE platform has included petabyte-scale remote sensing data with high-performance computing 10 capabilities and powerful algorithm libraries (Gorelick et al., 2017). Presently, many successful attempts have been conducted using the GEE platform, such as mapping forest dynamics (Xiong et al., 2020), terrace (Cao et al., 2021), and surface water (Pekel et al., 2016). It is convenient to obtain and process satellite data using the GEE platform. The combination of massive satellite observations and a flexible development environment makes it possible to derive annual dynamics of crop phenology with fine resolution in China. 15 In this study, we extracted spatial and temporal patterns of maize phenology indicators in China from Landsat observations using the GEE platform. The derived phenology indicators include v3 (the date when the third leaf is fully expanded) and maturity (i.e., when the dry weight of maize grains first reaches the maximum) phase. We mapped annual phenological indicators of maize at a fine resolution (30m) from 1985 to 2020, using full achieve of Landsat images. The remainder of this paper is organized as below: Section 2 introduces the study area and datasets, Section 3 presents the method used in this study, 20 Section 4 and Section 5 describes the results with discussion and the derived dataset, respectively, and an ending mark is provided in Section 6.

2 Study areas and datasets
We selected China's main maize producing area as our study area (Fig. 1). Maize is one of the major crops in China and is planted over a wide region, the sown area and production account for 36% and 39% of food crops in 2019 (China Statistical Yearbook, 2020), respectively. The planting pattern and phenology of maize are highly heterogeneous due to the influence of climate conditions, soil properties, and anthropogenic activities (e.g., sowing date) (Wu et al., 2010). The spring maize is mainly 5 distributed in Northeast China, dominated by the single cropping type. However, summer maize is mainly planted in the Huang-Huai-Hai Plain (Fig. 1b), where the double cropping system (rotation between winter wheat and summer maize) is commonly seen (Luo et al., 2020). In addition, there is also a certain amount of maize in other provinces (e.g., Xinjiang province). The growth period of summer maize spans roughly from June (after the harvest of winter wheat) to October compared with that of spring maize from May to October. Furthermore, the maize in Northeast China is mainly rain-fed. In 10 contrast, irrigation is needed for maize commonly exists in Huang-Huai-Hai Plain and Northwest China (arid and semi-arid areas) (Wu et al., 2010). Under these diverse cropping systems, phenology dates (such as v3 and maturity) of maize varied significantly over space. We used the Landsat satellite data as the primary data source to characterize the phenology changes of maize in China. We used all available Landsat surface reflectance data in this study, including images obtained from Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Operational Land Imager (OLI), from 1985 to 2020. The Landsat surface reflectance data have been corrected for the radiometric and topographic effects. The atmospheric effect has also been corrected using the Landsat ecosystem disturbance adaptive processing system (LEDAPS) (Masek et al., 2006). Clouds and shadows 5 were removed using the function of the mask procedure (Zhu and Woodcock, 2012). Therefore, all available clear-sky pixels of Landsat observations over the past three decades were used in our study.
Maize maps from multiple resources were adopted to constrain the region of crop phenology mapping. The distribution map of maize in Northeast China was derived using Sentinel-2 data (You et al., 2021), and the used maize map was the classification result in 2019. You et al. (2021) derived the crop maps in Northeast China using a random forest classifier, with optimized 10 features including spectral, temporal, and textural characteristics (gray-level co-occurrence matrix). In other provinces, the maize maps were obtained using the "temporal similarity assessment" approach proposed by Dong et al. (2020). The distribution of maize is mainly determined by comparing the similarity of the vegetation index series of unknown pixels with a referred curve derived from maize fields. The retrieved maize datasets have been validated with survey data with reliable performance (Fig. S1). The accuracy of the maize map in Northeast China is 0.85 (more than 8000 samples for cross-validation 15 in 2019), and that of maize maps in other provinces is 0.79 (about 2000 samples for validation). Given that the original resolutions of these two classification maps are 10m (i.e., Northeast China) and 30m (i.e., other provinces), we aggregated the 10m maize map to 30m in our study. It is worth note the maize distribution map is consistent across different years in our study due to the relatively stable planting situation as one of the major crops (Sun et al., 2007;Li et al., 2008). Of course, we also admitted that certain dynamics exist in maize distribution due to the changing maize price, climate conditions, and choice of 20 farmers across different years. Mapping the maize dynamics at the national scale in China is challenging because of the scarcity of massive ground samples. There is also no public available maize dynamic product with fine spatial resolution and a long temporal span. So we kept the maize distribution map consistent and derived dynamics of maize phenology indicators with tolerable errors.

7
In addition, we also collected massive datasets to validate our results, such as the agricultural meteorological stations (AMS), PhenoCam network, and the MODIS phenology product (MCD12Q2). First, the records in AMS include phenology information of major crops (such as maize, wheat, and rice) in China, with large spatial and temporal ranges (Luo et al., 2020;Huang et al., 2019a). Crucial phases during the maize growth, including v3 (i.e., the date when the third leaf is fully expanded), and maturity (when the dry weight of maize grains first reaches the maximum) phases, were recorded in the AMS. Thus, this 5 dataset can validate the mapped phenological indicators from remote sensing, and we collected AMS phenology records of the spring and summer maize (Fig.2). Second, the in-situ PhenoCam observation was derived from digital cameras using the green chromatic coordinate (GCC) index, which is composited by visible wavebands and able to characterize the dynamic greenness of vegetation. Third, the MODIS phenology product (MCD12Q2) was also employed in our study to validate the results derived from Landsat observations. The phenological indicators (e.g., dates of green-up and dormancy) in the MODIS product were 10 mainly derived from the two-band enhanced vegetation index (EVI2) time series data (Gray et al., 2019). The multiple cycles (up to two) of crop rotations were also recorded in the MODIS phenology product, which is suitable for validation with our phenology results of maize in this study.

Methodology
We extracted the phenology indicators of maize using the full archive of Landsat images in GEE. The adopted framework includes three components (Fig. 3). First, we collected all available Landsat images during 1985-2020 in our study area and used the collected maize map as a mask. After the cloud removal, we constructed the long-term time series data of EVI for each pixel. Second, we fitted the long-term mean EVI curve using the harmonic model, which can delineate multiple cycles of 5 crop rotations and identify the number of cycles (i.e., used to distinguish spring and summer maize). Thus, two phenological indicators, the v3 (the date when the third leaf is fully expanded) and the maturity (when the dry weight of maize grains first reaches the maximum) phase, were determined from the long-term mean curve of spring and summer maize. Finally, we identified the annual dynamics of these two phenology indicators by measuring the difference in dates when the vegetation index in a specific year reaches the same magnitude as its long-term mean. Details of each component were given in the 10 following sections. Fig. 3: The adopted framework for deriving annual phenology dynamics  from Landsat time series data, including data preprocessing (a), mapping the long-term mean phenology indicators (b), and identifying the annual dynamics (c).

Data preprocessing
We implemented the data preprocessing step in the GEE platform. First, we used the quality layer in the Landsat surface reflectance data to remove clouds and shadows. Thus, all available clear-sky pixels can be used to enrich the Landsat 5 observations. Second, we excluded non-maize areas using the maize map, which can significantly reduce the computational requirement. Third, we calculated the EVI indicator using Eq. (1) to minimize the impact of soil and clouds; meanwhile, it is sensitive to vegetation growth and dormancy (Huang et al., 2019a;Li et al., 2019).
where NIR, RED, and BLUE represent surface reflectance of the corresponding spectral bands in Landsat. Parameters of G, 10 L, C 1 , C 2 were used to correct the disturbance of aerosols and soil background, as suggested with values of 2.5 (G), 1(L), 6(C 1 ), and 7.5(C 2 ) in Huang et al. (2019a).

Long-term mean phenology indicators
We derived the long-term mean maize phenology indicators including v3 and maturity. First. we sorted all available EVI observations according to their day of the year (DOY) and fitted the annual crop cycle using the harmonic model (Eq. 2). 15 Compared with other fitting approaches, the harmonic model can well delineate multiple seasonal cycles of crops within one year, with clear physical meanings for each parameter (de Beurs and Henebry, 2010;Chen et al., 2018;Lee et al., 2020).

∑ cos sin
( 2) where f(t) is the fitted EVI value, t is the Julian date of a particular observation, and T is the maximum value of the time variable. bi and ci are coefficients for intra-annual change of the EVI time series data. a1 and a0 represent the slope and 20 intercept of EVI change among different seasonal cycles. n represents the maximum number of harmonic components, and it needs to be calibrated according to different situations. Considering the double-crop (winter wheat-summer maize rotation system) and the planting patterns of winter wheat (planted in autumn of the first year and harvested in the second year), we set n as six due to the good fitting performance in our study after trial and error test using multiple sites in different regions.
Then, we identified spring and summer maize according to the cycles of the fitted curve (Fig. 4). Spring and summer maize can be identified using the information of EVI cycles. For instance, summer maize always has two crop cycles, notably different from spring maize, with only one cycle. To identify maize with different cycles, we calculated the derivative of the fitted harmonic model and identified the peaks of these cycles. When the EVI peak value before the maize part exceeded 40% of the maximum EVI value of the maize cycle (Gray et al., 2019;Wu et al., 2010), we regarded it as double cropping, and the 5 second part of it is summer maize (Fig. 4b); otherwise, it is a single crop (i.e., spring maize) (Fig. 4a). 10 Finally, we adopted the dynamic threshold approach to derive the v3 and maturity dates from the green-up and green-down segments (shadowed blue areas in Fig. 4). These two segments were derived from the cycle of maize growth, separated by the point with peak values of EVI. Given that the EVI amplitudes of green-up and green-down are different for the spring (g in Fig. 4a) and summer maize (g1 and g2 in Fig. 4b), we determined the v3 and maturity dates according to their EVI amplitudes accordingly. For the spring maize, the v3 and maturity dates were defined as the dates with 10% EVI amplitude (g in Fig. 4a) 15 during the green-up segment and 50% EVI amplitude (g in Fig. 4a) during the green-down segment (Huang et al., 2019a), respectively. Similarly, for the summer maize, the EVI amplitude during the green-up and green-down segments was referred to g1 and g2 in Fig. 4b to determine v3 and maturity, respectively.

Annual dynamics of phenology indicators
We adopted a similar approach to Li et al., (2017) to estimate the annual dynamics of phenology indicators. First, we adopted a self-adjusting strategy to determine the rational range of EVIs during the green-up and green-down periods (shaded areas in Fig. 5). These ranges were determined using the derived long-term mean curve, and they can be used to filter outliers in individual years. Then, we measured the difference in dates when the EVI in a specific year reaches the same magnitude as its 5 long-term mean (Fig. 5). The mean value of the date difference between the observations and the long-term mean was adopted as the annual variability of phenological indicators.

10
The solid and empty circles are long-term and year-specific enhanced vegetation index (EVI) observations. The definitions of all the acronyms are as follows: EOP: end of peak; SOP: start of peak.

Performance of the harmonic model
The harmonic model can well delineate the seasonal dynamics of EVI for spring and summer maize. Spring maize belongs to 15 the single cropping type and has one rotation cycle. In contrast, summer maize is mainly distributed in the double cropping system area (i.e., the second rotation cycle is summer maize). These crop growth cycles can be well detected by the EVI time series data from Landsat observations (Fig. 6). The fitting performance from the harmonic model suggests the fitted line can well delineate the growth phase of crops across different regions and types.

Comparison with records from the AMS
The derived long-term mean maize phenology indicators from Landsat observations are consistent with the records from the AMS (Fig. 7). We compared results derived from Landsat and AMS from 2001 to 2010, during which period the AMS observations can be maximally used. Due to the lack of accurate locations of observed crops in AMS (i.e., only station 10 locations), we measured the uncertainties of phenological indicators of maize within the range of 5km to the station. The 13 adopted approach performed well in extracting summer maize phenology. The correlations of v3 and maturity dates of summer maize are 0.60 and 0.80, respectively (Fig.7 c-d). Besides, the RMSE of v3 and maturity dates are 5.20 and 6.38 days.
Nevertheless, the correlation of derived maturity dates of spring maize and corresponding records from AMS is relatively low, likely attributed to the discrepancies in the definitions between remotely sensed results and AMS. For instance, the v3 phase in AMS is defined as the date when the third leaf is exposed from the second leaf sheath, and the maturity is defined as the 5 date when the dry weight of maize grains first reaches the maximum, more than 80% of the outer bracts of the plants turn yellow, and the filaments become dry . These definitions in AMS are challenging to be measured from remote sensing, and they are slightly different in terms of their definitions. Annual dynamics of derived phenology indicators (i.e., v3 and maturity) also agree well with the AMS observations (Fig.8). 5 The comparison of annual results is similar to that of the long-term mean phenology. In general, the annual dynamics of phenological indicators in summer maize are better than that of spring maize (especially at maturity phases), and this finding 15 is consistent with previous studies (Huang et al., 2019a). The correlations of phenology indicators (i.e., v3 and maturity) of summer maize derived from Landsat and AMS are 0.34 and 0.59, respectively. And for spring maize, the correlations of v3 and maturity indicators from the two datasets are 0.51 and 0.16. The difference between these two datasets is mainly attributed to (1) lacking accurate locations of the crop in the AMS data; (2) the crop planting patterns may be altered over the years (Fig.   9); and (3) different definitions. 5

Comparison with PhenoCam data
Using the phenology mapping approach in this study, we observed a good agreement between Landsat derived and PhenoCam 5 derived phenological indicators (Fig.10). We chose the United States (US) because no PhenoCam data are available in China.
The same approach adopted in China for crop phenological indicator mapping was used in the US with agriculture sites where 18 PhenoCam data are accessible. Thus, the feasibility of our approach can be evaluated. The phenology dataset provided by Richardson et al., (2018) is extracted from continuous observations of vegetation growth collected by digital cameras.
PhenoCam sites in Fig. 10 are mainly distributed in agriculture ecosystems, with records spanning from 2015 to 2018.
Definitions of phenology indicators from Landsat and PhenoCam are consistent i.e., definitions of transition_10 and transition_50 date when VI series data crossed 10% and 50% EVI2 amplitude (Richardson et al., 2018). The correlations of 5 v3 and maturity dates from Landsat and PhenoCam are 0.74 and 0.63, respectively, with the root mean square error (RMSE) of 7.61 (v3) and 7.11 (maturity) days. Observations from these two datasets are located around the 1:1 line, suggesting the adopted mapping approach of phenology from satellite data can well match the in-situ observations. Possible reasons behind explaining their difference can be attributed to (1) different vegetation indices used (i.e., EVI and GCC) and (2) the scale effect caused by their data sources. 10

Comparison with MODIS phenology dataset
The derived phenology indicators from Landsat and MODIS have a consistent temporal trend (Fig.11). MODIS phenology 5 product (MCD12Q2) provides multiple phenology indicators (e.g., midgreendown). For areas with two vegetation cycles, we 20 selected the phenology indicators of the second cycle (summer maize) for comparison. In the green-down segment of each crop cycle, the MCD12Q2 product provides three phenology indicators, i.e., dormancy, midgreendown, and senescence, defined as 90%, 50%, and 10% of the segment EVI2 amplitude in a specific cycle. We selected the midgreendown indicator in the MODIS phenology product to compare in this study because it has the same definition as the maturity date in Landsatderived results. We aggregated the fine-resolution maize data to the same resolution as MODIS and only kept those relatively 5 pure pixels (maize pixels accounting for more than 50% of them) for comparison. We found the temporal trends of derived phenology indicators of spring and summer maize from Landsat images are consistent with those derived from MODIS data in most years (Fig.11b). Our approach can well capture the crop growth phases' dynamics (i.e., delay and advancement). The magnitude difference between maturity date derived from Landsat observations and midgreendown derived from MCD12Q2 is within three days in most years. Different data sources and fitting methods (i.e., spline fit was used in MCD12Q2) likely 10 cause discrepancies between the two phenology datasets. In addition, it is worthy to note that there is a high correlation of the maturity dates derived from Landsat and MODIS (Fig.11c). corrections (c). Cases 1-2 are the spring maize, and case 3 is the summer maize. Each scene represents a 1.5km x 1.5km square. Note that solid lines represent the mean phenology at the regional scale, and the shadowed areas represent the range from 25 th to 75 th quantiles of 5 maturity date derived from Landsat. Phenology indicators derived from Landsat observations also have a close spatial pattern to the MODIS phenology product but more spatial details (Fig.12). For example, there is a noticeable advancement of maturity in 2018 and a delay in 2015 (red boxes in Fig. 12), and these variations are captured by the two phenology datasets successfully. Besides, we can note that Landsat-derived phenology indicators (e.g., maturity) depict the difference in crop growth stages with more spatial details 10 compared to the MODIS phenology product.

Analyze with climate data
Summer maize has a higher requirement for hydrothermal conditions (especially for temperature) than spring maize (Fig. 13). 5 Northeast China Plain and Huang-Huai-Hai Plain (Fig. 1b) are the two largest maize-producing areas in China, about 60% of 23 spring maize is grown in the Northeast China Plain (Fig. 13b), and more than 80% of summer maize is distributed in Huang-Huai-Hai Plain (Fig. 13c). In addition, we used the monthly mean air temperature and the total precipitation from May to October (Peng et al., 2019) in the study area for analyses from 2011 to 2020. Overall, the mean total precipitation and mean temperature change range within the spring maize planting areas are larger than in the summer maize. Meanwhile, summer maize growing areas are mainly distributed in high-temperature areas (i.e., above 20℃). These results suggest that summer 5 maize has a higher requirement for hydrothermal conditions than spring maize. We observed a noticeable difference in the temporal trends of the derived maize phenology indicators before and after 2000 (Fig. 14). The temporal trends of derived phenology indicators, including v3 and maturity date of spring and summer maize 15 before and after 2000, are notably different. For climate variables, the temperature within maize planting areas has a steeper 24 upward trend (the slope is more than 0.5℃ per decade) before 2000 than after 2000. The mean total precipitation shows different trends before and after 2000. It is worth noting that the precipitation within the spring maize producing area has a diverse and sharper tendency compared with that of the spring maize grown area. For phenological indicators, the changes in spring maize phenology mainly concentrated in the segments after 2000. The v3 date is advanced (-0.37 days per year), and the maturity date is delayed (0.38 days per year). The v3 and maturity indicators of summer maize have an advanced tendency 5 before 2000, while the maturity date is delayed after 2000. The annual dynamics of maize phenology indicators may be partly attributed to the rising temperature and annual variations of total precipitation. In this research, we did not consider the impact of other factors (such as photoperiod and genotype of maize) on the variations of maize phenology and the response of maize phenology and growth season duration to climate change was also not comprehensively considered. The temporal trends of phenological indicators (i.e., v3, maturity) and climate variables (i.e., mean temperature and mean total precipitation) during the growing period (from May to October), from 1985 to 2020. Two segments (i.e., 1985~2000 and 2001~2020) were independently fitted due to their distinct difference in temporal trends. We provided the temporal trend of variables across the study area and the interannual variations within different major agricultural zones.

Data availability
This dataset provides the annual dynamics of maize phenological indicators with a fine spatial resolution (30m) and a long temporal span  in China. The extracted phenology indicators include v3 (the date when the third leaf is fully expanded) and maturity (when the dry weight of maize grains first reaches the maximum). The format of this dataset is GeoTiff, with a spatial reference of WGS84. Each file in this dataset is named based on phenology indicator, start year, end year, and 5 province. We divided the maize phenology into three parts : 1985-2000, 2001-2010, and 2011-2020 (Table 1). We included 17 provinces in our study, i.e., Beijing, Tianjin, Hebei, Henan, Shanxi, Shaanxi, Shandong, Hubei, Anhui, Jiangsu, Inner Mongolia, Ningxia, Gansu, Xinjiang, Heilongjiang, Jilin, and Liaoning. The derived annual maize phenology data in China from 1985 to 2020 are available at https://doi.org/10.6084/m9.figshare.16437054 (Niu et al., 2021).

Conclusions
In this study, we generated the first annual maize phenology product with a fine spatial resolution (30m) and a long temporal span  in China, using all available Landsat images on the GEE platform. First, we extracted long-term mean phenology indicators (including v3 and maturity) from multi-year Landsat observations using the harmonic model. Second, 15 we identified the annual dynamics of phenological indicators by measuring the difference of dates when the EVI in specific years equals the fitted value.
The maize phenology product derived from Landsat scenes agrees with the commonly used phenology dataset. Our derived maize phenology datasets consistently meet the in-situ observations from the AMS and the PhenoCam phenology network. In addition, the phenology dataset in this study has similar temporal trends and can provide more spatial details than the MODIS 20 27 phenology product. In addition, we observed a noticeable difference in the temporal trend of maize phenology before and after 2000, which is likely attributable to increasing temperature and annual variations of precipitation.
The extracted maize phenology dataset has great implications for crop field management and studies of the response of maize phenology to the changing environment. There are noticeable differences in crop growth due to diverse local climates, soil properties, and anthropogenic activities (such as sowing dates). The derived phenology product with a fine spatial resolution 5 can delineate the difference and provide corresponding information to improve the field management and yield estimation (Zeng et al., 2020;Bolton and Friedl, 2013). In addition, this phenology product can also be used to investigate the response of crop phenology to global warming (Badeck et al., 2004;Niu et al., 2021). However, this study does not consider land cover changes (e.g., urban expansion and planting system change), which needs to be further investigated. For example, the maize distribution was regarded as persistent in our study over the past decades. 10

Author contributions
Quandi Niu, Xuecao Li, and Jianxi Huang designed the research, performed the analysis, and wrote the paper; Hai Huang, Xianda Huang, Wei Su, and Wenping Yuan revised the manuscript.