Method and system for identifying crop growth cycle and type

Abstract

The present invention discloses a method and system for identifying crop growth cycles and types, including: obtaining remote sensing images of crop growth in a long time series in a research area and preprocessing them; calculating and obtaining multiple vegetation indices and filtering them to form a feature grid on the time series, performing random point sampling on the feature grid to form a data set, and dividing the data set into a training set and a test set; constructing an XGBoost model, training the XGBoost model with a training set, and then testing the trained XGBoost model with a test set, selecting an XGBoost model with accuracy that meets the requirements as an optimized XGBoost model; obtaining remote sensing images of crops to be predicted, calculating multiple vegetation indices and inputting them into the optimized XGBoost model for prediction. The present invention not only ensures the stability of data quality during cross-year classification, but also eliminates numerical anomalies caused by interference factors such as clouds, thereby improving the prediction accuracy of the model.

Description

Technical Field

[0001] The invention belongs to the technical field of agricultural research and management, and specifically relates to a method and system for identifying crop growth cycles and types. Background Art

[0002] Satellite remote sensing technology is a key tool for modern environmental monitoring and earth science research. Multispectral images have many applications in agriculture, hydrology, environmental protection and other fields due to their good temporal continuity, wide information coverage area and large amount of information. In the agricultural field, satellite remote sensing provides detailed crop growth data, soil moisture information, and pest and disease monitoring through long-distance monitoring of the surface, which greatly improves the accuracy and timeliness of agricultural resource management.

[0003] Remote sensing technology can identify crop types and estimate planting areas through satellite images. By analyzing multispectral images obtained at different time points, it is possible to monitor the growth of multiple stages in the crop growth cycle, including germination, growth and harvest. Typical examples include the MODIS (Moderate Resolution Imaging Spectroradiometer) and LANDSAT satellites of the National Aeronautics and Space Administration (NASA), which provide rich agricultural monitoring data. Through time series analysis of these data, agricultural inputs can be adjusted in a timely manner, such as the best time for irrigation and fertilization. Satellite remote sensing is also partially used in assessing soil moisture. The data obtained by the European Space Agency (ESA) SMOS and NASA's SMAP (Soil Moisture Active Passive) satellite projects can accurately reveal soil moisture conditions, which is of great significance for guiding agricultural irrigation management. Such information can also help assess drought conditions, predict the decline in grain production, and take corresponding countermeasures. Remote sensing data can also be used to detect and monitor the development of pests and diseases. Using the reflectivity of specific bands, it is possible to detect subtle changes in crop leaves, so as to detect potential problems in the early stages of the disease, so that farmers can take targeted treatment measures in time to prevent the spread of the disease.

[0004] Taking China's northern agriculture as an example, by analyzing years of satellite data, researchers can identify wheat and corn planting patterns, evaluate water utilization during different rice growing stages, provide reference for regional irrigation planning, and provide a basis for precise pesticide application through disease monitoring.

[0005] Due to weather factors, cloud areas often appear in multispectral remote sensing images. These cloud areas not only cause the loss of collected information, but also bring great difficulties to further processing tasks such as target detection and object recognition, thus producing erroneous results in subsequent analysis. At the same time, when we consider time series information as input, clouds will also destroy the integrity of time series in different regions, thereby restricting the application of remote sensing images. Summary of the invention

[0006] The purpose of the present invention is to provide a method for identifying crop growth cycles and types in view of the shortcomings of the prior art. The method can eliminate numerical anomalies caused by interfering factors such as clouds and improve the prediction accuracy of the model.

[0007] In order to solve the above technical problems, the present invention adopts the following technical solutions:

[0008] A method for identifying crop growth cycles and types comprises the following steps:

[0009] Step 1: Obtain a long-term series of crop remote sensing images of the study area and pre-process the remote sensing images;

[0010] Step 2: Calculate the preprocessed remote sensing image to obtain multiple vegetation indices, sort the calculated vegetation indices by time and filter them, extract the values ​​of each vegetation index of each pixel in multiple growth cycles from the filtered time series, form the extracted values ​​of each vegetation index of each growth cycle into a new band and superimpose them to form a feature grid, then perform random point sampling on the feature grid to obtain a data set, and divide the data set into a training set and a test set;

[0011] Step 3: Build an XGBoost model for predicting crop growth cycle and type, train the XGBoost model with the training set, test the trained XGBoost model with the test set, and select the XGBoost model with the required accuracy as the optimized XGBoost model;

[0012] Step 4: Obtain the remote sensing image of the crop to be predicted, calculate its multiple vegetation indices, and then input the multiple vegetation indices into the XGBoost model optimized in step 3 to predict the planting type and growth cycle of the crop.

[0013] Furthermore, the preprocessing in step 1 includes performing spatial correction and atmospheric correction on the remote sensing image.

[0014] Further, the multiple vegetation indices in step 2 include LSWI, GLI, RNDVI, MSAVI, IRECI and SAVI.

[0015] Furthermore, the calculation formulas of multiple vegetation indices are as follows:

[0016]

[0017]

[0018]

[0019]

[0020]

[0021]

[0022] In the formula, NIR represents the near infrared band; SWIR represents the short infrared band; G represents the green band; R represents the red band, B represents the blue band; EDGE1 represents the red edge 1 band; EDGE2 represents the red edge 2 band, EDGE3 represents the red edge 3 band; L is the adjustment parameter.

[0023] Furthermore, the method for filtering the vegetation index in step 2 includes:

[0024] The improved SG filtering method is used to smooth the vegetation index, and the formula of the filtering function is as follows:

[0025] g i = ;

[0026] In the formula, f i Represents the original data value in the time series; g i is the smoothed value, which is c n and f i The linear combination of n is the moving window width for filtering, nL and nR correspond to the left and right edges of the signal component respectively; c n is a quadratic polynomial fitting function, then by solving the least squares fit, ranging from nL to nR, we get c n , quadratic polynomial fitting function c n for:

[0027]

[0028] Here, t corresponds to the day of the current year in the NDVI time series.

[0029] Furthermore, before filtering the vegetation index, the Kriging interpolation method is used to interpolate the vegetation index. The specific formula is:

[0030]

[0031] Where: Z(x0, y0) is the weighted sum of the values ​​of n sampling points; Z(x i ,y i ) is the sampling point (x i ,y i ) value; i The weight coefficient is determined by the variation function, i=1,2,…,n.

[0032] Furthermore, in step 3, the test set labels are compared with the model prediction results, and the accuracy and Kappa coefficient are used to evaluate the classification accuracy of the XGBoost model in various scenarios.

[0033] Furthermore, the calculation formula of Kappa is as follows:

[0034]

[0035] In the formula, represents the proportion of correct classification; is the probability of a random prediction, which is the sum of the products of the actual and predicted numbers for all classes, divided by the square of the number of samples.

[0036] Another object of the present invention is to provide a system for implementing the above-mentioned method for identifying crop growth cycles and types, comprising:

[0037] The data acquisition module is used to obtain long-term series of remote sensing images of crops in the study area and pre-process the remote sensing images;

[0038] The vegetation index calculation module is used to calculate the remote sensing images preprocessed by the data acquisition module to obtain multiple vegetation indices;

[0039] The data set processing module is used to filter the vegetation index calculated by the vegetation index calculation module after sorting them by time, extract the values ​​of each vegetation index of each pixel in multiple growth cycles from the filtered time series, form the extracted values ​​of each vegetation index in each growth cycle into a new band and superimpose them to form a feature grid, then perform random point sampling on the feature grid to obtain a data set, and divide the data set into a training set and a test set;

[0040] The model building and training module is used to build an XGBoost model for predicting crop growth cycles and types. The XGBoost model is trained using a training set, and then the trained XGBoost model is tested using a test set. The XGBoost model with the accuracy that meets the requirements is selected as the optimized XGBoost model.

[0041] The prediction module is used to obtain the remote sensing image of the crop to be predicted, and input it into the vegetation index calculation module to calculate multiple vegetation indices, and then input the multiple vegetation indices into the optimized XGBoost model to predict the planting type and growth cycle of the crop.

[0042] Compared with the prior art, the beneficial effects of the present invention are as follows: compared with the existing classification task using single-date remote sensing images, the present invention collects remote sensing images within a longer time range in the study area and arranges them in chronological order, calculates the selected characteristic values ​​of each image, designs a filter function and initializes it, constructs an independent filter function for each pixel point, fits the filter function for each point, extracts node feature structures from the filter function of each point after fitting, and finally constructs a crop classification data set by a random sampling point method. This data set construction method can take into account time change information, namely the crop growth cycle, to ensure the stability of data quality during cross-year classification, and can also eliminate numerical anomalies caused by interference factors such as clouds, thereby improving the prediction accuracy of the model. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 A flowchart of a method for identifying a crop growth cycle and type according to an embodiment of the present invention;

[0044] Figure 2 It is a comparison diagram between the SG filtering method and the gas filtering method according to an embodiment of the present invention. DETAILED DESCRIPTION

[0045] The technical solutions in the embodiments of the present invention will be described clearly and completely below in combination with the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0046] It should be noted that, in the absence of conflict, the embodiments of the present invention and the features in the embodiments may be combined with each other.

[0047] The present invention will be further described below in conjunction with specific embodiments, but the present invention is not limited thereto.

[0048] like Figure 1 As shown, the embodiment of the present invention discloses a method for identifying the growth cycle and type of crops, comprising the following steps:

[0049] Step 1: Obtain a long-term series of crop growth remote sensing images of the study area and pre-process the remote sensing images;

[0050] The study area of ​​this embodiment is located in Heilongjiang Province. Heilongjiang Province is the northernmost and highest latitude province in China, and has a cold temperate and temperate continental monsoon climate. The average annual temperature is about 5°C. The average annual rainfall is about 550mm. The average annual sunshine time is between 2200 hours and 2500 hours. In the study area of ​​this embodiment, there are a total of 10 farms with a planting area of ​​about 5211893 hectares. The main large-scale crops are rice, corn, soybeans and wheat. The remote sensing image data of the study area uses Sentinel-2 L2A images as the data source. The Sentinel-2 L2A images are image products obtained after pre-processing such as atmospheric correction on the basis of spatially corrected L1C images. The image includes a total of 13 bands. The bands with a resolution of 10m include three visible light bands and an infrared band, with central wavelengths of 496.6nm (B2), 560nm (B3), 664.5nm (B4) and 835.1nm (NIR), respectively. The bands with a resolution of 20m include vegetation red edge (B5, B6, B7), narrow near infrared (B8A), and shortwave infrared (B11, B12), with central wavelengths of 703.9nm (B5), 740.2nm (B6), 782.5nm (B7), 864.8nm (B8A), 1613.7nm (B11) and 2202.4nm (B12), respectively. The coastal aerosol (B1) and bands with a resolution of 60m include water vapor (B9), with central wavelengths of 443.9nm (B1) and 945nm (B9), respectively. The Visible-light Sentinel-2 mission consists of two polar-orbiting satellites that are located in the same sun-synchronous orbit with a phase difference of 180° and have a 3-day revisit period in the area studied in this study, acquiring Sentinel-2 images of the study area from April 1 to September 1.

[0051] Step 2: Calculate the preprocessed remote sensing images to obtain multiple vegetation indices, and annotate them according to the crop types in the study area during the period to obtain a data set, and divide the data set into a training set and a test set;

[0052] In this step, based on the Sentinel-2 L2A remote sensing image in step 1, this embodiment calculates a total of 6 vegetation indices, namely LSWI, GLI, RNDVI, MSAVI, IRECI and SAVI. The specific calculation formula is as follows:

[0053]

[0054]

[0055]

[0056]

[0057]

[0058]

[0059] In the formula, NIR represents the near infrared band; SWIR represents the short infrared band; G represents the green band; R represents the red band, and B represents the blue band; EDGE1 represents the red edge 1 band; EDGE2 represents the red edge 2 band, and EDGE3 represents the red edge 3 band; L is an adjustment parameter, which is used to adjust the vegetation index to take into account the image caused by the degree of soil surface exposure. In this embodiment, the value is 0.5.

[0060] In this embodiment, the above vegetation indices are calculated for the Sentinel-2 images obtained from April 1 to September 1 of the study area. After the vegetation indices are sorted by time, the improved SG filter is used to perform a time series smoothing algorithm on multiple vegetation indices to eliminate numerical anomalies caused by interference factors such as clouds. After processing, the filtering results of each feature on the 15th of each month from April to August are extracted, and a total of 30 features are obtained. This process is implemented on Google Earth Engine (GEE). In this embodiment, the SG filtering method is improved, and the Savitzky-Golay (SG) smoothing method, also known as the least squares method or digital smoothing polynomial, can be used to smooth noisy signals. First, the filtering function is constructed:

[0061] g i = ;

[0062] In the formula, f i Represents the original data value in the time series; g i is the smoothed value, which is c n and f i linear combination of; n is the width of the moving window for filtering, nL and nR correspond to the left and right edges of the signal component respectively; if c n is a constant, defined as c n =1 / (nL+nR+1), then SG filtering becomes sliding average smoothing. The idea of ​​SG filtering is to find the filter coefficient c that retains the high-order moment n Therefore, when c n is not a constant, but a polynomial fitting function, usually a quadratic or quartic function, and then solve the least squares fit, ranging from nL to nR, to get c n In this embodiment, for a specific data set of time series within a moving window, the polynomial fitting function is defined as a quadratic polynomial to fit a specific range of fi:

[0063]

[0064] Here, t corresponds to the day of the current year in the NDVI time series.

[0065] However, SG filtering cannot process data at different time intervals, so proper interpolation of vegetation index data before filtering and reconstruction can not only facilitate further processing, but also greatly improve the credibility of medium, low quality and noise pixels. In this embodiment, Kriging interpolation is used for interpolation. Kriging interpolation believes that the change of attributes in space is very irregular and discontinuous. Using a simple smoothing function will result in more errors. A variation function that can reflect spatial variability should be used. The specific formula is shown as follows:

[0066]

[0067] Where: Z(x0, y0) is the weighted sum of the values ​​of n sampling points; Z(x i ,y i ) is the sampling point (x i ,y i ) value; i (i=1, 2, ···, n) is the weight coefficient determined by the variogram. Different variogram fitting models are suitable for different spatial environments. Common variograms include exponential, Gaussian, and spherical models. This study uses a spherical model. The variogram model has three parameters: range is the range of the model, that is, the effective range of the interpolation result; sill is the base value, that is, the constant variance after exceeding the range; nugget is the nugget value, that is, the semi-variance value when the distance is 0.

[0068] Figure 2 In the figure, the filtering diagrams of the improved SG algorithm, the SG algorithm and the HUNTS filtering method are shown respectively. It can be seen from the figure that the improved SG filtering method can best retain the time series information of the index while filtering out the disturbance information.

[0069] After filtering the six vegetation indices, the values ​​of each pixel on the 15th of April to August were extracted from the smoothed time series, and the values ​​of each index on each date were combined into new bands and superimposed to form a feature grid. The data set required for model input was produced by random point sampling of the feature grid.

[0070] Step 3: Build an XGBoost model for predicting crop growth cycle and type, train the XGBoost model with the training set, test the trained XGBoost model with the test set, and select the XGBoost model with the required accuracy as the optimized XGBoost model;

[0071] This embodiment constructs an XGBoost model for crop classification prediction. XGBoost is an optimization of a gradient boosting machine (GBM) model, which also belongs to a forward ensemble machine learning method, aiming to integrate multiple weak learners into an efficient model. The algorithm iteratively constructs new base learners that are highly correlated with the negative gradient of the loss function to maximize model performance. In XGBoost, the use of the squared loss function not only penalizes large deviations of the target output, but also ignores small residuals, making the calculation of pseudo-residuals and the direction of the optimal solution in each iteration more efficient. Through the clever use of the squared loss function, XGBoost pays more attention to the accurate handling of errors during training, thereby improving the robustness and generalization performance of the model.

[0072] After establishing the basic framework of the XGBoost model, the training set is used to train the XGBoost model. The training process is as follows:

[0073] Initialize model parameters: determine hyperparameters such as the number of trees (n_estimators), tree depth (max_depth), learning rate (learning_rate), and loss function (objective);

[0074] Calculate the initial prediction value: For regression problems, the average value of all training samples is generally used as the initial prediction value; for classification problems, the initial prediction value can also be determined based on the category distribution and other conditions.

[0075] Iterative training:

[0076] a. For each tree: Calculate the negative gradient: Calculate the negative gradient of the loss function based on the prediction value of the current model. Fit the residual: Use the training data to fit the regression tree (or classification tree) with the negative gradient. Update the model: Update the model's prediction results based on the fitted tree; b. Determine the optimal split point based on the loss function and update the tree structure;

[0077] Regularization: XGBoost includes regularization terms for both tree complexity and model complexity to prevent overfitting.

[0078] Termination condition: When the preset number of trees is reached, or a certain termination condition (such as the depth of the tree, the minimum number of samples in a node, etc.) is met, the iteration stops.

[0079] Finally, the accuracy of the XGBoost model is evaluated through the test set, and the model with the highest accuracy is selected as the optimized XGBoost model. When evaluating the accuracy, the test set labels are compared with the model prediction results, and the accuracy and Kappa coefficient are used to evaluate the classification accuracy in various scenarios. The overall accuracy represents the ratio of correctly classified pixels to the total number of validation pixels. Compared with the accuracy, the calculation of Kappa is based on the confusion matrix, which compares the model's results with the labels, taking into account the randomness and consistency between the predictions and labels, which also gives it a unique advantage in evaluating classification results. The calculation formula of Kappa is as follows:

[0080]

[0081] In the formula, represents the proportion of correct classification; is the probability of a random prediction, which is the sum of the products of the actual and predicted numbers for all classes, divided by the square of the number of samples.

[0082] The final optimized XGBoost hyperparameter settings are shown in Table 1

[0083] Table 1 XGBoost hyperparameters

[0084]

[0085] In order to evaluate the information integrity and stability of the input data set produced by the above method in this embodiment, this embodiment uses the crop planting type data of Yunshan Farm and Longzhen Farm in Heilongjiang Province in 2019 for model training, and the data of 10 farms in Heilongjiang Province from 2019 to 2021 for model testing, and the data shown in Table 2 are obtained.

[0086] Table 2 Classification prediction accuracy of each farm from 2019 to 2021

[0087]

[0088] It can be seen from Table 2 that the model is used to predict the remote sensing data from 2019 to 2021. The prediction accuracy of the model is high and meets the accuracy requirements, which shows that the method of preparing a data set using the method of this embodiment can significantly improve the prediction accuracy of the model.

[0089] Step 4: Obtain the remote sensing image of the crop to be predicted, calculate multiple vegetation indices based on the remote sensing image, and input the multiple vegetation indices into the XGBoost model optimized in step 3 to obtain the planting type and growth cycle of the crop.

[0090] Another object of the present invention is to provide a system for implementing the above-mentioned method for identifying the crop growth cycle and type, comprising:

[0091] The data acquisition module is used to obtain long-term series of remote sensing images of crops in the study area and pre-process the remote sensing images;

[0092] The vegetation index calculation module is used to calculate the remote sensing images preprocessed by the data acquisition module to obtain multiple vegetation indices;

[0093] The data set processing module is used to filter the vegetation index calculated by the vegetation index calculation module after sorting them by time, extract the values ​​of each vegetation index of each pixel in multiple growth cycles from the filtered time series, form the extracted values ​​of each vegetation index in each growth cycle into a new band and superimpose them to form a feature grid, then perform random point sampling on the feature grid to obtain a data set, and divide the data set into a training set and a test set;

[0094] The model building and training module is used to build an XGBoost model for predicting crop growth cycles and types. The XGBoost model is trained using a training set, and then the trained XGBoost model is tested using a test set. The XGBoost model with the accuracy that meets the requirements is selected as the optimized XGBoost model.

[0095] The prediction module is used to obtain the remote sensing image of the crop to be predicted, and input it into the vegetation index calculation module to calculate multiple vegetation indices, and then input the multiple vegetation indices into the optimized XGBoost model to predict the planting type and growth cycle of the crop.

[0096] The above are only preferred embodiments of the present invention, and are not intended to limit the implementation methods and protection scope of the present invention. Those skilled in the art should be aware that all solutions obtained by equivalent substitutions and obvious changes made using the contents of the specification of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for identifying crop growth cycle and type, characterized in that: The steps include: Step 1: Obtain a long-term series of remote sensing images of crops in the study area and pre-process the remote sensing images; Step 2: Calculate the preprocessed remote sensing image to obtain multiple vegetation indices, sort the calculated vegetation indices by time and filter them, extract the values ​​of each vegetation index of each pixel in multiple growth cycles from the filtered time series, form the extracted values ​​of each vegetation index of each growth cycle into a new band and superimpose them to form a feature grid, then perform random point sampling on the feature grid to obtain a data set, and divide the data set into a training set and a test set; Step 3: Build an XGBoost model for predicting crop growth cycle and type, train the XGBoost model with the training set, test the trained XGBoost model with the test set, and select the XGBoost model with the required accuracy as the optimized XGBoost model; Step 4: Obtain the remote sensing image of the crop to be predicted, calculate its multiple vegetation indices, and then input the multiple vegetation indices into the XGBoost model optimized in step 3 to predict the planting type and growth cycle of the crop; The method for filtering the vegetation index in step 2 includes: The improved SG filtering method is used to smooth the vegetation index, and the formula of the filtering function is as follows: g i = ； In the formula, f i Represents the original data value in the time series; g i is the smoothed value, which is c n and f i The linear combination of n is the width of the moving window for filtering, nL and nR correspond to the left and right edges of the signal component respectively; c n is a quadratic polynomial fitting function, then by solving the least squares fit, ranging from nL to nR, we get c n , quadratic polynomial fitting function c n for: ； Among them, t corresponds to the day of the current year in the NDVI time series; Before filtering the vegetation index, the Kriging interpolation method is used to interpolate the vegetation index. The specific formula is: ； Where: Z(x0, y0) is the weighted sum of the values ​​of n sampling points; Z(x i ,y i ) is the sampling point (x i ,y i ) value; i The weight coefficient is determined by the variation function, i=1,2,…,n.

2. The method for identifying crop growth cycle and type according to claim 1, characterized in that: The preprocessing in step 1 includes spatial correction and atmospheric correction of remote sensing images.

3. The method for identifying crop growth cycle and type according to claim 1, characterized in that: The multiple vegetation indices in step 2 include LSWI, GLI, RNDVI, MSAVI, IRECI and SAVI.

4. The method for identifying crop growth cycle and type according to claim 3, characterized in that: The calculation formulas for multiple vegetation indices are as follows: ； ； ； ； ； ； In the formula, NIR represents the near infrared band; SWIR represents the short infrared band; G represents the green band; R represents the red band, B represents the blue band; EDGE1 represents the red edge 1 band; EDGE2 represents the red edge 2 band, EDGE3 represents the red edge 3 band; L is the adjustment parameter.

5. The method for identifying crop growth cycle and type according to claim 1, characterized in that: In step 3, the test set labels are compared with the model prediction results, and the accuracy and Kappa coefficient are used to evaluate the classification accuracy of the XGBoost model in various scenarios.

6. The method for identifying crop growth cycle and type according to claim 5, characterized in that: The calculation formula of Kappa is as follows: ； In the formula, represents the proportion of correct classification; is the probability of a random prediction, which is the sum of the products of the actual and predicted numbers for all classes, divided by the square of the number of samples.

7. A system for implementing the method for identifying crop growth cycle and type according to any one of claims 1 to 6, characterized in that: include: The data acquisition module is used to obtain long-term series of remote sensing images of crops in the study area and pre-process the remote sensing images; The vegetation index calculation module is used to calculate the remote sensing images preprocessed by the data acquisition module to obtain multiple vegetation indices; The data set processing module is used to filter the vegetation index calculated by the vegetation index calculation module after sorting them by time, extract the values ​​of each vegetation index of each pixel in multiple growth cycles from the filtered time series, form the extracted values ​​of each vegetation index in each growth cycle into a new band and superimpose them to form a feature grid, then perform random point sampling on the feature grid to obtain a data set, and divide the data set into a training set and a test set; The model building and training module is used to build an XGBoost model for predicting crop growth cycles and types. The XGBoost model is trained using a training set, and then the trained XGBoost model is tested using a test set. The XGBoost model with the accuracy that meets the requirements is selected as the optimized XGBoost model. The prediction module is used to obtain the remote sensing image of the crop to be predicted, and input it into the vegetation index calculation module to calculate multiple vegetation indices, and then input the multiple vegetation indices into the optimized XGBoost model to predict the planting type and growth cycle of the crop.

Images