Mapping and Monitoring Rice Agriculture with Multisensor Temporal Mixture Models

:


Abstract:
Rice feeds more humans than any other crop on Earth.Accurate prediction of the timing and volume of rice harvests therefore has considerable global importance for food security and economic stability, especially in the developing world.Optical and thermal satellite imagery can provide critical constraints on the spatial extent of rice planting and the timing of rice phenology.We present a novel approach to the mapping & monitoring of rice agriculture using Temporal Mixture Models (TMMs) derived from parallel spatiotemporal analyses of coincident optical and thermal Landsat image time series.Using the Sacramento Valley of California as a test area, we characterize regional rice phenology in terms of both fractional vegetation abundance (Fv) and brightness temperature (Tb).We compare satellite Tb retrievals to station data and find uncorrected Tb to compare with the upper bound of the envelope of air temperature observations to within 3°C on average.Results from parallel spatiotemporal analyses of coincident Fv and Tb image time series over the 2016 & 2017 growing seasons suggest that TMMs based on single year image time series can provide simple and accurate maps of crop timing, while TMMs based on dual year image time series can provide comparable maps of year-to-year crop conversion.Fv time series data show particular promise for estimating crop timing, while Tb appears particularly well suited for discriminating between rice and other crops.We also build a sample model midway through the 2018 growing season to illustrate a potential near-realtime monitoring application.Field validation confirms that the mid-2018 monitoring model provides an accurate upper bound estimate of the spatial extent and relative timing of the rice crop, even under conditions of relative data scarcity.The implications of these results could have potential utility for further analyses of precision agriculture, pest management, evapotranspiration (ET) and cropping practice verification.

Introduction
Rice agriculture is critical for global food security.Much of the developing world relies on rice production for both subsistence and/or commercial purposes.Rice is the largest food source for the world's poor and more than half of the world's population relies on rice as its staple food (Muthayya et al., 2014).In terms of nutrition, rice provides 21% of global human per capita energy and 15% of per capita protein.More land area is used for rice production than any other agricultural crop (Global Rice Science Partnership, 2013).
Predictability of global rice production requires accurate forecasting of both planting and harvest.Mapping and monitoring efforts can improve prospective rice harvest estimates by constraining fundamental variables such as planted area and timing of key phenological transitions.In addition, concerns important to commercial growers like optimization of water and chemical additive use, assessment of extent and severity of weeds and pests, and independent verification of alternative cropping practices also rely on timely landscape scale biogeophysical observations.One cost effective tool to meet these needs is spatiotemporal analysis of satellite imagery.Free decamater (10 to 100 m) resolution multispectral optical Earth observing imagery already exists through the Landsat program, with data extending back to the early 1980s.The Landsat program has been supplemented since 2015 by the addition of the Sentinel-2 constellation (Drusch et al., 2012).Additional hectometer (100 to 1000 m) resolution imagery is available through the MODIS, Sentinel-3 and VIIRS programs.The global archive of meter to hectometer resolution optical satellite imagery is expected to continue to grow into the future as further missions are planned in both public and private sectors.This rich set of observationswith global coveragecan facilitate both retrospective and prospective analyses for a wide range of agricultural applications.
Using optical satellite imagery to accurately map and monitor agriculture requires rigorous characterization of the crop's life cycle and cropping practices in order to confidently relate time series of reflectance spectra to time series of bigeophysically relevant quantities.Fortunately, in recent decades significant effort has been devoted, and much success has been achieved, in studying and improving rice agriculture in the myriad contexts it is grown around the world (Khush, 2003).As a result, rice is one of the best-understood crops on Earth from an agronomic perspective.However, significantly less work has been done to characterize how the stages of rice agriculture are recorded by optical satellite imagery.To date, most applications of satellite imagery for rice agriculture have focused on discrete classification and mapping of rice extent, rather than on the physical properties and processes that characterize the phenological cycle of rice crops.
An obvious distinction between rice and most other crops is that rice is often (but not always) grown in paddies.This fundamental difference results in a distinctive physical signature for rice when imaged by satellites at both optical VSWIR (Visible through Shortwave Infrared, 0.4-2.5 m) and TIR (Thermal Infrared, 9-14 m) wavelengths.The evolution of the rice crop throughout a growing season therefore occupies a unique trajectory through the geophysical parameter space of temperature and optical reflectance.This distinction provides leverage for producing accurate cropspecific maps of rice, which are not generally possible with other crops.
In this work, we present a straightforward, physically-(as opposed to statistically-) based approach for the mapping and monitoring of rice agriculture using time series of coincident Landsat VSWIR and TIR observations.Using the Sacramento Valley of California as our example area, we apply the methodology of spatiotemporal analysis, described in detail by (Small, 2012), to the special case of a parallel analysis of the twin data streams from coincident Landsat 8 OLI (Operational Land Imager) and TIRS (Thermal Infrared Sensor) image time series.The resulting VSWIR and TIR data spaces each characterize complementary features of the biogeophysical evolution of the agricultural landscape.Using endmembers of distinct temporal trajectories, independently derived from each data stream, we demonstrate our method by mapping both the phenological timing and presence/absence of rice crops across 2 years.Finally, we also show the seasonal evolution of the landscape in Temperature-Vegetation space.This suggests more detailed investigation into the spatiotemporal analysis of evapotranspiration (ET) as an attractive avenue for future work by potentially directly addressing several of the 10 outstanding research questions identified by the recent review of (Fisher et al., 2017).
While a considerable body of literature exists focusing on mapping rice from optical satellite imagery, to our knowledge previous work has largely focused applying a variety of complex statistical algorithms to VSWIR data, usually without using thermal information.Some previous studies have included thermal into their algorithms, but rarely in a physically based way.Our approach has the benefits of both relative mathematical simplicity and explicit grounding in biogeophysical processes.

Rice in the Sacramento Valley
The Great Central Valley of California is comprised of two major units: the northerly Sacramento Valley and the southerly San Joaquin Valley.The climate and soil of the Sacramento Valley are well suited to rice production: clay-rich soils which are easily able to support paddy water dominate much of the valley floor and hot, dry summers provide a supportive growing environment for the rice plant with minimal risk of mid-or late-season rains damaging the crop.As a result, rice yields in the region are some of the highest in the world, with average yields now exceeding 9000 kg/ha (8000 lbs/acre).California is a major rice-growing region, producing 75% of the medium grain and 98% of the short grain rice grown in the United States.California grows more rice than any U.S. state other than Arkansas, which dominates in long grain rice production.Of the rice production in California, approximately 95% occurs in the Sacramento Valley, mainly in Colusa, Sutter, Butte, and Glenn counties.For a thorough yet concise background on the history and current state of rice production in California, see (Geisseler & Horwath, 2013).
The same clear, dry summer climate that supports rice production in the Sacramento Valley also benefits satellite imaging.This is reflected in the data archive, as over 10 cloud-free satellite imaging acquisitions were acquired for every growing season since the launch of Landsat 8 in 2013.When present, smoke from summer wildfires hinders satellite imaging.However, wildfires in California are most common in late summer, after the rice crop has been planted & established and high frequency imaging is less critical for crop monitoring.These factors, combined with the economic productivity of the California rice crop and the abundance of ground-based observational data, suggest that the Sacramento Valley may be a nearly ideal location to characterize the evolution of rice agriculture in the geophysical parameter space relevant to optical satellite imaging.

Spatiotemporal Analysis and Temporal Mixture Models
Our approach to the quantitative analysis of satellite image time series is described in detail by (Small, 2012).The method can be decomposed into two stages: characterization and modeling.The characterization stage involves performing an Empirical Orthogonal Function (EOF) analysis of the satellite image time series.This initial EOF analysis considers the image time series as a matrix with each geographic pixel occupying a row and each time step occupying a column.An image time series therefore generally contains a large (generally >10 6 ) number of rows and a much smaller (generally <10 3 ) number of columns.The EOF analysis is based on a Principal Component (PC) transform which computes the covariance (or correlation) matrix of this dataset and then performs an eigenvalue/eigenvector decomposition in order to determine the dominant modes of uncorrelated spatial and temporal variance.The transformed spatiotemporal data are then examined in this new, optimized space, referred to as the Temporal Feature Space (TFS) and depicted graphically as scatterplots of the spatial Principal Components in conjunction with the corresponding temporal EOFs.The most distinct temporal patterns (i.e.phenologies) occupy apexes of the data cluster in the TFS while (in the case of linear mixing) binary mixtures of these patterns occupy sharp edges between the apexes.The temporal patterns of the geographic pixels associated with these apexes are identified as temporal endmembers (EMs).
The temporal EMs are then used as the basis for a linear Temporal Mixture Model (TMM).This model treats every pixel time series as a linear combination of the EM time series identified during the characterization phase of the analysis.Inversion of the mixture model effectively unmixes each pixel time series to yield an estimate of the relative contribution (fraction) of each temporal EM to each pixel time series in the dataset.The result is a continuous map showing the relative contribution of each temporal EM pattern at each point in space.The physical meaning of the model result is therefore straightforward (quantitative description of similarity to each EM) and the model inversion is mathematically transparent (ordinary least squares inversion of an overdetermined system of linear equations).This approach also has the benefit of generality, as no portion of the analysis imposes assumptions specific to particular geophysical variables or sensing modalities.The only implicit assumption is that spatiotemporal variance in the data represents potentially useful information.

VSWIR versus Thermal Imaging
Multispectral optical satellite imaging generally leverages two electromagnetic domains: the reflected optical wavelengths and the emitted thermal wavelengths.VSWIR instruments measure Earth-leaving radiance in the 0.4 to 2.5 m wavelength range.In this portion of the spectrum, the majority of outgoing radiance is comprised of reflected solar illumination.Data from these instruments can be easily reduced to apparent reflectance, or the fraction of incident sunlight which is reflected by the Earth surface & atmosphere at each pixel, by normalizing the measured radiance (outgoing) by the known solar illumination (incoming).Reflectance can be used to discriminate between Earth materials on the basis of chemical absorptions and microscale physical structure.
In contrast, thermal satellite imaging measures Earth-leaving radiance in the 9 to 14 m range.Because the solar emission curve falls off precipitously at optical wavelengths, the radiance measured by thermal imaging sensors is dominated by the temperature and emissivity of the Earth surface and atmosphere.Because temperature and emissivity are measures of fundamentally different physical attributes than reflectance, VSWIR and TIR data are often analyzed separately.However, the combined use of VSWIR and TIR imagery remains an active area of researchparticularly for characterization of evapotranspiration and surface energy balance.
For some applications, atmospheric corrections are attempted for both VSWIR and thermal satellite imagery.Atmospheric correction models attempt to remove the effect of absorption and scattering along the two-way path (for VSWIR) or one-way path (for thermal) of radiation through the atmosphere.However, atmospheric correction is dependent upon the ancillary data and choice of model used.For many applications the magnitude of the signal far exceeds that of the correction, and atmospheric correction does not significantly change the analysis outcome.For this analysis, we do not attempt atmospheric correction because we rely on the large differences in reflectance and thermal emission between rice and other land cover types throughout the growing season.

Data & Preprocessing
All data used in this analysis are publicly available.MODIS data were downloaded using both the USGS EarthExplorer website (http://earthexplorer.usgs.gov/)and the MODIStools R package (Busetto & Ranghetti, 2016).Landsat data were downloaded using the USGS GloVis website (http://glovis.usgs.gov/).Weather station data were downloaded from the CIMIS website (https://cimis.water.ca.gov/).MODIS thermal data were downloaded as 8-day composite images (MODIS Product MOD11A2, Version 6) and converted to Daytime Estimated Land Surface Temperature by multiplying by the scaling factor 0.02.MODIS vegetation data were downloaded as 16-day vegetation index composite images (MOD13Q1).The Enhanced Vegetation Index (EVI) product was used and values were retrieved by multiplying by the scale factor 0.0001.
The Landsat 8 satellite platform hosts both the Operational Land Imager (OLI) and the Thermal Infrared Sensor (TIRS) instruments.The OLI instrument operates in the 0.4-2.3m optical wavelength range and the TIRS instrument operates at much longer (10 -13 m) thermal wavelengths.OLI data were calibrated from DN to exoatmospheric reflectance and TIRS data were calibrated to apparent brightness temperature using the standard calibration procedures described by the Landsat Data Users Handbook (USGS, 2016).All TIRS data shown are from Band 10 and OLI data are from bands 2 to 7.
The physical bases and foci of this study are the relationships between the annual phenological cycles of the rice crop and other vegetation (agricultural and indigenous), seasonal variations in soil moisture (from precipitation and irrigation), evaporation and transpiration of water and their combined impact on surface temperature.Understanding of these relationships is facilitated by conversion of the measured optical reflectances to fractional abundance of the four most distinctive and functionally important components of the landscape: vegetation, soil, water and shadow.In this study, we focus on the relationship between water and vegetation, but this approach also provides information on soil and shadow which are described in a separate study.The 6-band OLI multispectral reflectance images were used to derive estimates of Substrate (rock & soil), Vegetation, and Dark (water & shadow) subpixel fraction abundance by inverting a linear spectral mixture model using the standardized global endmembers of (Sousa & Small, 2017).The result is a set of coincident image time series of Substrate, Vegetation and Dark fractions and brightness temperatures for each 30 x 30 m pixel in the study area for the 22 cloud-free overpasses of the study area in 2016 and 2017.

Regional Thermal Setting
Figure 1 shows a tritemporal thermal composite image of the study area used for this analysis in regional context.The northern portion of the Great Central Valley of California (area within & around the white box) is one of the most productive rice growing regions on earth.The Central Valley lies between the Sierra Nevada mountains to the east and the Coast Ranges to the west.Further east beyond the Sierra Nevada lies the geophysical provinces of the Basin and Range and Mojave Desert.To the west of the Coast Ranges lies the Pacific Ocean.The area to the north of the valley is mountainous and includes the Cascades Volcanic Range.
Land use in the valley is dominated by intensive agricultural production of a wide range of commodity and specialty crops.The northern (Sacramento Valley) portion is cooler and wetter than the southern (San Joaquin Valley) portion.The Sacramento Valley is characterized by Köppen Climate Csa (Hot Summer Mediterranean).The vast majority of precipitation is delivered between the late fall and the early spring, with summers being characterized by dry, hot days with sparse cloud cover.A wide range of soil types exist in California and are distinguishable on the basis of VSWIR reflectance (Sousa & Small, 2018).Clay-rich soils are common in the northern portion of the valley, providing a nearly ideal environment to naturally support the standing water commonly used in intensive rice agriculture.Water in the region is highly controlled through a network of government-run water projects.Effectively all summer agriculture in the Central Valley is irrigated.
The colors on Figure 1 illustrate at a glance the typical thermal evolution of the region over the first half of the calendar year.In this figure, the Blue, Green, and Red channels of the image represent brightness temperatures in January, May, and August of 2016, respectively.No locations are bright blue because the landscape is universally colder in January than in May or August.The Mojave Desert appears as yellow (Red + Green) because its temperature generally mimics the annual cycle of solar insolation (bottom time series).In contrast, the Basin & Range further north takes longer to heat up and maps with a more reddish hue.High elevations in the Sierra Nevada and Coast Ranges appear as dark red because seasonal snowpack results in extremely low winter temperatures that then rapidly rise once the snowpack melts in late spring/early summer.Agriculture in the Central Valley and Snake River Plain appears as green because surface temperature in late spring/early summer generally increases before crops have reached maturity.The low thermal inertia of bare, dry fields therefore results in high brightness temperatures in May.Once crops have grown, however, the associated evapotranspiration (ET) significantly cools the landscape during the summer growing season.Typical annual thermal phenologies for Central Valley Rice from the MODIS sensor are shown.In comparison, fields in the Central Valley left fallow, and the grasslands on the foothills of the Sierra Nevada and Coast Ranges, show annual thermal cycles more similar to the desert than the adjacent agricultural fields.Coastal areas, dominated by the thermal inertia of the Pacific Ocean and inundated in fog for much of the year, have similar brightness temperatures year-round.Inland lakes show somewhat more of an annual cycle, but still are characterized by relatively cool, stable brightness temperatures, therefore mapping as dark.

Weather stations
The California Irrigation Management Information System (CIMIS) maintains an extensive network of field observations of meteorological data throughout the state of California.Four of these stations are relevant for the purposes of this study.Figure 2 shows a map of their locations within the area of the white box from Figure 1.The stations bound the region dominated by rice agriculture (darker green) to the north, south and west.A small mountain range, the Sutter Buttes, is located near the center of the rice zone, providing a convenient landmark for image interpretation.
Hourly time series of soil temperature (green) and air temperature (red) are plotted for each of the stations for both years in the study, along with Landsat-derived brightness temperatures for a each of the 9 pixels in a 3x3 box surrounding the station location.For stations Durham and Verona, the Landsat brightness temperature clearly captures the upper portion of the envelope of air temperatures with very little bias.This is also true for Biggs in 2016, but in 2017 the two records diverge somewhat at the air temperatures of the Biggs station increase but the brightness temperatures decrease.The opposite is observed with the Williams station (installed August 2016), where Landsat brightness temperatures are consistently higher than the station air temperatures.The Williams air temperatures are also consistently cooler than the soil temperatures for 2017 but not 2016.
The variability in agreement between the station data and brightness temperatures may be related to variations in microclimate related to the siting of each station.While all four stations are immediately surrounded by irrigated pasture, the Durham and Verona pastures are situated within a matrix of tree and row crops, while the Williams and Biggs pastures are surrounded by rice and/or wetlands.It is possible that the more saturated surface hydrology surrounding the Biggs and Williams sensors introduces atmospheric boundary layer effects that impact brightness temperature retrievals.Regardless of these complexities, it is clear that the shape and rough amplitude of the Landsat brightness temperature time series generally agree with the station data, in some cases remarkably precisely, especially given that no atmospheric correction of the Landsat thermal data was attempted.

Vegetation Phenology
Figure 3 demonstrates a representative annual phenological progression of the region using a series of false color composite images from the calendar year 2016.Each image is enhanced identically, so colors are directly comparable.In these images, the red channel corresponds to the Shortwave Infrared (SWIR, 2.2 m), the green channel to the Near Infrared (NIR, 0.86 m), and the blue channel to the visible wavelengths (Vis, 0.55 m).Soils thus generally appear as brown, photosynthetic vegetation as bright green, and standing water as black.As vegetation senesces, it grades from bright to dull green, through orange and yellow, and into light brown.Images are only shown from April through November because winter months are frequently covered in cloud and rice is not grown during winter, with many fields intentionally flooded to provide habitat for migrating birds.
Rice is generally grown in the clay rich soils to the east and west of the Sacramento and Feather River channels.The rice growing area can be visually identified as the area of standing water (black) in the May 26 image.The reason for this is due to the phenology of the rice crop in the area, summarized as follows: Many, but not all, rice fields are flooded during winter (November through February or March) to provide habitat for migrating birds.In April, the flooded fields are drained and all rice fields are plowed in preparation for the rice crop.Rice fields are then generally flooded in May and seeded by airplane into standing water.The rice then greens up during the summer months and begins to senesce in September, continuing to senescence and being harvested through September and October.By November the harvest is complete.After harvest, some fields are left bare and others are flooded for winter bird habitat, and the cycle begins again.
While rice dominates the landscape, other vegetation also exists.The foothills of the Sierra Nevada, Coast Ranges, and Sutter Buttes largely host rainfed grasslands used for grazing.The phenology of these regions is nearly opposite that of rice agriculture because the grasslands are green during winter and generally senesced and non-photosynthetic during the summer months.In addition, a wide range of other row crops and orchards are grown in the valley with a diversity of cropping systems and phenologies.In this portion of the Central Valley, these other crops are generally grown close to the Sacramento and Feather River channels.These regions also have more complex networks of field geometries and smaller field sizes.Human settlements also occupy a small portion of the landscape, hosting an even more complex mixture of evergreen and deciduous trees, shrubs, and grasses.
The time series in the lower right of Figure 3 show typical phenological progression of rice crops (solid) versus native grassland (dashed) in terms of both vegetation abundance (green) and temperature (red).Lines indicate data from MODIS composites and points indicate data from Landsat time series.The more frequent revisit time of the MODIS sensor provides context for the higher spatial resolution Landsat observations.For the MODIS sensor, vegetation is quantified as EVI but for Landsat it is quantified as more accurate Vegetation Fraction (Fv) from a 3-endmember spectral mixture model.Both quantities are scaled from 0 to 1 and the values have been shown to be comparable, with a small bias (Sousa & Small, 2017).
As expected from the image time series, the vegetation abundance curves for the rice and grassland pixels are nearly 180 out of phase.Rice has a strong peak in vegetation abundance during the summer growing season, and also a minor peak during the winter fallow season, presumably from non-rice vegetation growing in flooded fields.The grassland pixels have very low amounts of photosynthetic vegetation during the hot and dry summer months because they are rainfed, but reach broad peaks during the winter rainy season of comparable but somewhat smaller amplitude than the agriculture.The benefit of using both MODIS and Landsat for these purposes is evident in that there were no cloud free Landsat images available over a 3-month period during the wintertime, but sufficient usable MODIS images were acquired to produce a reasonable composite product.

Thermal Phenology
The fundamental physical process underlying the thermal phenology of the region is the seasonal insolation curve.Insolation is sinusoidal with a minimum at the winter solstice and a maximum at the summer solstice, explaining the winter minima for both rice and grassland.However, grassland achieves a much greater summer temperature than the rice because the dry soil and non-photosynthetic vegetation (or absence of vegetation) can support very little evapotranspirative cooling.As a result, the shape of the grassland thermal curve in Figure 3 is roughly sinusoidal, albeit with minor skewness.On the other hand, rice agriculture is not only irrigated but grown in standing water for much of the summer, resulting in substantial cooling both from evaporation of the paddy water below the canopy and transpiration of the respiring and photosynthesizing vegetation.This results in a relatively stable brightness temperature during the summer months, yielding a flatter top to the temperature curve.
Although pixel level phenological comparisons can be suggestive, these results must generalize in a spatially coherent way if they are to be used for the purposes of mapping.Figure 4 illustrates the spatiotemporal evolution of the field of brightness temperatures and subpixel fraction abundances for the same acquisitions used Figure 3.The left pair of panels shows a tricolor composite image where Red, Green, and Blue correspond to relative fractional abundances of substrate (soil or non-photosynthetic vegetation), photosynthetic vegetation, and dark materials (shadow or water).These relative abundances were computed using the 3-endmember spectral mixture model described in the Methods section.Comparison of these images to the false color composites from Figure 3 provides an intuitive illustration of how the 3 EM spectral mixture model represents the high dimensional reflectance information present in multispectral satellite imagery.
The right pair of panels corresponds to the time series of thermal imagery co-acquired with the optical imagery by the Landsat 8 TIRS sensor.The image time series is scaled from purple to red, with values over 50 C represented as white.While air temperatures of 50 C in the region are very uncommon, air temperatures in the 40 -45 C range do regularly occur throughout the region in summer months.Brightness temperatures of the land surface can also achieve much higher values than air temperatures.
The phenological progression of the thermal image time series complements that of the fraction image time series.In late April, the only prominent spatial patterns in the thermal field are due to elevation, rivers and lakes, and transient clouds.However, as the spring progresses into summer, clear differences emerge between rice and grassland areas.In the late August, September, and October images, differences within the region of rice agriculture become apparent, likely on the basis of evapotranspiration decreasing at variable rates as some fields senesce sooner than others, and some fields are drained before others to prepare for harvest.Once the fields are harvested, substantial differences between land covers are greatly reduced.In fall, the remaining vegetation can actually be somewhat warmer than its dry surroundings due to the greater heat capacity of its leaf water.
Finally, the scatterplots in the lower left of the thermal images show the relationship between brightness temperature (x-axis) and vegetation fraction (y-axis) for every pixel in the image.These scatterplots are thus an alternate coordinate system with which to understand the dataset.The scatterplots are color-coded so warmer colors represent more pixels and cooler colors represent fewer pixels.The grayscale background silhouette of the full multitemporal point cloud shows the composite space of all the images together.Because so many pixels in the image correspond to rice agriculture, they generally form the largest cluster in temperature-vegetation (TV) space and are represented by warmer colors.In April, for instance, rice fields are recently plowed so they are cool and totally unvegetated, so the warm colored cluster plots at the bottom center of the space.In May, the fields are flooded they become even cooler but their vegetation fraction does not change appreciably, so the cluster moves slightly to the left.In late June, fields vary widely from just planted (late crop) to nearly mature (early crop), resulting in a wide range of vegetation fraction.ET cools the crop so its temperature does not change appreciably, resulting in a vertically elongated cluster.For the rest of the summer the fields do not change appreciably in vegetation or temperature, so the cluster is stable, until September when many fields are drained and begin to senesce, resulting in the cluster elongating toward higher temperatures.Harvest of the rice crop results in a stepwise change of the TV properties of a rice field.The separation of the point cloud splitting into two distinct clusters by the October 1 image agrees with 2016 USDA estimates that 54% of the California rice crop was harvested by October 9 (Childs, 2016).By mid-November, all the rice has been harvested and the landscape shifts into its winter state.Interestingly, the few usable winter images that exist plot in a nearly disjoint portion of the vegetation temperature space, suggesting that the spatial relationship between vegetation and temperature during winter may be dominated by fundamentally different physical processes than during summer.

Spatiotemporal Characterization
The preceding figures can be summarized by the following set of observations: 1) The thermal phenology of rice agriculture is substantially different in amplitude and shape from other land cover types in the region; 2) The parallel evolution of both thermal and vegetation phenology can be explained in terms of the surface hydrologic cycle and growth cycle of the multiple phases of rice crops; 3) The spatiotemporal variations in brightness temperature have substantial differences from those of the vegetation abundance, despite their interdependence; and 4) The spatiotemporal variations in both brightness temperature and vegetation fraction can be explained using fundamental physical principles.
Put together, observations 1, 2 and 3 suggest that including Landsat brightness temperature in a phenological analysis could add information that is not present using vegetation abundance alone.Observation 4 suggests that such a phenological analysis could be based on straightforward physical principles with a bare minimum of model complexity.
A potential approach to this mapping problem is the spatiotemporal analysis method of (Small, 2012).This approach is rooted in spatiotemporal characterization using Empirical Orthogonal Function analysis of the image time series, followed by modeling using linear temporal mixture models.One primary benefit of this approach is that it imposes no assumptions about the functional form of the phenology, but rather characterizes the temporal patterns based on the data itself .Other benefits to the method are its simplicity, robustness, and generalizability, as well as its ability to generate results with straightforward physical interpretations.Because they are based on physical principles and easily identified temporal endmembers representing known phenologies, maps derived from linear model inversion provide a degree of uniqueness of solution that is almost never provided by discrete thematic classifications which are based on ad hoc selection of land cover classes and training data.
In order to use this approach with a dual data stream, two decisions must be made.The first decision is whether to incorporate both streams into a single analysis or to analyze each in parallel.We present the results of parallel characterization method here because it is conceptually simpler and achieves the objectives of the current study.The combined characterization and analysis of thermal and vegetation phenology lends itself naturally to the study of ET dynamics and is the focus of a separate study.
The second decision is whether to analyze the two years together or separately.Because significant benefits can result from either approach, we present results of both ways in Figures 5a   and 5b, respectively.Because the two single-year characterizations yielded similar results, we present only the 2016 results as characteristic of both years of the dataset.
For each case, we began by first conducting an EOF analysis of the relevant image time series.We chose unnormalized (covariance-based) rotations in each instance.Normalized (correlation-based) rotations were also investigated, but the resulting feature spaces were less informative.In each case, over 65% of the variance was represented in the first 3 dimensions of the data.While this number is low enough to suggest that informative structure may be present in the higher dimensions of the dataset, the three low-order dimensions capture the most relevant phenological patterns necessary to distinguish rice from other crops and to characterize its seasonal and interannual variability.
Figures 5a and 5b summarize the characterization stage of the parallel analyses for single year (2016) and dual year (2016+2017) time series, respectively.Characterization is based on the 24 x 24 km subset shown in the white box in Figure 4. We used this spatial subset for our characterization because it is dominated by rice agriculture with a wide range of crop timing.In every case, the loadings of the first 3 spatial PCs are shown as scatterplots.These scatterplots can be conceptualized as representing the location of each geographic pixel in the image time series in a 3-dimensional space, known as the Temporal Feature Space (TFS), in which the axes represent the uncorrelated temporal patterns of maximum variance.Clusters in this space correspond to sets of geographic pixels with similar temporal trajectories of the quantity of interest over the course of the two years of the study.Corners in the space represent endmember (EM) temporal patterns which represent pixels with the most distinctive temporal patterns in the image time series.Pixel time series laying inside a convex hull connecting the EMs can be represented as weighted linear combinations of the EM time series.Because these EMs represent observationally derived basis functions, they will form the basis of the forward model to produce fraction maps.
Figure 5a shows a comparison of the TFS for image time series of both vegetation fraction and brightness temperature for the year 2016.The TFS of the vegetation fraction image time series shows four distinct apexes representing the temporal EMs: Early Rice, Late Rice, Wetlands/Evergreen, and Water/Fallow.Time series of vegetation fraction corresponding to these EMs are shown in the lower right.All 4 EMs are clearly distinct on the PC 3 vs 2 scatterplot.The point cloud in this scatterplot forms a cross shape, with the axis between EMs 1 and 2 corresponding to phenological timing and the axis between EMs 3 and 4 corresponding to overall vegetation abundance.The Water/Fallow EM forms the sharpest corner, as expected given the very small amount of variability expected in the vegetation fraction of water bodies and fallow fields through time.Early and Late EMs are less sharp but still clearly defined, indicating substantial variability in the timing of the phenological signal.However, the Early and Late clusters are also visibly distinct, indicating two clear phenological groups.The Wetlands EM is the most diffuse, as expected given the wide range of vegetation types, hydrological regimes, and land management strategies for wetlands in the area.
The TFS of the brightness temperature image time series shows substantially more clustering than that of the vegetation fraction time series.This indicates that more pixels have more similar brightness temperature trajectories than vegetation fraction trajectories.At least 7 EMs are identifiable.The most distinct from each other are Water (EM 4; blue) and fallow fields (1; yellow).Water time series are particularly distinctive as their very low PC 2 and very high PC 1 and 3 values position them as disjoint from the remainder of the dataset.The remainder of the space partitions into a) differences in timing of rice phenology, and b) differences between rice, non-rice agriculture, and wetlands.Wetlands (EM6) occupy the corner of the point cloud closest to water, and also form an axis grading into fallow fields (EM1).The Fallow-Wetland axis is nearly orthogonal to the Early-Late axis of rice phenology, most clearly seen in the PC 1 vs 3 projection.Non-rice agriculture (EM5) plots on a continuum between the Early Rice and Fallow EMs.Finally, double cropping (D) is clearly identifiable in both the Fv and Tb feature spaces.Few fields in the spatial subset used for this rotation practice double cropping, so this EM is sparsely populated.
Figure 5b shows characterization of the dual year 2016+2017 image time series.In both time series, the early/late phase information that dominated the single year time series is suppressed (though still present), and the structure of the TFS is dominated by year-to-year differences in crop presence or absence.In the Fv space, both Wetlands (EM4; green) and Water/Fallow (EM5; blue) are still clear EMs.The largest cluster, however, is now that of fields cropped in both 2016 and 2017 (EM3; red).This cluster is clearly distinct from the small cluster of fields cropped in 2017 only (EM2; light magenta) and the much larger cluster of fields cropped in 2016 only (EM1; yellow).The fact that the 2016 only cluster is much larger than the 2017 only cluster is concordant with official USDA estimates of 42,000 fewer acres of California rice planted in 2017 than 2016 due to higher prices for competing commodities and severe early season flooding (Childs & Skorbiansky, 2017;NASS, 2017).
The dual year Tb space again shows more clustering and overall complexity than the corresponding Fv space.At least 6 EMs are identifiable.Again, the most distinct from each other are Water and fields that are fallow in both years (F/F).Water time series are again particularly distinctive as their very low PC 1 and PC 3 values position them as disjoint from the remainder of the dataset.The remainder of the space partitions into a) differences in timing of phenology for fields under rice cultivation in both years, and b) the presence and absence of a rice crop in each year.The axis corresponding to phenology is nearly orthogonal to the axis corresponding to presence/absence of crop in each year.The potential utility of this information is discussed below.
While information about the phase of the rice is clearly present in both single year datasets, the Fv space EMs captures more end-of-season variability than the Tb EMs.This could be due to gradual browning of the top canopy layer over weeks to months (resulting in progressive decline in Fv over the second half of the growing season) being accompanied by continual cooling by ET of the green vegetation and paddy water below (resulting in minimal change of Tb).However, the Tb dataset clearly shows enhanced ability to discriminate between rice and non-rice crops, likely due to early season paddy water producing a unique signature in Tb but not Fv.For the dual year characterization, both Fv and Tb clearly discriminate between areas cropped both years versus each individual year, but again Tb shows superior ability to discriminate between rice and non-rice crops.Clearly, analyzing both datasets in parallel may yield substantial benefit over only using one, especially given that the two are co-acquired in all standard Landsat image acquisitions.

Temporal Mixture Modeling
Figure 6 shows the result of the linear temporal mixture model inversions using sets of EMs from each characterization to map distinct spatiotemporal patterns of rice agriculture.The top row shows the full study area, and the bottom row shows the 24 x 24 km inset area within the white box.
The left column shows the result of a 3 EM temporal mixture model of vegetation fraction phenology to map timing of the 2016 rice crop using the Early and Late EMs obtained in Figure 5.In this map, green corresponds to Early, blue corresponds to Late, and red corresponds to complete absence of photosynthetic vegetation during the growing season.The saturation of each color corresponds to the similarity of each pixel to each of the EM time series.Greater saturation implies greater similarity to the corresponding EM (or binary mixtures for subtractive colors, while less saturation implies less similarity of the corresponding pixel to any of the temporal EMs.The latter is expected for pixels with phenologies not included in the model.Individual fields are clearly identifiable as either very early (pure green), very late (pure blue) or a mixture of the two (dark cyan).Intrafield heterogeneity in phenology is also clearly present in some cases, showing portions of individual fields growing faster than other portions.The potential utility of Landsat's ability to resolve intrafield spatiotemporal variability is discussed below.
The right column shows the result of a 3 EM temporal mixture model of 2016+2017 thermal phenology showing crops in both years (cyan), only 2016, magenta, or neither (dark yellow/orange).The unique temporal signature of rice thermal phenology allows it to be readily identified from other types of agriculture, which map as dull colors corresponding to combinations of EMs.Field-to-field variations in the similarity of each pixel time series to the temporal EMs are present, likely due to a combination of soil moisture during the fallow year and/or greenup phase during the cropped year.Some intrafield variability is also present, likely for similar reasons.
While the vegetation fraction TMM maps crop timing with high accuracy, it does not explicitly distinguish between rice and non-rice agriculture.As a result, non-rice agriculture maps in a variety of ways, potentially mimicking rice, depending on its phenological characteristics.Fortunately, the thermal image time series can explicitly capture the phenological signature of rice by leveraging the large early-season difference between cold, flooded, recently planted rice fields and hot, dry, recently planted non-rice fields.Interpreting the two of these phenology maps together therefore provides maximum information about both the location and timing of rice agriculture.

Near-Realtime Monitoring & Field Validation
To illustrate one potential application of this methodology, we present the results of a TMM generated in the middle of the current (2018) growing season.2018 presents an unusually challenging case for the model because only the first part of the growing season is available with only 2/3 as many usable images as the same period in 2016 and 2017.The prevalence of cloud cover precludes using many of the spring 2018 Landsat 8 thermal images and the smoke plumes from multiple severe wild fires contaminate many of the late summer images.Through the end of July, only 4 usable thermal images were collected, in comparison to 6 images through the end of July for each of the 2016 and 2017 growing seasons.
Figure 7 shows false color composite and thermal image time series for these 4 images.The false color composite images show a clear difference in crop timing between the rice growing areas on the eastern versus western portions of the valley, with the dividing line being approximately at the Sacramento River channel in the north and the Sutter Bypass in the south.The wide range of planting times in 2018, and the overall unusually late crop, is due to complex water management circumstances.Because of the relatively dry winter and associated low reservoir levels, uncertainty existed in early spring about the expected water allocations that rice farmers would receive.Heavy April rains then boosted allocations, but also forced farmers to wait for the clay-rich soils to dry before it was possible to bring tractors onto the fields to laser level them.The fields which had already been prepared could be flooded and planted on time, but those which had not were forced to plant significantly later than usual.A brief description of the situation is described by (Linquist, 2018) and expected impacts of this situation are broadly described by (Childs & Skorbiansky, 2018).
A TMM condenses the information from these 4 images of brightness temperature into a single map, show in Figure 8.In this model, red corresponds to grasslands and fallow fields, green corresponds to early rice, and blue corresponds to water and non-rice crops.Despite the limited data availability in 2018, and the fact that this map is produced mid-season with no data on senescence or harvest, the areas in which rice are grown are clearly identifiable.The broad east-west dichotomy in planting date described by (Linquist, 2018) is clearly evident in the discrepancy between bright and dull green map color.
This TMM is validated with over 1650 km of field transects collected between July 26 and 29 (white lines), including over 8500 field photos and 380 VNIR reflectance spectra.The locations of the field spectra are labeled in white and will be discussed further in Figure 10.Comparison of field photos to map output shows the model to produce <1% false negatives (i.e.fields planted as rice but not captured by the model).However, approximately 11% of fields in the model were identified as false positives (fields in which rice was not planted, but mapped by the model as rice).The reason for the prevalence of false positives in the 2018 model is straightforward, as only 1 image captures the rice in the phase in which it is dominated by standing water.Because the thermal signature of standing water can sometimes be similar to that transpiring vegetation, the thermal time series of any crop which was bare in mid-April but already fully green by June 1 would mimic that of the rice.This thermal mimicking could potentially be resolved by building a more complex model in which both vegetation fraction and brightness temperature are used.However, for the purposes of simplicity and illustration of the method outlined in this work, we defer a detailed examination of the strengths and weaknesses of more complex models to a future study.
VNIR field spectra were collected with an ASD HandHeld2 field spectrometer on July 29.Wildfires were active in the Coast Ranges to the west and north of the study area from mid July through late August, 2018.As a result, substantial atmospheric contamination was present during both during collection of the field spectra and satellite imaging immediately before, during, and after the field campaign.When compared to both same-day Sentinel-2 spectra (corrected with the ESA Sen2Cor algorithm) and Landsat 8 spectra from August 4, 2018 (corrected with the NASA/USGS LaSRC algorithm), the field spectra acquired were higher reflectance across all wavelengths than the satellite spectra.We interpret this as resulting from substantial indirect illumination by haze from the fires, view angle discrepancies due to the field spectrometer being held slightly off-nadir, and the unavoidable difference in measurement distance between the Spectralon white standard (used for calibration) and the distance from the rice canopy necessary to sample a representative mixture of foliage and exposed paddy water (or soil).Fortunately, however, the phenological progression depends on relative differences in reflectance with wavelength rather than absolute differences between sensors.In order to facilitate comparison of spectral shapes, the overall amplitudes of the field spectra were adjusted by simple scaling by constants in the range 0.5 to 0.8.
Figure 9 shows the corrected field spectra at full spectral resolution (thin gray) and the mean spectrum (thin black) for reference, as well as convolved with the broadband spectral response functions of the VNIR bands of the Landsat 8 OLI sensor (thick blue).These are shown in comparison to the actual observed Landsat 8 spectra from the August 4 Landsat 8 image (thick gold).In some cases (A-C), the agreement between the synthetic and observed Landsat spectra is remarkably good, even despite the significant atmospheric problems discussed above and the nearly 1 week between field and satellite data collection.In other cases, the spectral shape at visible wavelengths is distorted (D-F) or there is an overall loss of contrast (G-I) in the satellite spectra.These imperfections are likely due a combination of 1) changes in the state of the atmosphere in the 10s of seconds between field calibration to white standard and collection of reflectance spectra, and 2) the inability of the atmospheric correction model to capture the complex scattering and absorption effects of the smoke plumes.
In addition to the comparison with broadband observations, field spectra show subtle features which are not captured by multispectral instruments.Changes in the brightness of the visible red are likely due to Chlorophyll-b and other plant pigments.In addition, changes in the curvature at the base of the red edge (0.69 -0.71 m), as well as variations in slope at the top of the red edge (0.8-0.9 m) have the potential to be leveraged in future work for investigations into nutrient stress and phenological progression of plant functional traits.
Finally, sample field photos showing rice at a variety of stages are shown in Figure 10.In late July, 2018, rice in the Sacramento Valley could be observed in a range of stages of maturity (photos A-C).The prevalence of weeds widely varied among fields (D and E).Abundant waterfowl (F) could be observed feeding in the rice fields.Viewed from above, weedy rice fields (G) were clearly distinct from weed-free fields (H).Substantial standing water background was observed in nadir-looking photos.The importance of view angle can be clearly seen when comparing (H) and (I), as rice fields of this age are fully opaque when viewed obliquely but not when viewed from nadir.A wide diversity of other agricultural land covers were also present in the valley, including sunflowers (J), fruit and nut orchards (K), rice straw (L), bare soil (M), alfalfa (N) and tomatoes (O).

Discussion
The topological structure of the low order temporal feature spaces provides a detailed yet intuitive characterization of the thermal and vegetation phenology of the 2016 and 2017 rice crops in the Sacramento Valley.The clearly distinguished temporal endmembers spanning the feature space correspond to distinctive and easily verified land cover phenologies with straightforward interpretations in terms of vegetation abundance and brightness temperature at decameter scales.The clear distinction of the temporal endmembers, and the fact that they bound almost all the pixels in the feature space, allows them to be used as the basis for a linear mixture model that can be inverted to yield pixel scale similarity metrics for each temporal endmember.The resulting endmember fraction maps provide continuous field maps of both vegetation and thermal phenology of multiple rice crops as well as other land cover types.While the continuous field maps provide the most information on the spatiotemporal evolution of the rice crops, the clear presence of distinct clusters within the temporal feature spaces suggests that the continuous field maps could be discretized into thematic classifications if necessary.The spatiotemporal clustering highlights clear distinctions between the multiple phases of rice cropping, as well as from other land cover types with different phenologies.
The ability to rigorously map and monitor the state of rice agriculture with optical satellite imagery has a wide range of potential implications.Obvious implications include potential improvement in harvest date prediction for food security and commodities applications, weed management, and sub-field level crop health and nutrient stress.Less obvious implications include potential application to evapotranspiration (ET) studies, water use models, and pest management strategies based on spatial network approaches.We discuss each of these potential implications in turn.

Harvest prediction
In many regions of the world the prediction of rice harvest is hindered by lack of basic information.While providing accurate estimates of yield from satellite imagery can be difficult because of the range of factors that can contribute to the harvest, it is likely that fundamental variables such as planting date and total area under cultivation could be better constrained with the joint VSWIR-Thermal temporal mixture model approach presented here.In particular, TMMs based on single-year time series show the potential to map crop timing and dual-year time series show the potential to identify changes in the location and overall abundance of planted fields.This could improve one set of inputs to complex yield models which also rely on a significant amount of additional ancillary information, as well as work synergistically with other rice crop monitoring systems such as the methodologies proposed by (Boschetti & Stroppiana, 2009), (Torbick et al., 2011) and (Nelson et al., 2014).
2. Intra-field variability: weed and nutrient management Examination of decameter image time series from sensors such as Landsat and Sentinel reveals considerable intra-field heterogeneity in the time series of green-up and senescence.These within-field variations likely arise as a result of imperfections in field leveling, flooding, and drainage; variability in weed and nutrient management; and heterogeneity in planting density.The information present in these variations could be leveraged using active management strategies for precision agriculture.In conjunction with unmanned autonomous vehicles, localized application of fertilizer and herbicide could be targeted to areas of concern with a potentially significant reduction in environmental impact and cost to the grower.The growing prevalence of VNIR hyperspectral imaging sensors has the potential to contribute significantly to these problems by discriminating between individual absorption features such as those visible in the field spectra of Figure 10.
3. Pest management Some pests require proximity in order to transmit from one field to another.Large, interconnected regions growing similar crops may provide transmission opportunities for these pests.Transmission could be minimized by a basin-scale management approach based upon strategic disruption of spatial networks of agriculture.Accurate crop-specific maps could provide a useful input into a network-based approach to understanding agricultural landscape evolution such as explored in (Sousa & Small, 2016).
4. Evapotranspiration and water use ET has been shown to be remarkably consistent for rice fields throughout the growing season (Linquist et al., 2015).Indeed, the relative stability of the brightness temperature and area of illuminated photosynthetic vegetation of rice fields are a part of the reason why the joint VSWIRthermal mapping approach is so well suited to rice in particular.This is shown explicitly in the spatiotemporal composite Brightness Temperature versus Vegetation Fraction (TV) space for all 22 images in the 2016-2017 time series, shown in Figure 11.
The composite TV space shown in Figure 11 is particularly illuminating when viewed in relation to the 8 of its subsets shown earlier in Figure 4.The wide range of Tb values shown at nearzero Fv is primarily due to the substantial discrepancy between the relatively high heat capacity of water and the much lower heat capacity of dry soil (and synthetic surfaces like pavement).In contrast, the much more constrained range of values at high Fv values is likely a reflection of the thermal stability of dense, well-watered, respiring vegetation.The warm-colored clusters in this space represent groups of geographic pixels with similar TV properties.At low values of Fv these are likely bare fields with similar moisture contents and tilling practices.In contrast, the clustering at high Fv values likely corresponds to rice fields imaged at different dates, with differences in water temperature resulting from fluctuations in air temperature due to weather within the growing season.
At least two additional features of this plot suggest potential avenues for future work in application to ET, following the Triangle Method introduced by (Carlson et al., 1994) and reviewed by (Carlson, 2007).First, the remarkably straight edges suggest the Sacramento Valley may be a particularly good location for this approach to ET mapping.This could be due to a combination of the large field size, prevalent irrigation, and dominance of agriculture in the landscape.Second, a clear division between the growing season and the winter is also evident, suggesting the relationship between temperature and vegetation can be described as alternating between (at least) two physical regimes over the course of the year.Moving forward, spatiotemporal analysis of satellite-derived ET estimates has the potential to make progress toward all 5 of the requirements of the Path Forward outlined by (Fisher et al., 2017).Detailed comparison of ET estimates from this simple method of estimation versus more complex ET models present an attract avenue for future work.

View Angle & Flooding Presence
Because rice is an erectophile plant, small differences in off-nadir viewing geometry can produce significant differences in reflectance.Specifically, the background of paddy water or soil is most pronounced when nadir-looking and rapidly becomes less significant when viewed off-nadir.Sufficient data now exists to quantify this effect at various stages throughout the rice life cycle.Understanding the size of this effect has the potential to determine at what stages of development the flooding or drying of the substrate below the rice can be accurately determined, with potential implications for verification of alternative cropping practices such as alternate flooding & draining of fields.
6. Sharpening the thermal field Thermal imagery is generally collected at coarser spatial resolution than VSWIR imagery because the amount of radiance being measured at thermal wavelengths is much less than that being measured at VSWIR wavelengths.This allows for the use of higher resolution VSWIR imagery to sharpen lower resolution TIR imagery collected simultaneously (see (Zhan et al., 2013) for a review of thermal disaggregation approaches).The consistency of the relationship between vegetation fraction and brightness temperature could be leveraged to develop a method to sharpen thermal imagery to VSWIR spatial resolution specific to rice-dominated landscapes.While not particularly necessary in a landscape with large fields sizes like the Sacramento Valley, this could provide substantial benefit for ET monitoring in areas dominated by smallholder rice agriculture which are present throughout the developing world.

Integration of new data streams
While retrospective analysis of agriculture from satellite imagery must rely on whatever data has already been collected, in the future increased data availability has the potential to significantly advance the ability of optical remote sensing to constrain dates of phenological transition.As revisit times shorten and spatial resolution increases, the potential for near-realtime monitoring and discrimination between subtle signals improves considerably.In particular, incorporation of nowactive Sentinel-2 VSWIR imagery and soon-to-be-available Sentinel-3 thermal imagery will be an obvious next step in the near future.
In addition, every daytime satellite overpass has a corresponding nighttime overpass.This time is often used for optical Earth imaging satellites to conduct housekeeping procedures like data downlinking, cooling, and calibration, but recent improvements in sensor and platform design allow for more frequent nighttime imaging.The potential to leverage day/night pairs of thermal imagery for ET purposes is another avenue for future work.
Finally, future sensors providing hyperspectral capability in the VSWIR and multispectral capability in thermal wavelengths may provide significantly enhanced ability to discriminate between differences in specific absorptions due to leaf chemistry as well as subtle differences in emissivity and brightness temperature.The recent placement of the ECOsystem Spaceborne Thermal Radiometer Experiment (ECOSTRESS) sensor onto the International Space Station, and the planned launch of hyperspectral missions such as NASA's Surface Biology and Geology (SBG) and DLR's EnMap, have the potential to substantially enhance the ability of the scientific community to map and monitor rice agriculture.

Conclusion
The combined spatiotemporal analysis of decameter scale VSWIR and thermal imaging shows considerable promise for improvements to the mapping and monitoring of rice agriculture.Using the Sacramento Valley of California as a test area, we characterize the phenology of rice agriculture in the area in terms of both vegetation abundance and brightness temperature.Characterization of single-year time series readily yields information about early/late crop timing, while dual-year characterization provides information about year-to-year crop conversion.Surprisingly, parallel spatiotemporal analyses of coincident vegetation fraction (retrieved from VSWIR reflectance) and thermal imagery show considerably more clustering in the thermal image time series than in the vegetation fraction in both single-and dual-year time series.Investigation of the temporal feature spaces suggests that vegetation fraction time series data show particular promise for estimating crop timing, while brightness temperature appears particularly well suited to distinguish rice from other crops.Sample monitoring results from the 2018 growing season show considerable promise for near-realtime mapping of rice phenology using partial-season TMMs.Field validation confirms that mid-season monitoring models can provide an accurate upper bound estimate of the spatial extent and relative timing of the rice crop, even under conditions of relative data scarcity.These results could have potential implications for food security, precision agriculture, pest management, and ET.