Crop Monitoring Using Sentinel-1 Data: A Case Study from The Netherlands

Abstract

Agriculture is of huge economic significance in The Netherlands where the provision of real-time, reliable information on crop development is essential to support the transition towards precision agriculture. Optical imagery can provide invaluable insights into crop growth and development but is severely hampered by cloud cover. This case study in the Flevopolder illustrates the potential value of Sentinel-1 for monitoring five key crops in The Netherlands, namely sugar beet, potato, maize, wheat and English rye grass. Time series of radar backscatter from the European Space Agency’s Sentinel-1 Mission are analyzed and compared to ground measurements including phenological stage and height. Temporal variations in backscatter data reflect changes in water content and structure associated with phenological development. Emergence and closure dates are estimated from the backscatter time series and validated against a photo archive. Coherence data are compared to Normalized Difference Vegetation Index (NDVI) and ground data, illustrating that the sudden increase in coherence is a useful indicator of harvest. The results presented here demonstrate that Sentinel-1 data have significant potential value to monitor growth and development of key Dutch crops. Furthermore, the guaranteed availability of Sentinel-1 imagery in clouded conditions ensures the reliability of data to meet the monitoring needs of farmers, food producers and regulatory bodies.

1. Introduction

The Netherlands is the second largest exporter of food and agricultural products in the world, exporting 65 billion euro of agricultural produce annually. These exports represent 17.7% of total Dutch exports [1] and the sector is expected to grow further in the next 15 years [2]. The provision of timely, reliable, accurate and high resolution satellite remote sensing data are essential to facilitate a transition from parcel-level decision-making towards precision agriculture [3,4]. This transition towards precision agriculture is expected to yield increased productivity, lower environmental impact, transparent production and smarter production methods. Furthermore, the commercialization of remote sensing data processing and added-value product generation has the potential to become a valuable export commodity.

The current abundance of high-resolution optical data offers unprecedented opportunities for agricultural applications. Satellite remote sensing provides valuable information for many users from individual farmers to food producers, as well as national and international governmental agencies. Optical imagery can be used to predict yields, delineate management zones, support variable rate application and to monitor inter-field, intra-field and interannual variability [5,6,7,8].

The European Union aims to reform the Common Agricultural Policy (CAP) in 2020. The Integrated Administration and Control System (IACS) implements area-based, direct and rural development payments for farmers using high, and very high resolution satellite imagery. Given the abundance of open and free satellite data available through the European Commission’s Earth Observation programme (Copernicus), the EU is encouraging member states to change the IACS system to use Sentinel imagery for 100% monitoring of agricultural parcels. Therefore, many EU member states are actively investigating the feasibility of exploiting Sentinel data for monitoring crops [9].

The Netherlands Space Office (NSO) makes satellite imagery freely available for The Netherlands with a view to stimulating the development of innovative applications in agriculture and other domains [10]. In addition to freely available imagery from MODIS, Sentinel-2, etc., NSO acquires satellite imagery from commercial providers and makes it freely available for The Netherlands. However, the reliability of optical imagery in The Netherlands is severely undermined by cloud cover. Van der Wal et al. [11] used 20 years of data from a Royal Netherlands Meteorological Institute (KNMI). weather station data at Eelde (the Netherlands) to highlight the influence of cloud cover on the availability of optical imagery in The Netherlands. They showed that there is about a 20% chance of obtaining a clear (<2 Oktas) satellite acquisition during the growing season. While UAVs have been embraced by some growers, large-scale monitoring still depends on satellite remote sensing.

Low frequency (1–10 GHz) radar penetrates cloud cover and does not depend on solar illumination, so its potential to provide the timely and reliable observations essential for agricultural applications has been recognised for many decades [12]. More importantly, radar observations in this frequency range are sensitive to the moisture content, size, shape and orientation and roughness of the vegetation constituents as well as the moisture content, texture and roughness of the underlying soil [12,13,14,15,16]. In other words, low frequency radar is well suited to sensing soil moisture as well as water content and structural changes in agricultural crops. Several decades of ground-based experiments, and campaigns based on airborne and spaceborne sensors have demonstrated the suitability of low frequency radar for agricultural applications, particularly soil moisture monitoring, crop classification and crop monitoring [17]. The European Space Agency’s (ESA) Sentinel-1 Mission offers unique opportunities in terms of the operational use of radar observations for agricultural monitoring for two reasons: (1) its revisit time is unprecedented and (2) the imagery is freely distributed. The two satellites—the C-band Sentinel-1 Mission (Sentinel-1A and 1B) was launched in 2014 and 2015, respectively. They are in the same orbital plane, providing an average revisit time of two days above ${45}^{°}$ N/S, and achieve global exact repeat coverage every two weeks. The default acquisition mode over land in Europe is the Interferometric Wide-swath (IW) mode that provides both VV and VH data [18]. Temporal coverage varies across the globe, but combining ascending and descending tracks from both satellites yields observations every 1–2 days in The Netherlands.

The high temporal revisit of Sentinel-1 provides insight into the temporal variability of vegetation, which has been exploited in many recent classification studies [19,20,21,22,23]. The current study, however, is motivated by recent research highlighting the potential of Sentinel-1 for crop monitoring during the growing season. Veloso et al. [24] compared Sentinel-1 data to Normalized Difference Vegetation Index (NDVI) estimates from optical data and ground observations of precipitation, temperature, green area index and fresh biomass. They demonstrated that Sentinel-1 data, particularly the VH/VV ratio, could yield useful information on crop development. In particular, they highlighted the potential to distinguish between crops based on the temporal variation of backscatter. They also showed that, for barley and maize, the NDVI and VH/VV agreed well with Green Area Index (GAI) and fresh biomass. More recently, Vreugdenhil et al. [25] analyzed time series of Sentinel-1 backscatter and polarization ratio (VH/VV), and explored their relationship to vegetation water content (VWC), height, biomass and leaf area index (LAI). They showed that Random Forest modelling could be used to estimate VWC from Sentinel-1 imagery. The relationship between Sentinel-1 observables and VWC varied during the growing season, however, with variations in backscatter in some periods being dominated by structural changes. The studies, conducted in France and Austria, illustrate the potential to use Sentinel-1 imagery for crop monitoring.

The current study is focused on one of the most productive agricultural areas in The Netherlands to explore the value of Sentinel-1 in monitoring regionally important crops, namely sugar beet, potato, maize, winter wheat and English rye grass. By comparing Sentinel-1 imagery to hydrometeorological data and ground measurements of phenological stage, and vegetation height, it is shown that the time series of Sentinel-1 backscatter data reflects moisture and structural changes associated with phenological development of crops during the growing season in this region. It is shown that key dates of interest (emergence and closure dates) can be mapped using Sentinel-1 backscatter data. Finally, it is shown that, in addition to backscatter data, the coherence between consecutive Sentinel-1 images is influenced by the structural changes associated with (potato) haulming and harvest.

2. Data and Methods

2.1. Study Area

The study was conducted in the Flevopolder, a region of reclaimed land with an area of 970 km${}^2$ that was drained in 1957 with the main purpose of agriculture. Figure 1 shows the location of the Flevopolder in The Netherlands, as well as the spatial distribution of the five crops considered. Parcel boundaries and crop types were determined from the Basisregistratie Gewaspercelen (BRP) [26]. The land surface in this region is flat and lies ±3 m below sea level. Soil at the surface is clay overlaying a sand layer at about 2 m depth. Capillary rise from the shallow groundwater is a major stabilizing control on soil moisture. The average minimum temperature during winter is $-3.3{}^{°}$C, the average maximum temperature during summer is 22 ${}^{°}$C and the mean annual precipitation is 797 mm per year [27]. Weekly crop growth stage, crop height and soil moisture data were collected in 24 agricultural parcels on the Flevopolder. The crop types were sugar beet (five parcels), potato (four parcels), maize (five parcels), wheat (five parcels) and English rye grass (five parcels).

2.2. Hydrometeorological Ground Data

Hydrometeorological sensors were installed at a silage maize parcel at the Aeres Practijkcentrum in Dronten (52.53${}^{°}$ N, 5.66${}^{°}$ E) for continuous ground measurements. The growing season of this parcel lasted from May to the end of September 2017. A weather station was installed on the 24th of May, adjacent to the maize parcel. A Decagon ECH2O Rain Model ECRN-100 tipping bucket was mounted on the station; however, it failed and was replaced on 23 June. Precipitation data from the KNMI station at Dronten (Station Number 364) were used in the interim [28]. The weather station also included an Apogee SP-212 pyranometer (solar radiation), a Davis Cup anemometer (wind and gust speed, and wind direction), and a HOBO Temperature/RH Smart Sensor Model S-THB-M008 (temperature and relative humidity).

Radar backscatter, particularly at C-band, is sensitive to the presence of water on the canopy [29,30,31,32,33]. Therefore, three Decagon Dielectric Leaf Wetness Sensors were installed in the maize parcel in Dronten to measure the presence and duration of interception or dew on the leaf. To investigate the wetting and drying at different levels, the sensors were installed on three heights in the canopy: 27.5, 42.5, and 145 cm from the ground. Precipitation data were used to determine whether water droplets on the sensor were likely to be interception or dew.

Root zone soil moisture was monitored by installing five Decagon ECH2O EC-5 soil moisture sensors between two rows in the middle of the parcel, at five different depths: $-5,-10,-20,-40,-80$ cm. The accuracy of the factory generic calibration is approximately 3 to 4% [34]. Results from the sensor at 20 cm were excluded from the analysis because the sensor produced spurious values throughout the study period.

2.3. Ground Data at 24 Parcels

In each of the 24 parcels, two sampling locations were identified and marked at the start of the growing season. The sampling locations were located 20 m from the parcel boundary to avoid edge effects. These sites were visited approximately once per week, weather permitting.

Soil surface roughness parameters (root mean square (rms) height and correlation length (L)) were determined using digital photos of a 1 m long grid board during the bare soil period for the maize, sugar beet, potato and wheat parcels. At each sampling location, the grid board was moved twice, so that three adjacent meters were photographed. The images were digitally combined to produce a synthetic 3 m data series. The images were digitized manually, sampling every 1 cm along the soil surface. Surface soil moisture at each sampling point was measured using ML3 ThetaProbe Soil Moisture Sensors [35]. The parcel’s soil moisture was estimated as the average of eight measurements, four at each sampling location. In potato parcels, soil moisture measurements were made at the crest and trough of the soil mounds. Crop growth stage was determined by visual inspection, based on the BBCH scale [36,37]. Where the BBCH stage differed across parcels, the “average” value is indicated in the following sections. Crop height was measured, and a photo archive was generated to document the development stage.

2.4. Sentinel-1 Data

Sentinel-1 includes C-band Synthetic Aperture Radar (SAR) operating at a centre frequency of 5.405 GHz. The observation mode over (non-polar) land is the Interferometric Wide (IW) mode providing dual polarization (VV and VH) imagery over a 250 km swath at a 5 × 20 m spatial resolution.

Combining data from Sentinel-1A and 1B, the Flevopolder is covered by 4–5 tracks (See

Table 1. Sentinel-1 IW data available over the study area.

Relative Orbit	Pass	Local Time	Min. Inc. Angle (°)	Max. Inc. Angle (°)
37	DESC	06:49	38.9	41.9
161	ASC	18.32	44.7	46.1
88	ASC	18:24	36.6	40.4
15	ASC	18:15	30.0	31.5
110	DESC	06:58	30.0	33.7

). Track 110 was not considered because it only covers part of the domain. This potentially provides an average of 20–25 acquisitions per month. This study uses 60 images from Sentinel-1A and 1B from relative orbit 88 with a temporal resolution of six days between 4 January 2017 to 30 December 2017. An evening overpass time was chosen to minimize the influence of dew on the vegetation (see Section 3.1.2). Ground Range Detected (GRD) data were obtained via Google Earth Engine, and Single Look Complex (SLC) products for the coherence analysis were downloaded from the Sentinel Hub [38,39].

Time series of copolarized (VV) and cross-polarized (VH) normalized cross section (${\sigma }^{°}$) were extracted from the GRD products for all parcels shown in Figure 1. The number of parcels for each crop type is shown in

Table 2. Number of parcels for each of five crop types in the study area.

Crop Type	Number of Parcels
Maize	335
Potato	886
Sugar beet	763
Wheat	1048
English Rye Grass	1286

. Pre-processing steps including radiometric calibration, removal of thermal noise offset and orthorectification with radiometric correction for residual slope effects had been done by Google Earth Engine using a Sentinel-1 Toolbox [40].

Spatial multilooking was performed by averaging across all pixels within each parcel polygon. Hence, the radiometric resolution, or precision, for a single parcel depends on the parcel area. Typically, about 100 independent looks are available per hectare, resulting in a radiometric precision n of 0.5 dB for a parcel of 1 ha.

Interferometric coherence can be defined as the amplitude of the complex correlation coefficient between two radar acquisitions. Given two interferometric complex SAR images $s_1$ and $s_2$ (Sentinel-1 SLC products), the coherence amplitude is described as:

$$\gamma =\frac{\mid \langle s_1{s}_2^{*}\rangle \mid }{\sqrt{\langle s_1{s}_1^{*}\rangle {\langle s_2{s}_2^{*}\rangle }^{`}}},$$

(1)

where |..| indicates the absolute value, $\langle ..\rangle$ indicates an ensemble averaging operation, and ${}^{*}$ denotes the complex conjugate product.

A time series of interferometric coherence was obtained by processing Sentinel-1 IW mode images with a temporal base line of six days using the SNAP tool (version 6.0.0) [41,42], and the SNAP Graph Processing Tool (GPT) for batch processing. A flow chart of the processing steps for each pair of images is shown in Figure 2. The relevant sub-swath, bursts and polarization (VV) were extracted from the original Sentinel-1 SLC product. This was done separately, but with the same parameters, for each pair of images. The Sentinel-1 precise orbit files were then applied to both images separately. In the Back-Geocoding step, the orbital and DEM information was used to geometrically co-register the SLCs. The Enhanced Spectral Diversity (ESD) function was used to correct for both range and azimuth shifts over the burst overlap regions of the TOPS data. After the interferogram was generated, the TOPSAR-Deburst function was used to generate a spatially continuous product, the Flat Earth phase term was removed using the TopoPhaseRemoval function, and phase adaptive filtering (Goldstein Phase filtering) was applied to improve visualization and aid the subsequent unwrapping step. Finally, terrain-based geocoding was performed to produce the final product in geographical coordinates so that it could be combined with the parcel boundary information.

Detecting Emergence, Closure and Harvest Date

Emergence refers to the appearance of the first shoots and leaves above the ground. After emergence, VH backscatter is expected to increase as leaf development results in an increase in volumetric scattering. Therefore, emergence will be estimated as the time at which the slope of the VH backscatter becomes positive. Closure is defined as the moment when leaves from adjacent rows meet. It is a key input for crop growth models currently used for harvest forecasting and planning. For example, the SUikerbieten Model (SUMO) has been in operational use in The Netherlands since 1996 to provide bi-weekly forecasts of sugar beet yields, sugar content, and refined sugar yield [43]. These estimates are essential to provide timely advice to growers and to plan for the transport and processing of harvested beets to processing plants. Later closure dates generally correspond to lower yields [44]. A previous proof-of-concept study found the slowing down of leaf production around the closure date caused the VH backscatter to stabilize [27]. Hence, the closure date could be estimated as the time at which the slope in the VH backscatter becomes zero.

To estimate the emergence and closure dates, a curve is first fit to the time series of VH backscatter from 1 April to 1 July. This period is assumed to encompass the vegetative stages for sugar beet, potato and maize. The minimum and maximum points of this curve are determined and assumed to correspond to the emergence and closure dates. The estimated dates were validated visually using photos from the monitored parcels.

In Section 3.2, it will become clear that backscatter alone is not always sufficient to detect harvest in some crops. It will also be demonstrated that the disturbance to soil and vegetation properties and structure during harvest produce a change in coherence that can be used to detect harvest. In particular, it will be illustrated that, when the crop is harvested, the scene becomes more static and the coherence increases. A rule-based classification will be applied to identify harvest at a parcel level.

3. Results and Discussion

3.1. Hydrometeorological Data

3.1.1. Weather Station Data

Daily precipitation, average daily temperature and solar radiation, and 15-min relative humidity and wind speed are shown in Figure 3. Total precipitation in the period between 27 May and 19 September was 468.6 mm. The start of this period was characterized by two dry intervals of about two weeks with maximum daily average temperatures reaching up to 24 ${}^{°}$C. More than one third of all rain fell between 8 and 19 September. The highest peaks of solar radiation and temperature were found in June and July, while August was relatively moderate. In September, both solar radiation and temperature decreased substantially.

3.1.2. Interception and Dew

Figure 4 shows the presence of water on the sensors at three different heights in the canopy. Most periods during which water was detected coincide with rain events, and are indicated as ‘interception’. Water detection in periods without rain are indicated as ‘suspected dew’. The lower, middle, and upper sensor were wet 20.5, 15.4, and 25.7 percent of their measuring time, respectively. The upper sensor (Figure 4a) detected water more frequently than the lower sensors, particularly during periods of low-intensity rainfall, e.g., 22–25 July and 17–20 August. During high-intensity rainfall, e.g., on 8 September, some interception was detected on all sensors.

Table 3. Percentage of days with water on canopy at morning and evening overpass times.

Sensor	A.M. (%)	P.M. (%)
upper	45.6	19.1
middle	26.1	8.0
lower	35.2	11.4

shows the percentages of days with water on the canopy at times coinciding with satellite overpasses (Table 1). Moisture was detected on the canopy almost three times more often during the morning (descending) than the evening (ascending) passes. In the following analysis, only data from the ascending pass (Relative Orbit 88) will be used to avoid the confounding effect of vegetation surface water on backscatter.

3.1.3. Soil Moisture

Figure 5 shows the volumetric soil moisture ($\theta$) measured at four depths at the Aeres Practijkcentrum from 13 June to 20 September. In general, the season was characterized by well-watered growing conditions. From Figure 5, the soil moisture at 40 cm and 80 cm saturate close to 0.5 and 0.52 m${}^3$/m${}^3$, respectively. Porosity values higher than 0.5 are not uncommon in this clay soil, with higher values likely at depth [45]. The soil moisture values at 80 cm remain close to saturation throughout the study period due to precipitation and the availability of moisture from shallow groundwater. The frequent rain events in July and August allow the soil moisture at 5 cm and 10 cm to rise continuously, reaching near saturation by 9 September.

Figure 6 shows the surface soil moisture measured in each crop type throughout the growing season. Warm temperatures and limited precipitation in late-May and June are reflected in the relatively dry surface soil moisture values before 30 June. After that date, regular precipitation ensures that the surface soil moisture in all crop types remained relatively high. Temporal variations are similar across the different crop types due to the limited spatial extent and soil texture homogeneity of the study domain. Spatial variations are largest in the potato parcels due to differences between the crests (drier) and troughs (wetter) of the ridges. This is particularly clear during the drier period in late-May. Surface soil moisture in the grass parcels exhibits limited spatial variability. In the remaining crops, the variability is limited following large precipitation events or persistent precipitation (e.g., 30 June and 15 September).

3.2. Sentinel-1 Time Series

3.2.1. Maize

The time series of Sentinel-1 backscatter (VV, VH and VH/VV), as well as phenological development stages and measured crop height for maize parcels are shown in Figure 7. The crop height measurements indicate that the maize plants emerged around mid-May. In this region, maize is typically left to ripen in the field and harvested from mid-September onwards. Weekly phenological stage monitoring therefore covered the development from the early vegetative stages (BBCH = 12) to fully ripe (BBCH = 89).

From 1 January to mid-May, VH backscatter is around 12 dB lower than VV backscatter, but the temporal variations in both are quite similar. These are dominated by the surface soil moisture. Variations in space can be attributed to variations in roughness, soil texture and row orientation/geometry. The three sharp decreases in VV and VH on 22 January, 9 February and 12 December 2017 are the result of frost and snow cover, respectively [46] and are observed in the other crops too (Figure 8, Figure 9, Figure 10 and Figure 11). The increase of VV and VH before 5 March is caused by melting snow and precipitation on 26 February and 5 March. A heavy rain event on 16 April caused an abrupt increase in both VV and VH. It is noteworthy that these events have little or no influence on the VH/VV ratio. This reflects its reduced sensitivity to soil moisture variations compared to the individual VV and VH backscatter [25].

During the leaf development stage (BBCH 12 to 30) from 17 May to 8 June, the VH backscatter and VH/VV ratio increased by 7 and 3.5 dB, respectively. This is mainly due to an increase in volume scattering as the newly formed leaves unfold. From stem elongation (BBCH 31 to 39), through to the beginning of tassel emergence (BBCH 51), the maize height increased from 40 cm to 2 m. This rapid accumulation of above ground biomass, from 8 June to 10 July, resulted in increases in the VH backscatter and VH/VV ratio of 4.0 and 2.43 dB.

In general, VV backscatter also increases slightly during the leaf development and stem elongation phases. This can be explained by an increase in the double bounce between the vertical stalks and the soil [47,48]. While some studies report an increase in VV backscatter associated with maize biomass [48,49], others have reported no correlation or even an inverse relation between maize biomass and VV backscatter intensity [24,50]. The VV and VH backscatter start to plateau as the maize reaches its maximum height and the grain begins to develop (BBCH 71). This is consistent with previous studies showing that radar backscatter from maize parcels tends to saturate around the tasseling stage (BBCH 51) as the LAI reaches values of 2–3 [25,51,52].

During the late fruit development (BBCH 75) and ripening stages (BBCH 89) from 23 August until mid-September (BBCH 85 to 89), VH backscatter is relatively constant with an average of $-16.2$ dB. At the same time, the VV backscatter increased slightly from $-9.5$ to $-7.5$ dB. The combined effect is a decrease in VH/VV ratio consistent with the conclusion of Vreugdenhil et al. that the dynamics VH/VV reflect changes in vegetation water content [25].

Abrupt decreases in backscatter are observed in the individual parcels on various dates between mid-September and mid-October as maize is harvested. Harvesting reduces the backscatter, particularly in VH, and increases the spread in VV, VH, and VH/VV backscatter among the parcels. Harvesting the vegetation removes the dense maize canopy which reduced the sensitivity of both VV and VH backscatter to soil moisture signal during the growing season [24,49,53]. However, one reason why the harvest of the maize, sugar beet and potato is more difficult to detect in the autumn months is the elevated soil moisture. The expected reduction in backscatter from the vegetation is not observable against the high and variable backscatter from the moist soil as precipitation events occur throughout the usual harvest period. From mid-October onwards, the variability between the parcels remains relatively constant and fluctuations due to precipitation (25 of October and 12 of November 2017) and snowfall (12 of December 2017) are observed.

3.2.2. Potato

Figure 8 shows the Sentinel-1 and in situ data from January 2017 until January 2018 from 886 potato parcels in the Flevopolder. From November to March, the soil is bare so the dynamics in VV and VH backscatter are due to surface soil moisture variations and land management practices (tillage). Similar to the maize time series, the drops in backscatter associated with snow/frozen soil events on 22 January, 9 February and 12 December 2017 are observed in the potato parcels. The increase in standard deviation of backscatter in both polarizations from mid-March is due to the formation of ridges for potato cultivation. Differences in the orientation of the ridges with respect to the radar look direction results in variations in the backscatter from bare soil. In the potato parcels, VV backscatter in parcel DB is always higher than the others. This may be due to the influence of ridge orientation with respect to Sentinel-1 range direction. To illustrate this effect, Figure S6 (Supplementary Materials) shows the VV and VH backscatter from relative orbits 37 and 88 for field DB, which has a comparable incidence angle but is in a descending rather than an ascending track. Co-polarized backscatter (VV) is higher in orbit 88 from end of March (start of soil work) until first week of July (closure date). This suggests that differences in backscatter among the potato parcels is in part due to the influence of ridge geometry. After closure, the difference is reduced considerably as fractional cover of the underlying soil is high. In contrast, the backscatter values in VH polarization from the two orbits are quite similar, consistent with the idea that cross-polarized backscatter is less sensitive to ridge geometry. It is worth noting that the variability in VH backscatter is less than that in VV, suggesting that cross-polarized backscatter is less sensitive to directional scattering due to ridge geometry. Similar observations were made by Wegmuller et al. [54] in their analysis of roughness anisotropy in the same study area. This difference in sensitivity of co- and cross-polarized backscatter in different relative orbits should be explored as a potential source of information on the ridge geometry or tillage practices.

By the start of the phenological monitoring, the potatoes had already reached BBCH 61 (the beginning of flowering). The rapid increase in both VH and VV backscatter from 22 May onwards corresponds to the period of leaf formation and stem elongation. During this time, the potato plants enter a stage of exponential growth until the ground cover reaches a maximum in the middle of flowering stage at BBCH 63. Above ground biomass stabilizes during the flowering stage (BBCH 63, 27 June). The backscatter stays relatively constant, with only slight variations in VH due to heavy rain on 3 and 15 July and 8 August.

During the fruit develpment (BBCH 71 to 79) and ripening of fruit and seeds (BBCH 81 to 90) stages, the average height of the plants decreases. Variations in planting date, and field management practices (e.g., ridge structure and nutrient application) mean that the amount of above ground biomass can vary between parcels. In addition, the rate at which the plants reach the yellowing and brownish stages varies. Moreover, the reduction in vegetation water content increases sensitivity of the backscatter to the underlying soil layer. Therefore, variability in radar backscatter in VV and VH polarizations increase during this time. This is clear from the the increasing standard deviation in backscatter from 10 August to mid-September, particularly in the VH backscatter.

The abrupt decrease in VH backscatter from 13 to 28 September is related to haulming, the process of destroying the haulms (i.e., stalks or stems) prior to harvesting. In The Netherlands, chemical haulming is typically carried out 10–20 days before the intended harvest date to facilitate an equally ripening crop so that all tubers form resilient skin. Destruction of the above-ground plant reduces volumetric scattering and increases sensitivity to surface soil moisture and roughness. Hence, the subsequent increase in variability among the parcels. The rise in radar backscatter from 1 to 10 of October was caused by a rainy period from 29 September until 11 of October during which 67.4 mm of precipitation was observed at the KNMI station at Dronten [28].

The reduction in standard deviation in the VV, VH and VH/VV data after mid-October suggests that harvesting has been completed by that time. The main effect of mechanical harvesting is a perturbation in surface roughness after which the backscatter dynamics are primarily sensitive to soil moisture. This could explain the decrease in backscatter around this time. Similar to the case for maize, the combined influence of the changing roughness and increased in soil moisture level as a result of many rain events make it difficult to pinpoint the exact harvest time in the backscatter data. In Section 3.3.4, it will be shown that harvest is more readily identifiable using coherence data.

3.2.3. Sugar Beet

Figure 9 shows the Sentinel-1 and in situ data from January 2017 until January 2018 based on 763 sugar beet parcels in the Flevopolder. There are just three BBCH stages noted for sugar beets. The first BBCH stage measurement was during the rosette growth period and the last was when the crop cover was complete (BBCH 39), so the difference observed was primarily a change in ground cover from 60% to 90%. The later stages concern the development of harvestable beet root, and cannot be evaluated from visual inspection of the above-ground plant. Crop height was measured through to September 2017.

During March and April, the standard deviation across the parcels is high due to variations in sowing date, soil roughness and row orientation. The backscatter values start to converge in late-May during the leaf development stages. This coincides with a rapid increase in backscatter, particularly in VH as the development of new leaves increases the plant’s above-ground biomass. The sugar plant produces broad leaves in both horizontal and vertical direction which leads to an increase in volume scattering. The backscatter in both VV and VH stabilizes in mid-June and remains relatively constant until mid-August, with the exception of the influence of precipitation on 9 June and 8 August. Leaves cover 90% of the ground by 10 July (BBCH 39), so this stability in backscatter is due to stability in the structure and moisture content of the above-ground vegetation. There is limited sensitivity to increases in soil moisture due to smaller precipitation events during this period as the canopy is attenuating the signal from the soil completely. The influence of the larger events on 9 June and 8 August is clear in the VV and VH backscatter but is barely discernible in the VH/VV ratio. From August to mid-September, there is a noticeable decrease in the average VH/VV, and the standard deviation across parcels increases. This is due to the combined effect of very slight increases in backscatter, particularly in VV. This increase may be due to the increase in soil moisture (Figure 6), and an increased sensitivity to the surface soil moisture as the leaves lose moisture at the end of the summer. The standard deviation across parcels increases considerably from September to December coinciding with the harvest period. Sugar beet is harvested throughout this period. Hence, a sudden drop in the mean backscatter across all parcels is not expected. At an individual parcel scale, backscatter may decrease due to harvest. However, in Figure 9, these are difficult to distinguish from fluctuations due to increase in soil moisture level during the harvest period.

3.2.4. Winter Wheat

The time series of Sentinel-1 backscatter (VV, VH and VH/VV) and in situ measurements of crop height and phenological growth stages for 1048 wheat parcels are shown in Figure 10. The wheat was sown from mid-October to mid-November 2016 and harvested in August 2017. Four distinct periods can be identified in the backscatter time series.

From mid-March to mid-May, both VV and VH decrease due to the increase in biomass during the tillering and stem elongation period. VV backscatter is dominated by the direct ground and canopy contributions. The decrease in VV during the stem elongation is therefore attributed to the increasing attenuation due to the predominantly vertical structure of the wheat [24,25,55,56,57,58]. The decrease observed in VH is less severe than in VV. VH backscatter is influenced by double-bounce between the wheat stems and the ground [55,59] as well as a volume scattering component [48]. During the tillering and stem elongation periods, the increased attenuation of the double-bounce is countered by an increase in volume scattering from the new biomass [24,48]. The net effect is an increase in the VH/VV ratio, affirming its suitability as a measure of fresh biomass [25].

From mid-May to 1 June, during the booting and heading periods, VH backscatter remains stable while the VV backscatter starts to increase as the flag leaf opens and the inflorescences emerges. Brown et al. [55] found that the C-band VV backscatter during this period was primarily generated by the flag leaves and/or ears when they conducted detailed 3D imaging of wheat during this growth stage. The booting stage is therefore identifiable as the time at which the increase in backscatter due to stem elongation tapers off, but the VV backscatter increases due to the moisture content of the flag leaves and emerging ears.

From 1 June to 8 August, both VV and VH backscatter increase steadily by 8–9 dB as the winter wheat goes through the flowering, fruit development and ripening stages. This is due to the increased sensitivity to the ground contribution (both direct and double-bounce) as the vegetation loses internal moisture. Mattia et al. [56] emphasized the significance of the heading period as a turning point at which the radar backscatter becomes primarily sensitive to soil moisture rather than above ground biomass variations. The VH/VV increases by just 2 dB during this period.

Finally, the sudden drop in both VV and VH backscatter between 8 and 14 August is due to the harvest period. After this date, the standard deviation in backscatter between parcels increases. This is due to the diversity in post-harvest land management practices and the sensitivity to surface soil moisture and roughness which varies between parcels.

3.2.5. English Rye Grass

Figure 11 shows the Sentinel-1 and in situ data from January 2017 until January 2018 for 1286 English Rye Grass parcels in the study area. The snow/frozen soil events on 22 January, 9 February and 12 December 2017 result in dips in the VV and VH backscatter in the grass parcels. The heavy rain events on 9 June, 3 July, 8 August, 13 September and 7 October result in peaks of similar magnitudes in both VV and VH backscatter. The effect of precipitation/soil moisture on the VH/VV ratio is therefore limited, confirming that VH/VV is a good indicator of fresh biomass and vegetation water content [24,25].

The VH/VV ratio has a limited seasonal cycle due to the presence of grass cover throughout the year. A slight increase is observed in the summer months as the grass responds to increased solar radiation and warmer temperatures. This increase is due to a decrease in both VV and VH backscatter as the increasing fresh biomass attenuates the backscatter from the soil. The decrease in VV is larger than that observed in VH resulting in an increase in VH/VV during this time. Conversely, the lower temperatures from August onwards result in a decrease in biomass and an increase in sensitivity to soil moisture.

Dips in the grass height observations suggest that mowing occurred around 23 May, 15 June, 13 July, 17 August and 7 September. Slight decreases in VH backscatter are detected after these events. However, backscatter during this period is also influenced by precipitation and soil moisture variations, which makes it difficult to distinguish mowing events from soil moisture changes.

3.3. Mapping Key Dates

3.3.1. Emergence Date

The sugar beet emergence date was estimated as the point at which the emergence of the first leaves resulted in an increase in volumetric scattering from the vegetation, and hence an increase in VH backscatter. Several relationships were fit to the 6-day Sentinel-1 VH backscatter data from 1 April to 1 July, a period assumed to contain the vegetative stages. The highest goodness of fit (coefficient of determination) was obtained using a third order polynomial. The average goodness of fit across all parcels was $R^2=0.97$. For each parcel, the emergence date was determined as the minimum of the fitted curve. Results for the entire study area are mapped in Figure 12a. The estimated emergence date for 93% of parcels occurred between day 115 (25 April) and day 130 (10 May) with the remainder emerging either one week earlier or one week later.

The first available photos of the sugar beet parcels are shown in Figure 12b–e, along with the estimated emergence date. For the RBW parcel in Figure 12b, the first pair of leaves is visible (BBCH 12). Therefore, it is plausible that shoot emergence (BBCH 9) occurred four days previously. For the other three monitored parcels, the first photo is 10 to 14 days after the estimated emergence date. However, visual inspection of Figure 12c–e indicates that these sugar beets were already at BBCH 17 to 19 by the time there were first photographed. The consistency in estimated emergence date across all parcels, together with this limited photo validation, suggest that the estimated emergence dates are credible. Similar results obtained for potato and maize are provided in the Supplementary Materials.

3.3.2. Closure Date

The closure date was estimated as the date of the maximum of the third order polynomial fit to the data. Results for sugar beet are mapped in Figure 13a. The closure date for 93% percent of parcels (710 parcels) was estimated to be between Day 166 (15 June 2017) and Day 171 (20 June 2017). Photos taken from five monitored parcels around the estimated closure dates are shown in Figure 13b–k. For each of the five monitored parcels, it is clear that closure date occurred between the dates of the photos in the top and bottom rows. The leaves from adjacent rows are several centimetres apart in Figure 13b–d, but are clearly touching in Figure 13g,i. In Figure 13e, the leaves appear to be almost touching and the estimated closure date is two days later. In the next photo at that location (Figure 13j), the leaves are providing almost full ground cover, so it is plausible that closure is closer in time to the first photograph.

In practice, closure date is also determined by visual inspection and is an inherently subjective measure. Additional photos around the expected closure date, and a larger number of photos distributed across the study domain are essential to draw a more robust conclusion. However, the limited range of estimated closure dates, their timing and the agreement with photographs from the field suggest that this method provides a reasonable estimate of the closure date.

3.3.3. Sugar Beet Harvest Date

Figure 14a shows an example of parcel-averaged NDVI from 6 August to 30 November 2017 estimated from Sentinel-2 imagery. This parcel is chosen because it is known that it was harvested before 15 September. The parcel-averaged coherence between pairs of Sentinel-1 VV images is shown in blue for the same period. The coherence reported for a given date indicates the coherence between the image on that date and its predecessor.

The NDVI is greater than 0.7 until 3 September and decreases to around 0.1 by 15 September. This is consistent with the parcel being harvested between these dates. The coherence between images is low (<0.3) in the period when the NDVI is high. This is to be expected in a vegetated field where movement of leaves, fluctuations in water content and structure, and variations in soil moisture or leaf surface water result in phase differences between consecutive images. Conversely, in October and November, when the NDVI is close to zero, the coherence is generally above 0.6 as the scene is comparatively stable. Therefore, harvest is associated with a significant increase in coherence. Occasionally, the coherence decreases when soil moisture variations result in phase differences between consecutive images.

The high coherence on 19 September suggests that the phase differences were small between 13 and 19 September, i.e., that harvest occurred between 7 and 13 September.

Figure 14b shows a box plot of the parcel-averaged coherence for all 763 sugar beet parcels in the study area during the same period. The coherence is generally below 0.4 until 13 September, suggesting that all parcels are covered in vegetation. From 7 October onwards, the median slowly increases until it reaches almost 0.9 by 30 November and the interquartile range (IQR) gradually shifts towards higher values. From 19 September to 25 October, there are so few parcels with high coherence that they are considered outliers of the distribution. By 30 November, the entire IQR has a coherence greater than 0.6 suggesting that most of the parcels have been harvested.

Figure 14c shows the results from a simple rule-based classification formulated based on the results in Figure 14a,b. An increase in coherence is considered “significant” if it is greater than 1.5 times the standard deviation of coherence during the reference vegetated period (6 August–6 September). Coherence values of 0.5 are considered “high”. In Figure 14c, a parcel is considered harvested when the coherence undergoes a significant increase and reaches a “high” coherence value.

The red curve is the percentage of harvested parcel observed estimated by field agents/SuikerUnie. The two estimates converge after 18 November, and both suggest that 85% of parcels have been harvested by 30 November. Prior to that date, the rule-based approach estimates a higher percentage of harvested parcels. One contribution to the mismatch in Figure 14c is that the coherence is averaged to parcel level before the rule-based classification. This leads to ambiguities if the parcel is not all harvested at once, something quite common in this study area.

3.3.4. Potato Haulming and Harvest Date

Figure 15 shows a box plot similar to Figure 14b for all potato parcels in the study area. The median coherence is less than 0.3 throughout August during the ripening and senescence stages as the loss in moisture content and associated changes in structure cause phase differences between consecutive images. The median coherence between images on 13 and 19 September is 0.5, suggesting that a transition to a more stable scene occurred between 7 and 13 September. Recall from Figure 8 that there was a decrease in mean backscatter from mid-September and an increase in the standard deviation of backscatter across the potato parcels. Together, Figure 8 and Figure 15 suggest that haulming occurred during this period. Harvesting can typically be expected around 10–20 days later. However, frequent rainfall in September and October produced some decreases in coherence during these months. This precipitation also possibly delayed harvest to longer than 20 days after haulming. Nonetheless, increases in coherence generally increases during October. Coherence values greater than 0.6 from 31 October onwards suggest that most parcels were harvested before 25 October.

4. Conclusions

Results presented here illustrate the potential of Sentinel-1 SAR backscatter and interferometric coherence for crop monitoring and the detection of key dates for important agricultural crops in The Netherlands. Furthermore, the prevalence of rainy (and cloudy) conditions during this growing season underscore the need to use Sentinel-1 data by itself, or combined with optical data to provide timely and reliable information on crop monitoring in The Netherlands.

Time series analyses showed that, for each of the crops considered, structural and biomass changes associated with crop development influenced the backscatter throughout the season. Similar to the results of Vreugdenhil (2018) and Veloso (2017), the VH/VV ratio proved to be particularly useful as it reduces the influence of soil moisture. It is particularly sensitive to the increase in fresh biomass during the vegetative stages, and decreases during senescence as the vegetation water content decreases.

Key dates such as emergence and closure dates were estimated by fitting polynomials to the time series of backscatter, and were validated using field photos. Harvest detection proved to be more difficult in all crops apart from wheat. For sugar beet and potato, coherence data were needed to detect the harvest date.

This study used data from a single beam mode resulting in a temporal resolution of six days, replicating the more common Sentinel-1 acquisition strategy. Combining all beam modes would yield almost daily data in The Netherlands, improving the temporal resolution of emergence, closure and harvest date estimates. The sensitivity of interferometric coherence to structural and fresh biomass change is under-used at present, particularly in crop monitoring. The temporal density of Sentinel-1 data, particularly in Europe, offers an opportunity to capitalize on the synergy between backscatter and coherence data for agricultural monitoring.