European Rice Cropland Mapping with Sentinel-1 Data : The Mediterranean Region Case Study

Rice farming is one of the most important activities in the agriculture sector, producing staple food for the majority of the world's growing population. Accurate and up-to-date assessment of the spatial distribution of rice cultivated area is a key information requirement of all stakeholders including policy makers, rice farmers and consumers. Timely assessment with high precision is, e.g., crucial for water resource management, market prices control and during humanitarian food crisis. Recently, two Sentinel-1 (S-1) satellites carrying a C-band Synthetic Aperture Radar (SAR) sensor were launched by the European Space Agency (ESA) within the homework of the Copernicus program. The advanced data acquisition capabilities of S-1 provide a unique opportunity to monitor different land cover types at high spatial (20 m) and temporal (twice-weekly to biweekly) resolution. The objective of this research is to evaluate the applicability of an existing phenology-based classification method for continental-scale rice cropland mapping using S-1 backscatter time series. In this study, the S-1 images were collected during the rice growing season of 2015 covering eight selected European test sites situated in six Mediterranean countries. Due to the better rice classification capabilities of SAR cross-polarized measurement as compared to co-polarized data, S-1 cross-polarized (VH) data were used. Phenological parameters derived from the S-1 VH backscatter time series were used as an input to a knowledge-based decision-rule classifier in order to classify the input data into rice and non-rice areas. The classification results were evaluated using multiple regions of interest (ROIs) drawn from high-resolution optical remote sensing (SPOT 5) data and the European CORINE land cover (CLC 2012) product. An overall accuracy of more than 70% for all eight study sites was achieved. The S-1 based classification maps reveal much more details compared to the rice field class contained in the CLC 2012 product. These findings demonstrate the potential and feasibility of using S-1 VH data to develop an operational rice crop monitoring framework at the continental scale.


Introduction
Europe is the fourth largest importer of rice in the world (Figure 1a).Over the last five years, Europe's annual rice imports were on average about 1 million tons (milled basis), ranging between 873 million tons in 2011 and 1190 million tons in 2015 [1].Meanwhile, the size of rice cultivated areas has either reduced or remained stable around 420,000 hectares (Figure 1b).Information about how the area of rice croplands varies from year to year is an important piece of information for the European economy being relevant for: risk management for the insurance industry [2], environmental reporting, Delta to a continental scale (Mediterranean region).To achieve this objective, we used a dense time series stack of S-1 backscatter data as input to map rice paddy area at fine spatial scale over eight study sites in six European countries.

Study Sites Characteristics
In the European Union (EU) the total area of cultivated rice is about 430,000 hectares.The growing areas are mostly located in the Mediterranean countries [1] (Figure 2), where the summer seasons are warm and dry.The normal growing season is either from April/May to September/October or from May to October/November depending on the temperature.Italy and Spain are the top rice producer EU member countries followed by Greece, Portugal and France [37].[1]; and (b) the EU rice balance for 2008 to 2013 [37].

Study Sites Characteristics
In the European Union (EU) the total area of cultivated rice is about 430,000 hectares.The growing areas are mostly located in the Mediterranean countries [1] (Figure 2), where the summer seasons are warm and dry.The normal growing season is either from April/May to September/October or from May to October/November depending on the temperature.Italy and Spain are the top rice producer EU member countries followed by Greece, Portugal and France [37].

Study Sites Characteristics
In the European Union (EU) the total area of cultivated rice is about 430,000 hectares.The growing areas are mostly located in the Mediterranean countries [1] (Figure 2), where the summer seasons are warm and dry.The normal growing season is either from April/May to September/October or from May to October/November depending on the temperature.Italy and Spain are the top rice producer EU member countries followed by Greece, Portugal and France [37].Rice varieties grown in Europe mostly belong to the Japonica (70%) and Indica (30%) species group.The recent evolution of Japonica and Indica areas (in hectares) in the EU member States is shown in Figure 3. Rice is mostly grown in congregated areas such as in the Po valley in Italy, the Rhône delta in France, and the Thessaloniki area in Greece.In Spain, rice cultivation is more scattered; rice growing areas are found in the Aragon region, the Guadalquivir valley, the Ebro delta and Valencia Albufera [38].In Portugal, rice cultivation area is concentrated mainly in three regions: the Tagus and Sorraia valleys, the Mondego, and the Sado and Caia river valleys [39,40].In Turkey, the most productive regions are Thrace and Marmara, which are producing 10-15% of the total national rice production.Rice varieties grown in Europe mostly belong to the Japonica (70%) and Indica (30%) species group.The recent evolution of Japonica and Indica areas (in hectares) in the EU member States is shown in Figure 3. Rice is mostly grown in congregated areas such as in the Po valley in Italy, the Rhône delta in France, and the Thessaloniki area in Greece.In Spain, rice cultivation is more scattered; rice growing areas are found in the Aragon region, the Guadalquivir valley, the Ebro delta and Valencia Albufera [38].In Portugal, rice cultivation area is concentrated mainly in three regions: the Tagus and Sorraia valleys, the Mondego, and the Sado and Caia river valleys [39,40].In Turkey, the most productive regions are Thrace and Marmara, which are producing 10-15% of the total national rice production.Rice planting in Europe involves direct seeding into flooded soil (water seeding) or dry soil (dry seeding).With both the farming methods, the floodwater is maintained until the harvest season.Water is drained several times from the field prior to harvesting so that the fields can be dried and harvesting equipment can pass through.Fields are also drained for early foliar herbicide treatments, then re-flooded within few days.The drainage period allows certain weed species to germinate in the aerobic environment (Figure 4).

Sentinel 1A Data
For this study, we accessed archived Interferometric Wide Swath (IW) mode S-1 data acquired during the rice crop growing cycle in 2015, from April to early November to completely cover every test site.The S-1A SAR IW mode acquisitions come from eight different tracks, whereas each test site is covered by the same track (Figure 2).The details for each test site regarding data acquisition, location, date, path and range of incidence angle are shown in Figure 4.All the images were obtained from ESA as standard Level 1 GRD (ground-range detected) high resolution images.Only the S-1 IW acquisitions with VH polarization were selected for rice cropland classification.

Optical Data and Ancillary Data
For the study sites in France, and Italy, multi-temporal Spot 5 (10 m spatial resolution) for the year 2015 were downloaded from the SPOT website [41].These optical datasets were used to produce reference classification maps for validating the S-1 rice maps.Rice cropland area-a vector dataset with a minimum mapping unit (MMU) of 1 hectare-over Seville in Spain for the year 2015 has been Rice planting in Europe involves direct seeding into flooded soil (water seeding) or dry soil (dry seeding).With both the farming methods, the floodwater is maintained until the harvest season.Water is drained several times from the field prior to harvesting so that the fields can be dried and harvesting equipment can pass through.Fields are also drained for early foliar herbicide treatments, then re-flooded within few days.The drainage period allows certain weed species to germinate in the aerobic environment (Figure 4).

Sentinel 1A Data
For this study, we accessed archived Interferometric Wide Swath (IW) mode S-1 data acquired during the rice crop growing cycle in 2015, from April to early November to completely cover every test site.The S-1A SAR IW mode acquisitions come from eight different tracks, whereas each test site is covered by the same track (Figure 2).The details for each test site regarding data acquisition, location, date, path and range of incidence angle are shown in Figure 4.All the images were obtained from ESA as standard Level 1 GRD (ground-range detected) high resolution images.Only the S-1 IW acquisitions with VH polarization were selected for rice cropland classification.

Optical Data and Ancillary Data
For the study sites in France, and Italy, multi-temporal Spot 5 (10 m spatial resolution) for the year 2015 were downloaded from the SPOT website [41].These optical datasets were used to produce reference classification maps for validating the S-1 rice maps.Rice cropland area-a vector dataset with a minimum mapping unit (MMU) of 1 hectare-over Seville in Spain for the year 2015 has been obtained from the Institute of Statistics and Cartography of Andalusia's website [42] Two others vector dataset of rice cropland area-over Valencia in Spain and Thessaloniki in Greece for the year 2015 were retrieved from the ERMES (An Earth Observation Model Based Rice Information Service)'s website [43] Google Earth imagery and Sentinel-2 [44] optical data were used to produce reference data for the rest of study areas.In addition to this, the European CORINE land cover (CLC 2012) product was also used for comparing classification results across all European test sites in a consistent manner.The data set is not as detailed and accurate as the other two datasets.However, it is the only available reference datasets covering all our study regions.Figure 5 shows the land cover composition over the eight selected test sites in 2012.Rice cropland is the most common land cover in Lombardia, Italy with 33.6%, followed in Ebro, Spain (18.8%) in Valencia, Spain (12.4%), and only 2.3% of rice cropland in Mondego, Portugal.Inland water makes up over 10% of the total land covers in Camargue, France and below 1.5% for the others.It is the same for wetland class, e.g., 7% in Camargue, France, and below 2% for the others.
obtained from the Institute of Statistics and Cartography of Andalusia's website [42] Two others vector dataset of rice cropland area-over Valencia in Spain and Thessaloniki in Greece for the year 2015 were retrieved from the ERMES (An Earth Observation Model Based Rice Information Service)'s website [43] Google Earth imagery and Sentinel-2 [44] optical data were used to produce reference data for the rest of study areas.In addition to this, the European CORINE land cover (CLC 2012) product was also used for comparing classification results across all European test sites in a consistent manner.The data set is not as detailed and accurate as the other two datasets.However, it is the only available reference datasets covering all our study regions.Figure 5 shows the land cover composition over the eight selected test sites in 2012.Rice cropland is the most common land cover in Lombardia, Italy with 33.6%, followed in Ebro, Spain (18.8%) in Valencia, Spain (12.4%), and only 2.3% of rice cropland in Mondego, Portugal.Inland water makes up over 10% of the total land covers in Camargue, France and below 1.5% for the others.It is the same for wetland class, e.g., 7% in Camargue, France, and below 2% for the others.

Methodology
The time series algorithm used in this study was introduced by Nguyen et al. [36].It consists of six steps as illustrated in Figure 6 5) classification using a knowledge-based decision-tree approach; and (6) accuracy assessment based on reference data.

Pre-Processing
All the selected S-1 SAR IW images were pre-processed (orbit correction, radiometric calibration, resampling and geocoding) using ESA's Sentinel-1 Toolbox.The geocoding step involved a Range Doppler Terrain correction algorithm that uses the elevation data from the 1 arc-second DEM product from the Shuttle Radar Topography Mission (SRTM) and POD orbit state vectors provided by ESA.Notice that precise orbits files (POD) are produced few weeks after the acquisition and they are automatically downloaded from ESA website [45].In this process, data are resampled and geo-coded to a grid of 10 m spacing preserving the 20 m spatial resolution.

Methodology
The time series algorithm used in this study was introduced by Nguyen et al. [36].It consists of six steps as illustrated in Figure 6 5) classification using a knowledge-based decision-tree approach; and (6) accuracy assessment based on reference data.

Pre-Processing
All the selected S-1 SAR IW images were pre-processed (orbit correction, radiometric calibration, resampling and geocoding) using ESA's Sentinel-1 Toolbox.The geocoding step involved a Range Doppler Terrain correction algorithm that uses the elevation data from the 1 arc-second DEM product from the Shuttle Radar Topography Mission (SRTM) and POD orbit state vectors provided by ESA.Notice that precise orbits files (POD) are produced few weeks after the acquisition and they are automatically downloaded from ESA website [45].In this process, data are resampled and geo-coded to a grid of 10 m spacing preserving the 20 m spatial resolution.

Identification of Potential Rice Pixels
For the identification of potential rice growing areas, our approach is to threshold the dynamic range backscatter image to identify image pixels that change more than the defined threshold value (dB).Threshold value selection depends on the nature and expected changes in the magnitude of VH backscatter and the SAR geometry (e.g., incidence angle).A generalized threshold for the rice fields can only be determined if the optimum SAR data acquisition is guaranteed (e.g., SAR observations are available in the flooded and vegetative stages).Otherwise, the threshold must be optimized considering the data acquisition and the constraints of local crop calendar.S-1A provides at least one acquisition after every 12 days over the selected study sites (see Figure 4); after the launch of S-1B the temporal sampling was reduced to 5-6 days.Based on the visual interpretation of optical imagery for the selected period and expert knowledge acquired from the ancillary data, a threshold of 8.5 dB was used to extract the potential rice pixels referred to as _ .Figure 7 (first column) shows sensitivity images based on dry reference (P05) and wet reference (P95) images for all study areas, and the results of applying thresholds of 8, 8.5 and 9 dB to these images are illustrated in following columns of the Figure 7.The CLC 2012 map is used to indicate the compartment boundaries, and pixels exhibiting change below and upper the threshold (8.5 dB) are classified as potential rice cropland areas.Combining Figure 6 with Figure 7 suggests a threshold on sensitivity of 8.5 dB to extract the potential rice cropland pixels.Lowering the threshold to 8.0 dB increases area outside the potential rice cropland areas in all the study sites leads to computation time increases.Raising the threshold to 9.0 dB decreases the potential rice cropland areas may leads to numerous rice pixels are omitted.Thus 8.5 dB seems the better choice as long as we accept that physically, processing time, as well as classification precision.With this threshold, the majority of the rice cultivated areas are correctly picked out in all data sets.From the basic concept of the approach, it is clear that classification error will occur for the rice pixels whose variation is less than the selected threshold.Different of farming activities during the growing season (rice varieties, water level in the fields, density of rice plants in the fields), SAR acquisitions period may give rise to this type of misclassification.

Identification of Potential Rice Pixels
For the identification of potential rice growing areas, our approach is to threshold the dynamic range backscatter image to identify image pixels that change more than the defined threshold value (dB).Threshold value selection depends on the nature and expected changes in the magnitude of VH backscatter and the SAR geometry (e.g., incidence angle).A generalized threshold for the rice fields can only be determined if the optimum SAR data acquisition is guaranteed (e.g., SAR observations are available in the flooded and vegetative stages).Otherwise, the threshold must be optimized considering the data acquisition and the constraints of local crop calendar.S-1A provides at least one acquisition after every 12 days over the selected study sites (see Figure 4); after the launch of S-1B the temporal sampling was reduced to 5-6 days.Based on the visual interpretation of optical imagery for the selected period and expert knowledge acquired from the ancillary data, a threshold of 8.5 dB was used to extract the potential rice pixels referred to as σ o V H_potential .Figure 7 (first column) shows sensitivity images based on dry reference (P05) and wet reference (P95) images for all study areas, and the results of applying thresholds of 8, 8.5 and 9 dB to these images are illustrated in following columns of the Figure 7.The CLC 2012 map is used to indicate the compartment boundaries, and pixels exhibiting change below and upper the threshold (8.5 dB) are classified as potential rice cropland areas.Combining Figure 6 with Figure 7 suggests a threshold on sensitivity of 8.5 dB to extract the potential rice cropland pixels.Lowering the threshold to 8.0 dB increases area outside the potential rice cropland areas in all the study sites leads to computation time increases.Raising the threshold to 9.0 dB decreases the potential rice cropland areas may leads to numerous rice pixels are omitted.Thus 8.5 dB seems the better choice as long as we accept that physically, processing time, as well as classification precision.With this threshold, the majority of the rice cultivated areas are correctly picked out in all data sets.From the basic concept of the approach, it is clear that classification error will occur for the rice pixels whose variation is less than the selected threshold.Different of farming activities during the growing season (rice varieties, water level in the fields, density of rice plants in the fields), SAR acquisitions period may give rise to this type of misclassification.

Time Series Filtering
Filtering the backscatter time series has the purpose of reducing the short-term influence of environmental conditions and noise inherent in the S-1 data due to speckle and other noise-like influences.The processed output is a smoothed backscatter signal ( _ ), which will be used for the extraction of different phenological stages of rice (e.g., a start of the season, heading time, and length of season).For the temporal filtering of the _ time series a Gaussian smoothing filter (with the standard deviation of 3 dB for the kernel) was used.A detailed investigation on the selection of kernels with different standard deviation was already reported by Nguyen, Wagner et al. 2015 [16].Moreover, in order to discriminate the rice pixels from the other land cover classes we have empirically defined a list of three more static thresholds based on the _ values.

Time Series Filtering
Filtering the backscatter time series has the purpose of reducing the short-term influence of environmental conditions and noise inherent in the S-1 data due to speckle and other noise-like influences.The processed output is a smoothed backscatter signal (σ o V H potential_smooth ), which will be used for the extraction of different phenological stages of rice (e.g., a start of the season, heading time, and length of season).For the temporal filtering of the σ o V H_potential time series a Gaussian smoothing filter (with the standard deviation of 3 dB for the kernel) was used.A detailed investigation on the selection of kernels with different standard deviation was already reported by Nguyen, Wagner et al. 2015 [16].Moreover, in order to discriminate the rice pixels from the other land cover classes we have empirically defined a list of three more static thresholds based on the σ o V H potential_smooth values.

Extraction of Vegetation Phenology Parameters
For this step, we calculated three phenological indicators, namely, date of beginning of season (DoS), date of maximum backscatter (DoM), and Length of the season (LoS).This allows to delineating the areas as "rice paddy", and all other areas were placed in a generalized "non-rice" class.These parameters, introduced by [37], are summarized in Table 1.
Table 1.Phenological parameters for the rule-based classification.

DoS
During the growing season, the date of the beginning of season is defined as the first local minima in σ o V H potentia_smooth time-series.

DoM
During the growing season, the date when backscatter reaches a maximum value is defined as the local maxima in σ o V H potentia_smooth time-series.This date must come after the date of the beginning of season, where it reaches its local minimum backscatter value.

LoS
The length of the season is defined as the number of days difference between DoM and DoS.

Rice Paddy Identification
Due to the high temporal variability in the SAR backscatter signal across the different study sites, the raw output from the thresholding of phenological profiles contained some noisy pixels.This implies that most fields were not fully classified as rice and non-rice class at the pixel level.Therefore, we constrained the minimum mapping unit to the average farm size in the Mediterranean region, which is 1/4 hectares.This means that no polygon was composed of fewer than 25 S-1 pixels (20 m spatial resolution, 10 m pixel spacing).We implemented this through a post-classification processing step, whereas a majority/minority analysis with the window size of 5 × 5 pixels was applied to remove the small pixel groups in order to obtain refined classification results.

Accuracy Assessment
For validation and evaluation of classification results, standard accuracy assessment measures were used, i.e., kappa coefficient, overall accuracy, omission error, and commission error.Rice cropland vector layer 2015 for study site in Spain (Seville), rice cropland maps created through interpretation of SPOT-5 data for three study sites in Spain (Valencia), Italy, and France; rice cropland raster layers 2006 from CLC 2006 for Marmara-Thrace, Turkey; and rice cropland raster layers from CLC 2012 for all study sites were used.

Temporal Rice Backscatter Signature from Sentinel-1 SAR Data
Seasonal VH time series are shown in Figure 8a-h, for Mondego (Portugal), Seville (Spain), Valencia (Spain), Ebro delta (Spain), Camargue (France), Lombaria (Italy), Thessaloniki (Greece) and Marmma-Thrace (Turkey), respectively.In addition, Figure 8a-h also shows the results of smooth backscatter profiles (magenta color), and the date of maximum and minimum backscatter values for the selected rice fields (represented by the red and blue dots, respectively).To complement this analysis, Figure 9 shows the VH backscattering coefficients and false color composites of the study sites.
Figure 9 (column 4) shows the color composites, which are created by using the multi-temporal SAR acquisitions in order to highlight the temporal characteristics within the rice fields.The red, blue, and green colors in these figures (Figure 9, column 4) correspond to the images acquired during the flooded/seeding (April/May), heading (August/September) and post harvested (October/November) period, respectively.The blue and dark green color regions in these figures (Figure 9, column 4) indicate the rice cropland areas.At the beginning of the growing season (April-May) the σ o V H values from the rice fields were very low due to flooding after sowing.Thus, in the early stage of rice crop growth the fields appeared as dark areas in SAR images.This condition corresponds to those scenes that were acquired during the months April and May (see black areas in Figure 9, the first column).
In general, during this period the rice fields show low backscatter values, e.g., less than −20 dB (for reference, see Figure 8).However, there are several conditions related to the soil preparation and wind speed/direction that must also be considered.Some rice fields require a special tillage, such as furrows, or the presence of low water level in the field for a shorter period due to the diversity in farming activities among different parts of the region.If these furrows are shallow and under light wind conditions, the surface is smooth irrespective of the furrow direction.However, if the furrows are deep, the furrow direction becomes an important factor regarding SAR image acquisition geometry.In case that the deep furrows are nearly parallel to the radar viewing direction, the surface seems smooth in radar image.On the other hand, if the deep furrows are perpendicular to the radar viewing direction, the surface becomes strongly rough, and the signal backscatter becomes very strong.These results are in accordance with Brisco et al. (1991) when evaluating the effect of tillage row direction in relation to the radar's look direction using radar backscattering coefficient from three different radar frequencies [45].The Sentinel-1 data are also influenced by the incidence angle, whereas the strength of the radar's backscatter signal gradually decreases with the increase of incidence angle.Therefore, backscatter is in general higher over the Marmma-Thrace (Turkey) and Mondego (Portugal) sites, where the incidence angles were 35 • and 36 • compared to the other regions where the incidence angle was about 40 • (Figures 4 and 8).
In the second stage (vegetative stage), backscatter value increases as the vegetation grows (e.g., plant size increases), and eventually the SAR images show no significant difference between rice fields and other agricultural fields or vegetated areas (see bright areas in Figure 9, second column).One month to 45 days after the start date of the growing season it reaches the first peak in June/July.Specifically, a simple visual inspection to the σ o V H images which were acquired at the end of July (where DoY is around 210, Figure 8a-h)) reveals that the backscatter value has dropped suddenly, despites the fact that vegetation is fully developed.This anomaly is observed for all study sites, except Mondego (Portugal) and Lombaria (Italy).It is clearly illustrated in the case of Seville (Spain) and Valencia (Spain) in the Figure 8.The scale of these changes is quite different for different study sites due to different agricultural practices.Like in other European regions, rice crop in Seville and Valencia (Spain) begins in mid-April with the deep placement of fertilizers under dry conditions.Flooding starts during the first week of May and then the seeds are sown.At the end of June there is a short dry period of ten days [46].Furthermore, this period is characterized by the increase of the rice plant height and the number and size of tillers that lead to make free space among narrowed or blocked rice stems.As a result, the absence of water from the rice fields will minimize the double bounce effect of SAR signals and this explains the decrease in backscattering values at this stage.For the other regions, throughout the rice cultivation period, water is commonly kept at a depth of 4-8 cm, and drained away 2-3 times during the growing season to improve the crop rooting, reduce the algae growth and to allow application of herbicides.For the reproductive stage (August/September), backscatter values continuously keep increasing until they achieve the maximum value in September.During this period the values of backscattering coefficient vary between the range of −17 dB and −13 dB, which might get influenced by the variations in incidence angle, water level in the fields, cultivation activities or the rice varieties.However, during the ripening phase a slight decrease in SAR backscatter signal is observed.One potential reason for this could be due to the fact that the plants will dry before the harvesting.Normally, rice fields are drained towards the end of August to allow harvesting [47].After harvest, the fields can have diverse conditions, either bare and dry fields or covered with weeds in wet conditions.Fields may also be flooded due to local farming activities, e.g., in some areas fields are flooded again until January for duck hunting [46].This event is clearly visible at the study sites in Spain (Figures 8b-d and 9).However, the levels of change are quite different among the study sites due to the differences in farming activities.Meanwhile, in other regions, some small crops or weeds cover is possible right after the rice crop is cultivated.The radar backscatter is thus variable and in most cases can have high backscatter values.Moreover, as a result of rain which is typically After harvest, the fields can have diverse conditions, either bare and dry fields or covered with weeds in wet conditions.Fields may also be flooded due to local farming activities, e.g., in some areas fields are flooded again until January for duck hunting [46].This event is clearly visible at the study sites in Spain (Figures 8b-d and 9).However, the levels of change are quite different among the study sites due to the differences in farming activities.Meanwhile, in other regions, some small crops or weeds cover is possible right after the rice crop is cultivated.The radar backscatter is thus variable and in most cases can have high backscatter values.Moreover, as a result of rain which is typically high in the winter season in the Mediterranean region, the water content in bare soil increases, and this may also explain the increase of backscattering.Consequently, there is little possibility to interpret this last stage of the cultivation.

Thresholds Selection
Optical data sources (Spot 5, Sentinel-2, and Google Earth imagery), vector dataset and expert knowledge were used to create polygons for threshold selection.These polygons were carefully selected and digitized across eight study sites.Training areas were considered in terms of distribution and variety of rice fields (e.g., size, density, rice variety) and SAR geometry characteristic (e.g., incidence angle) to build robust training data set.These polygons used for threshold determination were then excluded from the set of polygons used for accuracy assessment.Time series analysis was carried out using all available S-1 IW Mode data between April and November 2015 (Figure 8 and Section 4.1).The overall behavior of σ o V H is comparable to our previous investigations [37].However, the dynamic range over different regions varies due to different incidence angles and farming activities.In Figure 8, it is important to note that the response signal patterns from all sample parcels are very consistent.Therefore, we can conclude that these temporal signatures are associated with the rice crop type, and the empirical thresholds of three phenological parameters based on the σ o V H potentia_smooth values can be used for the identification of rice cropland.This study, we followed the approach of Nguyen et al. [15] who used phenological parameters (derived from σ o V H profile) to map the rice areas of eight sites in the Mediterranean region.
For rice areas mapping, the first two key parameters, σ o DoM -the peak (maximum) of VH backscatter and ∆σ o the amplitude backscatter (the difference between σ o DoM -backscatter at date of maximum backscatter and σ o DoS -backscatter at the begin of the growing season) are very critical.The peak and valley of VH backscatter (i.e., σ o DoM , σ o DoS ) within the one rice-growing cycle can be identified by a local extrema algorithm.To eliminate unrealistic peaks, a threshold for VH backscatter is required.These two parameter thresholds (i.e., σ o DoM ≥ −19 dB and ∆σ o ≥ 2.5 dB) were selected based on a conservative deduction from training data, limited census statistics, and expert knowledge [36].The corresponding positions (date) and backscatter values of these two points showed in the Figure 8 where big-red dots are DoS(s) and big-blue dots are DoM(s) respectively.All areas that practice rice crops start the growing season between May-June where the lowest backscatter is observed.Rice flowering period then comes after more than one and a haft month (July-August) where S-1 respond signal in the temporal signature reaches the peak.
The third parameter, length of the season (LoS) which is defined as the temporal distance (number of days difference) between the date of beginning of season and the date when maximum backscatter value is recorded.This temporal distance has to be greater than the shortest possible rice growing cycle and smaller than the longest possible rice growing cycle.From crop calendar and our expert knowledge about rice growing season in the study areas, the threshold of temporal distance which is about 50 and 120 days, respectively.If all three conditions are met then the pixel is classified as rice, otherwise as non-rice area.

Spatial Distribution and Comparison of S-1 Derived Rice Area with Reference Data
The temporal backscatter signatures is shown in Figure 8 (and Figure S1 in the supplementary material) served to define specific thresholds which were applied to the S-1 time series in order to generate rice area maps for the growing season of 2015.The classification results of rice fields over eight selected sites in the Mediterranean region are shown in Figure 10.The classification accuracy has been assessed by comparing the classification results to the 2015 rice cropland vector layer for four study sites in Seville (Spain), Valencia (Spain), Lombardia (Italy), and Camargue (France); and CLC (2012) for all study sites.The confusion matrices are presented in Tables 1 and 2, respectively.To provide direct comparison, the classification results in all study sites: Seville and Valencia in Spain, and Thessaloniki in Greece were evaluated by using the reference vector dataset from the same year 2015; Lombardia (Italy), and Camargue (France) were evaluated by using higher spatial resolution (SPOT 5, 10 m); and Earth Imagery and Sentinel-2 data were used to produce reference data to validate the classification results for the rest of study areas (Table 2).The Kappa coefficients, rice user accuracy (commission error), and rice producer accuracy (omission error) were, respectively, 0.76, 70.2%, and 82.7% for Mondego, Portugal; 0.87, 86.8%, and 89.2% for Seville, Spain; 0.85, 81.7%, and 93.3% for Valencia, Spain; 0.79, 70.3 and 95.3% for Ebro delta, Spain ; 0.85, 80.2%, and 94.0% for Camargue, France; 0.82, 79.8%, and 86.2% for Lombardia, Italy; 0.79, 74.1%, and 94.0% for Thessaloniki, Greece; and 0.76, 71.1%, and 84.1% for Marmara-Thrace, Turkey (Table 1).The lower accuracy level was observed in Mondego, Portugal, Ebro delta, Spain, Thessaloniki, Greece and Marmara-Thrace, Turkey, because of most rice fields in these regions associated with the small farms, and scattered throughout the mixed agricultural landscape and, thus, were easily omitted with non-rice classes.One potential source of error lies in the time difference between the reference data and the Sentinel-1 image acquisition.For example, the reference data were produced by digitizing homogenous sites of rice fields based on the image visualization and interpretation of existing high-resolution Google Earth imagery and Sentinel-2 optical data, while the classification maps were produced from 2015 S-1A data.Moreover, a number of different of farming activities during the growing season (rice varieties, water level in the fields, and density of rice plants in the fields) and SAR acquisitions period could also cause an increase in mapping errors.Although the proposed approach can properly detect most rice areas, some land cover types with high variation in backscatter values (e.g., wetland or seasonal water bodies areas) can cause commission errors.To provide direct comparison, the classification results in all study sites: Seville and Valencia in Spain, and Thessaloniki in Greece were evaluated by using the reference vector dataset from the same year 2015; Lombardia (Italy), and Camargue (France) were evaluated by using higher spatial resolution (SPOT 5, 10 m); and Earth Imagery and Sentinel-2 data were used to produce reference data to validate the classification results for the rest of study areas (Table 2).The Kappa coefficients, rice user accuracy (commission error), and rice producer accuracy (omission error) were, respectively, 0.76, 70.2%, and 82.7% for Mondego, Portugal; 0.87, 86.8%, and 89.2% for Seville, Spain; 0.85, 81.7%, and 93.3% for Valencia, Spain; 0.79, 70.3 and 95.3% for Ebro delta, Spain ; 0.85, 80.2%, and 94.0% for Camargue, France; 0.82, 79.8%, and 86.2% for Lombardia, Italy; 0.79, 74.1%, and 94.0% for Thessaloniki, Greece; and 0.76, 71.1%, and 84.1% for Marmara-Thrace, Turkey (Table 1).The lower accuracy level was observed in Mondego, Portugal, Ebro delta, Spain, Thessaloniki, Greece and Marmara-Thrace, Turkey, because of most rice fields in these regions associated with the small farms, and scattered throughout the mixed agricultural landscape and, thus, were easily omitted with non-rice classes.One potential source of error lies in the time difference between the reference data and the Sentinel-1 image acquisition.For example, the reference data were produced by digitizing homogenous sites of rice fields based on the image visualization and interpretation of existing high-resolution Google Earth imagery and Sentinel-2 optical data, while the classification maps were produced from 2015 S-1A data.Moreover, a number of different of farming activities during the growing season (rice varieties, water level in the fields, and density of rice plants in the fields) and SAR acquisitions period could also cause an increase in mapping errors.Although the proposed approach can properly detect most rice areas, some land cover types with high variation in backscatter values (e.g., wetland or seasonal water bodies areas) can cause commission errors.For all the study sites, CLC from 2012 were also used for comparison and checking the consistency of classification results over all study areas.Despite the limitations of the CLC data (e.g., there was a three years difference between CLC and S-1 datasets) the results suggest that the application of S-1 time series data for rice area mapping has produced consistent results for all the test sites with overall accuracies (in quotation {}) ranging from 77.7% to 98.9% (kappa average at 0.53).The average accuracies at Camargue (France), Thessaloniki (Greece) and Mondego (Portugal) are lower than that of Seville and Valencia (Spain).The results showed relatively high error rate for both commission and omission measure (see Table 3 for details).The best results were achieved over Seville and Valencia (Spain) with kappa of 0.70 (Rice omission 31.9%) and 0.82, respectively (Rice omission 22.6%), which is

Figure 2 .
Figure 2. Study sites (red) and spatial extent of the used S-1 scenes (blue).

Figure 2 .
Figure 2. Study sites (red) and spatial extent of the used S-1 scenes (blue).

Figure 2 .
Figure 2. Study sites (red) and spatial extent of the used S-1 scenes (blue).

Figure 3 .
Figure 3. Rice varieties areas in the EU member states [37].

Figure 3 .
Figure 3. Rice varieties areas in the EU member states [37].

Figure 4 .
Figure 4. Growing season of paddy rice crop and the available S-1A SAR scenes over eight selected test sites.

Figure 4 . 21 Figure 5 .
Figure 4. Growing season of paddy rice crop and the available S-1A SAR scenes over eight selected test sites.

Figure 5 .
Figure 5. Relative contribution of harmonized land-cover categories in eight selected test sites.

Figure 6 .
Figure 6.Systematic workflow of S-1 SAR (VH polarization) data processing, rice area classification and validation.

Figure 6 .
Figure 6.Systematic workflow of S-1 SAR (VH polarization) data processing, rice area classification and validation.

Figure 7 .
Figure 7. Sensitivity images and potential rice cropland area extent results using different thresholds.

Figure 8 .
Figure 8. Temporal evolution of the backscattering coefficients derived from VH polarization (where, −20 dB is base line).

Figure 9 .
Figure 9. Columns 1-3: images of σ o V H acquired at three dates (from left to right) and; Column 4: images of false color composite over the part of study sites.

Table 2 .
Confusion matrices of the accuracy assessment.