A Modiﬁed Shape Model Incorporating Continuous Accumulated Growing Degree Days for Phenology Detection of Early Rice

: Using a shape model (SM) is a typical method to determine the phenological phases of crops with long-time-series satellite remote sensing data. The average AGDD-based shape model (AAGDD-SM) takes temperature into account compared to SM, however, the commonly used daily average temperature is not sufﬁcient to determine the exact AGDD owing to the possibly signiﬁcant changes in temperatures throughout the day. In this paper, a modiﬁed shape model was proposed for the better estimation of phenological dates and it is incorporated into the continuous AGDD (CAGDD) which was calculated based on temperatures from a continuous 24 h within a day, different from the calendar day or the average AGDD indicators. In this study, the CAGDD replaced the abscissa of the NDVI growth curve over a 5-year period (2014 to 2018, excluding 2015) for a test site of early rice in Jiangxi province of China. Four key phenological phases, including the reviving, tillering, heading and anthesis phases, were selected and determined with reference to the ﬁeld-observed phenological data. The results show that compared with the AAGDD-SM, the method proposed in this paper has basically improved the prediction of each phenological period. For those cases where the average temperature is lower than the minimum temperatures (K1) but the effective accumulated temperature is not zero, more accurate AGDD can be calculated according to the method in this paper.


Introduction
The phenology of plants refers to a series of seasonal regular changes in the growth and development that adapt to the climate, such as germination, branching out, leaf development, flowering, fruiting, defoliation, and dormancy, and the corresponding dynamic phases of plants, which are called the phenological phases [1]. Phenology was widely applied in the agricultural field; accurate phenological information can help to develop reasonable agricultural production plans and facilitate crop cultivation and pest control [2][3][4][5][6][7][8].
Generally, there are three methods for observing the phenological phases [16]: manual observation, near-ground sensors and satellite-based remote sensing. Manual observation is the most common way for phenological phase determination over the past centuries and several phenological observation networks have been established across the world, such as the National Phenology Network (NPN) in the United States and the China Phenological Observation Network (CPON) in China [17][18][19]. Although the phenological network can provide accurate phenological information, it is time-consuming and lacks the capability of In this paper, a modified shape model was proposed for the better estimation of phenological dates, and it incorporated the continuous AGDD (namely CAGDD) which was calculated based on the temperatures of a continuous 24 h within a day, which is different from the calendar day or the average AGDD indicators. The advantage of this method is to use the whole-day duration combined with the sine method to calculate a more accurate AGDD. It is believed that the CAGDD can depict the contribution of temperature factors appropriately and reflect the growth of crops accurately. Compared with CD-SM and AAGDD-SM, the accuracy of CAGDD-SM is obviously improved, and this method can extract phenology more accurately. In addition, in order to prove the rationality of the combination of SM and the sine method in calculating AGDD, this paper also made two comparative experiments, which fully demonstrated the reasonableness of the proposed CAGDD-SM.

Study Area
The research area (117°4′43.38″E, 28°57′24.96″N) is located in Leping city, Jiangxi Province of China, covering an area of about 13 hectares (Figure 1). Rice is planted all year round in rotation of early rice and late rice.  Taking early rice as an example, the transplanting period of early rice usually starts from the end of April to the beginning of May, and early rice matures at the end of July (due to climate reasons, early rice is transplanted while not sown directly). The local early rice (Jiangxi Zhongzao 33) did not change in variety from 2013-2018. The reason for choosing this area was that this paddy field is the closest to the phenological observation station set by the Leping city Meteorological Bureau, which could therefore provide accurate ground observation phenological information of early rice. The observation station is located southeast of the paddy field, about 1.5 km away. Compared with previous research, the most prominent characteristic of this study area is that the selected area is not a standard and unified experimental field, but an ordinary piece of cultivated land, which is more practical. In addition, the Leping city Meteorological Bureau also provided daily hourly temperature data from 2014 to 2018. These years were consistent with the observation data of phenology that they could provide. The reliability of the data has been validated by professional departments.

Remote Sensing Data
The remote sensing data used in this paper was MODIS daily products-MOD09GQ (collection6, title: h28v06, 250 m) and MOD09GA (collection 6, title: h28v06, 500 m) from April to July 2014-2018. The main purpose of the MOD09GA products was to generate cloud masks to eliminate MOD09GQ products in which the study area was often covered by clouds. The mask file was generated by using the quality evaluation band in MOD09GA and then resampled to the same resolution as MOD09GQ. Finally, the cloud-containing images of MOD09GQ products in the study area were filtered by the mask file. Finally, the cloudless image of the study area could be obtained to generate the growth curve. The NDVI [45][46][47] values were calculated using band 1 (red band) and band 2 (near the red band) in MOD09GQ (Equation (1)). Taking the images of the study area in 2014 as an example, Figure 2 gives the NDVI images at different imaging times and it is shown that that the brighter pixels, which indicate higher NDVI values, can be observed in July compared to that in April.
(due to climate reasons, early rice is transplanted while not sown directly). The local early rice (Jiangxi Zhongzao 33) did not change in variety from 2013-2018. The reason for choosing this area was that this paddy field is the closest to the phenological observation station set by the Leping city Meteorological Bureau, which could therefore provide accurate ground observation phenological information of early rice. The observation station is located southeast of the paddy field, about 1.5 km away. Compared with previous research, the most prominent characteristic of this study area is that the selected area is not a standard and unified experimental field, but an ordinary piece of cultivated land, which is more practical. In addition, the Leping city Meteorological Bureau also provided daily hourly temperature data from 2014 to 2018. These years were consistent with the observation data of phenology that they could provide. The reliability of the data has been validated by professional departments.

Remote Sensing Data
The remote sensing data used in this paper was MODIS daily products-MOD09GQ (collection6, title: h28v06, 250 m) and MOD09GA (collection 6, title: h28v06, 500 m) from April to July 2014-2018. The main purpose of the MOD09GA products was to generate cloud masks to eliminate MOD09GQ products in which the study area was often covered by clouds. The mask file was generated by using the quality evaluation band in MOD09GA and then resampled to the same resolution as MOD09GQ. Finally, the cloudcontaining images of MOD09GQ products in the study area were filtered by the mask file. Finally, the cloudless image of the study area could be obtained to generate the growth curve. The NDVI [45][46][47] values were calculated using band 1 (red band) and band 2 (near the red band) in MOD09GQ (Equation (1)). Taking the images of the study area in 2014 as an example, Figure 2 gives the NDVI images at different imaging times and it is shown that that the brighter pixels, which indicate higher NDVI values, can be observed in July compared to that in April.

Methods
Cloud detection was conducted on the daily MOD09GQ product obtained from April to July between 2014 and 2018, and seventy-seven images were left (2015 was excluded due to the existence of too many cloudy images). The framework of the proposed method was shown in Figure 3. Firstly, the AGDD value for a certain day was obtained by the sine

Methods
Cloud detection was conducted on the daily MOD09GQ product obtained from April to July between 2014 and 2018, and seventy-seven images were left (2015 was excluded due to the existence of too many cloudy images). The framework of the proposed method was shown in Figure 3. Firstly, the AGDD value for a certain day was obtained by the sine fitting method from the continuous 24 h temperatures, and the NDVI-AGDD growth curve was obtained in which the AGDD value was the independent variable and the NDVI value was the dependent variable. Then this NDVI-AGDD growth curve was smoothed using the double-sine fitting method [48]. These smoothed curves from different years were then further averaged to obtain the improved shape model. Finally, the generated shape model was used to match the smoothed growth curve of another year to obtain the optimal matching parameters and was then combined with the predefined phenological parameters of the shape model to estimate the phenological period of the area. curve was obtained in which the AGDD value was the independent variable and the NDVI value was the dependent variable. Then this NDVI-AGDD growth curve was smoothed using the double-sine fitting method [48]. These smoothed curves from different years were then further averaged to obtain the improved shape model. Finally, the generated shape model was used to match the smoothed growth curve of another year to obtain the optimal matching parameters and was then combined with the predefined phenological parameters of the shape model to estimate the phenological period of the area.

Calendar Day-Based Shape Model
The calendar day-based shape model [38] assumes that the growth of the crop is correlated to the growth time, and therefore utilizes the calendar day as the independent variable to fit a growth trajectory curve of the crop for each year. Furthermore, the shape model can be built from the crop's multi-year growth curves. After that, the optimal matching parameters between the generated shape model and the growth curve for another year can be calculated, and these parameters are used to estimate the phenological dates. This kind of method mainly consists of the following steps: (1) Given the temporal NDVI values obtained from image sequences of the study area in year i, where ∈ 1, and M is the number of years to generate the shape model, the growth trajectory curve of year i can be fitted and further smoothed, as illustrated in Figure S1.
(2) The fitted curves for M years are then averaged to obtain the shape model, which can effectively suppress the inter-annual differences between the growth curves, as Figure  4 shows. With the constructed shape model, another growth curve from the target area for year j, where ∉ , , can be matched with the shape model to generate the optimal matching parameters, as Equation (2) shows:

Calendar Day-Based Shape Model
The calendar day-based shape model [38] assumes that the growth of the crop is correlated to the growth time, and therefore utilizes the calendar day as the independent variable to fit a growth trajectory curve of the crop for each year. Furthermore, the shape model can be built from the crop's multi-year growth curves. After that, the optimal matching parameters between the generated shape model and the growth curve for another year can be calculated, and these parameters are used to estimate the phenological dates. This kind of method mainly consists of the following steps: (1) Given the temporal NDVI values obtained from image sequences of the study area in year i, where i ∈ [1, M] and M is the number of years to generate the shape model, the growth trajectory curve of year i can be fitted and further smoothed, as illustrated in Figure S1.
(2) The fitted curves for M years are then averaged to obtain the shape model, which can effectively suppress the inter-annual differences between the growth curves, as Figure 4 shows. With the constructed shape model, another growth curve from the target area for year j, where j / ∈ [1, M], can be matched with the shape model to generate the optimal matching parameters, as Equation (2) shows: where x is the calendar day which is temporally resampled at equally spaced intervals and N is the corresponding number of sample points. f(x) is the growth curve from the target area for year j, and g(x) is a function of the generated shape model, as illustrated in Equation (3): where h(·) is the obtained shape model, xscale, yscale and tshift are the matching parameters to be calculated by minimizing Equation (2). The search ranges for each parameter were empirically determined as follows: 0.3 < xscale < 1.5, 0.3 < yscale < 1.5 and −100 < tshift < 100.
where x is the calendar day which is temporally resampled at equally spaced intervals and N is the corresponding number of sample points. f(x) is the growth curve from the target area for year j, and g(x) is a function of the generated shape model, as illustrated in Equation (3): where h(•) is the obtained shape model, xscale, yscale and tshift are the matching parameters to be calculated by minimizing Equation (2). The search ranges for each parameter were empirically determined as follows: 0.3 < xscale < 1.5, 0.3 < yscale < 1.5 and −100 < tshift < 100. An example of the matching procedure is demonstrated in Figure S2. (3) With the generated optimal matching parameters, the phenological phases of year j can be determined as Equation (4).
where is an estimated phenological date of year j, and is the preliminary defined phenological date of the shape model which can be obtained by averaging field phenological data of M years in the study area.

The Continuous Accumulated Growing Degree Days
According to the law of AGDD [39], the growth rate of a crop has a strong correlation with the temperature, and the crop can only continue its development state within a specific temperature range. Generally, a three-cardinal temperature system was commonly used to estimate the growth of crops in which three different temperature thresholds (K1-K3) were considered [44]. Given that the current temperature at any time is T, and are the minimum and maximum temperatures, respectively, allowing for the growth of the crop, which means the crop will not grow if or . is the optimum temperature which means the growth rate of the crop will increase with the current temperature if , but remains unchanged if . Table 1 gives the threecardinal temperatures of the early rice which are provided by the local Meteorological Bureau. An example of the matching procedure is demonstrated in Figure S2.
(3) With the generated optimal matching parameters, the phenological phases of year j can be determined as Equation (4).
where x est is an estimated phenological date of year j, and X s 0 is the preliminary defined phenological date of the shape model which can be obtained by averaging field phenological data of M years in the study area.

The Continuous Accumulated Growing Degree Days
According to the law of AGDD [39], the growth rate of a crop has a strong correlation with the temperature, and the crop can only continue its development state within a specific temperature range. Generally, a three-cardinal temperature system was commonly used to estimate the growth of crops in which three different temperature thresholds (K 1 -K 3 ) were considered [44]. Given that the current temperature at any time is T, K 1 and K 3 are the minimum and maximum temperatures, respectively, allowing for the growth of the crop, which means the crop will not grow if T < K 1 or T > K 3 . K 2 is the optimum temperature which means the growth rate of the crop will increase with the current temperature if Table 1 gives the three-cardinal temperatures of the early rice which are provided by the local Meteorological Bureau.  30 40 Considering the influence of temperature, Zhou et al. [43] introduced the AGDD value, which was calculated from the daily average temperature, to replace the calendar day as the independent variable of the shape model to reflect the contribution of temperature to the growth of the crops (Equation (5)). This method is named as AAGDD-SM in this paper.
where GDD is the calculated growing degree day value for day k, and AGDD is the sum of GDD values from day 1 to day n. T avg is the daily average temperature, T min and T max are the lowest and highest temperatures of day k, respectively. However, the average temperature used often underestimates the possible huge difference in the temperatures within a day. Moreover, the given assumption of Equation (5) is also flawed. For example, if the temperatures in one day are all greater than K 3 , it is deduced by Equation (5) that the crop will continue to grow, which is not supported by the AGDD law or common sense. In this paper, another criterion was designed for the proposed CAGDD-SM method according to the three-cardinal temperature system strictly, as shown in Equation (6).
where T is the temperature at any moment of the continuous 24 h. Therefore, a modified GDD value can be calculated by fitting the temperature curve with the hourly temperatures in one day (02:00 a.m. of day k to 02:00 a.m. of day k + 1). The sine fitting method was used to generate the temperature curve by determining three key points, as explained in Figure S3. The red curve is the fitted daily temperature curve. The general procedure was as follows: Given the fitted daily temperature curve, eight cases were considered giving the threecardinal temperature thresholds listed in Table 1, as shown in Figure S4. The red curve is the fitted temperature curve and the daily GDD value was represented by the shaded area which can be calculated, as Equation (7) shows, taking case F as an example: The generated shape models with the proposed continuous AGDD and the average AGDD (namely CAGDD-SM and AAGDD-SM) are shown in Figure 5.
where is the fitted temperature curve and is the intersection point of with the cardinal temperature lines.
Finally, the AGDD value for day k can be calculated by accumulating the GDD values, and it was further used to replace the calendar day in the shape model as the independent variable, to better reflect the correlation between the shape model and the development of the crop, thereby obtaining more accurate information about the phenological phases. The generated shape models with the proposed continuous AGDD and the average AGDD (namely CAGDD-SM and AAGDD-SM) are shown in Figure 5. Furthermore, the newly proposed CAGDD-SM method was used to match the growth curve of each year, and the optimal matching parameters of the corresponding year were obtained. These matching parameters can therefore be used to estimate the AGDD value for each phenological period of a certain year, as Equation (8) shows: in which is obtained by averaging the AGDD values corresponding to the ground observation dates of each phenological period in multiple years. Then, we can estimate the phenology in calendar date from AGDD. The specific conversion between them is described in Section 3.1. With the estimated phenological periods of different years, the accuracies of different phenological extraction methods for the years 2014-2018 (excluding 2015) can be calculated, as Equation (9) shows:  Furthermore, the newly proposed CAGDD-SM method was used to match the growth curve of each year, and the optimal matching parameters of the corresponding year were obtained. These matching parameters can therefore be used to estimate the AGDD value for each phenological period of a certain year, as Equation (8) shows: in which AGDDX S 0 is obtained by averaging the AGDD values corresponding to the ground observation dates of each phenological period in multiple years. Then, we can estimate the phenology in calendar date from AGDD. The specific conversion between them is described in Section 3.1. With the estimated phenological periods of different years, the accuracies of different phenological extraction methods for the years 2014-2018 (excluding 2015) can be calculated, as Equation (9) shows: where i represents the number of observation years, X i 0 refers to the actual observation value of the phenological period, and X i est refers to the predicted value of the phenology.

Supplement of Continuous Accumulated Growing Degree Days Model
According to the sine method, the daily temperature curve was fitted. Due to the influence of some abnormal weather, it seems that the fitted curve was not very reasonable, as shown in Figure S5a. Therefore, when fitting the temperature curve, it should be as close as possible to the real temperature curve. Then, according to the rule of calculating the growing degree day by the sine method, the integral of the temperature curve under this fitting method was calculated in combination with the cardinal temperatures, which was the growing degree day of this method. This fitting method was called the curve fitting method. As shown in Figure S5c,d, the curve fitting method was the sum of three sine functions, which was the closest fitting to the original curve.
In addition, because the daily hourly temperature in these years was known, the hourly GDD can be calculated by combining the cardinal temperatures, and the GDD of the day can be obtained by adding up the hourly value for 24 h and dividing it by 24 (Equation (10)). This method was called the time interval control method. GDD = ( 24 ∑ j=1 HGDD j )/24 (10) in which HGDD represents the hourly GDD of the day.

Phenological Extraction
In this experimental part, the observed phenological periods of years 2014-2018 (excluding 2015) were used to generate the averaging phenological parameters of the shape model, and also to evaluate different phenological extraction methods by matching their growth curves with the shape model to obtain the estimated phenological periods. In this way, the observed phenological periods of the years 2014-2018 (excluding 2015) can be used as the ground truth. Two shape model-based methods, CD-SM and AAGDD-SM were introduced for comparison with the proposed CAGDD-SM method. The generated average phenological periods, X s 0 , and their corresponding AGDDX s 0 values for the CD-SM and CAGDD-SM methods are listed in Table 2, respectively. Given the obtained optimal matching parameters for the different shape models with the annual growth curve, the estimated phenological periods of CD-SM can be obtained directly with Equation (4), while that of CAGDD-SM and AAGDD-SM can be obtained, as Tables 3 and 4 illustrates. Based on the GDD calculation, the relationship between the calendar day of different phenological dates and its corresponding AGDD value can be built, as shown in Tables 3 and 4. For example, the estimated AGDD value of the reviving phase of the year 2018 using the shape model is 133.18, and it is calculated from the daily 24-h temperatures that this AGDD value is between days 11-12. Therefore, the estimated calendar day of the Reviving phase for the year 2018 can be determined as the 11th day.
It can be seen from Table 5 that the total accuracy error of method CAGDD-SM (RMSE = 23.64) is better than that of method AAGDD-SM (RMSE = 25.31) and CD-SM (RMSE = 27.87); this has obvious advantages in the reviving and tillering phases. The error of tillering obtained by CAGDD-SM is 2.27 times smaller than that obtained by CD-SM. In any year, the maximum value, minimum value and average value of daily growing degree day in the first two phenological periods were much higher than those in the latter two phenological periods, especially in 2016, where the minimum value of growing degree day in the first two phenological periods was three times that in the latter two phenological periods ( Figure 6). In other words, the difference in AGDD obtained by observation between two adjacent days is much larger in the first two phenological periods than in the last two phenological periods. Because the GDD of the last two phenological phases is smaller than that of the first two phenological phases, it leads to the greater floating range of the predicted calendar day after the conversion of the last two phenological periods. Therefore, the estimation errors of the heading date and anthesis date will be relatively high. By observing the daily growing degree day in each phenological period, it was found that the daily growing degree day in the vegetative growth period was generally high, but the growing degree day in the reproductive growth period was relatively low. Therefore, when the estimated error of accumulated growing degree days was the same, the error of calendar time after conversion in the vegetative growth period will be smaller. It can be seen from Figure 7a that the error range of the phenological period estimated by the shape model based on calendar time was 0-10 days, and the error was about 10 days from the reviving period to the anthesis period. Compared with the former, most of the estimation errors of each phenological period in Figure 7b tended to be closer to the centerline, that is, the errors were smaller. The result shown in Figure 7b is illustrated in Figure 6, and it is also the conversion error caused by the difference between the GDD in different stages. Compared with the other two phenological periods, the daily GDD of the heading and anthesis are smaller, so when it is converted from AGDD to the calendar date, its floating days are larger, resulting in larger errors. Therefore, in the last two phenological periods in Figure 7b, there are some cases with large deviation days.

Verification of Rationality of Extracting Phenological Period by CAGDD-SM Method
In Figure 8, the red curve is the temperature curve fitted by the CAGDD-SM method, the blue curve is the temperature curve fitted by the curve fitting method (CFAGDD-SM), and the three horizontal lines refer to the cardinal temperatures.
It can be seen from Figure 8 that for some abnormal temperature curves, the curve fitted by the sine method was not as close to the original temperature curve as the curve fitting method. Therefore, in this paper, the curve fitting method (CFAGDD-SM) and time interval method (HAGDD-SM) were used to calculate the accumulated growing degree days, and the other steps were the same as the CAGDD-SM, and the phenological errors of these two methods were obtained (Table 6).   It can be seen from Table 6 that the effect of CAGDD-SM was the best, which was better than CFAGDD-SM. This was contrary to the expected result of CFAGDD-SM. The reason was that the law of accumulated growing degree days assumes that the temperature during the day and at night has the same effect on crop development, which was unreasonable. According to the daily hourly temperature data, the lowest temperature is at around 5:00 to 6:00 in the morning and the highest temperature is at around 2:00 to 4:00 in the afternoon. It can be seen from Figure 8 that the sinusoidal method enlarged the weight of the daytime and weakened the influence of night on the GDD (it was assumed that daytime is 6-18 o'clock). Therefore, CAGDD-SM was more effective than CFAGDD-SM. The date in Figure 8 was just an example. Any date can be used to explain that CAGDD-SM is better than CFAGDD-SM. For example, the two unusual date curves in Figure S5 can also be explained in this way. Similarly, this explains why HAGDD-SM was less effective than AAGDD-SM, because AAGDD-SM calculated GDD by the lowest and highest average values, and the average temperature was in the daytime range, which enlarged the daytime weight. While HAGDD-SM was calculated based on the same influence of day and night on GDD. Therefore, CAGDD-SM had the best effect.

Verification of Rationality of Extracting Phenological Period by CAGDD-SM Method
In Figure 8, the red curve is the temperature curve fitted by the CAGDD-SM method, the blue curve is the temperature curve fitted by the curve fitting method (CFAGDD-SM), and the three horizontal lines refer to the cardinal temperatures. It can be seen from Figure 8 that for some abnormal temperature curves, the curve fitted by the sine method was not as close to the original temperature curve as the curve fitting method. Therefore, in this paper, the curve fitting method (CFAGDD-SM) and time interval method (HAGDD-SM) were used to calculate the accumulated growing degree days, and the other steps were the same as the CAGDD-SM, and the phenological errors of these two methods were obtained (Table 6). It can be seen from Table 6 that the effect of CAGDD-SM was the best, which was better than CFAGDD-SM. This was contrary to the expected result of CFAGDD-SM. The reason was that the law of accumulated growing degree days assumes that the temperature during the day and at night has the same effect on crop development, which was

Conclusions
Temperature is the most important factor, which is used to extract phenology instead of calendar day as a dependent variable. As for the index of temperature-accumulated growing degree days, it is assumed that the growth rate of crops is positively correlated with temperature, and only when the active growing degree day exceeds the minimum temperature of crops will the crops begin to grow. The active growing degree day defined by accumulated temperature theory is the average value of daily minimum temperature and maximum temperature, which ignores the difference of daily temperature and misses the growing degree day for some days; although the average temperature is lower than the minimum temperature of crops, there are some moments when the temperature is higher than the minimum temperature of crops. Therefore, the sine method is used to calculate the growing degree day every day, which can avoid the above situations.
The daily GDD can be calculated more accurately using the sine method, and the accumulated growing degree days calculated using this method are combined with the shape model (CAGDD-SM). Compared with the accumulated growing degree days calculated using the traditional and common average temperature method, the estimation errors of each phenological period extracted by CAGDD-SM are better than that of the AAGDD-SM method, and the overall error is also better than that of the shape model with calendar time as the dependent variable. CAGDD-SM reduces the situation that some AGDD in AAGDD-SM are calculated as 0 but plants are still growing, so the prediction of the phenology of the former is improved compared to that of the latter. Although the precision of the heading date and anthesis date extracted by CAGDD-SM were worse than that without temperature, this phenomenon appears not only in this experiment; Zhou (2020) also showed that the precision of the reproductive growth period extracted by AAGDD-SM was worse than that of the vegetative growth period.
To verify the rationality of the CAGDD-SM method, this paper also uses the HAGDD-SM method and the CFAGDD-SM method to compare with CAGDD-SM, and the results show that the CAGDD-SM method is still the best. In addition, different from previous studies, the farmland in this study area is an area of ordinary cultivated land, which has neither advanced irrigation equipment nor unified and strict planting density planning. Agricultural producers rely on their own experience to plant and irrigate this cultivated land. Therefore, the research in this paper has more practical significance and achieved good results. Although the method in this paper needs daily hourly temperature data, which may be more difficult to obtain than the average temperature method, the method provided in this paper is based on a more accurate and reasonable AGDD calculation method to predict, and the obtained phenological period is more accurate. For example, for a specific date, the average temperature is less than K1. According to the calculation of the AAGDD-SM, the GDD on that day is 0, but the plants have not stopped growing. CAGDD-SM can greatly reduce these situations. Taking early rice as an example, the transplanting period of early rice usually starts from the end of April to the beginning of May, and the early rice matures at the end of July. The daily temperature difference during this period is not as obvious as that in winter, so the advantages of method A are not so prominent. The emphasis of this paper is to provide a more accurate phenological prediction method.
Therefore, the research in this paper has great potential in crop phenology estimation at a regional scale and provides a lot of useful information for agricultural production. Accurate phenological information can help to make a reasonable agricultural production plan, and facilitate work such as crop cultivation, pest control, and crop yield estimation. The research in this paper has the potential to be applied to the phenological period estimation of other crops, such as cotton, soybean, and corn. Future work is based on the existing research, considering the difference of accumulated growing degree days between day and night, to further improve the estimation of crop phenological period. In addition, the emphasis of this paper is to consider more reasonable and accurate temperatures based on SM, and more factors can be considered in the future, such as the influence of water, soil and different vegetation indexes. This paper adopts a more reasonable method to calculate AGDD and combines it with the shape model to get a more accurate estimation of the phenological period. Finally, the rationality of this proposed method (CAGDD-SM) is also verified and can provide useful help for agricultural production planning.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs14215337/s1, Figure S1: The fitted growth curve from MODIS NDVI data for year i; Figure S2: Matching of shape model with growth curve; Figure S3: An example of 24-h temperature fitting; Figure S4: Eight cases to calculate the GDD value. (A,B) refer to the case where the daily maximum temperature is greater than K1 and less than K2. (C-E) refer to the case where the daily maximum temperature is greater than K2 and less than K3. (F-H) refer to the case where the daily maximum temperature is greater than K3; Figure S5: The temperature curves fitted by the sine method and the curve fitting method. (a,b) are temperature curves fitted by the sine method; (c,d) are temperature curves fitted by curve fitting method.  Data Availability Statement: Ground observation phenological information and daily hourly temperature data are not available. Data were obtained from the Leping city Meteorological Bureau but are confidential. Remote sensing data can be accessed and downloaded from the United States Geological Survey (USGS). This data can be found here: https://ladsweb.modaps.eosdis.nasa.gov/search (accessed on 1 July 2022).