Multi-Season Phenology Mapping of Nile Delta Croplands Using Time Series of Sentinel-2 and Landsat 8 Green LAI

Abstract

Space-based cropland phenology monitoring substantially assists agricultural managing practices and plays an important role in crop yield predictions. Multitemporal satellite observations allow analyzing vegetation seasonal dynamics over large areas by using vegetation indices or by deriving biophysical variables. The Nile Delta represents about half of all agricultural lands of Egypt. In this region, intensifying farming systems are predominant and multi-cropping rotations schemes are increasing, requiring a high temporal and spatial resolution monitoring for capturing successive crop growth cycles. This study presents a workflow for cropland phenology characterization and mapping based on time series of green Leaf Area Index (LAI) generated from NASA’s Harmonized Landsat 8 (L8) and Sentinel-2 (S2) surface reflectance dataset from 2016 to 2019. LAI time series were processed for each satellite dataset, which were used separately and combined to identify seasonal dynamics for a selection of crop types (wheat, clover, maize and rice). For the combination of L8 with S2 LAI products, we proposed two time series smoothing and fitting methods: (1) the Savitzky–Golay (SG) filter and (2) the Gaussian Processes Regression (GPR) fitting function. Single-sensor and L8-S2 combined LAI time series were used for the calculation of key crop Land Surface Phenology (LSP) metrics (start of season, end of season, length of season), whereby the detection of cropland growing seasons was based on two established threshold methods, i.e., a seasonal or a relative amplitude value. Overall, the developed phenology extraction scheme enabled identifying up to two successive crop cycles within a year, with a superior performance observed for the seasonal than for the relative threshold method, in terms of consistency and cropland season detection capability. Differences between the time series collections were analyzed by comparing the phenology metrics per crop type and year. Results suggest that L8-S2 combined LAI data streams with GPR led to a more precise detection of the start and end of growing seasons for most crop types, reaching an overall detection of 74% over the total planted crops versus 69% with S2 and 63% with L8 alone. Finally, the phenology mapping allowed us to evaluate the spatial and temporal evolution of the croplands over the agroecosystem in the Nile Delta.

1. Introduction

Satellite Earth observation data are increasingly used for analyzing dynamic land surface processes at different spatial and temporal scales [1,2,3]. Their processing into easy interpretable products become a valuable tool in relevant contexts, such as food security, assisting agricultural practices or studying the vegetation response to environmental changes [4,5,6]. Current agriculture challenges are related to the growing global population aiming to yield increase and sustainable intensification of cropland production [7,8,9], which requires regular crop growth monitoring in agroecosystems. Based on satellite data, changes in vegetation dynamics and crop growth evolution can be studied through the estimation of Land Surface Phenology (LSP) metrics from vegetation indices or biophysical variables time series [10]. Among the LSP metrics relevant for croplands, the most important ones are start-of-season (SOS) and end-of-season (EOS), marking the green-up and senescing of the crop growth cycle [11]. LSP products have been traditionally developed from coarse to moderate spatial resolution satellites, such as Moderate Resolution Imaging Spectroradiometer (MODIS), Advanced Very High Resolution Radiometer (AVHRR) and Visible Infrared Imaging Radiometer Suite (VIIRS) [12,13,14,15,16,17]. Thanks to the frequent image acquisitions of these missions, derived LSP products allowed regional to global long-term phenology detection of vegetation across diverse land cover types, but with low spatial detail and with the occurrence of mixed pixels over heterogeneous landscapes. For cropland growth cycle monitoring, however, both fine spatial and temporal resolutions are required to enable characterizing individual crop fields at their different phenological stages [18,19]. Available imagery from Landsat 8 (L8) at 30 m, with a single orbit and a 16-day repeat cycle, favoured time series and data fusion studies for identifying phenology patterns at the field scale [20,21,22,23,24]. On the other hand, ESA’s Sentinel-2 (S2) constellation (S2A/S2B) provides imagery at an even shorter revisit frequency, up to 5 days, and with similar spatial and spectral characteristics as L8. The combination of both satellites, therefore, provides optimal systematic land imaging coverage adequate for timely cropland mapping applications [25]. NASA’s initiative to create the Harmonized Landsat and Sentinel-2 (HLS) dataset brings forth freely available time series of spectral and radiometric consistent surface reflectance imagery with a uniform gridding from both satellites [26]. The HLS dataset has already been used in vegetation monitoring studies, showing a good performance and sufficient capability to capture seasonal dynamics of vegetation phenology, with high agreement compared to coarse-resolution satellites estimations, but at a finer spatial resolution of 30 m [21,27]. The HLS dataset was also positively evaluated in its feasible shorter revisit time to monitor grasslands [28], and proved to be able to achieve high overall accuracies of cropland mapping over different regions across the world [29]. Hence, the HLS dataset emerged as a potentially ideal data source for cropland phenology monitoring.

When it comes to LSP metrics estimation, widely used products derived from satellite data streams include spectral bands combination of vegetation indices (VIs) and biophysical variables obtained from either physical or empirical retrieval models [10]. Similarly to the broadly used VIs, i.e., normalized difference vegetation index (NDVI), enhanced vegetation index (EVI) and two-band EVI (EVI2) [11], which used alone can present a limited sensitivity to changes of vegetation photosynthesis dynamics [30] and high biomass growth stages [31,32], biophysical variables such as Leaf Area Index (LAI) can reflect the seasonal evolution of photosynthetic active vegetation in addition to serve as indicator of canopy structure condition [33]. Moreover, while established VIs represent an arithmetic formulation based on a few bands to maximize the spectral response to vegetative components, an LAI multiband retrieval model can exploit the full spectral information available from the satellite [34,35,36]. LAI is an essential input for land surface and atmospheric energy exchange modeling [37] and has been used for global scale phenology monitoring as well as to map cropping systems [38,39]. Despite the lower proportion of LAI-based phenology studies published, it was previously shown as a biophysical variable capable to improve the detection performance over croplands with respect to the EVI [40].

The timing of cropland phenology seasonality and growing seasons length are typically estimated as LSP metrics corresponding to key dates of the crop growth cycle by time series processing and fitting functions approaches [41,42,43]. As pre-processing for later extraction of LSP metrics, data smoothing methods are usually applied in order to reduce time series disturbances and possible artifacts [44,45]. Established methods include double logistic and asymmetric Gaussian functions and the moving Savitzky–Golay filter for noise mitigation of remote sensing time series [10]. More recently, machine learning fitting techniques emerged as attractive alternatives for bringing the possibility of further perform time series gap-filling with promising results [46]. In particular, the Gaussian Processes Regression (GPR) algorithm proved to be a competitive gap-filling technique, being able to additionally provide an approximation of the estimates uncertainty [47,48].

Compared to other vegetation canopies, time series analysis over croplands can be particularly challenging due to highly dynamic seasonal growth cycles, management practices and the complex spatial patterns in areas with high cropping frequency, with more than a single growing season per year [49]. LSP metrics calculation allow to map the cropland extent and characterize seasonality of different crop types. The simplest way to estimate cropping frequency is by counting the number of growing seasons as the number of peaks corresponding to a complete growth cycle in the time series [12,50]; while a single cropping time series presents one growth cycle per year, the identification of two to three crop cycles within a single year is more difficult. Moreover, multi-cropping patterns can span over consecutive years with different crop types rotations [22,51]. Standard phenology detection methods essentially extract transition dates representative of crop growth stages from the time series curve, either based on identifying curvature changes or defining a magnitude thresholds [44,46]. Particularly, threshold-based methods are generally adapted by the user accounting for the variability of cropland seasonal dynamics of the study site by setting temporal constraints and maximum magnitude thresholds to dismiss noisy peaks or minimize the non-vegetated spectral influence [10].

Altogether, the aim of this study is to demonstrate the added value of using harmonized S2 and L8 data streams, i.e., HLS dataset, with respect to single-sensor based approaches for characterizing the phenology of multi-cropping systems. To this end, we present a workflow using green LAI time series derived from the HLS data stream. The following objectives were defined: (1) to adapt an S2 GPR retrieval model to L8 imagery in order to obtain consistent estimations from both sensors data; (2) to derive, respectively, distinct green LAI time series and subsequently combine them; (3) to extract phenology indicators of different crop types and assess the consistency of the crop seasonal dynamics retrieved from different time series datasets during the same time period; and (4) to generate final crop phenology descriptor maps. The remainder of this paper is structured as follows. Section 2 describes the generation process of green LAI time series used for LSP metrics estimation. Section 3 presents the results of crop phenology characterization and mapping, while major findings and limitations of this work are discussed in Section 4. Finally, conclusions and future research lines are provided in Section 5.

2. Materials and Methods

With the ambition to characterize multi-season cropland phenology with HLS data [26], we propose the workflow sketched in Figure 1. First, HLS tile time series were separately processed for green LAI estimation applying a retrieval model originally developed for S2 [52], and here adapted to the spectral configuration of L8. As such, we obtained two single-sensor green LAI time series, which were separately processed with a Savitzky–Golay (SG) [53] smoothing filter based on a local average moving window. Second, both single-sensor LAI data streams were combined using two time series processing methods: (1) SG and (2) Gaussian Process Regression (GPR) [54], a machine learning fitting technique able to perform temporal gap-filling interpolation [48]. Next, the four green LAI data collections generated were analyzed on their ability to identify multi-season growth cycles by extracting a set of key LSP metrics (start of season, end of season, length of season) per crop type at parcel level. To do so, two different threshold strategies were applied defining either a seasonal or a relative amplitude value [46,55]. The performances of both methods were assessed by quantifying the number of growth seasons detected per crop type during the period of time studied, and then compared against crop calendar data. Finally, the phenology extraction process was applied to the combined GPR LAI time series for mapping the annual evolution of cropland phenology and cropping frequency over the study site. The complete data processing chain and analysis are detailed in the following sections.

2.1. Study Area

The Nile Delta accounts for nearly half of all agricultural lands of Egypt. This region is characterized by an arid climate, with relatively moderate temperatures and low precipitation occurring mostly during winter. Annual rainfall ranges from 100 to 150 mm and annual mean temperature is around 21.5 °C. A variety of crop types are cultivated by irrigation with water supply coming from the Nile. Major summer crops are rice and maize. Winter wheat is the third most important cash crop in terms of production area, while the Egyptian clover is a major winter forage crop. These crops are alternatively cultivated, allowing to identify multiple crop cycles within a single year. Extended perennials annual crops include grape and citrus species. The absence of remarkable variations in climate conditions allow to establish a reference crop calendar, with regular plant and harvest dates (

Table 1. Established crop calendar over the study site (Tanta, El Gharbia governorate, Egypt).

Crop	Planting Date	Harvest Date
rice	15 May	10 October
maize	1 May	10 August
wheat	10 November	1 May
clover	1 October	1 May

). It should be noted that these are the dates usually employed for crop management over all years, and here provide a reference framework to enable evaluation of the identified crop phenology evolution.

Information on crop type planting was collected from 2016 to 2019 within the study site (Figure 2), which is an agricultural area covering approximately 326 ha in the region of Tanta, within El Gharbia Governorate. From 162 parcels available in total, 55 are under crop rotation schemes with more than 1 crop type per year, planted and harvested twice annually. The remaining 107 parcels correspond to single-cropping systems with the same annual crop type. We focus on major cash crops, namely, rice, maize, wheat and clover, excluding fruit and tree crops. Parcels’ delimitation into a geospatial vector file for carrying out the later processing per parcel crop type was provided by the University of Tanta (Egypt).

2.2. Harmonized Landsat and Sentinel-2 (HLS)

NASA’s Harmonized Landsat and Sentinel-2 (HLS) surface reflectance imagery data streams were used to conduct this study. The HLS dataset is routinely produced, over a limited number of regions across the world, from observations of the Operational Land Imager (OLI) aboard Landsat 8 (L8) and the Multi-spectral Instrument (MSI) aboard Sentinel-2A/B (S2), providing seamless products from both sensors, based on applying a common atmospheric correction algorithm [56] and performing radiometric and geometric adjustments [26]. HLS provides separately BRDF-normalized surface reflectance, respectively, derived from MSI and OLI, using a fixed solar angle and a common nadir view. For the bands common to both sensors, the MSI spectral bands are adjusted to OLI’s spectral bandpasses, which are used as reference. The HLS projection and gridding correspond to the tilling system used by Sentinel-2, which is UTM/WGS84 projection [26]. We used the current version of the HLS dataset (v1.4) corresponding to tile 36RUV (T36RUV) at 30 m resolution, from the beginning of 2016 to the end of 2020. Only images with less than 50% of cloud coverage were selected to further processing. This led to a collection of 221 S2 images and 173 L8 images.

2.3. Green LAI Retrieval

The HLS L8 and S2 datasets were separately processed into independent data streams of green LAI based on an earlier S2 green LAI model that was developed within the framework of the H2020 project SENSAGRI (Sentinels Synergy for Agriculture (http://sensagri.eu/, accessed on 24 January 2022) [52]. Gaussian Processes Regression (GPR) was chosen as retrieval modeling technique as a competitive machine learning regression algorithm in vegetation properties mapping applications in terms of prediction accuracy [36,57], and noteworthy for providing uncertainty estimates of the predictions [58], as well as insight on the bands relevance in the trained model [59]. In brief, a GPR model was trained from data pairs of green LAI ground-based measurements and their corresponding S2 reflectance spectra, which was either synthesized from hyperspectral data or extracted from simultaneous S2 imagery, covering the S2 spectral bands of 10 and 20 m spatial resolution, i.e., 10 bands in total [60]. The training dataset comprises a variety of crop types over diverse agricultural sites across Europe. Additionally, training samples of forested and non-vegetation areas were added to improve the performance of the retrieval model and to make the model widely applicable. A total of 218 data samples were used to train the green LAI retrieval model, which was validated based on an independent dataset (R${}^2$ = 0.7, RMSE = 0.67 m${}^2$/m${}^2$). Detailed information about data sampling and validation procedures can be found in [52]. Given the spectral configuration compatibility between S2 and L8 (

Table 2. S2 and L8 spectral bands.

	Visible & NIR								SWIR
S2 band	2	3	4	5	6	7	8	8A	11	12
Wavelength (nm)	490	560	665	705	740	783	842	865	1610	2190
Spatial resolution (m)	10	10	10	20	20	20	10	20	20	20
L8 band	2	3	4	-	-	-	-	5	6	7
Wavelength (nm)	480	560	655	-	-	-	-	865	1610	2200
Spatial resolution (m)	30	30	30	-	-	-	-	30	30	30

), a new model adapted to L8 was developed. The spectral data of the described training dataset were resampled according to the band settings of L8; 6 bands ranging from the visible to the shortwave infrared (SWIR). Subsequently, these data were used to train a new GPR model. A general introduction into GPR can be found in [61], and a summary of the generic formulation is given in the next section.

2.4. Gaussian Processes Regression

In general, the GPR technique models the relation between input samples $\mathbf{x}\in {\mathbb{R}}^D$ and output observations $y\in \mathbb{R}$ as $y=f\left(\mathbf{x}\right)+ϵ$, where $ϵ$ is an additive Gaussian noise with zero mean and variance ${\sigma }_n^2$, and $f\left(\mathbf{x}\right)$ is a Gaussian-distributed random vector with zero-mean and covariance matrix $\mathbf{K}(\mathbf{x}, \mathbf{x})$, i.e., $f\left(\mathbf{x}\right)\sim \mathcal{N}(\mathbf{0},\mathbf{K})$. The role of the Covariance matrix is to encode the similarity between each combination of the input samples ${\mathbf{x}}_i$ and ${\mathbf{x}}_j$ using a kernel function $k({\mathbf{x}}_i, {\mathbf{x}}_j)$. The covariance design is of paramount relevance, as it must take into account the main properties of the variable to be modeled. Concerning the retrieval of vegetation properties such as LAI from Earth observation data, the asymmetric Square Exponential (SE) kernel is usually preferred due to its capability to (1) successfully approximate smoothly varying functions and (2) take into account asymmetries in the feature space [54].

The asymmetric SE defines the covariance kernel function as:

$$k({\mathbf{x}}_i,{\mathbf{x}}_j)={\sigma }_s^2\exp \left(-\frac{1}{2}\sum _{b=1}^D{\left[{\displaystyle \frac{{\mathbf{x}}_i\left(b\right)-{\mathbf{x}}_j\left(b\right)}{{\sigma }_b}}\right]}^2\right),$$

(1)

where ${\sigma }_s^2>0$ represents the output variance and ${\sigma }_b$ is related to the spread of the training information along the input dimension b in a way that the inverse of ${\sigma }_b$ describes the relevance of band b in the prediction process: the higher ${\sigma }_b$ the lower informative content of b. The covariance matrix is completely defined once the kernel’s free parameters and the noise variance ${\sigma }_n^2$ are set. These terms, usually referred to as GPR model’s hyperparameters, can be collectively denoted as $\mathit{\theta }=\{{\sigma }_s^2,{\mathit{\sigma }}^2,{\sigma }_n^2\}$, where $\mathit{\sigma }=[{\sigma }_1,\cdots ,{\sigma }_D]$.

The Bayesian framework of GPR allows estimating the distribution of f at any test point $x_{\ast }$ conditioned on the information carried by the training data. According to its formulation, $f\left(x_{\ast }\right)$ is normally distributed with mean and variance given by:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill f\left(x_{\ast }\right)=k_{\ast }^T{(\mathit{K}+{\sigma }_n^2{\mathit{I}}_N)}^{-1}y\\ \hfill {\sigma }_f^2\left(x_{\ast }\right)=c_{\ast }-k_{\ast }^T{(\mathit{K}+{\sigma }_n^2{\mathit{I}}_N)}^{-1}{k}_{\ast }\end{array}\end{array}$$

(2)

where N is the number of training samples, $k_{\ast }={[k(x_{\ast },x_1),\dots ,k(x_{\ast },x_N)]}^T$ is an $N\times 1$ vector containing the similarity between $x_{\ast }$ and the training input information, $y={[y_1,\cdots ,y_N]}^T$ is the training output and $c_{\ast }=k(x_{\ast },x_{\ast })+{\sigma }_n^2$.

The probability of the observations given the model’s hyperparameters $p\left(y\right|x,\mathit{\theta })$ is provided by the marginal likelihood over the function values f [61], whose logarithmic expression is:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \log p\left(y\right|x,f)=-\frac{1}{2}{y}^T{\left(\mathit{K}+{\sigma }_n^2{\mathit{I}}_N\right)}^{-1}y-\frac{1}{2}\log |\mathit{K}+{\sigma }_n^2{\mathit{I}}_N|-\frac{n}{2}\log 2\pi \end{array}\end{array}$$

(3)

Equation (3) is made up of three terms: the first one is essentially a data-fit term, the second one represents a complexity penalty and the last one is just a normalizing constant. The maximization of the marginal likelihood, i.e., the minimization of Equation (3), provides directly the optimum value of $\mathit{\theta }$. This optimization procedure is usually referred to as training the GPR [61,62]. Once $\mathit{\theta }$ values have been estimated, the prediction of y for a new input vector $x_{\ast }$ is given along with its uncertainty by Equation (2).

Finally, the formulation for the time series domain is straightforward. By imposing the dimension of input data be unitary, i.e., $D=1$, and substituting $\mathbf{x}$ for t, the hyperparameters of the GPR model to be employed to gapfill the collection of a generic surface property $P_S$ such as LAI become ${\mathit{\theta }}_{\mathit{t}}=\{{\sigma }_{st}^2,{{\mathit{\sigma }}_{\mathit{t}}}^2,{\sigma }_{nt}^2\}$, and the corresponding SE kernel for covariance estimation is given by:

$$k_t(t_i,t_j)={\sigma }_{st}^2\exp \left(-\frac{1}{2}{\left[{\displaystyle \frac{t_i-t_j}{{\sigma }_t}}\right]}^2\right),$$

(4)

where $t_i$ and $t_j$ denotes two generic acquisition dates. The estimation of LAI at any dates is still given by Equation (2), with the new covariance matrix ${\mathbf{K}}_{\mathbf{t}}$ being now calculated using the optimized ${\mathit{\theta }}_{\mathit{t}}$ among the elements of the input time series.

2.5. Smoothing and Gap-Filling

Smoothing filters allow to mitigate the noise effect to produce enhanced time series for later feature extraction [10,44]. In particular, the Savitzky–Golay (SG) filter minimizes the least-squares error in fitting a polynomial function to noisy data within a moving window centered at a point [53]. In this study, the smoothing window size (span) was set to 7 and the polynomial function to a quadratic. The resulting smoothed value ($g_i$) is obtained from a linear combination according to the following expression:

$$g_i=\frac{{\sum }_{n=-nL}^{nR}{c}_n{f}_{i+n}}{n}$$

(5)

where $f_i$ represents the original data point value, n is the width of the moving window and $nL$ and $nR$ are the left and right edges of the defined window frame. Then, each data value is replaced by a locally weighted average of nearby data points, and thus, the greater the value of the span is, the smoother is the fitted curve. SG is considered essentially a low-pass filtering method, which tends to preserve high frequency signal components, and responds less effectively in case of high noise levels [53], although it is commonly applied for data smoothing in crop phenology monitoring studies [39,63,64,65].

Optical remote sensing data streams are affected by cloud cover and unfavorable atmospheric conditions, reducing the data temporal continuity [66,67]. Compared to smoothing filters, gap-filling techniques (i.e., fitting functions) allow to reconstruct complete time series with missing data. Given the good precedent of the GPR algorithm as fitting function in reconstructing optically derived LAI time series due to cloud-cover [47,48], the GPR formulation was here adapted as fitting function for gap-filling purposes.

2.6. Generation of Green LAI Time Series Collections

Both S2 and L8 GPR LAI models were implemented in an HLS image processing chain for generating LAI time series, which were analyzed for different configurations: S30, L30 and S30+L30. The former two were derived from HLS S2 and L8 surface reflectance dataset, denoted as S30 and L30, respectively. For the latter, previous to the cropland phenology detection, the single-sensor LAI time series derived from each sensor were merged in a single time series by adding up both time series and applying independently and pixel-wise: (1) the local moving Savitzky–Golay smoothing filter and (2) GPR fitting for interpolation and gap-filling processing. These methods were run over all image acquisition dates, leading to continuous time series data streams, hereafter termed as SL30${}_{SG}$ and SL30${}_{GPR}$, respectively.

2.7. Crop Phenology Estimation

LAI time series smoothing and reconstruction are previous steps to facilitate calculation of phenology metrics [45,68]. The phenology extraction process was conducted for each green LAI time series configuration, over the five years of available time series, focusing on multi-crop parcels, whose major crop types are maize, rice, wheat and clover. Each of the LAI time series collections (S30, L30, SL30${}_{SG}$, SL30${}_{GPR}$) were first spatially averaged per parcel and subsequently analyzed to characterize their seasonal dynamics over time by calculating the following key LSP metrics corresponding to crucial stages of crops growth cycle [44]: (1) start-of-season (SOS) or green-up, (2) end-of-season (EOS) or dormancy, (3) length-of-season (LOS)—as the time span in days between SOS and EOS and (4) the seasonal integral area under the curve between SOS and EOS (Area).

The procedure for identifying the number of seasons from a time series is according to previous conventional approaches [69,70]. It makes use of consecutive local minimum, maximum and minimum of the curve to extract a single growing season, when the time series surpasses predefined thresholds. These specific constraints (i.e., thresholds) allow to reduce contributions from undesired artifacts at low frequencies. Consequently, the phenological cycles were evaluated using two criteria [46]: (1) peak prominence and (2) minimum separation thresholds. The former one corresponds to the minimum vertical distance from the season maximum to any side of the peak either reaching a local minimum or an endpoint, i.e., it measures how much the peak stands out due to its intrinsic height and its location relative to other peaks. The latter criteria determines the minimum separation in time between two consecutive potential peaks, so that only the largest maximum within the predefined threshold is selected. Here, the minimum prominence applied ranged from 10% to 30% the difference between the absolute maximum and minimum of a given parcel time series. The minimum separation was set to 90 days. Peaks not meeting the described criteria were discarded.

Regarding the determination of start and end of cropland growth cycles across time, two established strategies were evaluated on their seasonality detection capability, i.e., based on a seasonal or a relative amplitude threshold [44,46,71]:

• Seasonal. Each individual growing season extracted is analyzed to identify key phenological dates (SOS and EOS), when the upward and downward part of the LAI curve defining the growing season reaches a certain percentage fraction of the seasonal amplitude (difference between the maximum and the average of the two local minima per season), respectively.
• Relative. A fixed relative amplitude value is calculated as the difference between the mean maximum and mean minimum of the whole time series and set for all seasons detected, so that the start and end of season occur when the LAI curve reaches a certain percentage fraction of this relative amplitude.

After testing and comparing a range of threshold values on their seasonality detection capability, we defined an optimum threshold of 30% for both methods. Regardless of the seasonality detection approach applied, LOS is calculated as the difference between EOS and SOS. Finally, each season area refers to the integral encompassed between the LAI curve describing a season and the curve minima corresponding to SOS and EOS.

2.8. Analysis Setup

First, for each crop parcel and LAI time series collection, the key LSP metrics were extracted from valid growing seasons detected by independently applying the seasonal and relative methods. Then, the corresponding crop type was assigned based on the SOS and EOS dates retrieved throughout the studied years. LSP metrics were statistically analyzed per crop type and year. The availability of a crop planting record for each parcel allowed to assess the growing seasons detection accuracy for each crop type. Finally, spatio-temporal crop seasonality patterns were evaluated over the four years by mapping phenology descriptors over the study site. The smoothing, gap-filling and phenology extraction steps were carried out based on the time series modeling processor imported from the in-house developed scientific toolbox Decomposition and Analysis of Time Series software (DATimeS, v1.12) [46]. DATimeS incorporates conventional interpolation techniques and advanced machine learning fitting algorithms for time series analysis, as well as vegetation phenology modeling methods, with the advantage that unevenly spaced satellite time series can be processed.

3. Results

3.1. Green LAI Time Series

The HLS green LAI time series are well able to capture seasonal dynamics from 2016 to the end of 2020 over croplands practicing different cropping systems. See Figure 3 for an illustration. In this example, multi-cropping rotation parcels are characterized by more than one growth cycle per year, with generally two green-up and dormancy events, while in single crop type parcels, the growing season is annually completed. Both cases reflect a good agreement between the two single-sensor LAI estimations, which reveal a relatively high density of cloud-free observations, without large time gaps of missing data, throughout crops growing season. Reconstructed time series derived from the two single-sensor LAI time series provide a more smoothed and temporally continuous LAI profile. This effect is clearly visible over single cropping parcels, which initially present nosier LAI estimations, indicating high spatial heterogeneity and temporal variability.

3.2. Crop Phenology Characterization and Evaluation

The temporal evolution of the targeted phenology metrics was computed for the four LAI time series configurations and four crop types. Figure 4 shows the boxplot results of day-of-year (DOY) of the crop phenology metrics along with the planting and harvest dates from the reference crop calendar (Table 1). For later convenience and brevity, only the results obtained from the seasonal method are shown.

Rice and maize SOS were usually detected about 30 days after the planting date. Whereas rice EOS was estimated prior to the harvest date, maize EOS usually fell about 30 days after the harvest. For wheat, SOS was detected more closely to the planting date, and occasionally even before, while EOS was mostly estimated before the harvest date. Summer and winter crops can be clearly differentiated by the LOS values. For the former, the mean value is approximately 90 days, while for the latter it oscillates around 190 days. More than for LOS, a more significant dispersion is shown for Area, suggesting a higher variability in the estimation of this parameter. In general, all time series collections provided similar estimations of the key phenological dates over the four years. Inter-annual variations in phenology are reflected by all LAI time series collections, for instance, the sooner dormancy of rice and maize in 2018, matching the sooner green-up of wheat and clover also reported in 2018.

Table 3. Mean value and standard deviation (SD) in days of the key LSP metrics (SOS (DOY), EOS (DOY) and LOS (days)) per crop type and time series collection, comparing the two phenology detection methods (sesonal and relative). An overall standard deviation average in days was calculated from all crop types and each time series collection.

Time Series	Rice
	$\overline{\mathit{S}\mathit{O}\mathit{S}}\pm$SD		$\overline{\mathit{E}\mathit{O}\mathit{S}}\pm$SD		$\overline{\mathit{L}\mathit{O}\mathit{S}}\pm$SD
	Seasonal	Relative	Seasonal	Relative	Seasonal	Relative
S30	171 ± 16	173 ± 17	265 ± 13	263 ± 14	94 ± 17	90 ± 21
L30	165 ± 10	161 ± 21	264 ± 13	263 ± 14	99 ± 14	102 ± 24
SL30${}_{SG}$	166 ± 12	166 ± 19	264 ± 15	265 ± 13	98 ± 17	98 ± 22
SL30${}_{GPR}$	166 ± 16	170 ± 19	260 ± 13	260 ± 15	94 ± 18	90 ± 24
	Maize
	$\overline{\mathit{S}\mathit{O}\mathit{S}}\pm$SD		$\overline{\mathit{E}\mathit{O}\mathit{S}}\pm$SD		$\overline{\mathit{L}\mathit{O}\mathit{S}}\pm$SD
	Seasonal	Relative	Seasonal	Relative	Seasonal	Relative
S30	169 ± 17	169 ± 15	259 ± 15	256 ± 19	90 ± 17	87 ± 18
L30	165 ± 12	165 ± 12	254 ± 14	252 ± 16	89 ± 14	87 ± 18
SL30${}_{SG}$	166 ± 13	167 ± 13	258 ± 15	255 ± 18	92 ± 15	89 ± 17
SL30${}_{GPR}$	167 ± 11	168 ± 11	254 ± 14	253 ± 16	87 ± 10	85 ± 15
	Wheat
	$\overline{\mathit{S}\mathit{O}\mathit{S}}\pm$SD		$\overline{\mathit{E}\mathit{O}\mathit{S}}\pm$SD		$\overline{\mathit{L}\mathit{O}\mathit{S}}\pm$SD
	Seasonal	Relative	Seasonal	Relative	Seasonal	Relative
S30	294 ± 23	290 ± 30	110 ± 24	116 ± 29	181 ± 34	192 ± 45
L30	294 ± 32	290 ± 39	108 ± 18	109 ± 27	179 ± 40	184 ± 53
SL30${}_{SG}$	298 ± 17	295 ± 19	113 ± 21	114 ± 23	180 ± 33	184 ± 36
SL30${}_{GPR}$	303 ± 17	296 ± 28	106 ± 22	107 ± 24	168 ± 32	177 ± 40
	Clover
	$\overline{\mathit{S}\mathit{O}\mathit{S}}\pm$SD		$\overline{\mathit{E}\mathit{O}\mathit{S}}\pm$SD		$\overline{\mathit{L}\mathit{O}\mathit{S}}\pm$SD
	Seasonal	Relative	Seasonal	Relative	Seasonal	Relative
S30	302 ± 19	298 ± 26	110 ± 17	113 ± 21	173 ± 24	180 ± 34
L30	307 ± 18	301 ± 24	112 ± 25	117 ± 26	170 ± 33	181 ± 37
SL30${}_{SG}$	303 ± 18	300 ± 16	113 ± 21	118 ± 19	176 ± 30	184 ± 28
SL30${}_{GPR}$	310 ± 17	305 ± 18	110 ± 16	114 ± 16	165 ± 27	174 ± 25
	Total
	${\overline{\mathit{S}\mathit{D}}}_{\mathit{S}\mathit{O}\mathit{S}}$		${\overline{\mathit{S}\mathit{D}}}_{\mathit{E}\mathit{O}\mathit{S}}$		${\overline{\mathit{S}\mathit{D}}}_{\mathit{L}\mathit{O}\mathit{S}}$
	Seasonal	Relative	Seasonal	Relative	Seasonal	Relative
S30	19	22	17	21	23	29
L30	18	24	17	21	25	33
SL30${}_{SG}$	15	17	18	18	24	26
SL30${}_{GPR}$	15	19	16	18	22	26

summarizes the global mean and standard deviation of SOS, EOS and LOS for each time series and crop type, comparing the seasonal and relative threshold methods. For the seasonal method, standard deviations for SOS range from 10 to 17 days for summer crops and from 17 to 32 days for wheat and clover. Standard deviation values were analogously registered for the EOS parameter, with a larger variability around the phenological dates detected for winter crops. The range of variation of LOS was also comparable to the previous metrics for summer crops, but not for winter crops, for which the standard deviations were especially higher. Particularly for wheat, a maximum value of 40 days (L30) was obtained.

Results obtained from the relative threshold method show for summer crops similar SOS, EOS and LOS values as the ones from the seasonal method with differences between 1 and 4 days. In the case of winter crops, differences are between 1 and 7 days for SOS and EOS and from 4 to 11 days for LOS. The relative method provided the higher LOS associated standard deviation, up to 53 days (L30), for wheat. In general, for winter crops, the seasonal method estimates green-ups later and dormancies sooner, consequently leading to shorter crop seasons. Considering all crop types together, the seasonal method proved to be able to reduce the overall average standard deviation of SOS and EOS calculated for each of the time series collections. In particular, SL30${}_{GPR}$ with the seasonal method detected growing seasons across all crop types with a lower overall average standard deviation than the other time series collections, i.e., 15 days for SOS, 16 days for EOS and 22 days for LOS.

The number of growing seasons detected for each crop type and time series configuration is shown in

Table 4. Number (#) and accuracy detection (%) of crop seasons with respect of the number of crops planted per time series collection and crop type from 2016 to 2019 over all available parcels. The highest accuracy (%) per crop type and method columns is bolded.

Time Series	Rice				Maize				Wheat				Clover				Total
	Seasonal		Relative		Seasonal		Relative		Seasonal		Relative		Seasonal		Relative		Seasonal		Relative
	#	%	#	%	#	%	#	%	#	%	#	%	#	%	#	%	#	%	#	%
S30	36	60	33	55	124	71	105	60	60	67	52	58	73	72	66	65	293	69	256	60
L30	36	60	36	60	110	63	99	57	55	61	57	63	69	68	63	62	270	63	255	60
SL30${}_{SG}$	42	70	37	62	131	75	118	67	58	64	52	58	75	74	71	70	306	72	278	65
SL30${}_{GPR}$	42	70	39	65	131	75	115	66	63	70	56	62	79	77	70	69	315	74	280	66

. As an accuracy assessment, the detection capability was quantified for each crop type by comparing the total number of seasons detected with respect to the reference data of crop fields planted throughout the whole time period studied, expressed as a percentage. From both single-sensor-derived LAI time series, a similar number of phenology cycles were detected for rice and clover using the seasonal method, while for maize and wheat, S30 was able to capture a higher number than L30. The relative method led to a lower number of detected crop seasons for all time series collections and crop types, thus achieving lower detection rates. Particularly with this method, L30 identified more crop seasons over wheat and rice parcels than S30, while S30 identified more crop seasons over maize and clover. Overall, both seasonal and relative detection methods applied to multi-sensor combined time series captured the largest number of growing seasons for all crop types, achieving with GPR-gapfilled LAI time series the highest accuracy with a correct identification of 74% and 66%, respectively. The second best performance of the two methods were obtained with SG filtered time series, with an accuracy of 72% and 65% with the seasonal and relative method, respectively. The seasonal method and GPR are hereafter used for illustrating the subsequent phenology characterization analyzes and cropland mapping carried out.

Distinctions in crop seasonal patterns can be observed more precisely from the distribution of phenology metrics over the major crop types parcels shown in the histograms in Figure 5. This figure refers to the four years of SL30${}_{GPR}$ time series, which were grouped into 15 days bins for SOS, EOS and LOS, while differences of phenology timing between winter and summer crops is apparent, differences in green-up and dormancy dates were not as obvious within crops of the same season. Rice and maize SOS dates follow a close distribution; while rice reaches dormancy later, although over the same time range of distribution as maize. LOS number of days are slightly higher for rice than maize. Likewise, the Area values of both crops present a uniform overlap in the distribution, with relatively higher values for rice than maize. As for winter crops, wheat shows an earlier green-up and a slightly later dormancy than clover, which was already reported in Table 3. Particularly, wheat green-up tends to take place before clover by 15 days ahead. Wheat season length distribution is more irregular than clover, reaching higher values beyond 200 days. Differences in Area are not prominent between the two crop types.

A comparison of the LOS against the Area metric is provided in Figure 6. Both parameters show a strong linear correlation with a slope close to 1 (0.98), with higher data pair dispersion for winter crops, with a few LOS outliers values lower than 100 days. In this case, shorter crop seasons are due to an earlier retrieved EOS (Figure 5), while for summer crops, very high LOS values are a consequence of an earlier SOS estimation (Figure 5). All crop types considered, the Area parameter shows a relatively proportional increase regarding the season length. Particularly for summer crops, a remarkable Area variability occurs when approaching the maximum LOS value of 100 days when the highest Area values are also estimated.

3.3. Cropping Frequency and Phenology Mapping

The complete HLS LAI time series dataset (S2 and L8) were processed with the GPR interpolation technique with a 10-day regular spacing, as the temporal sampling used for mapping phenology metrics using the seasonal method over the study site. Cropping frequency was mapped as the number of crop growth cycles detected within the period of time studied, allowing to clearly observe the spatial distribution of the two main cropping patterns, for each year from 2016 to 2019 (Figure 7). Multi-cropping parcels being predominant, single-cropping parcels are mainly located in the central region of the study site. The extension covered by each cropping scheme varies locally over the years. Notably, a substantial proportion of pixels located in the center of the study region were not mapped in 2016, as no crop season were detected. Furthermore, when considering all years, the least crop seasons were detected within the multi-cropping area in 2018. This was partly due to the phenology detection method omission inferred from the time series, but mostly because of a delay on cropland seasonality, resulting in a growing winter season computed for the next year, in 2019.

Cropland phenology temporal evolution of the study site is presented next by distinguishing between summer (Figure 8) and winter crop growing cycles (Figure 9) for the multi-cropping pixels where two crop seasons are identified within a year. It should be noted that while SOS of both summer and winter season crops takes place during the same year, EOS dates corresponding to the latter ones are detected in the next year. Single-cropping pixels with only one growing season identified within a given year are mapped with the same value in the two figures. For better visualization, non-vegetated pixels are masked out applying a maximum value threshold: those pixels whose maximum LAI time series value does not exceed 1 are masked.

LSP metrics of summer season crops (Figure 8) were in general similarly retrieved across the four years. A sooner dormancy occurred in 2018 in accordance with a shorter season length estimation, as well as a notable decrease for Area, shown in Figure 4. Regarding the winter season, although the LSP metrics reveals a similar phenology spatial pattern across the four years, especially in 2016 and 2017, differences can be observed because of an early green-up detected in 2018, shown in Figure 4, consistent with the early dormancy of summer season crops within the same year. Compared to Figure 8, the Area metric exhibits a noisier pattern, though with a regular evolution in time, except in 2017 and 2018, where a general increase and decrease with respect to other years is, respectively, observed.

In both Figure 8 and Figure 9, some pixels were not mapped, since no crop growth cycle was detected within that season year, corresponding to those pixels with a lower number of seasons identified (Figure 7). This was particularly visible in the maps of 2016 and to a lesser extent of 2018. This is caused by intrinsic disturbances in the LAI time series and due to the fact that, occasionally, the time series does not contain the complete growing cycle (see Figure 3), which prevents the phenology method to compute a season.

4. Discussion

Having developed a workflow for LSP estimation and cropping patterns mapping from HLS data, the following steps were accomplished: (1) obtaining comparable LAI time series from HLS double-source data streams based on an S2 retrieval model; (2) contrasting the performance of two different threshold-based phenology detection methods; (3) assessing the LSP metrics estimation; and (4) analyzing the time series smoothing and gap-filling methods applied. These steps are further discussed below, followed by encountered limitations and future opportunities of the workflow.

4.1. Adaptation of S2 LAI Model to HLS

The availability of harmonized moderate spatial resolution (30 m) satellite products provides a high temporal resolution time series dataset advantageous for cropland monitoring purposes. As opposed to single-source satellite data streams, integration of S2 and L8 increases the temporal frequency of observations, usually leading to improved crop type classification and cropping intensity mapping [50,72]. Moreover, it was previously determined that data from both satellites can provide similar accuracies in LAI retrieval [73]. Furthermore, the HLS two-sensor data collection has shown a superior capability for LAI retrieval as opposed to single-sensor data streams separately in addition to not underestimating high LAI values [74]. The HLS dataset has been previously used for vegetation monitoring, revealing great improvements at detecting phenology at a finer scale from vegetation indices time series, commonly NDVI and EVI2 [21,27,65]. HLS is capable of capturing cropland phenology similarly to higher spatial resolution imagery such as PlanetScope 3 m spatial resolution [75], and is adequate for spatially scaling phenology field observations from GPR-based models [76]. Other multi-source based studies fused HLS data with other satellite data (e.g., Sentinel-1, GOES-16 ABI, VIIRS), allowing high-frequency, quality observations for an improved vegetation phenology estimation [77,78,79]. Here, four years of an HLS surface reflectance dataset was separately processed into a green LAI time series; while the GPR model was originally trained for estimating LAI from S2 data, the initial GPR model was successfully adapted to the spectral configuration of L8, and showed a compatible and consistent LAI evolution between the different harmonized satellites collections (S30 and L30). Therefore, despite the different original satellites spatial resolution and spectral configuration, our results demonstrated the benefits of the HLS data stream, likely thanks to a common radiometric adjustment and atmospheric correction [26,80].

4.2. Comparison of Two Threshold-Based Phenology Detection Methods

The phenology metrics extraction procedure, based on pre-defined amplitude threshold approaches as implemented in DATimeS, allowed to consistently capture cropland seasonal evolution over four consecutive years. Although threshold methods represent the simplest approach to extract phenology metrics from remote sensing time series [10], these methods are easy to tune according to the study location [81] and represent the most common method used to study vegetation dynamics [11]. This study compared two types of thresholds-based approaches, i.e., relative and seasonal amplitude; while the relative amplitude method makes use of a fixed amplitude threshold for the entire time series, the seasonal amplitude method, defined in relative terms, implies an adaptive calculation of phenology timing to individual growing seasons. This approach proved to be more flexible in determining vegetation dynamics over the study area, where typically different crop types are rotated in consecutive years following a periodic scheme. Therefore, each season’s amplitude value varies over time and do not always meet the threshold defined in absolute terms. The relative method has been widely applied for detecting green-up and dormancy events of single growing seasons within a year [16,20,41,82]. Yet, our findings suggest that this threshold is not suitable to be applied for longer term cropland seasonality detection over regions with different inter-annual crop species rotations. However, a quantitative comparison analysis with other well-established phenology detection methods may be necessary to further confirm these findings.

It must also be remarked that SOS and EOS dates, and consequently, the LOS and the Area metrics, can vary depending on the used threshold value, i.e., lowering the threshold generates earlier green-up dates and at the same time the estimation can be affected by noise and small variations in the time series [20]. Usual applied thresholds range from 10% to 50% [11], and may need to be adjusted for each crop type [16,20,71]. With the aim of establishing a standard procedure capable of phenology characterization over single and multi cropping parcels, we proposed a common threshold value of 30%, which allowed to capture as many growing seasons as possible over both types of cropping schemes. Correspondingly, other thresholds (prominence and minimum separation days) were defined as the optimum values in view of the results obtained trying to maximize the number of detected seasons. However, the presented workflow application to monitor larger regions might need another thresholds settings to account for different agroecosystems. Thus, a global parametrization of the threshold-based phenology detection methods should be further investigated.

4.3. Evaluation of Crop Phenology Detection

In the absence of in situ crop growth stages data, fixed crop calendar data were used as a reference to evaluate the detection of key phenological events in a general time frame and to additionally assess the suitability of the defined detection threshold. SOS dates were detected after planting and varied within 2–4 weeks, in accordance with previous studies [16,20,43,65], while EOS dates were generally retrieved closer to the harvesting, around 1–3 weeks ahead, except in the case of maize, for which dates were estimated later. Consistency in determining crop phenology was further proven, with a reliable correlation found between two key phenology metrics, such as LOS and Area. It opens opportunities for modeling new agricultural applications in future research, considering, for instance, that the Area parameter can be used to estimate crop yield [83,84]. Furthermore, phenology-based estimation of cropping frequency permitted to monitor the distribution of cropping patterns over the study site. When applied in an operational framework, the derived yearly maps of cropland dynamics can, among others, help assessing if specific crop species rotations affect soil quality [85,86] or the pressure over water resources [51]. Furthermore, characterizing the spatial distribution of different phenological patterns can be useful for distinguishing between vegetation types across heterogeneous landscapes [13,22] in addition to evaluate the impact of climatic changing conditions [50].

4.4. High Temporal Resolution for Phenology Detection Improvement

Combined multi-sensor time series were generated using smoothing and gap-filling techniques, obtaining a considerable noise reduction and data continuity, which eventually benefited the automated phenology detection. Each of the analyzed data streams was able to consistently monitor crop progress, with identified similar yearly patterns of crop growth for the same crop type. These patterns varied from year to year, likely due to weather conditions or changes in management practices. No substantial discrepancies in the mean estimates of the LSP metrics were observed between the different LAI time series configurations, mainly due to the initially available number of cloud-free observations from S2 and L8 single-sensor collections over the study site. As such, it permitted to automatically extract crop phenology seasons with a similar timing as the combined reconstructed LAI time series. The results also revealed that combining time series can reduce the dispersion around the estimation of the key crop phenology dates analyzed. More significant is the improvement in the detection of the number of growing seasons, providing S2 an overall accuracy of 69% in comparison to 63% obtained by L8. These results confirm the superior performance of S2 as compared to L8 as also demonstrated in previous studies (see review [49]). The combined time series data flow led to the highest accuracies, with GPR gap-filling being more effective than SG, i.e., 74% and 72%, respectively, in addition to a reduction in the variability of the number of days associated with the key phenological dates.

Given these promising results, the workflow should, however, be extended over other regions where cloud contamination represents a more limiting factor for crop monitoring, so that more substantial improvements derived from multi-sensor data combination can be expected. Furthermore, expected in near future, the recent launch of Landsat 9 (L9) will further contribute to reducing the time gap of data collection, since the combined L8 and L9 revisit time will be each 8 days [87]. In anticipation of a higher frequency of satellite acquisitions, preferences tend to move towards faster time series processing smoothing methods with flexibility to preserve intra-seasonal variations (e.g., Savitzky–Golay [53] and Whittaker [88,89]) [90].

4.5. Limitations and Future Opportunities

Cropland monitoring through phenology mapping requires a long-range time series data stream to proper capture crop growing stages, by means of continuous and smooth data. GPR represents an adequate technique for generating continuous time series that can be widely applied without using any ancillary data but with an important inherent disadvantage, i.e., the computational cost. In this respect, some efforts have been made to overcome the computational inconvenience by optimizing the calculation of the GPR hyperparameters per crop type, obtaining a considerable processing time reduction and hardly affecting the fitting accuracy retrieving phenology parameters per pixel, although it is necessary to have a land cover map at disposal [47]. Concerning this constraint, the advent of cloud computing platforms such as Google Earth Engine (GEE) allow to process satellite data for large-scale vegetation mapping with a high-performance computing capacity [91]. With ambition to expand the applicability of this study, the here developed workflow could be potentially implemented into GEE, where a similar GPR LAI retrieval model was recently successfully imported for green LAI mapping from S2 over the Iberian peninsula [48]. Furthermore, image fusion of different satellites (Landsat and MODIS) for time series smoothing and gap-filling has already been investigated [23,92] and a framework for cropland phenology monitoring in GEE using GPR has been recently devised [93]. These works suggest that our proposed workflow can be almost directly integrated into GEE, implying that the phenology detection can be realized anywhere in the world and for any time window. However, some additional efforts will be required to enable integrating the here presented flexible phenology detection method, i.e., with a seasonal-based amplitude threshold. GEE-based phenology studies currently rely on a double logistic model, whereby the SOS and EOS metrics are derived with a half-maximum criterion method, which is unable to distinguish between two or more seasonal cycles [93,94].

5. Conclusions

A workflow was developed for cropland phenology estimation based on a green LAI time series derived from NASA’s Harmonized Landsat 8 (L8) and Sentinel-2 (S2) surface reflectance dataset. It was applied over an agricultural region in the Nile Delta from 2016 to 2020. Starting from a green LAI retrieval model previously trained with GPR for S2, an adapted GPR model was trained for L8. Both models were implemented in an image-based retrieval chain, obtaining two different LAI time series data streams, which were used separately and combined to extract LSP metrics for characterizing crop rotation of major cash crops over the study area (maize, rice, wheat, clover). Two time series processing methods were proposed to combine both single-sensor LAI data streams: the (1) the Savitzky–Golay (SG) smoothing filter and (2) the GPR gap-filling fitting technique. Differences in detecting cropland phenology seasonality were compared for the four resulting LAI data streams (i.e., L8, S2 and combined SG and GPR). The evolution of cropland seasonal dynamics were observed by generating yearly phenology maps. Overall, the combined LAI time series data by means of GPR, followed by SG, allowed to detect correctly a higher number of within-year growing seasons per crop type within the same time period, as opposed to the single-sensor LAI time series. Summarizing, this study demonstrated that distinct multiband green LAI retrieval models applied to HLS data can retrieve compatible estimations that were successfully used to reconstruct enhanced time series. A dense temporal image data stream, with a higher cloud-free observations availability provided by L8 with respect to the individually use of S2, proved to be crucial for proper monitoring crop rotations of dynamic and heterogeneous agricultural lands. Although our study focused on the Nile Delta agroecosystem, the proposed workflow is neither site- nor crop-specific, and therefore, can be easily applied over other regions where both S2 and L8 images are available.