In-season Mapping of Crop Type with Optical and X-band Sar Data: a Classification Tree Approach Using Synoptic Seasonal Features

Supplementary Information Figure S1. Distribution of fields sampled for the reference dataset derived from CUAA 2013 map and used for training and validation of the crop mapping approach implemented. Image background is mid-July EVI (16 July 2013) layer, represented in grey tones.


Introduction
The increasing demand for information on crop acreage for agricultural monitoring in support of private and public decision makers requires the production of reliable crop maps [1,2].Up-to-date information on agricultural land use is necessary for crop planning and management: e.g., for estimating biomass and yield, analyzing agronomic practices, assessing soil productivity, monitoring crop phenology and stress.Earth Observation (EO) techniques have been widely exploited in agriculture and agronomy for the advantages offered when compared to in situ and statistical surveys: frequency of acquisitions, synoptic view, and multi-dimensional content.Satellite remote sensing also constitutes the only source of consistent historical data for long-term analysis over large areas, e.g., for the identification of anomalous conditions in vegetation development driven by climatic variability [3].Moreover, EO data are available already during the growing season, whereas official statistics on crop acreages are often provided at the end of the season or later, thus being not useful for supporting in-season crop management.Since crop productivity quickly responds to unfavorable growing conditions, timeliness in delivering information on crop status is an important operational requirement [4,5], e.g., for mitigating the impact of crop stress conditions, especially for summer crops, which are prone to water stress in the dry summer months [6,7].
EO satellite data have been used for agricultural monitoring since the launch of the Landsat-1 system in the early 1970s and their potential for distinguishing different crops has been shown across various environmental conditions, and with many different data sources and methodologies, e.g., in [8,9].The Landsat archive constitutes the longest record of multi-spectral data available at medium spatial resolution, and has been used for crop mapping purposes at regional scale [10][11][12], using either spectral response and/or vegetation indices [13][14][15][16].The opening of the Landsat archives in 2008 has pushed forward the implementation of data analysis and image classification techniques based on multi-temporal features and time series analysis [17].Multi-temporal analysis techniques have been applied as well to coarser resolution data such as NOAA-AVHRR [18] and NASA-MODIS data [1,19], taking advantage of high revisit time for these sensors [20].Other satellite data too, with spectral and spatial features similar to Landsat, have been used for crop mapping achieving satisfactory results, e.g., IRS LISS data [21,22].
Data acquired by Synthetic Aperture Radar (SAR) active sensors have also been exploited for crop mapping and monitoring, especially during the last two decades.C-band data have been used for mapping rice [23][24][25][26], wheat [27], and multiple crops [28][29][30][31][32]. L-band data have been used too, although with generally poorer performance [33,34].More recently, with the launch of the TerraSAR-X and COSMO-SkyMed satellites, the use of X-band SAR data has largely expanded, mainly thanks to the higher spatial and temporal resolutions and theoretical flexibility of these platforms [35,36].Concerning X-band SAR data, different polarimetric configurations have been tested for crop mapping, from vertical-based, e.g., in [37], to horizontal-based polarization, e.g., in [38]; comparative studies using multiple polarizations have been carried out as well, e.g., in [39].However, to our knowledge no agreement has been reached so far on the best polarization configuration for crop mapping.
The integration of SAR and optical sensors for agricultural applications is a recent topic, aiming at reducing the impact of optical and SAR specific limitations (i.e., dependence on solar and clear sky conditions for the former, and on signal noise and stability for the latter).SAR and optical data integration takes advantage of their complementarity in terms of sensitivity to vegetation and soil characteristics [40][41][42]: plant biomass, soil moisture and surface texture for SAR, spectral response of canopy-background system and photosynthetic features for optical sensors [31,32,43].Successful integration examples are the work of Michelson et al. [44], and more recently some large scale studies using multi-source data from RADARSAT-1, Envisat ASAR, SPOT and Landsat sensors [29,37,45,46].The algorithms used for land cover mapping with both optical and SAR data range from maximum likelihood [46] and neural network ensembles [46,47] for crop classification, to maximum likelihood with iterated conditional modes [48] and Random Forest [49] for regional and urban land cover targets, up to fuzzy scores aggregation for burned area mapping [50].With the advent of new generation satellites, e.g., Landsat 8, Sentinel-2, WorldView-3, as well as COSMO-SkyMed, TerraSAR-X, RADARSAT-2, and Sentinel-1, crop mapping applications can be more timely and reliable, in particular during the early growth stages, and the operational use of such techniques will be promoted.
This paper describes a classification tree approach for in-season crop mapping, which exploits features derived from multi-temporal optical, Landsat 8 Operational Land Imager (OLI), and X-band SAR, COSMO-SkyMed, data for producing reliable in-season crop maps over temperate areas.The proposed approach builds on the analysis of separability between different crops to identify the best performing proxy combinations and synoptic seasonal features as crop mapping inputs.Classification tree approaches can handle input features of different types and derived from different sources, and are directly interpretable and adaptable, being structured as a set of simple rules [51][52][53].These characteristics make classification tree approaches both efficient and effective, especially for delivering mapping algorithms which are to be used in operational contexts.
The main objective is to define a classification tree approach for producing a crop map early in the summer season, i.e., around mid-July [54], to support agricultural management in Northern Italy.Spectral features for the winter and summer crop seasons (named synoptic seasonal features) are extracted from the temporal profiles of a set of proxies derived from optical and SAR data.Different proxies were used, sensitive to vegetation and soil conditions and able to characterize the dynamics of different crop types throughout the growing season: Spectral Indices (SIs) from optical data and/or backscatter and interferometric coherence information from X-band SAR data.Most of the literature using multi-temporal information for crop mapping focus on the use of temporal profiles of spectral indices derived from optical data [55][56][57][58], and only recently some authors successfully included SAR backscatter profiles for rice mapping [59,60].The novelty of our work are the use of synoptic seasonal features integrating optical and SAR data, and the delivery of crop type mapping already during the early stage of growth; eventually, the transferability to a growing season different from the one used for developing the approach is tested.

Study Area
The study area is located in Lombardy region, northern Italy (Figure 1); it lies south of Milan and it is bordered by the Po river.The area, covering around 1100 Km 2 , is mostly flat and intensively cultivated.Climate is continental, with annual temperature changes between January and July up to 20 °C and average precipitations of 850 mm/year.The most economically valuable crop types cultivated in Lombardy are (in per cent of total cropland area): maize (38.5%), temporary and permanent meadows used for forages (34.1%), rice (10.1%), winter cereals (wheat and barley, 7.5%), soybean and other legumes (2.7%), and vegetables (1.4%) [61].The two major crop seasons run from October to June and from April to October for winter and summer crops, respectively.Barley is the prevailing winter crop, typically flowering in April-May and harvested in May-June.The main summer crops, covering most of the cropland area and consuming most of irrigation resources, are maize and rice.Maize, sown between April and early May, reaches the peak of the vegetative phase in July and is harvested from the end of August through September.Maize is often also sown in double cropping practices for fresh forages or silage, after meadows or winter cereals, in integrated crop-livestock systems.Rice is usually sown later than maize, from the second half of April to late May, reaching the flowering stage in late July or early August, and it is harvested from late September onwards.

Satellite Data
The remotely sensed dataset is composed of 13 COSMO-SkyMed (CSK) and 14 Landsat 8 OLI images covering the spring/summer seasons of the years 2013 (18 April-23 July) and 2014 (05 April-19 July), as shown in Table 1.CSK data have been consistently acquired by the same satellite of the COSMO SkyMed constellation (CSK 1) in single polarization (HH) and interferometric mode.CSK dataset of both years is acquired from the same flight track and with constant orbital configuration.This configuration allowed us to exploit the dataset not only for X-band intensity calculation, but also for extracting repeat-pass interferometric coherence, with absolute values of perpendicular baselines ranging from 190 to 479 m for 2013 CSK dataset, and from 363 to 1020 m for 2014 CSK dataset.The CSK product was acquired as StripMap HIMAGE in Single Look Complex (SLC) format, HH polarization, descending pass, with look angle of 24.1 degrees.The OLI dataset was collected by taking advantage of two overlapping WRS-2 paths (193 and 194, rows 28-29) that guarantees a theoretical revisiting frequency of 7-9 days over the study area.Landsat 8 OLI scenes were retained only when overall cloud cover was less than 10%.

Reference Data
Three thematic levels of crop types were considered for covering the variability of crops cultivated in the study area: a detailed level (level 2-Lev2), an intermediate level (level 1-Lev1), and a generic level (level 0-Lev0), as summarized in Table 2. Lev2 is composed of 12 classes: early maize (Ma1), medium maize (Ma2), late maize (Ma3), early rice (R1), late rice (R2), dry seeded rice (R3), early soybean (Sb1), late soybean (Sb2), winter crop (WC), double crop (i.e., winter crop followed by a summer crop; DC), forages (i.e., permanent and temporary meadows used as fodder; Fo), and forestry-woodland (either natural or man-made; F-W).Lev1 groups the subclasses of maize, rice and soybean crops into three mono-type classes: maize (Ma), rice (R), and soybean (Sb), thus delivering a total of 7 crop cover classes (Ma, R, Sb, WC, DC, Fo, F-W).Indeed, Lev1 classes represent the target crop types for operational use of the early crop map, but since some Lev1 classes showed multimodal SIs temporal profiles due to different sowing dates of various cultivars (e.g., long and short cycle maize), we split some summer crop classes into different sub-classes in order to effectively calculating class separability, thus composing Lev2 classes.Finally, Lev0 was derived by further aggregating summer crop types (Ma, R, Sb) into a unique summer crop (SC) class, leading to a total of 5 land cover classes (SC, WC, DC, Fo, F-W).For building the reference dataset to be used for crop type classification implementation, a set of 570 crop fields (almost 9000 pixels), belonging to 12 different crop classes, have been identified for the 2013 spring-summer season based on the Annual Agricultural Land Use Map of Lombardy region (Carta Uso Agricolo Annuale, CUAA); this map is produced and distributed by the "Ente Regionale per i Servizi all'Agricoltura e alle Foreste" (i.e., the regional agency for agriculture and forest services of the Lombardy regional government) on an annual basis and relies on farmers' declarations as the primary source of information [62].
The crop categories of the CUAA legend are not consistent with crop classes defined on the basis of spectral response from remotely sensed data.For example, the CUAA crop category "maize" includes both single and double crop cultivations, which are characterized by different temporal profiles of the Spectral Indices derived from OLI data, thus leading to two distinct classes.Furthermore, the CUAA category "forages" includes all crops cultivated for animal consumption, i.e., some winter cereals, fodder and managed grasslands (alfalfa and similar), which are grown and mowed several times per season.Finally, no official figures are provided for the accuracy of the CUAA product.Therefore, the CUAA 2013 map was not used as direct reference data source, but it was used as base information for extracting sample fields, which have been confirmed by visual assessment of high resolution satellite photos covering part of the study area (acquired on 22 March and 10 August 2013, from Google Earth), in situ observations for a limited number of fields (survey along main roads using camera and GPS) and interpretation of multi-temporal profiles of optical scenes (i.e., for extracting Lev2 subclasses from Lev1 by assessing season timing of EVI peak, and for delineating double cropped fields).The reference dataset for 2013 growing season is described in Table 2, including the number of fields sampled for each class and level.The sampling was done on the basis of random selection of spatially distributed points from CUAA 2013 map within each crop class, followed by checking for correct class labelling consistently with semantic crop classes included into our target legend.Finally, a subdivision of the reference set was made on a per-field basis, with 2/3 of the fields used for training and 1/3 used for validation (see Table 2).This procedure allowed the attribution of crop type labels consistent with target crop classes in 2013 season, and it is used together with the satellite development set.For 2014 season, a validation set was constructed starting from CUAA 2014 map and using the same checking procedure described for 2013 reference set.2014 validation set is composed of a total of 3759 pixels, with class cardinality ranging from 289 to 753 pixels, and it is consistently used for assessing mapping accuracy derived from satellite transferability set.
The satellite and reference data have been divided into two separate datasets: (i) a development set, used for best input feature selection, training of the classification algorithm and accuracy assessment (i.e., for developing the crop mapping approach), made of the satellite data from the year 2013 (7 CSK, 6 OLI) and the training and validation samples extracted and checked from CUAA 2013 (Table 2) and (ii) a transferability set, used for validation of the crop mapping approach implemented for a different growing season (i.e., testing the transferability of the approach), composed of the satellite data acquired in 2014 (6 CSK, 8 OLI) and the validation sample extracted and checked from CUAA 2014.

NIR(b5OLI)
RED(b4OLI) RED(b4OLI) EVI was developed as an enhanced version of NDVI, including correction for background and atmospheric disturbances; the spectral bands of near infrared (NIR) and visible red (RED) are supplemented by information from the visible blue (BLUE) band, by using optimal weighting [66].EVI provides information about vegetation vigor, linked to biomass and fractional cover.As a complement to EVI, we included NDFI, which is an index developed for the detection of surface water in flooded rice areas; NDFI, originally introduced as NDFI2 or NDSIB2B7 [67], is the normalized difference of the RED and the short wave (SWIR, 2.1-2.2µm) spectral bands.NDFI provides information about soil moisture and flooding conditions, especially relevant for paddy rice fields.The ratio of the RED and GREEN reflectance values (RGRI) was included due to its sensitivity to photosynthetic efficiency and leaf pigments [68].
Figure 2 shows EVI (a) and NDFI (b) multi-temporal profiles extracted from the 2014 OLI dataset (covering the whole growing season: 16 March-23 October), together with acquisition dates of CSK and OLI for the development dataset (2013); the grey bar highlights the temporal range adopted as early crop map production deadline (mid-July).The profiles qualitatively well describe the cycles of the major winter and summer crops of the study area, thus promoting the SIs selected as suitable candidates for monitoring crop dynamics during the season.
CSK images were pre-processed with MAPscape-RICE software [59] for (i) mosaicking single frames into slant range continuous strips and (ii) co-registration of images using orbital information and automatic spatial matching based on cross-correlation.
Two different proxies of crop conditions were derived from CSK data: X-band backscattering coefficient sigma nought (σ°), related to plant biomass and morphology and soil conditions (moisture and roughness), already used for crop mapping by Fontanelli et al. [69], and repeat-pass interferometric coherence (γ), related to the evolution of surface scattering properties of canopy/background system (plant height and density) during the season [70,71].The σ° was derived through three processing steps: (a) multi-temporal speckle filtering according to the approach developed by De Grandi et al. [59,60,72], to balance differences in reflectivity between images at different times, (b) geocoding and radiometric calibration, using a Digital Elevation Model (SRTM DEM, at 90m equivalent ground resolution) and the radar equation, in which scattering area, antenna gain patterns and range spread loss were considered, and finally (c) normalization on local incidence angle, according to the cosine law.The interferometric coherence γ maps were produced using the complex data of image pairs of consecutive acquisitions [73], with a temporal baseline of 16 days (32 days for the 08 June-10 July 2014 pair).The multi-temporal σ° and γ maps were finally geocoded to UTM 32N WGS84 and spatially resampled to the same spatial resolution of L8 OLI images (30 m), by aggregating the average value of 10 × 10 pixels at the original resolution of 3 m.

Multi-Temporal Proxies Test
In order to capture the distinct seasonal patterns of different crops, we divided the multi-temporal dataset into two periods: (i) from April to the beginning of June, when winter crops are harvested; and (ii) from May-June, when summer crops emergence starts, to mid-July.Synoptic seasonal features over the two periods were computed for EVI, RGRI, σ° and γ.For NDFI, only pre-summer features were derived, being it related to flooding in rice cultivation.Seasonal proxies used for the development set (2013) are summarized in Table 3.Each seasonal proxy is made a series of values corresponding to dates falling into the Day Of Year (DOY) range representative of winter crop and early summer crop seasons.
For each proxy, the following synoptic seasonal features were extracted for the winter and summer periods: minimum value over the seasonal range (min), maximum value over the seasonal range (max), mean value over the seasonal range (ave), standard deviation over the seasonal range (std) (Figure 3), and the asymmetry index of the seasonal proxy scores histogram, or skewness (ske), not included in Figure 3.

In-Season Crop Type Classification
The classification scheme and the input features were selected to satisfy as much as possible crop mapping pre-operational requirements, thus (i) providing a product with high thematic mapping accuracy; (ii) being transferable to different years and (iii) building on rules both simple and interpretable, even to those who are non-experts in pattern recognition and remote sensing.Given these constraints, we implemented a rule-based classification tree, which grants both flexibility and robustness, and support the use of multi-source data [51][52][53].The scheme is implemented using the Classification and regression Tree (CT) algorithm, in the extension of J48 java routine, programmed in WEKA 3.6 [74].J48 CT routine is shaped on C4.5 [75] and consists of a recursive algorithm, that generates a classification tree through iterative partitioning of the feature space by using the information gain (computed from the entropy function) of each attribute for a set of cases [76].Each node in a tree is associated to a set of two or more cases.The attribute with the highest information gain is selected for each node, and the optimal threshold for continuous attributes is computed.For avoiding too complicated tree structures and over-fitting, embedded pruning capabilities are implemented into J48 CT algorithm, according to a given confidence level.For our J48 CT implementation we allowed only binary splits for each node, and used online pruning, with confidence factor of 0.25 and sub-tree raising option.For dealing with possible over-fitting and minimizing the tree size, we set the minimum number of classified instances per each final node equal to 200, which is approximately half of the size of the smallest Lev1 crop class in the training dataset (forages, 447 pixels, see Table 2).CT outputs a set of hierarchical rules with optimized decision boundaries in form of thresholds, which can be straightforwardly implemented for image classification.An additional output of CT is the assessment of class attribution error for each tree node, which is a useful metric for ex-post tree re-structuring, in case of high accumulation of misclassified instances in some branches.
The performance of the CT schemes implemented was tested with different combination of proxies and features, by computing the confusion matrix and derived accuracy metrics [77].During this phase, results achieved with the CT approach were compared to the ones achieved by Random Forest (RF) classification [78], which is currently acknowledged as the upper limit reachable using state of the art classification tree algorithms and multi-source data [79][80][81].

Selection of Input Features
Descriptive statistics were extracted for each Lev2 class and different synoptic seasonal features combinations: (i) σ° features; (ii) σ° and γ features; (iii) EVI features; (iv) EVI and NDFI features; (v) EVI, NDFI and RGRI features; (vi) EVI, NDFI, RGRI, and σ° features; (vii) EVI, NDFI, RGRI, σ° and γ features.Lev2 class-by-class separability was computed for these combinations using the Jeffries-Matusita Distance (J-MDIST) [82], and aggregated as per-class separability by averaging all possible pairings comprising a specific class.As a rule of the thumb, good separability is generally set at J-MDIST higher than 1.9 [83].Lev2 classes used for extracting class-by-class separability were further summarized into average class separability scores at Lev1.The best performing combinations of seasonal proxies were then selected by maximizing separability, with J-MDIST~2.Following the selection of best proxy combinations, the feasibility of reducing the synoptic seasonal features input set was assessed based on the overall performance of preliminary crop classification tests using CT and RF, with three different sets of features: (i) minimum, maximum, mean, standard deviation, and skewness (min-max-ave-std-ske); (ii) minimum, maximum, and mean (min-max-ave); and(iii) minimum and maximum (min-max).The rationale behind this choice is to assess the effect in terms of mapping accuracy, when discarding features mostly affected by year-to-year variability of satellite image frequency: i.e., we first excluded standard deviation and skewness, and then the average, which in case of few scenes available can be more biased than extreme values (min-max).
At the end of the test phase, we retained the input set (i.e., the combinations of proxies and synoptic seasonal features) that granted as much reduction as possible in the number of features, granting the lowest difference of overall accuracy between CT and RF results, at Lev1.

Validation
The CT scheme derived for the selected input set was applied and validated over 2013 and 2014 seasonal data, development and transferability sets, by calculating confusion matrices and derived metrics: Overall accuracy (OA), Kappa coefficient of agreement (κ) and per-class Commission (CE) and Omission (OE) Errors [77].Validation was carried out on samples independent from the training set, with per-class cardinality either proportional to the one of the training set (case t), or of actual cropland cover calculated from CUAA reference information (case r).Crop mapping performance was assessed at two thematic levels: Lev1, and Lev0.
An overview of the whole methodological approach described in Section 3 is given in Figure 4.

Selection of Synoptic Seasonal Features
Since the separability scores for Lev2 classes were generally not high enough for expecting good performance in early mapping, Lev1 (i.e., 7 crop classes) was chosen as the target level for our crop mapping approach.Class separability scores aggregated at Lev1 (Table 4) show some interesting response to the different input proxy combinations (derived from optical and/or X-band SAR data).
Using only X-band SAR information, J-MDIST increases consistently for all classes when adding γ to σ°: from +0.041 for better separated classes (rice and forages), to +0.150 and up to +0.204 for crop classes difficult to separate using only σ° (maize, soybean, and double crop).Using optical proxies, a slighter but consistent increment in J-MDIST (0.037 to 0.053) is observed by adding NDFI to EVI for summer crops (maize, rice, and soybean).Further addition of RGRI brings a small, yet consistent, increment in J-MDIST (0.009 to 0.023) over summer crops.As regards separability achieved by using optical and SAR proxies together, an increment is observed adding σ° to the full optical feature set; since separability is already close to the maximum value (J-MDIST~2) the increment is lower for maize and soybean classes (+0.007 to +0.010).No significant increment (+0.000 to +0.002) is granted by further adding γ.
Based on separability scores shown in Table 4, we kept only the best performing combinations in terms of overall separability (minimum JMDIST > 1.98): EVI+NDFI+RGRI (ERN, J-MDIST > 1.983) and EVI+NDFI+RGRI+σ° (ERN+s, J-MDIST > 1.997).Since the scores achieved using EVI+NDFI+RGRI+σ° and EVI+NDFI+RGRI+σ°+γ synoptic features are not significantly different, we decided to discard the combination including γ to keep the feature set as simple as possible.These two combinations were used as input for the classification approach development.Table 5 shows OA of preliminary crop classification tests using synoptic seasonal features extracted from ERN and ERN+s combinations at Lev1 and Lev0 as input for CT, with RF scores as reference.Results achieved with the complete features set (min-max-ave-std-ske) show that: (i) at Lev1, OA achieved with CT increases from 85.3% (ERN) to 89.0% (ERN+s), with a gap towards RF of 5.8%-8.9%;(ii) at Lev0, very high OA is scored by CT (96.7-97.8%),reducing the gap towards RF to 1.3%-2.5%.Reducing the input features to min-max-ave set did not produce a sensible decrement in OA, with maximum decrement in OA of −0.8% across different proxy combinations and thematic levels, while increments in OA for CT are observed using ERN: +1.2% (at Lev0) and +1.7% (at Lev1).Further reducing input features to min-max did not significantly change OA at Lev0, but a decrement up to 2.1% was observed for Lev1.
Figure 5 shows per-class omission and commission errors (OE and CE) at Lev1 for CT and RF.When the min-max-ave set is used, no significant increase of per-class errors is observed, compared to the use of a full set of features (min-max-ave-std-ske; Figure 5b-e).Instead, the use of min-max-ave set and CT fed with optical only features (ERN) contributes to a reduction of OE for forages (25%, Figure 5b), and of CE for forestry-woodland (9%, Figure 5e).When the input feature set is further reduced to min-max (Figure 5c-f), higher errors for double crop (+11% OE, +5% CE), and some overestimation of winter crop (+8% CE) are observed.The best performing synoptic seasonal feature set was therefore identified as min-max-ave, which was therefore selected as best option input for early crop mapping.(a) Table 5. Overall Accuracy (OA) for CT at Lev1 and Lev0, using ERN and ERN+s proxy combinations and min-max-ave-std-ske, min-max-ave and min-max seasonal features.ERN = EVI+NDFI+RGRI; ERN = EVI+NDFI+RGRI+σ° = ERN+s.Results achievable using RF are shown as reference.

In-Season Crop Type Classification
Figure 6 shows the two CT schemes implemented using the J48 algorithm for to the 2013 training dataset with the min-max-ave synoptic seasonal feature set, applied to the ERN (CTmin-max-ave(ERN), Figure 6a) and ERN+s (CTmin-max-ave(ERN+s), Figure 6b).
The CT scheme developed using only optical features (CTmin-max-ave(ERN)) shows the first split of tree nodes based on winter season RGRI maximum (RGRIW max ≥ 1.361), separating pixels which are vegetated in May-June (right side) from non-vegetated ones (left side), thus distinguishing summer crops from winter crops, double crops and other agricultural land cover.The tree branches on the right side of Figure 6a further split winter and double crops from non-sown land covers (forages, forestry-woodland) based on two combinations of EVI mean in early summer (EVIS ave ), since the latter classes show high green fractional cover already during spring, when summer crops are not yet sown.All these nodes show crop type attribution errors lower than 2.4%, except for winter crop class (9.6% cumulated node error).Left branches in Figure 6a are populated by summer crop classes, characterized by low RGRIW max .Below these branches, the main splits are based on NDFI maximum (>0.089) to separate flooded rice fields (very accurate, with node error of 0.1%), and EVI maximum in winter-spring season (EVIW max > 0.298) to identify early-cycle crops (mostly maize).Crop type detection in lower level branches are due to a combination of EVI and RGRI features from both winter-spring and early summer features, and are meant to separate a mixture of rice, maize and soybean; these branches are characterized by cumulative node error above 20% (high misclassification rate).The classification scheme developed using integrated optical and X-band SAR features (CTmin-max-ave(ERN+s)) shows a main split for RGRIW max ≥ 1.361, consistently with CTmin-max-ave(ERN).The branches on the right side highlight the contribution of X-band minimum backscatter in early summer for the discrimination of forages from forestry-woodland (σ°S min > 0.055).In the left branches, the major rice class (flooded rice) is first identified based on maximum NDFI (as in CTmin-max-ave(ERN)) while maize or rice pixels showing very high X-band backscatter in early summer (σ°S max > 0.682) are separated from a mixture of rice, maize and soybean.This mixture of summer crops is further untangled by a combination of winter season EVI, mean summer backscatter and RGRI peak scores.As previously noted for Figure 6a scheme, these leftmost branches are characterized by high class attribution errors (>20%).Figure 7 shows the in-season early crop maps produced using CTmin-max-ave(ERN+s) scheme applied to the 2013 and 2014 datasets.

Validation
Table 6 summarizes accuracy metrics (OA and κ) for the crop maps shown in Figure 7 and computed for two different validation sets, with crop class cardinality proportional to either the training set (case t), or the actual cropland coverage of the area calculated from CUAA (case r).For 2013, the accuracy scores retrieved using the two sets are highly consistent, with case r giving slightly higher scores.OA achieved at Lev1 with 2013 ERN+s dataset, peaking at 91.8% over case r set (κ = 0.897), are 1.7% higher than when using only optical data (ERN).At Lev0, the performance of ERN and ERN+s data are nearly the same, with OA (κ) around 98% (0.960) for both case r and case t validation sets.Classification performance over the transferability set (2014) shows very good accuracy scores at Lev0, yet lower than for 2013: OA = 85.6% (κ = 0.740) using optical features and OA = 92.4% (κ = 0.861) using optical and σ° features.An increment of 6.8% in OA is achieved at Lev0 by integrating σ° for 2014, while for 2013 no enhancement was observed.At Lev1, less consistent results are observed: using ERN input an OA = 66.9% was reached, 23.2% lower than 2013 results, while when ERN+s input are used, a rebound of +19.7% in OA is achieved, reaching 86.6% (κ = 0.826).This result highlights the significant contribution of X-band SAR backscattering in terms of transferability of the approach, i.e., when the classification scheme is applied to a seasonal dataset different from the one used for algorithm development.The additional information brought by X-band SAR increases the robustness of the mapping approach at Lev1.Per-class errors were analyzed at Lev1, providing some insights into the disaggregation of global accuracy results (Figure 8 and Table 7).Using the 2013 validation dataset (case r), depicted in light and dark green bars of Figure 8 and in the two upper matrices of Table 7, OE and CE are consistently lower than 25% for all classes, with the exception of soybean (OE = 25.3%-40.7%,CE = 50.7%-55.0%,with either ERN or ERN+s input), which is mainly misclassified as maize.For the soybean class, the use of X-band σ° results in a OE reduction of 15.4%, due to less confusion with other summer crops, and in a CE reduction of 4.3%, due to less confusion with rice; a reduction of 4.9% CE for rice is also registered using ERN+s input, while errors over other crop types are stable.The best performances, with class errors not exceeding 10% and balanced between omission and commission, are achieved for maize (OE = 9.3%, CE = 10.6%),rice (OE = 9.3%, CE = 1.6%), forages (OE = 0.1%, CE = 0.5%), and forestry-woodland (OE = 0.6%, CE = 0.1%).Tendencies to overestimation for winter crop (CE = 16.8%) and to underestimation for double crop (OE = 23.7%) are observed, due to mutual confusion between these two classes.When the CTmin-max-ave scheme is applied to 2014 data (transferability set) in case r (Figure 8, light and dark orange bars, Table 7, two lower confusion matrices), some different patterns emerge: both CE and OE are greater than 2013 over most of the classes.A remarkable case is represented by soybean, with OE = 100.0%,meaning that this class is not represented in the classified pixels belonging to 2014 validation set.These pixels are mistakenly classified as rice (78%) and maize (22%).Other major discrepancies between 2013 and 2014 per-class accuracy occur for winter crop and double crop classes: for the former, the 2014 crop map is strongly underestimating (OE = 58%), while for the latter the rather conservative 2013 performance (CE = 0%) is not repeated for 2014 (CE > 40%, OE > 30%).The overestimation of winter crop is due to confusion with double crop, while the errors observed for double crop come from misclassification not only with winter crop, but also with forages.This could be due to the timing of 2014 acquisition dates, which less effectively capture the single and double crop dynamics compared to 2013 data.The high OE for forestry-woodland (56%) is instead due to class confusion with forages.
The integration of σ° for the 2014 dataset (ERN+s) results in very small changes of OE for soybean, winter crop and double crop, while it brings an improvement for maize (OE = 7.0%, CE = 4.9%), rice (OE = 0.2%, CE = 12.4%), forages (OE = 14.9%,CE = 20.0%), and forestry-woodland (OE = 0.3%, CE = 11.7%); this is mainly due to the reduced confusion between rice and maize and between forages and forestry-woodland.Still, per-class performances using ERN+s over 2014 are generally worse than for 2013, especially over crop classes more sensitive to the seasonal climatic conditions (e.g., double crop and soybean, usually sown late in summer season in Lombardy and thus prone to meteorological fluctuations).A different behavior is shown by rice, with a slight overestimation of class extent for 2014 dataset (<OE, >CE), compared to 2013 results.In summary, CTmin-max-ave performance using ERN+s input is better than with ERN for most of the crop types at Lev1 and outputs acceptable class accuracies ranging from 62.8% for double crop, to >90% for forestry-woodland (94.0%), rice (93.8%), and maize (94.6%).Important misclassification errors are observed for soybean class across the interannual dataset.

Error Reduction Strategy
Since Lev1 map assessment showed some misclassification for summer crop classes (especially on 2014 dataset), we tested an expert-based ex-post pruning of the CTmin-max-ave scheme, implemented by restructuring the set of rules for branches with higher class attribution cumulated error (>20%).The rationale is to test the performance of a classifier which reduces the overall classification error at the expenses of the detail of the thematic level.This way, we generate a crop map which is only partially at Lev1 detail (see Table 7 for Lev1 areal coverage percentage), by grouping summer crops, which are in the left branches of the schemes shown in Figure 6 together into a generic crop type label (generic summer crop, SCg) (Figure 9).These re-structured crop mapping schemes, implemented for both ERN and ERN+s are respectively named: CT'min-max-ave(ERN), shown in Figure 9a, and CT'min-max-ave(ERN+s) scheme, shown in Figure 9b.
Validation was carried out using the same validation sets used for assessing CTmin-max-ave results, by excluding the areas labelled as SCg, which are classified now at Lev0.As a consequence, CT'min-max-ave early crop maps do not cover the whole study area at Lev1 and some of the cropland is mapped at Lev0.Table 8 summarizes accuracy metrics of the CT'min-max-ave scheme calculated for two validation datasets: case t and case r.The Lev1 classified coverage ranges from a minimum of 84% (obtained for the ERN 2013 map), to a maximum of 94% of the total study site cropland area (for the ERN 2014 map).As expected, the global accuracy scores achieved with the ex-post pruned schemes are higher  Per-class error analysis (Figure 10) shows that the re-structured schemes provide an improvement only to CE of the summer crop classes, maize and rice (being soybean excluded as target here) with a significant reduction for rice CE (4.4%) over 2014 using ERN+s input set.In summary, the CT'min-max-ave showed slight but consistently better performance of the original scheme, and could be adopted as error reduction strategy when the proposed crop mapping approach is applied to different conditions; this would lead to a crop map with different thematic levels: Lev1, i.e., distinguishing the majority of rice and maize fields , and Lev0.

Conclusions
This study describes a classification tree approach for in-season crop mapping over Northern Italy during the early summer season (mid-July) from the integration of optical (Landsat 8 OLI) and X-band SAR (COSMO-SkyMed) synoptic seasonal features.A rule-based approach offers the advantage of being interpretable through rules/conditions applied to input features, which are representative of crop conditions and development.Results described could be applied to Northern Italy and, with minimal check and tuning by local experts, also to areas with similar environmental and agricultural characteristics (i.e., European temperate to Mediterranean areas).
Key findings and conclusions relevant for crop mapping applications are listed below:  The proposed approach produces early in-season (mid-July) crop type maps at two levels of thematic detail with the greatest accuracy obtained when both optical and SAR features (ERN+s set) are used as input: overall accuracy is 91.8% for the 2013 season and 86.6% for the 2014 season;  Best performing input features for effectively distinguishing 7 crop types (maize, rice, soybean, winter crop, double crop, forages, forestry-woodland) can be extracted from synoptic seasonal features calculated for winter and summer crops and derived from two combinations of remote sensing proxies for vegetation and soil conditions: i) EVI, NDFI and RGRI from OLI data (ERN set), and ii) the integration of OLI-derived proxies with CSK backscattering (ERN+s set). The contribution of X-band σ° (HH polarization) is relevant for promoting the transferability of the approach over a season (2014) different from the one used for developing the classification rules, with an increment of 19.7% in OA compared to crop maps produced using only optical input features;  The integration of CSK σ° reduces class errors (omission and commission) especially for crop types characterized by more seasonally stable agronomic patterns (i.e., maize, rice, forages, forestry-woodland);  Expert-based tuning and ex-post pruning are key assets when dealing with operational monitoring and can be used as error reduction strategy, delivering a modified early crop mapping scheme with hybrid thematic level output, and higher overall accuracy (90.3% using ERN+s input for 2014 season);  In the framework of agriculture management, the achieved overall accuracy at mid-July is considered satisfactory given the fact that the information on crops is provided early during the growing season as management requires, i.e., 2-3 months before end of season and harvesting;  The proposed the approach is interpretable and flexible enough for being exploited for mapping crops at different levels of detail and possibly exploiting different input data with similar spectral bands (e.g., Sentinel-2 MSI).

Figure 2 .
Figure 2. EVI (a) and NDFI (b) multi-temporal profiles extracted from the 2014 OLI dataset for major crops in Lombardy study site (covering the whole growing season from winter crop stem elongation phase to summer crops harvesting: March to October).Satellite acquisition dates of OLI and CSK images used for the development set (2013) are superimposed on the graphs.Grey color box represent the target temporal range for producing the early in-season crop map (i.e., mid-July).

Table 3 .Figure 3 .
Figure 3.An example of synoptic seasonal features (min, max, ave, std) extracted for the winter and summer periods from EVI 2014 multitemporal profiles (April to mid-July) over winter wheat (red) and maize (yellow) sample fields.

Figure 4 .
Figure 4. Flow chart of the overall methodological approach.

Figure 7 .Figure 8 .
Figure 7. Early in-season crop maps produced using the CTmin-max-ave(ERN+s) scheme over the study area, for the year 2013 (a) and 2014 (b).Crop maps are overlaid on mid-July EVI (16 July 2013 for panel a, 19 July 2014 for panel b), represented in grey tones.

Table 1 .
The X-band SAR and optical satellite acquisitions divided into the development (2013) and the transferability dataset (2014).

Table 2 .
Composition of the reference sample for 2013 development set, including training and validation sub-samples.

Table 4 .
Mean J-MDIST for each crop type class at Lev1, as a function of the combination of OLI and SAR proxies.Maximum separability corresponds to J-MDIST = 2.

Table 6 .
OA and κ computed for crop maps derived with the CT schemes of Figure 6 applied to the development (2013) and transferability (2014) sets, with ERN and ERN+s input features.case r = class distribution proportional to real case crop acreage (%); case t = same class distribution of training set; ERN = EVI+NDFI+RGRI; ERN+s = EVI+NDFI+RGRI+σ°.