Daily Retrieval of NDVI and LAI at 3 m Resolution via the Fusion of CubeSat, Landsat, and MODIS Data

Abstract

Constellations of CubeSats are emerging as a novel observational resource with the potential to overcome the spatiotemporal constraints of conventional single-sensor satellite missions. With a constellation of more than 170 active CubeSats, Planet has realized daily global imaging in the RGB and near-infrared (NIR) at ~3 m resolution. While superior in terms of spatiotemporal resolution, the radiometric quality is not equivalent to that of larger conventional satellites. Variations in orbital configuration and sensor-specific spectral response functions represent an additional limitation. Here, we exploit a Cubesat Enabled Spatio-Temporal Enhancement Method (CESTEM) to optimize the utility and quality of very high-resolution CubeSat imaging. CESTEM represents a multipurpose data-driven scheme for radiometric normalization, phenology reconstruction, and spatiotemporal enhancement of biophysical properties via synergistic use of CubeSat, Landsat 8, and MODIS observations. Phenological reconstruction, based on original CubeSat Normalized Difference Vegetation Index (NDVI) data derived from top of atmosphere or surface reflectances, is shown to be susceptible to large uncertainties. In comparison, a CESTEM-corrected NDVI time series is able to clearly resolve several consecutive multicut alfalfa growing seasons over a six-month period, in addition to providing precise timing of key phenological transitions. CESTEM adopts a random forest machine-learning approach for producing Landsat-consistent leaf area index (LAI) at the CubeSat scale with a relative mean absolute difference on the order of 4–6%. The CubeSat-based LAI estimates highlight the spatial resolution advantage and capability to provide temporally consistent and time-critical insights into within-field vegetation dynamics, the rate of vegetation green-up, and the timing of harvesting events that are otherwise missed by 8- to 16-day Landsat imagery.

1. Introduction

Satellite remote sensing offers unparalleled opportunities for land surface monitoring and characterization over space and time domains, with far-reaching socioeconomic benefits [1,2,3]. Growing global food security concerns [4] necessitate space-borne observing systems that provide time-critical sensing of agricultural systems together with actionable intelligence for optimizing sustainable management practices, agricultural efficiencies, and crop production [5,6,7]. In the context of precision agriculture, timely and repeatable information on intrafield variability of crop conditions is key to achieving sustainable resource allocation and utilization [8,9]. Similarly, agricultural management decisions (e.g., timing of irrigation, fertilizer application, and harvesting) could benefit tremendously from continuous and near real-time observations of crop phenological development [10,11].

Space-borne observing capacity has evolved significantly in terms of spectral, spatial, and temporal resolution since the launch of the first Land Remote Sensing Satellite (Landsat) in 1972 [12,13]. However, the inevitable trade-off between resolvable spatial detail and frequency of observations [14] continues to affect the capacity to rapidly detect evolving vegetation dynamics at spatial scales fine enough for effective resource management. These spatiotemporal restrictions have been partly overcome via multisensor data fusion algorithms designed to synergistically exploit sensor-specific attributes, such as the high spatial resolution (i.e., 30 m) of Landsat and high revisit frequency (i.e., daily) of MODIS [15]. As an example, the spatial and temporal adaptive reflectance fusion model (STARFM) [16] has been used to predict surface reflectance [17,18], leaf area index [19], land surface temperature [20], and evapotranspiration [21,22] at both high spatial and temporal resolution. However, the lack of actual observations to inform upon intrinsic (i.e., at the scale of interest) spatial features at the time of prediction can introduce uncertainties that may affect the ability to meet often stringent product accuracy criteria [23]. Additionally, the conditional requirement of a sufficient number of homogeneous pixel samples at the coarse scale may not be met over heterogeneous landscapes associated with a relatively small characteristic field size [16].

The recent deployment of two identical Sentinel-2 satellites [24] represents an important advance towards resolving the classical spatiotemporal divide. The combination of high spatial resolution (<30 m) Landsat 8 and Sentinel-2 observations can facilitate a global median average revisit interval of about three days [25], which highlights the key importance of multisensor integration and interproduct harmonization efforts [26,27]. However, the achievable frequency is still suboptimal for detecting rapid vegetation dynamics and the exact timing of abrupt change events, particularly in light of obfuscating atmospheric effects (e.g., clouds, dust) which limit visibility and reduce the achievable repeat rate. While daily imaging at very high spatial resolution (<5 m) may be facilitated by commercial operators of programmable sensor constellations (e.g., WorldView, Pleiades), the satellite tasking is area-limited and prohibitively expensive [12].

One solution may lie in the cost-effective observation strategy offered by emerging constellations of small satellites, which is on a fast-track to revolutionize the earth observation potential [14]. This new operational paradigm has been facilitated by the standardized CubeSat concept [28] that uses compact and light-weight (<1.33 kg) single-unit (1U; 10 × 10 × 11.35 cm) building blocks as the modular basis for forming larger satellite configurations (i.e., 3U, 6U, 12U) for a variety of application fields [29]. As CubeSats are comprised largely of commercial off-the-shelf components and deployed as secondary payload on reusable rockets, the changing economics of earth observation have made it feasible to launch CubeSats in flocks at an unprecedented rate [30]. Planet [31] is leading the way in exploiting this new resource and its constellation of 3U CubeSats deployed in low-earth orbits has realized daily global imaging over terrestrial surfaces with an orthorectified pixel size of 3.125 m.

The number of explorative investigations into the potential benefits afforded by CubeSat data for land surface characterization is still limited and based on relatively short and temporally sparse records available at the time of study. The opportunities of CubeSat detectable intrafield-scale NDVI variability for precision agriculture were first demonstrated by Houborg and McCabe [32], while McCabe et al. [30] showcased their potential for providing valuable temporal insights into vegetation dynamics and evaporation. However, the relatively inexpensive and less advanced sensor designs do not come without limitations. Spectral bandwidth and response limitations, together with low radiometric quality relative to rigorously calibrated conventional satellites (e.g., Landsat and Sentinel), result in cross-sensor inconsistencies, which can severely affect the interpretability and utility of CubeSat-based time series records [33]. The Cubesat Enabled Spatio-Temporal Enhancement Method (CESTEM) was designed to optimize the benefit afforded by CubeSats via radiometric normalization against Landsat 8 observations [33]. Application of CESTEM to a limited time series of CubeSat imagery demonstrated accurate reproduction of Landsat-consistent surface reflectances at the scale and frequency of CubeSat observations, highlighting potentially far reaching implications via high quality, timely, and temporally consistent sensing at subfield scales [33].

In this study, we demonstrate an expanded application of CESTEM to a temporally dense record of CubeSat scenes acquired over an agricultural farming area in Saudi Arabia. Key advancements include an improved technique for smooth phenology reconstruction and an extension to map Landsat-based leaf area index (LAI) at the spatial and temporal resolution of the CubeSat data (i.e., CESTEM-LAI). The overall objective is to demonstrate the potential of the +170 CubeSat constellation for producing Landsat-consistent vegetation metrics at the highest spatiotemporal resolution ever retrieved from space-based sensors. The specific objectives of the study are: (1) to evaluate and intercompare the robustness and consistency of NDVI time series derived from “raw” and CESTEM-corrected CubeSat data over a highly dynamic multicut alfalfa field; and (2) to demonstrate and evaluate the capacity of the CESTEM framework for spatiotemporal enhancement of Landsat-based LAI.

2. Materials and Methods

This study used CubeSat (Section 2.1), Landsat 8 (Section 2.2), and MODIS (Section 2.3) data acquired over an irrigated farming area (~2400 ha) in Saudi Arabia (Figure 1) over a six-month period from 24 April to 25 October 2017. Alfalfa and corn were the predominant crops under pivot (~800 m in diameter) irrigation during this period. While the corn cultivated at the farm is characterized by two distinct growing seasons, with planting dates in the spring and the fall, alfalfa is grown continuously as a multicut crop with a 30–40 day harvesting frequency [34]. During this time, the average daytime air temperature high varied from ~45 °C during the summer months to ~35 °C in October. Sky conditions were predominantly free of clouds although prone to occasionally high dust aerosol loadings [35].

2.1. CubeSats

CubeSat data were acquired from the PlanetScope (PS) constellation that is commercially operated by Planet [31]. The constellation currently consists of >170 active 3U (i.e., 10 × 10 × 30 cm) CubeSats (“Doves”). At the time of writing, the PS constellation is capturing daily RGB and near-infrared (NIR) imaging of the earth’s terrestrial surfaces with a 3–4 m ground sampling distance (GSD). The vast majority of the CubeSats are deployed into sun-synchronous orbits (~475 km orbit altitude) with a midmorning equatorial overpass time (9:30–11:30 a.m., local solar time) and a nadir GSD between 3.5–4 m. Smaller flocks of PS satellites still operate in the International Space Station orbit (400 km altitude), which results in a variable equatorial overpass time and a GSD of ~3 m. Each CubeSat is associated with unique spectral response functions, and in addition to being significantly broader than the discrete bands of Landsat 8, the bands in the visible domains are characterized by substantial overlap [33].

Here, we used the Level 3A PlanetScope Ortho Tile Product that merges a consecutive series of individual scenes (typically four to five) from a single satellite into 25 × 25 km tiles referenced to a UTM map projection and with an orthorectified pixel size of 3.125 m [33,36]. The products have been geometrically corrected with a positional accuracy typically better than 10 m root-mean-square-error. Planet recently implemented a fairly rigorous radiometric correction approach based on a combination of prelaunch calibration, lunar calibration, and on-orbit cross-sensor calibration based on crossovers with rigorously calibrated sensor systems such as RapidEye and Landsat 8. The updated calibration approach has been verified against a limited number of vicarious collections (i.e., 12 samples covering a narrow reflectance range) over the Railroad Valley desert site in Nevada, which suggested a radiometric accuracy within ±10% of the top of atmosphere (TOA) reflectance reference and an average constellation uncertainty on the order of 5–6% [36,37]. This suggests a significant improvement over the first generations of Planet CubeSats [32,33].

The Planet data were converted to TOA reflectance using the calibration coefficients associated with each scene. In order to investigate the impacts of atmospheric effects, conversion of the at-sensor radiance to surface reflectance was undertaken using the Second Simulation of the Satellite Signal in the Solar Spectrum (6SV) atmospheric correction approach using sensor-specific spectral response functions, scene-specific view and illumination geometry, and satellite-based atmospheric state data (Section 2.2).

2.2. Landsat 8

Landsat 8 (L8) data were acquired over the study region with an 8–16-day revisit period using data from overlapping 165/43 and 164/43 path/row acquisitions. Whereas the 165/43 acquisitions provide full farm coverage, the 164/43 scene swath edge is positioned on the western half of the farm, with the exact location dependent on the day of the overpass. As a result, the study region is divided into sections with an 8-day, 8–16-day, and 16-day repeat coverage (Figure 1). L8 is positioned in a sun synchronous orbit (~705 km orbit altitude) with an equatorial overpass time between 9:56 a.m. and 10:26 a.m. (local solar time). The Operational Land Imager (OLI) surveys the earth in a 185-km swath (with a maximum view angle of ±7.5°) acquiring data in nine spectral bands distributed over the visible to near-infrared (VNIR) and shortwave-infrared domain. However, only the blue (band 2), green (band 3), red (band 4), and NIR (band 5) bands were considered here. The radiometric reflectance accuracy of the OLI is generally within the 3% sensor design specifications, supported by extensive vicarious calibration campaigns [23]. The Level-1 Terrain and Precision (L1TP) corrected dataset [38] was used in this study.

At-sensor radiance was converted into surface reflectance using an implementation of the 6SV [39,40], specifically adapted to handle the impacts of dust aerosols and adjacency effects over desert environments [23]. Inputs to 6SV include scene-specific sun-sensor geometry, total ozone from the Ozone Monitoring Instrument (OMI), atmospheric profiles of water and temperature from the Atmospheric Infrared Sounder (AIRS), and aerosol optical thickness at 550 nm (AOT₅₅₀) from MODIS deep blue algorithm retrievals [23]. Further details on the atmospheric correction approach can be found in Houborg and McCabe [23,41].

Landsat-based leaf area index (LAI) was mapped over the study domain using prediction models previously developed over the study site on the basis of a hybrid machine-learning approach [41]. The technique utilizes an implementation [42] of the Cubist rule-based regression tree approach [43,44,45] to define the relationship between a suite of atmospherically corrected L8 vegetation indices (VI) and the target property (i.e., LAI). In Cubist, each terminal node of the learned tree is associated with a multivariate linear regression model and defined by a set of spectral threshold conditions (i.e., rules). A multiday hybrid LAI target training dataset was used, consisting of in situ measurements coupled with physically based retrievals independently derived via an implementation of the regularized canopy reflectance model (REGFLEC) [46]. The addition of REGFLEC-based LAI targets significantly improved model predictability and portability, yielding a cross-validated relative mean absolute difference (rMAD) of ~11% relative to an rMAD of ~18% for the in situ only training case. The optimal LAI prediction model was formulated as two multivariate regression equations combining VIs derived from visible, near-infrared, and shortwave-infrared L8 bands (see Equations (1a) and (1b) in Houborg and McCabe [41]). Interestingly, models relying only on standard indices formed from red and NIR band reflectances (e.g., NDVI) were associated with a substantially higher retrieval uncertainty (i.e., rMAD ~ 21%).

2.3. MODIS

MODIS surface reflectance data were acquired on days coincident with CubeSat acquisitions. The version 6 combined (i.e., Terra and Aqua) MCD43A4 product [47] was used, which provides daily 500 m surface reflectances in seven bands normalized to a nadir view direction and local solar noon using a semiempirical bidirectional reflectance distribution function (BRDF) [48]. The data are provided in a sinusoidal projection grid on 1200 × 1200 km tiles. Only the four bands in the VNIR were retained for further analyses. MODIS has a wide field of view (±55°) and corrections for reflectance anisotropy are generally needed to produce temporally consistent time-series data. The BRDF utilizes the best observations from both Terra and Aqua sensors collected over a 16-day period centered on the day of interest. Observations at the day of interest are emphasized in the daily retrievals. The associated quality flags were used to isolate high-quality retrievals that have an accuracy of ~5% [49].

2.4. CESTEM

This section describes the framework of the CubeSat Enabled Spatio-Temporal Enhancement Method (CESTEM) [33], including an extension (CESTEM-LAI) that maps LAI to the spatial and temporal resolution of the CubeSat data (Figure 2). At its core, CESTEM serves as a radiometric normalization technique to translate broadband (RGB + NIR) PS data acquired from multiple CubeSats into L8-consistent surface reflectances (i.e., blue, green, red, and NIR). The approach relies on the availability of L8 reference data on the CubeSat prediction day (P_day) to learn band-specific translational associations on a scene-specific basis using the Cubist decision tree regression approach, as detailed in Houborg and McCabe [33]. L8 reference data are drawn from the available pool of Landsat images, using MODIS-based information to quantify spectral changes over given CubeSat–Landsat acquisition timespans (Section 2.4.2). CESTEM-LAI (Section 2.4.4) uses the normalized L8-consistent PS data together with L8-based reference LAI in a similar fashion to generate maps of LAI at the spatial and temporal resolution of the CubeSat data. All of these processes have been designed for complete automation. A detailed description of the generic methodology and processing steps schematized in Figure 2 is provided in subsequent sections.

2.4.1. Data Preprocessing

PlanetScope Ortho Tile Product (Section 2.1), L8 L1TP (Section 2.2), and MODIS MCD43A4 (Section 2.3) images are acquired programmatically using Planet’s API, Google Cloud Platform, and the data pool portion of the Land Processes Distributed Active Archive Center (LP DAAC), respectively. In the case of PS acquisitions, tiles are only downloaded if the percentage of unusable data (i.e., as specified by the “unusableDataPercentage” field in the associated metadata file) is <60%. Additionally, if multiple tiles are available over the study domain on any given day, only the tile with the lowest percentage of unusable data is selected for download. All PS tiles are subsetted to conform with the study scene domain (Figure 1).

In order to reduce cross-sensor coregistration errors, successive PS scenes are co-aligned by optimizing their spatial autocorrelation, as described in Houborg and McCabe [33]. The automated technique determines pixel shifts in x and y directions in reference to a base image by maximizing the goodness of fit (i.e., R²) between estimates of NDVI extracted from corresponding high-NDVI contrast zones (e.g., such as around field boundaries). The pixel-shifted scene is then used as the basis for coregistering the next scene in the successive time-series record. While a small acquisition gap is important for effectively coregistering scene-specific spatial features, a secondary (already aligned) base image is incorporated to avoid error propagation in the pixel shifting. This reference is continuously updated, lagging the primary base image with a few days or more, and only involves clear-sky scenes as reflected in a high R² value (see below). As erroneous pixel-shifting is often the result of scene contaminations (as reflected in a lower R² value), the secondary base image will take preference over the primary base image in those situations.

The technique is also used to simultaneously detect and remove cloud contaminated PS scenes based on thresholding the R² values established between successive scenes. Since a near-perfect agreement (i.e., in terms of R²) is expected between corresponding spatial features in two scenes less than a few days apart (i.e., given minimal change in cover conditions), any reduction in the R² value below the given threshold (R²_thres) is likely attributed to scene contaminations such as clouds or high aerosol content. The R²_thres for scenes <5 days apart was initially set to 0.95 and continuously updated using a weighted average of R² values observed for successive scenes. Any scene producing an R² below the specified threshold is assumed to be cloud contaminated and removed from further analyses. For successive scenes >5 days apart, R²_thres defaulted to 0.7.

The MODIS products are resized to the domain of the PS imagery and reprojected into UTM using the MODIS Conversion Toolkit (MCTK) in ENVI batch mode and pixel-shifted to optimize the correlation between scale-consistent MODIS and PS NDVI (i.e., after resampling to 500 m resolution). The same pixel shift is used for all MODIS scenes. L8 data is also resized to the domain of the PS imagery, reprojected into the UTM zone of the PS imagery, and pixel-shifted (if needed) as described above to properly align with the closest (typically day-coincident) PS scene data (after resampling to 30 m resolution). The UTM zone harmonization is necessary when the L8 data is acquired from two different paths (Section 2.2). The Landsat data is atmospherically corrected as described in Section 2.2 and clouds are detected using a Python implementation [50] of the Fmask algorithm [51].

2.4.2. VNIR Reference Sampling

L8 reference data (${\mathbf{L}\mathbf{8}}_{\mathrm{Bi}}^{\mathrm{REF}}$) for each spectral band (B_i) at 30 m resolution are needed on a given CubeSat acquisition day to train the machine-learning model (Section 2.4.3) and to correct PS TOA reflectance data (${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{TOA}}$) into L8-consistent surface reflectances (${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{L}8}$) (Figure 2). The reference data are sampled from up to 10 day-coincident Landsat–PS scene pairs. Currently, the search for available Landsat scenes is conducted both backward and forward in time from the given prediction day (P_day), but that can be changed to backward-only to facilitate near real-time predictions. The reference sampling approach does not require spectrally invariant pixel targets and is able to account for temporal changes in the spectral signatures based on multitemporal MODIS data. With a daily imaging capability and Landsat-like VNIR band configurations, MODIS is well suited for quantifying the relative spectral changes in Landsat-consistent bands at the MODIS resolution (i.e., 500 m) over the various timespans between P_day and each of the day-coincident Landsat–CubeSat acquisition days. Downscaled (i.e., from 500 m to 30 m) relative surface reflectance changes are subsequently used to correct the Landsat pixel references for moderate (as defined by a threshold percentage, C_th) reflectance changes occurring over the acquisition timespans, as described in detail below.

Assessments of the relative change in surface reflectance at the Landsat 30 m scale requires downscaling via ${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{TOA}}$. To this end, histogram matching is first applied to day-coincident PS and MODIS scenes after aggregation of the PS data to 500 m resolution. In this way, the pixel values in the PS image are manipulated to mirror the distribution of the MODIS BRDF-adjusted surface reflectances (${\mathbf{MOD}}_{\mathrm{Bi}}^{\mathrm{SR}}$). The Cubist decision tree regression algorithm is then used to learn the relationship between the “raw” (i.e., TOA reflectances) and corrected (i.e., MODIS-consistent surface reflectances) 500 m PS pixel values, using the implementation of Cubist described in Houborg and McCabe [33]. Disaggregation to the 30 m scale is facilitated by applying the derived band-specific regression models to ${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{TOA}}$resampled to 30 m. Finally, a conservative three-point temporal smoothing algorithm is applied to the resulting time-series images (${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{MOD}}$). This algorithm goes through each scene sequentially and flags pixels that have a lower NDVI value than corresponding pixel retrievals from two scenes acquired immediately before and after (i.e., within five days) the scene in question. Since anomalous NDVI behavior of this kind is unlikely over such a short period in time, the B_i pixel value is replaced by an average of the two pixel retrievals weighted by the difference of the acquisition timespans (Δday) (i.e., 50:50 at Δday = 0 up to 80:20 at Δday = 5). The algorithm is then reapplied in a second iteration, this time flagging high NDVI pixels bounded by two low NDVI pixels.

While this produces MODIS-consistent PS spectral data (${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{MOD}}$) at 30 m, reflectance magnitudes may differ from the L8 surface reflectance distribution due to scale issues and substantial adjacency effects at high spatial resolutions over this specific environment [23]. This in turn will affect the utility of the data for reliably assessing the magnitude of Landsat-consistent reflectance changes. An NDVI-specific bias correction technique is implemented to make the change detection data more consistent with the Landsat surface reflectance data. This approach uses all data during day-coincident Landsat–CubeSat acquisitions to build band-specific bias correction models for discrete NDVI intervals (i.e., 0–0.1, 0.1–0.2, 0.2–0.3 … 0.9–1.0). The NDVI-specific correction models are developed for each spectral band (B_i) using a one-rule implementation of Cubist and ${\mathbf{PS}}_{\mathrm{NDVI}}^{\mathrm{MOD}}$ and ${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{MOD}}$ as explanatory input variables and ${\mathbf{L}\mathbf{8}}_{\mathrm{Bi}}^{\mathrm{SR}}$ divided by ${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{MOD}}$ as a target variable. The trained models are then applied to the full record of ${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{MOD}}$ scene data and the results multiplied by ${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{MOD}}$ to produce estimates (${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{MOD}\_\mathrm{L}8}$) more in line with the L8 based estimates (${\mathbf{L}\mathbf{8}}_{\mathrm{Bi}}^{\mathrm{SR}}$).

The ${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{MOD}\_\mathrm{L}8}$ dataset is used to calculate band-specific relative surface reflectance (SR) changes (dSR_k = [SR_j − SR_k]/SR_k) between any given CubeSat acquisition (j) on P_day and each pair (k) of day-coincident Landsat–CubeSat acquisition [33]. Reference VNIR samples are only drawn from L8 scenes associated with a band-combined dSR_k less than a given threshold percentage (C_th). Qualified reference samples are multiplied with the band-specific correction factors (1 − dSR_k) to account for reflectance changes occurring over the acquisition timespans. In the case where reference data can be drawn from multiple L8 scenes, weights are applied as a function of the acquisition timespan and the magnitude of surface reflectance change. This will serve to produce a smoother result by minimizing potential uncertainties (e.g., atmospheric contamination) associated with the L8 data. C_th is initially set to 15% but increases to 30% for pixels not selected as references under the initial threshold. These temporally adjusted reference samples (${\mathbf{L}\mathbf{8}}_{\mathrm{BI}}^{\mathrm{REF}}$) are used for model training as described next.

2.4.3. VNIR Model Training and Prediction

The 30 m resolution VNIR reference data (${\mathbf{L}\mathbf{8}}_{\mathrm{Bi}}^{\mathrm{REF}}$) specific to each CubeSat acquisition (Section 2.4.2) are used to train four band-specific prediction models relating ${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{TOA}}$ to L8 surface reflectance data. Cubist is used to produce the set of multivariate linear regression models [33]. For this purpose, the blue, green, red, and NIR prediction models use the corresponding PS band as the input explanatory variable after pixel aggregation to the Landsat spatial scale. The maximum number of rules (nrul) is set to 50 in order to effectively reproduce the spatial distribution of the reference data. Resulting band-specific models are then applied to full scene ${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{TOA}}$ data after pixel aggregation (i.e., 30 m) (Figure 2).

These multirule (nrul = 50) models are not directly transferable to the full resolution scale (i.e., 3.125 m) due in part to overfitting issues [33]. To remedy this, simpler two-rule (nrul = 2) regression models are constructed (using exactly the same 30 m inputs) and applied at the finer scale while enforcing scale consistency with the 30 m resolution estimates (nrul = 50) (i.e., pixel aggregation to 30 m will produce near-identical results). Consistency at the 30 m resolution scale is ensured by only utilizing the subpixel (i.e., 3.125 m) variability predicted by the two-rule models relative to the 30 m resolution pixel estimates. In this way, the predicted ~3 m data (${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{L}8}$) retain the spatial features and gradients present in the original 3.125 m resolution data (${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{TOA}}$).

2.4.4. CESTEM-LAI

This section describes the additional steps required to map LAI at the spatial and temporal resolution of the CubeSat data (${\mathbf{PS}}_{\mathrm{LAI}}^{\mathrm{L}8}$) (Figure 2). For this purpose, Cubist is combined with a random forest (RF) regression approach to optimize retrieval accuracies. While RF is powerful at reproducing features in seen observations, the learning process can be computationally demanding and resulting models may not transfer well to “unseen” space and time domains [41]. In comparison, Cubist is computationally efficient and models are generally more portable and better at handling extrapolation beyond the conditions of the training data [41].

Landsat-based LAI scenes (${\mathbf{L}\mathbf{8}}_{\mathrm{LAI}}$) that are day-coincident with a CubeSat acquisition are used for reference sampling in accordance with the approach outlined for VNIR data (Section 2.4.2). The change in LAI over given Landsat–CubeSat acquisition timespans is quantified based on initial PS LAI time series maps (30 m) (${\mathbf{PS}}_{\mathrm{LAI}}^{\mathrm{ini}}$) predicted on the basis of multirule (nrul = 20) Cubist regression models. These are constructed using a multi-day training dataset of ${\mathbf{L}\mathbf{8}}_{\mathrm{LAI}}$, paired with explanatory spectral variables retrieved from the CESTEM corrected surface reflectances (${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{L}8}$) acquired on the same day as the Landsat data. The spectral variables consist of six vegetation indices (VI) with known sensitivity to LAI and include the Normalized Difference Vegetation Index (NDVI) [52], the Optimized Soil Adjusted Vegetation Index (OSAVI) [53], the Green Normalized Difference Vegetation Index (GNDVI) [54], the Modified Triangular Vegetation Index (MTVI2) [55], and the Green Simple Ratio (GSR) [56]. The constructed models are applied to the full record of ${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{L}8}$ images (after VI conversion) for assessing the change in LAI over given acquisition timespans. Subsequently, LAI reference samples on P_day are drawn from the available pool of Landsat scenes (both before and after P_day), but only if the absolute change in LAI over the given timespan is less than the assigned threshold (arbitrarily set to 1.5 m³·m⁻³). Reference samples drawn from multiple scenes are combined and corrected for any change in magnitude occurring over the acquisition timespans to provide a weighted estimate of pixelwise reference LAI (${\mathbf{L}\mathbf{8}}_{\mathrm{LAI}}^{\mathrm{REF}}$) at 30 m, characteristic of conditions during each PS scene (Figure 2).

CESTEM-LAI integrates the RF ensemble-based decision tree approach [57] to learn the relationships between ${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{L}8}$-based VI and ${\mathbf{L}\mathbf{8}}_{\mathrm{LAI}}^{\mathrm{REF}}$ (30 m) on a scene-specific basis. In RF, a large ensemble of decision trees is grown, each using a random subsample of the training data and a random subset of VI variables derived from ${\mathbf{PS}}_{\mathrm{Bi}}^{\mathrm{L}8}$ (see above). Each tree consists of a hierarchy of connected nodes that are defined based on spectral threshold conditions and with each terminal node associated with a unique target prediction value. The final prediction is taken as the average of the individual decision tree outputs. RF is relatively insensitive to parameter specifications [41] that include the number of trees in the forest (ntree), the number of variables tried at each split (mtry), and the minimum number of samples required at terminal nodes (nodesize). In this study, default values for regression were used for mtry (i.e., number of variables divided by 3) and nodesize (i.e., 5), whereas ntree was set to 200. The distributed data structures (ddr) version [58] of the “randomForest” package in R [59] is adopted, as it uses parallel computing to speed up the training and prediction process, which makes it more adaptable to big data analytics [60]. With an assignment of 10 dedicated computing cores (Intel Xeon W-2125 [email protected] GHz), the time for training the RF model (using 27,000 samples) decreased from 63 s (i.e., with one core) to 14 s.

While RF is highly useful for learning potentially complex associations at the 30 m scale based on seen data, model application to unseen data beyond the conditions of the training dataset can be problematic due to inherent transferability and extrapolation issues [41]. In addition, RF model application at the ~3 m scale constitutes a substantial and time consuming computational load. To address this, a pair (nrul = 2) of localized multivariate regression relationships is constructed using Cubist and applied at the ~3 m scale while ensuring scale consistency with the 30 m resolution RF predictions (i.e., pixel aggregation to 30 m will produce near-identical results). This ensures a realistic distribution of predicted LAI at the 3 m scale (${\mathbf{PS}}_{\mathrm{LAI}}^{\mathrm{L}8}$) without impacting LAI magnitudes and variabilities predicted by the RF models at the 30 m resolution scale.

3. Results

A total of 103 PlanetScope images were retained for analyses over the six-month study period after automatic screening for cloud contaminations (Figure 3). The normalized frequency distribution of the number of days between successive acquisitions (Figure 3a insert) indicates a daily clear-sky imaging capacity 62% of the time (i.e., 64 out of 103 scenes were acquired one day after another acquisition) with an enhanced observing potential evident after 5 June (DOY 156). Impressively, 90% of the time the acquisition gap was ≤3 days. Eight scenes were identified as contaminated to some extent based on thresholding the spatial autocorrelation measure (i.e., R²) between successive scenes in the time series, as described in Section 2.4.1. Visual inspection of the scenes demonstrated the utility of this automated approach for retaining only clear-sky scenes in CESTEM processing. While cloud cover is less of an issue in this environment, aerosol loadings can be substantial and may vary significantly from day to day [46]. MODIS deep blue algorithm retrievals over the study domain during CubeSat acquisitions suggest an average AOT₅₅₀ of 0.395 ± 0.15 with a range of 0.118–0.854 (Figure 3b).

In terms of scene coregistration, established pixel shifts were less than or equal to one 3.125 m pixel in 98% of the cases, which indicates well-aligned images over this specific region and an improvement in scene coalignment compared to previous product distributions [33]. The absolute geolocation uncertainty was not directly assessed in this study but is reported to be 10 m root-mean-square error or better (Section 2.1).

The vast majority of scenes (n = 97) were acquired by CubeSats operating in sun synchronous orbits. Most of these (n = 94) had a consistent local overpass time between 9:37 a.m. and 9:45 a.m. (UTC + 3 h), about one hour earlier than the remaining three (10:49 a.m. to 10:56 a.m.). The rest (n = 6) were associated with CubeSats operating in the ISS orbit, with overpass times varying from 6:40 a.m. to 4:18 p.m. (UTC + 3 h) (Figure 3a). Sensor view zenith angles were generally low and less than 4° in 85% of the cases (Figure 3a). In addition to scene-specific sensing geometry, sensor-specific spectral response functions and radiometric quality may affect the consistency of the CubeSat data over time as shown in Houborg and McCabe [33] and explored further below.

3.1. NDVI Time Series Dynamics

The NDVI time series over a multicut alfalfa plot unit (180 × 180 m), calculated from both PS TOA and atmospherically corrected surface reflectances, demonstrate evolving growing cycles from low to higher NDVI, despite evidence of significant day-to-day variability (Figure 4a). The noise in the time series makes it challenging to infer the duration of the various alfalfa development stages and exact timing of key phenological transitions points (Figure 4a). The dynamic range in PS NDVI increases from 0.11–0.55 to 0.13–0.81 after atmospheric correction, but generally remains limited compared to the atmospherically corrected L8 data (0.29–0.92), available approximately every eight days over this section of the farm (Figure 1 and Figure 4a). Large NDVI differences are evident between day-coincident PS and L8 retrievals, particularly during high NDVI, which is partly attributed to inconsistent spectral bandwidths and spectral response functions [33]. However, a consistent bias is not evident (Figure 4a), making systematic correction difficult. While atmospheric contamination is a concern over this desert landscape, with significant impacts of dust aerosols (Figure 3b) and adjacency effects [23], the atmospheric correction does little in terms of reducing the unrealistic short-term fluctuations in NDVI. Overall, using the Planet data at this processing stage for precision agricultural applications would be challenging.

CESTEM was applied to the full record of PS TOA reflectance images using L8 surface reflectance images as reference. Figure 4b demonstrates the production of L8-consistent NDVI during each CubeSat acquisition based on a standard 16-day Landsat revisit. This resulted in a total of eight Landsat scenes available for reference sampling (i.e., dependent scenes) over the six-month period. After applying CESTEM, PS NDVI aligns closely with the dependent Landsat estimates and also tracks features in the independent Landsat estimates quite well (Figure 4b). Most importantly, well-defined and clearly discernable growing cycle characteristics of a multicut alfalfa field appear. A total of five individual growing cycles with a 30–40 day harvesting frequency can be clearly identified, although the frequency of observations is nonoptimal during the first two cycles (i.e., DOY 110–160). The vegetative stages are characterized by a rapid progression in NDVI followed by a period of minimal NDVI change towards plant maturity. Harvesting is clearly indicated by abrupt declines in NDVI, such as on DOY 159, 197, 239, and 274 (Figure 4b). The dynamic phenology and short growing seasons are realistic given the ideal growth conditions facilitated by high day time temperatures (35–45 °C) in combination with ample nutrient and water availability (i.e., pivot irrigation). In addition, the detected phenology is consistent with the MODIS-based information (Figure 4b). What is evident is that Landsat data alone, even with an eight-day frequency, is unable to detect the phenological transitions with sufficient accuracy.

CESTEM was also evaluated using the full eight-day Landsat record for reference sampling (i.e., 14 dependent scenes) in order to demonstrate the impact of data frequency on the CESTEM-corrected time series results (Figure 4c). While the overall NDVI phenology is quite similar between the two simulations (Figure 4b,c), the use of additional reference points has a significant impact on NDVI during the vegetative stages (e.g., DOY 164–174 and 201–210). Since Landsat reference samples during the greening-up phases (in the intermediate NDVI range) are lacking during the 16-day simulation (Figure 4b), reference data during these periods will either be completely missing (i.e., if C_th < 0.3) or inferred from high or low NDVI reference samples via band-specific correction factors derived from the MODIS-consistent surface reflectance change information (Section 2.4.2). While the downscaled MODIS-based NDVI time series depict growing cycles consistent with the CESTEM results, the green-up is typically less rapid and the timing of peak NDVI arrives later, which may be attributed to the 16-day compositing adopted for the BRDF normalization (Section 2.3). As a result, the change information established during these dynamic periods may be associated with greater uncertainty. This is best mitigated by having Landsat reference samples over the full range in NDVI, as that will reduce the impact of the spectral corrections (i.e., due to lower dSR values; Section 2.4.2). The lack of consistency between the downscaled MODIS (${\mathbf{PS}}_{\mathrm{NDVI}}^{\mathrm{MOD}}$) and Landsat NDVI magnitudes is noteworthy (Figure 4c), and is in large part attributable to scale issues and the impact of high-resolution adjacency effects that have shown to increase NDVI by up to 50% [23]. The adopted bias correction technique is effective at ensuring that the downscaled MODIS data (${\mathbf{PS}}_{\mathrm{NDVI}}^{\mathrm{MOD}\_\mathrm{L}8}$) is more Landsat consistent (Figure 4c), which is critical for reliable change detection and reference sampling within the context of CESTEM.

3.2. Spatiotemporal LAI Dynamics

The utility of CESTEM for the spatiotemporal enhancement of biophysical properties ultimately relies on the capacity of the adopted machine-learning technique to establish accurate and spatially representative associations between the CubeSat spectral data and the reference biophysical parameter information. Figure 5 showcases the agreement between L8 reference LAI and CubeSat LAI predicted via random forest machine-learning on DOY 293. The CESTEM-corrected VNIR data were used to construct the explanatory VIs for model training. Visually, predictions accurately reproduce the spatial features and magnitudes in the reference LAI map (Figure 5a). The associated density scatter plot confirms a strong interconsistency and clustering around the 1:1 line, with an R² of 0.991 and a relative mean absolute difference of 4.57%, excluding desert pixels (Figure 5b). The box-and-whisker plot resolves the mean absolute difference (MAD) statistic as a function of LAI, which demonstrates that the median MAD (i.e., the red line) is always lower than 0.09 m³·m⁻³ (Figure 5c). Similar evaluations based on the remaining day-coincident Landsat and CubeSat acquisitions reveal consistently good LAI retrieval performances with the RF models, yielding relative MADs ranging from 3.5% to 6.4% (

Table 1. Statistical metrics of the agreement between Landsat 8 and CubeSat-based LAI during day-coincident acquisitions. Evaluations were performed at the 30 m scale (excluding desert pixels) and include the coefficient of determination (R²), the mean absolute difference (MAD), the relative MAD (rMAD), and the relative mean bias difference (rMBD). $\mathrm{MAD}=\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{PS}}_{\mathrm{i}}{-\mathrm{L}8}_{\mathrm{i}}\right|$, $\mathrm{rMAD}=100\text{×}\left(\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{PS}}_{\mathrm{i}}{-\mathrm{L}8}_{\mathrm{i}}\right|\right)/\overline{\mathrm{L}8}$, $\mathrm{rMBD}=100\times \left(\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{PS}}_{\mathrm{i}}{-\mathrm{L}8}_{\mathrm{i}}\right)\right)/\overline{\mathrm{L}8}$.

DOY	149	158	165	174	181	197	206	213	222	245	270	277	286	293
R2	0.989	0.986	0.993	0.995	0.992	0.989	0.995	0.989	0.984	0.996	0.993	0.994	0.994	0.991
MAD	0.118	0.098	0.057	0.086	0.066	0.082	0.039	0.081	0.092	0.064	0.074	0.071	0.072	0.097
rMAD [%]	4.68	4.69	3.48	3.86	4.35	5.81	4.88	6.30	6.40	4.07	3.92	4.39	3.90	4.67
rMBD [%]	0.46	−0.93	0.13	0.22	0.53	1.02	0.82	1.57	−1.20	0.52	0.12	0.57	0.31	0.68

).

CESTEM produced LAI time series over the multicut alfalfa study-unit track the general features of the about eight-day Landsat-based LAI (dependent or independent) and add close to daily information (i.e., after DOY 156) on LAI dynamics (Figure 6). The value of the high-frequency LAI data is particularly evident during rapid transition events such as harvesting, which are often completely missed during the eight-day Landsat revisit period (Figure 6). The match between CESTEM and Landsat LAI during day-coincident acquisitions is not always perfect, as the scene-specific reference data constitute a weighted estimate of samples drawn from multiple Landsat scenes (Section 2.4.4). This smoothing approach was implemented in an effort to reduce uncertainties associated with the Landsat LAI retrievals [41]. In this case, the range in CESTEM predicted peak LAI over the course of five alfalfa growing seasons is constrained to ~3.9–4.7 m³·m⁻³, whereas Landsat records a peak LAI of 5.3 m³·m⁻³ (Figure 6). As for the spectral CESTEM corrections, pixelwise change information over given CubeSat–Landsat acquisition timespans represents an important prerequisite. While positively biased, the initial LAI (see Section 2.4.4) plotted in Figure 6 is highly correlated with the final CESTEM LAI and as such serves as an effective mechanism to quantify the change in LAI over relevant acquisition timespans (Figure 6).

Figure 7 displays a daily sequence of full resolution (i.e., 3.125 m) CubeSat LAI over the alfalfa pivot (~800 m diameter) highlighted in Figure 1 and Figure 6 during an eight-day green-up period (DOY 206–213) that is bounded by Landsat LAI retrievals (30 m) on DOY 206 and 213. The potential of the very high-resolution (i.e., daily, 3 m) CubeSat LAI retrievals for providing spatial and temporal insights into growth dynamics is clearly evident. With a factor of 10 increase in spatial resolution, intrafield-scale variability that is largely absent in the Landsat LAI becomes resolvable. Importantly, this allows for an improved delineation of field sections exhibiting optimal or suboptimal growth performance and may also provide direct insight into the early detection of disease and other crop health indicators. The benefit of daily image capture is particularly pronounced given the rapid crop development over such a short period of time.

The spatial and temporal advantages afforded by daily CubeSat LAI imagery are further emphasized in a three-day sequence capturing the timing and progression of a harvesting event commenced on DOY 181 and completed on DOY 183 (Figure 8a). In this instance, Landsat LAI is only available on DOY 181, missing critical day-to-day dynamics in vegetation response. Similarly, Figure 8b demonstrates that the exact timing of the harvesting event on DOY 222 can only be inferred from a daily sequence of CubeSat images. While the single Landsat acquisition on DOY 222 captures this specific harvest event, the exact transition day would be impossible to determine based on imagery captured 8–16 days earlier. In addition, the spatial resolution of the Landsat LAI on DOY 222 is generally too coarse for delineating the harvested areas, or progress, with high accuracy. Figure 8c demonstrates the potential of detail-rich CubeSat data for capturing intrafield scale variability in LAI during a green-up phase (DOY 254–257) from an alfalfa field that exhibits significant fine-scale heterogeneity. The closely spaced (i.e., intervals of ~50 m) circular wheel tracks from the rotating center-pivot irrigation system add further complexity to the interpretation of vegetation response and plant development, but these features are effectively resolved by the CubeSat data. In comparison, Landsat LAI represent a mixed signal significantly confounded by bare soil contributions from the wheel tracks. In addition, the temporal insights afforded by 8–16-day Landsat data present a significant limitation during dynamic vegetative periods such as this.

4. Discussion

Constellations of CubeSats in low-earth orbits constitute an unprecedented opportunity to monitor the terrestrial surface with enhanced spatial detail more frequently than ever before. This study demonstrated the capacity of Planet’s constellation of CubeSats, equipped with basic optical imagers (i.e., RGB + NIR), to monitor vegetation dynamics at high spatial (~3 m) and temporal (~daily) resolution. While daily imaging was available during the majority of the six-month study period, original (i.e., not CESTEM corrected) NDVI time series over a multicut alfalfa field indicated significant day-to-day variations not directly attributable to vegetation dynamics or phenology. While atmospheric influences due to dust aerosols were significant, atmospheric correction did not reduce the high-frequency noise. Cross-sensor inconsistencies in radiometric quality, spectral responses, and orbital configuration are the likely contributors to the noise observed in this time-series record. Even though Planet’s constellation of CubeSats undergoes fairly rigorous radiometric calibrations, as outlined in Section 2.1, the residual noise, as demonstrated here, represents a notable limitation for high-fidelity land surface characterization. Additionally, the limited spectral information content contained in broad and partially overlapping bands is inconsistent with the narrower VNIR band configuration typical of conventional multispectral satellite systems such as Landsat and Sentinel-2 [33,61]. While Planet have recently implemented their own 6S-based atmospheric correction protocols, the product was not in wide distribution until October 2017, which is at the end of our study period. Importantly, Planet’s approach does not include adaptations to account for localized aerosol and adjacency effects, which are key factors in our (and many other) study environments. In any case, atmospheric correction of the CubeSat data is not a requirement for running CESTEM that can accept raw digital counts or TOA reflectances as demonstrated in this study.

CESTEM was designed to realize the full potential afforded by CubeSat resources and to allow high-quality, very high-resolution spatiotemporal insights. At its core, CESTEM serves as a radiometric normalization mechanism that has a demonstrated capacity to translate CubeSat observations into L8-consistent imagery, while retaining the radiometric quality of actual L8 observations [33]. In the current study, CESTEM effectively corrected the “raw” multitemporal CubeSat scene data into L8-consistent surface reflectance estimates using an 8–16-day Landsat record as reference data. With this normalized dataset, key phenometrics, such as the start, end, and length of five well-defined and rapidly evolving alfalfa growing seasons, could be accurately retrieved. What is pressingly evident is that Landsat data alone, even with the eight-day frequency available over parts of the study domain (Figure 1), are incapable of detecting the high-frequency vegetation dynamics characteristic of multicut alfalfa. More generally, the time-critical and spatially explicit information needs of precision agricultural management cannot typically be satisfied by conventional satellite resources alone. On the other hand, integrative multisource approaches based on CubeSat imaging hold the potential to overcome traditional spatiotemporal constraints.

Temporal smoothing filters can be highly useful for assessing vegetation phenology and calculating phenological transition dates from satellite time-series data [11]. However, the robustness of the phenology reconstruction can be very sensitive to the choice of smoothing algorithm and parameter calibration [62], and ultimately depends on the availability of a sufficient number of high quality observations during phenological transition events [63]. In this case, common smoothing filters for phenology reconstruction [64] would have limited utility given the real risk of mistaking some of the high-frequency NDVI dynamics (i.e., such as the dips in NDVI between alfalfa growing seasons) as signal contamination. The addition of five-day Sentinel-2 data would certainly increase the capacity of conventional satellite resources for high-resolution land surface characterization, particularly if used synergistically in the manner envisioned by the Harmonized Landsat Sentinel-2 initiative [27]. When combined together, a global median average revisit of 2.9 days can be achieved [25], which constitutes a major advance towards obtaining clear-sky imagery in a timely fashion. Still, the temporal frequency is in many cases suboptimal for monitoring rapid land cover and land use changes in near real time [30], which requires clear-sky image acquisitions immediately before and after the change event. Even disregarding the enhanced spatial resolution, with a daily imaging frequency, the potential of Planet’s constellation of CubeSats for monitoring day-to-day variations in vegetation dynamics and characteristics (after careful corrections for cross-sensor inconsistencies) is clearly demonstrated in this study. Importantly, CESTEM can also be adapted to only exploit satellite data acquired prior to the CubeSat prediction day (see Section 2.4.2), which would enable near real-time predictions needed in the context of precision agriculture.

Unavoidable impacts of cloud contamination on data availability may be partly mitigated by exploiting the expanding capability of the constellation to survey the same location more than once a day. In this manner, the frequency of observations will approximate the revisit capability of the two-satellite (i.e., Terra and Aqua) MODIS configuration, with a truly impressive factor of 80–160 increase in spatial resolution (i.e., 250/500 m to 3.125 m). In addition to disrupting the temporal continuity of the observational CubeSat record, atmospheric contamination may affect the robustness of CESTEM retrievals, as extended timespans between clear-sky L8 overpasses will reduce the number of representative reference samples. Over cloud-prone regions, an effective near real-time phenology reconstruction via CESTEM would probably require a longer multiyear record of searchable L8 scenes together with harmonized Sentinel-2 data continuously updated with new reference scenes as they become available. In this context, the weighting of reference samples drawn from multiple scenes serves to effectively smooth out noise in the time-series reference data. As such, CESTEM also functions as a superior phenological smoothing filter that draws from actual high-frequency (~daily) CubeSat observations without the uncertainties associated with traditional smoothing and gap-filling approximations [65].

In addition to enhanced temporal insights, the spatial detail afforded by the CubeSat data within the framework of CESTEM can be highly advantageous for agricultural management at the precision scale via actionable intrafield information on vegetation dynamics and disturbances. This potential was demonstrated on the basis of CubeSat-derived LAI using CESTEM as the mechanism to enhance Landsat-based LAI both spatially and temporally. The CubeSat resolution (~3 m) was able to resolve intrafield-scale spatial features largely absent in the Landsat retrievals, and daily sequences of CubeSat LAI highlighted rapid and detail-rich crop dynamics completely missed by the 8–16-day Landsat imagery. The adopted RF machine-learning approach was shown to reproduce the L8 reference LAI with high accuracy (i.e., relative MAD ranging from 3.5–6.4%) based on a relatively limited suite of spectral variables derivable from the CESTEM corrected VNIR data. This was despite the fact that a more spectrally enhanced dataset, including critical bands in the shortwave infrared domain, was needed for the Landsat-based LAI retrieval in order to reproduce observed multiday features with sufficient accuracy [41]. Robust CubeSat-based LAI retrievals were facilitated by strong associations between LAI and standard indices in the VNIR [66] in combination with a repeatable redefinition of RF prediction models based on scene-specific reference data. This approach exploits the powerful capacity of data-driven machine-learning techniques to reproduce spatially explicit features in “seen” observations (e.g., reference LAI) while reducing the impact of intrinsic portability and transferability issues [41,67]. As such, CESTEM-LAI takes advantage of high quality and spectrally enhanced sensor data (e.g., Landsat 8) for reference data generation, while reserving the more spectrally limited CubeSat data exclusively for spatiotemporal enhancement purposes.

The downscaling from 30 m to 3.125 m resolution can result in errors due to overfitting and poor generalization [68] and was constrained by employing simpler Cubist-based regression functions at the full resolution scale. In addition to improving computational speed, this approach ensures realistic distributions at the ~3 m scale without impacting magnitudes predicted by the more complex models at the 30 m resolution scale. Although the scene-specific machine-learning approach still represents a limitation when processing large amounts of satellite scenes, continued developments in cloud computing and improvements in data-driven techniques [14] could offer enhanced mapping efficiency and a pathway towards a larger-scale cloud implementation.

A future integration of Sentinel-2 data is likely to add an expanded capacity to other vegetation biophysical properties (i.e., chlorophyll content or carotenoids) which rely on the additional bands in the red-edge for well-constrained retrievals [46,69,70]. In this case, correlative associations between red-edge and combined green and NIR signals, via pigment and vegetation biomass interdependences [71,72], could ensure effective CubeSat-based parameter enhancement on a scene-specific basis [33]. Extendibility to hydrological process variables may also be envisioned, including evapotranspiration (ET), given its strong link to optical remote sensing data and derivable metrics of vegetation density [30]. However, the capacity of CESTEM to effectively resolve spatial variability in water use, considering the complex impacts of irrigation, water availability, and vegetation stress [73], remains to be investigated further. In that context, a comparative analysis against traditional data-fusion algorithms for enhancing observing capacity [16,21], necessitated by the long-standing trade-off between the spatial and temporal resolution of conventional satellites (e.g., Landsat versus MODIS), would be highly relevant. Since CubeSat constellations overcome these traditional spatiotemporal constraints (i.e., by providing Landsat-like data at MODIS frequency), the CESTEM approach offers a unique avenue towards spatiotemporal enhancement with potentially superior retrieval accuracy.

With the key focus being directed towards demonstrating the temporal and spatial advantages of emerging CubeSat resources for vegetation monitoring, a rigorous validation of CESTEM results against independent data was not performed. Elaborate evaluations of temporal LAI dynamics and retrieval robustness should involve continuous records of ground-based LAI automatically observed based on downward-looking digital photos [74,75], upward-looking digital cover photography [76], or indirectly via tower-mounted measurements of spectral vegetation index data [77]. Similarly, evaluations of CESTEM-derived NDVI phenology should include comparisons against continuous NDVI records observed via near-surface remote sensing with digital cameras mounted on towers [78,79,80], as embodied by the PhenoCam network [81]. The integration of ancillary high-frequency and very high spatial resolution (i.e., <5 m) taskable satellite (e.g., WorldView, GeoEye, VeNUS) or UAV resources may be envisioned as an independent mechanism towards a time and scale-relevant verification of spectral and vegetation biophysical features in the CESTEM-derived product suite.

5. Conclusions

Here, we demonstrated the potential of a constellation of CubeSats for providing daily insights into vegetation dynamics at ~3 m resolution. When using the original Planet data, cross-sensor inconsistencies due to variations in orbital configurations, spectral responses, and radiometric quality resulted in a relatively noisy six-month NDVI time series, making phenology reconstruction difficult and susceptible to large uncertainties. The CubeSat Enabled Spatio-Temporal Enhancement Method (CESTEM), based on a multisource (i.e., CubeSat, Landsat 8, MODIS) and multiscale (i.e., 3 m, 30 m, 500 m) machine-learning technique, facilitated radiometric normalization of CubeSat top of atmosphere (TOA) RGB + NIR data into Landsat-8-consistent surface reflectances. The CESTEM-corrected NDVI time series revealed clearly discernable and smoothly evolving growing cycles with a crop-growth period on the order of 30–40 days. These results highlighted the impressive capacity of near-daily CubeSat imaging for resolving high-frequency vegetation dynamics and determining the exact timing of key phenological transitions.

In addition to providing temporally consistent phenology information via NDVI, CESTEM demonstrated its potential for enhancing Landsat-based LAI both spatially and temporally, while maintaining consistency with Landsat-scale LAI. Daily sequences of CubeSat-based LAI highlighted the spatial resolution advantage and provided critical temporal insights into within-field variations in vegetation health and condition, the rate of vegetation green-up, and the timing and progress of harvesting events, features that were largely uncaptured by the 8–16-day Landsat imagery.

While CubeSat constellations are on path to reshape our earth-observing potential, the realization of their full potential will likely require synergistic exploration with conventional satellite systems. By serving as a multipurpose data-fusion scheme for radiometric normalization, time-series smoothing, and spatiotemporal enhancement of biophysical properties, CESTEM provides this needed framework and represents a novel paradigm with application potential across diverse aspects of earth observation. Optimization of mapping efficiency (e.g., via cloud computing) and rigorous evaluation of retrieval robustness and observing capacity over a wide range of ecosystems and environmental conditions constitute key objectives for future research.