Next Article in Journal
Sustainable Land Use Diagnosis Based on the Perspective of Coupling Socioeconomy and Ecology in the Xiongan New Area, China
Next Article in Special Issue
Integrated Flood Impact and Vulnerability Assessment Using a Multi-Sensor Earth Observation Mission with the Perspective of an Operational Service in Lombardy, Italy
Previous Article in Journal
Evaluate Human Perception of the Built Environment in the Metro Station Area
Previous Article in Special Issue
Spatial Distribution Changes in Nature-Based Recreation Service Supply from 2008 to 2018 in Shanghai, China
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Development of a Pre-Automatized Processing Chain for Agricultural Monitoring Using a Multi-Sensor and Multi-Temporal Approach

1
Institute of Polar Sciences of the National Research Council of Italy (ISP CNR), Montelibretti, 00015 Rome, Italy
2
Institute for Advanced Studies of Pavia (IUSS), 27100 Pavia, Italy
3
Institute for Environmental Protection and Research (ISPRA), 00144 Rome, Italy
*
Author to whom correspondence should be addressed.
Submission received: 5 December 2023 / Revised: 29 December 2023 / Accepted: 4 January 2024 / Published: 12 January 2024

Abstract

:
Understanding crop types and their annual cycles is key to managing natural resources, especially when the pressures on these resources are attributable to climate change and social, environmental, and economic policies. In recent years, the space sector’s development, with programs such as Copernicus, has enabled a greater availability of satellite data. This study uses a multi-sensor approach to retrieve crop information by developing a Proof of Concept for the integration of high-resolution SAR imagery and optical data. The main goal is to develop a pre-automatized processing chain that explores the temporal dimension of different crop. Results are related to the advantage of using a multi-sensor approach to retrieve vegetation biomass and vertical structure for the identification of phenological stages and different crops. The novelty consists of investigating the multi-temporal pattern of radiometric indices and radar backscatter to detect the different phenological stages of each crop, identifying the Day of the Year (DoY) in which the classes showed greater separability. The current study could be considered a benchmark for the exploitation of future multi-sensor missions in downstream services for the agricultural sector, strengthening the evolution of Copernicus services.

1. Introduction

The agricultural sector is of fundamental importance, and it is one of sectors that is the most vulnerable to climate change [1]. At the same time, intensive and unsustainable agricultural practices deplete and degrade natural resources and key ecosystem services. In fact, in European countries, extreme events are responsible for reducing crop yields by 9% and 3.1% [2]. The European Union (EU) is significantly promoting sustainable agricultural practices [3]. One of the most important instruments provided by the EU is the Common Agricultural Policy (CAP), which aims to improve environmental protection and climate sustainable practices. The management of the EU’s CAP involves Member States through national paying agencies, which are responsible for monitoring the use of funding by the farming community [4]. The availability of information on agriculture and crops is crucial for successful agronomic planning and sustainable agriculture management [5]. However, official statistics on crop areas are often provided at the end of the season or later and are therefore not useful for supporting in-season crop management [6]. CAP regulations encourage the development of regular and systematic monitoring services [7].
Earth Observation (EO) has been used for agricultural monitoring since the early 1970s. The Copernicus Sentinel constellation allows the integration of SAR and optical sensors for agricultural applications, overcoming their specific limitations (i.e., dependence on the sun and clouds for optical and signal noise and stability for SAR) [8]. Optical sensors are sensitive to canopy response and photosynthetic characteristics, while SAR sensors are more sensitive to plant biomass, soil moisture, and surface texture [6].
Based on the above, Schiavon et al. (2021) asked institutional users to indicate the most important environmental and climate-related requirements in terms of current agro-environmental legislation and assessed the potential of current European EO technology to meet their monitoring needs. The results of this work show that users request the development of consolidated agricultural EO products, the provisioning of seasonal and annual information on crop production, and in-season indicators of crop development and status. Furthermore, they express the need to monitor the degree of diversification of agri-environmental measures in the agricultural sector in order to support green direct payments [7]. Thus, an effective multi-temporal cropland mapping is fundamental for the correct management of natural resources and for the prevention of risks due to climatic conditions such as drought. In fact, 14% and 12% of users indicated an interest in agricultural and water resources domains, respectively, and in the development of new EO products in these domains [7].
Optical sensors have demonstrated their ability to distinguish different crops [9,10], and in recent years, Synthetic Aperture Radar (SAR) sensors have also been used for crop mapping and monitoring [11,12,13]. Several authors have investigated the capabilities of SAR data in crop mapping. Successful results have been obtained with the C band, which is particularly valuable for rice phenology analysis [14]. Additionally, techniques based on X-band sensors showed encouraging results, especially dual-pol imagery [15]. Other types of crops mapped using SAR data include winter wheat, irrigated grassland, and summer cereals [11,16,17]. However, in addition to backscatter values, optical data can provide relevant and complementary information at the field and crop level [12,18]. Several authors investigated the synergic use of SAR and optical data, demonstrating that radar is very sensitive to soil moisture, especially in VV polarization, while optical data are very useful for determining vegetation characteristics [8,19,20]. However, all these case studies are based on the use of field data to train implemented algorithms, but these data are not always available, and the dependence on this information does not allow the development of an operational crop mapping service. Crop mapping remains a complex task for operational activities, is still highly dependent on field data, and relies on local knowledge of management practices [21]. It is of paramount importance to develop a methodology for crop mapping that is mostly based on the intrinsic properties of the crops and their phenology rather than on field data. In addition, the most used SAR data are in the C band, with a lower spatial resolution than the X band’s data resolution. Some improvement can be provided by the recent development of a new constellation of X-band sensors, significantly contributing to crop mapping thanks to higher spatial and temporal resolutions [6]. From this perspective, the development of an automated processing chain is necessary, especially to fulfill users’ need for targeted agricultural EO products with recurrent and frequent data to produce seasonal and annual information. This evaluation has provided important insights for the development of new sensors, also considering the investments planned by the European Union for further satellite missions [22]. Contextually, the Copernicus program delivers environmental information largely based on EO satellite data in the form of Copernicus Services, addressing six thematic areas. Copernicus operational services are not static but need to evolve with emerging user requirements and state-of-the-art methodologies for developing and testing algorithms for pre-operational and operational product prototypes and improving and developing future-specific Copernicus services with the potential for global applications [23].
It is paramount to bridge the gap between the user requirements and the technical limitations of existing sensors by developing applications for the new generation of satellites. The NOCTUA project was a pilot project for the development of a commercial service for collecting, processing, and distributing radar data designed for the Lombardy region (Northen Italy). This project foresaw the development of a new Low Earth Orbit (LEO) SAR satellite in the X band, with a revisit time of 11 days and a ground resolution of less than 1 m (in Spotlight mode).
NOCTUA has been considered a Proof of Concept (POC) of the IRIDE NEXT-GENERATION EU investment, which will develop an Italian LEO satellite constellation called IRIDE. This constellation will provide a downstream service at the national level to supply data from different types of sensors (multi-, hyperspectral, infrared, SAR), also considering existing missions, and to develop a series of operational services in several application domains: coastal, air quality, water resource, land cover, etc.
Several methodologies for crop mapping and monitoring have been developed in recent years. However, the shortcomings for operational use have not been broadly explored. A valuable use of crop mapping and monitoring algorithms exploiting EO data from the perspective of consistent operational product development and demonstrating a clear configuration of technical and operational users’ requirements of the relevant aspects of an operational service provision is still lacking. Thus, the main objective of this study is to develop a Proof of Concept (POC) of a pre-automatized processing chain based on a multi-sensor approach for crop mapping with the purpose of being potentially integrated into an operational service architecture, boosting the possibility of addressing users’ needs, with respect to their urgency, the closeness to the operational delivery process, and the availability of capacities. The integration of high-resolution SAR in X-band imagery and multispectral optical data has been explored to analyze the phenological stages of the crops, using an agricultural area located in the Lombardy region as a test site. Existing X-band SAR data (i.e., TerraSAR-X) have been used to simulate the NOCTUA signal and to integrate the vegetation biomass and vertical structure in the temporal dimension of different crop cycles.
The article is organized as follows: after the introduction of the current remote sensing-based crop monitoring techniques and applications, Section 2 describes the study area and addresses the input data used and the methodology implemented to pre-process and then process optical and SAR data. Section 3 describes the intermediate and final results. In Section 4, the final results are discussed, pointing out the advantages, limitations, and future developments, and finally, Section 5 summarizes the derived conclusions.

2. Materials and Methods

2.1. Study Area

The study area is located in the Lombardy region, Northern Italy (Figure 1), and lies south of Milan, encompassing parts of the Milan and Pavia provinces. It covers around 444 km2, is crossed by the Ticino River, and includes four main cities, namely Pavia, Gambolò, Vigevano, and Abbiategrasso, with populations of 70,636, 9730, 6076, and 32,425 people, respectively [24].
The Area of Interest (AOI) is in the Po plain, the largest irrigation basin in the region. Thirty-five percent of the national agricultural production occurs here in the Po River valley, which is a strategic area for the Italian economy. It generates the 40% of the country’s gross domestic product and has a population of more than 16 million [25]. In this area, the intense use of water for farming leads to water scarcity, which climate change will worsen [26]. This site has experienced several drought episodes (2003, 2005, 2006) [27], which make the area particularly vulnerable to the impacts of climate change. Flooding and droughts have a stronger impact on the delta region of the river basin [28]. However, the water supply faces both quantitative and qualitative pressures. Firstly, there is a high water demand, which originates from the agricultural sector. Additionally, the most commonly utilized irrigation method is surface irrigation, which generally lacks efficiency in water usage [29].
The test site is characterized by a temperate climate, with January being the coldest month, with a mean temperature of 2.9 °C, and July being the hottest at 25 °C. The area has an average annual precipitation of 750 mm [30]. Thirty-four percent of the total area is cultivated. The main types of crops are (in percentage of total cropland area) rice (50%), grain maize (10%), mixed herbage (8%), and tree crop (7%) [31,32]. The cultivation calendar of the main crops is illustrated in Figure 2 [9,33]. In this region, two major cropping seasons are distinguished, summer and winter; the main crop production takes place from April to the end of October.

2.2. Reference Data

To produce accurate crop maps and ensure the reliability of the results, declarations from farmers obtained from the Sistema Informativo Agricolo della Regione Lombardia (SIARL), referred to 2016, and provinces of Milan and Pavia have been requested [31,32]. This dataset consists of the location of agricultural parcels with information regarding the crop cultivated in each parcel. This information is provided by farmers and is mandatory for CAP funding applications. It contains information about the crops cultivated in each parcel, the total area of the parcel, and the area of the parcel dedicated to each crop. The data have been divided into two sets to obtain training and validation datasets.
A land use map has been obtained from the project “Destinazione d’Uso dei Suoli Agricoli e Forestali” (DUSAF) of the Lombardy region from 2015 [34]. From the Copernicus Land Cover Monitoring Service (CLCMS), a series of products have been freely downloaded: Imperviousness map (2018), Riparian zones product (2013), and EU-Hydro database (2013) [35]. These products have been used to derive a mask from excluding non-agricultural areas, such as urban areas, forests, riparian zones, and water bodies.

2.3. Satellite Data

The remote sensing dataset (Table 1) is composed of 21 images from TerraSAR-X (TSX) and 13 images from Sentinel-2 (S2) for 2016. The list of the images used is reported in Table 2, with the acquisition dates and the relative Day of the Year (DoY). Cloud cover for each image in the S2 dataset is reported, with the average value between the two tiles. The TSX sensor was chosen to simulate NOCTUA data because of the similarity of their technical characteristics.
The images were selected to include at least one acquisition per month for the winter season and two for the spring–summer period to have complete data coverage throughout the entire year and a better temporal resolution in the growing period of summer crops.
TSX images have been acquired in Stripmap mode, Single Look Complex format, single polarization (HH), and with the same flight track direction. S2 images have been acquired at Level 2A from Theia Catalogue (https://www.theia-land.fr/, accessed on 20 October 2021) [36]. The AOI is not completely included within a single tile of S2; for this reason, two tiles (T32TMR and T32TNR) for each image have been acquired to create a mosaic (Figure 3).
The implemented methodology is shown in Figure 4 and explained in detail in the following paragraphs.
Pre-processing and processing steps for the satellite data were carried out to build an automatized model to obtain the results directly. This model is fundamental for the realization of an automatized processing chain.

2.4. Pre-Processing

2.4.1. Filtering, Grouping and Class Selection

Based on the Declaration of Farmers dataset, no-crop areas have been excluded. The resulting dataset included more than 100 types of crop. Thus, these were grouped and listed in 24 classes according to characteristics of their growing cycle (e.g., silo maize and maize for the production of energy are harvested at the same stage of the growing cycle) and spectral similarity (e.g., barley and wheat were renamed “Winter cereals”). The fields with an area lower than 1 ha and the classes with less than 10 fields were discarded. Finally, the first seven classes in order of frequency (Figure 5) were selected and randomly sampled to create training and validation datasets using the centroid of each polygon. Seventy percent of the fields for each class were used to create the training dataset; the remaining thirty-percent composed the validation dataset [28].

2.4.2. Subsetting, Masking, and Mosaicking

S2 images were pre-processed applying a cloud mask, obtained from the Scene Classification raster at 20 m of resolution.
All the S2 images were resampled at 10 m with the nearest neighbor interpolation and clipped on the AOI. Although Sentinel-2 provides a range of multispectral bands with different spatial resolutions (from 10 to 60 m), and the lack of panchromatic band disables the direct production of a set of fine-resolution (10 m) bands. In this pre-automatized processing chain, the methodology for detecting crop phenology mainly exploits the original 10 m spatial resolution of S2 bands n. 2-4-8.
The tiles have been mosaicked and masked by applying a no-crop mask and water mask obtained from the declaration of farmers, DUSAF maps, and Copernicus products mentioned above.

2.4.3. Calibration and Geocoding

Using the associated orbital information, TSX images were calibrated and geocoded in WGS84—UTM 32N. Then, intensity images were obtained. This pre-processing was carried out using ®SARScape Analytics Toolbox (version 5.5). All the pre-processing steps were integrated into the Time Series tool (version 5.5). Also, this dataset was masked on AOI, water, and no-crop areas.

2.5. Processing

2.5.1. Multiple Endmember Spectral Mixture Analysis (MESMA)

The Multiple Endmember Spectral Mixture Analysis (hereafter MESMA) has been applied on a single S2 image to obtain a preliminary crop classification. To capture the greatest spatial variability in the phenology of agricultural patches, the S2 image of the maize and soybeans in the growing season, freshly harvested winter cereals, and freshly flooded rice was used. This algorithm consists of a spectral mixture analysis that considers a larger number of pure spectra, called Endmembers (EMs), for each pixel, as an extension of the most used Linear Spectral Mixture Analysis (LSMA) that does not consider the contrast between the materials within the pixel [22,37,38]. Thus, the MESMA overcomes the limitations of LSMA by requiring a model to meet the minimum fit, fraction, and residual constraints while testing multiple models for each image pixel. Therefore, this approach allows mapping more than four materials across an image [39].
The MESMA application procedure consists of four steps: (1) the selection of the EMs starting from a random sampling of the training dataset, (2) the creation of a spectral library, (3) the application of the spectral unmixing to obtain Fractional Abundance Maps (FAMs), and (4) shade normalization. At the end of the process, the FAMs were used to obtain an intermediate crop classification map. A majority analysis was carried out using a kernel size of 3 × 3 pixels to remove the noise and the spurious pixels within agricultural parcels. This analysis assigned the central pixel to the class to which most kernel pixels belonged [38,40]. The classification map was validated with a confusion matrix, using the validation dataset to calculate its accuracy.

2.5.2. Radiometric Indices

Three different radiometric indices were extracted from each S2 image as features for the crop classification: the Normalized Difference Vegetation Index (NDVI), the Leaf Area Index (LAI), and the Bare Soil Index (BSI).
The first two indices are widely used in studies on vegetation and are calculated using Blue (Band 2), Red (Band4), and NIR (Band8) channels of S2 (Equations (1) and (2)):
NDVI = (NIR − Red)/(NIR Red)
LAI = 3.618 × EVI − 0.118 where EVI = 2.5 × (NIR − Red)/(NIR 6 × Red − 7.5 × Blue 1)
In particular, the LAI, defined as the total one-sided area of leaf tissue per unit ground, is a biophysical indicator used to represent the dimension of the crop canopy and its variation over time. LAI measurements have been widely adopted for crop monitoring and modeling applications, being a key state variable associated with processes including light interception and the soil–crop water balance [41].
The BSI [42] is a specific index to identify bare soil in agricultural contexts. It combines NDVI with the Normalized Difference Build-Up Index and allows better identification of exposed soil surfaces and uncultivated areas based on soil characteristics. It is calculated using Blue (Band2), Red (Band4), NIR (Band8), and SWIR (Band11) channels of S2 (Equation (3)):
BSI = ((SWIR1 Red) − (NIR Blue))/((SWIR1 Red) + (NIR Blue))
NDVI and BSI provide a good way to distinguish bare soil from vegetation, while LAI and NDVI can better analyze the crop’s phenological stages to distinguish various crop types.
Then, for each crop, multi-temporal profiles of indices’ mean values were included in the decision tree to properly analyze crops’ seasonal behaviors.

2.5.3. Time Series Analysis

TSX images were used to build a Time Series Analysis using the SARScape Analytics® toolbox, provided by ®NV5 Geospatial Solutions, Inc (Broomfield, CO, USA). A geocoded dataset calculating multi-temporal statistics from intensity images was created to detect and extract temporal changes [43].
In bare soils, the radar backscatter depends on the roughness and the dielectric constant of the first centimeters of soil (1 cm for X band) [19]. During the year, any changes in the agricultural soil are due to soil tillage, the humidity condition, and the land cover variation (e.g., presence of vegetation). Plowing, seeding, and irrigation operations result in a change in surface roughness. For these reasons, it is possible to associate SAR signal variation with the variation of these features over time.
The main statistics considered, such as (1) minimum value, (2) minimum date associated with low roughness values, (3) maximum value, (4) maximum date associated with high roughness values, (5) coefficient of variation (CoV, ratio between standard deviation and the mean), (6) mean value, (7) standard deviation were finally mapped, and a multi-temporal profile of mean backscatter values for each crop has been obtained.

2.5.4. Decision Trees

Decision trees were built based on the intermediate products obtained from the previous analysis. We focused on two different approaches: the first is based on the use of optical data only (Figure 6a), and the second uses a synergic approach including both optical and SAR data (Figure 6b). These two approaches verified the hypothesis that integration radar and optical sensors have an advantage in the accuracy of crop classification, analyzing the multi-temporal behavior of optical radiometric indices and the backscatter signal.
Figure 6 shows the decision trees. For both approaches, the first branch of the decision trees was chosen to distinguish summer crops from perennial (tree crops) and winter crops. For the subsequent tree nodes, thresholds were selected based on the values of the indices. The selection was made by looking at multi-temporal profiles and choosing dates when the classes showed greater separability. In addition, for the S2-based approach, the FAMs obtained from the MESMA were used to classify vetch and soybean. To maximize the accuracy, each node was evaluated using the training dataset extracted from the reference data, selecting the variables and thresholds that would give better results.
The resulting maps were filtered to remove the noise and the spurious pixels, as explained in Section 2.5.1.
The final classification maps were validated using confusion matrices.

3. Results

3.1. MESMA-Based Crop Classification

The intermediate crop classification map obtained by MESMA and the confusion matrix are shown in Figure 7. This classification has an overall accuracy (OA) of 37%. Considering the producer accuracy, the best-classified crops are soybean and maize with 54% and 44% accuracy, respectively. On the other hand, lowest values of accuracy occur for tree crops (17%) and pasture (10%).
The decision tree uses fractions to distinguish the maize class and separate the vetch from the other crops.

3.2. Radiometric Indices

The training data were divided into crop classes, and the mean value of each radiometric index was calculated for each class. The temporal profile of each radiometric index was obtained using these mean values. NDVI and LAI temporal profiles (Figure 8) show a clear pattern for the summer crop (e.g., maize), with the peak values in the summer season. All the profiles are coherent with the crop calendar shown in Figure 2. The only exception is represented by winter cereal. This class shows a maximum peak of NDVI in August (DOY 223) and minimum values in November and December (DOYs 313 and 363), while according to the crop calendar, these peaks should occur in April–May and August, respectively. This behavior can be explained by the local agricultural practice that is well identifiable in the NDVI multi-temporal profile (Figure 9): winter cereals (WC), mainly ray and durum wheat, are sown in October/November, and the plants sprout in January. The cereals are harvested in May–June, and forage is sown immediately after. Finally, they are harvested in August, and the cycle begins again.
LAI temporal profiles showed a slightly different trend than NDVI, mainly due to the structure of the plants. Using this index, we obtained added-value information on crop typologies that maximized class (Figure 8b).
The temporal BSI trend provided information on the presence of bare soil. The profile in Figure 10 shows that summer crops (e.g., maize) have the minimum value of this index in summer months, consistently with their growing season according to the crop calendar. Moreover, tree crops (TCs) showed a different behavior, with a significantly lower presence of bare soil.

3.3. Time Series Analysis

From the time series analysis, we obtained information regarding the temporal trend of backscatter values, such as (1) the minimum date associated with low roughness values, (2) the maximum date associated with high roughness values, and (3) the coefficient of variation—hereafter called “CoV” (the ratio between standard deviation and the mean). In particular, the latter one provided information on the variation in land cover classes over time.
The analysis of CoV calculated for each crop has been used as proxy to identify the stability of the signal strictly correlated with the stability of the crop over time. As is shown in Figure 11, this proxy assumes higher values for the crops that show a higher variability within the year. The tree crops class is an exception, showing a lower CoV value because it is more stable over time.

3.4. Decision-Tree-Based Crop Classifications

The crop classification maps obtained using decision trees are shown in Figure 12, with the confusion matrices calculated using the validation dataset. Although the overall accuracy (OA) is slightly different between the two classifications (44% and 42% for TSX + S2 and S2, respectively), the producer accuracy for each crop showed some considerable differences.
In particular, the crop-classification map obtained using TSX + S2 data (Figure 12b) shows a greater accuracy than the classification based only on S2 data (Figure 12a). Tree crops showed a producer accuracy of 70% for TSX + S2 and of 57% for S2, with an increase of 13%. Also, the Vetch class showed an increase of the accuracy using the TSX + S2 approach, from 34% to 39% for S2 and TSX + S2, respectively.
In general, four out of seven classes showed an increase in the producer accuracy with the TSX + S2 classification approach, with the exception of winter cereals, maize, and soybean.

4. Discussion

The present work explores a multi-sensor and multi-temporal approach to develop a Proof of Concept (POC) for an automatized processing chain for monitoring crop cycles. It can be applied to monitor crops in different seasons, providing a reliable tool for identifying agricultural practices such as crop rotation. The results can support the Common Agricultural Policy (CAP) and the Italian institutional users, which encourage the development of mapping and monitoring through semi-automatic algorithms for identifying crops and agricultural activities throughout the year [7]. The added value of this POC is related to the advantage of retrieving vegetation biomass and a vertical structure, supporting intra-annual crop-rotation mapping, which influences both the sustainability of agricultural practices and the payments related to the European CAP.
The added value of a multi-sensor approach has been investigated using TerraSAR-X (TSX) and Sentinel 2 (S2) data in synergy. As demonstrated by several authors, SAR data are widely used for rice and grassland monitoring [11,14,21], but the results of this work are related to the advantage of using a multi-sensor approach to retrieve many classes of crops. The increase in the OA of comparing a single-sensor approach (S2) with a multi-sensor approach (S2 + TSX) becomes remarkable considering the producer accuracy of specific crops (e.g., tree crops), with an improvement of 13% due to the use of the coefficient of variation (CoV) from the radar backscatter. The CoV was one of the most relevant parameters in supporting the separation of tree crops from the other classes.
A fundamental step was the analysis of the separability of crop classes, analyzing the time trends of the main vegetation and soil descriptors (radiometric indices) and the backscatter signal. If we consider only the S2-based variables, pasture and soybean have similar responses, but by integrating the TSX multi temporal backscatter, it is possible to identify a peak for soybean on three dates (4 August, 26 August, and 6 September).
The limited separability of some classes can explain the low values of OA. For example, rice and vetch’s multi-temporal patterns are similar for all variables. The limitation related to the separating certain crops that show a similar multi-temporal pattern can be overcome by integrating the processing chain variables derived from SAR data (e.g., the flooded period of rice). Furthermore, the farmers’ declaration used to train the processing chain is affected by a bias: it is an annual declaration, so only one crop on each field is reported without considering an intra-annual crop rotation.
Thus, the analytical approach developed in this study is independent of the availability of field data since the parameters used (i.e., radiometric indices and backscatter signal) are mostly related to the intrinsic properties of the crops and their phenology (i.e., vegetation biomass and vertical structure). These characteristics remain almost unchanged under similar climatic conditions. This approach has been developed for the Lombardy Region, and it is geographically limited to the specific study area until the agroclimatic conditions do not change enough to determine different patterns in crop phenology. To apply the processing chain to different agroclimatic scenarios, it would require new training datasets regarding the crop typologies and the crop typologies, and this could open new issues related to the overlapping of radiometric signals among crop types. The MESMA decision tree should also be re-adapted in terms of thresholds, considering that this part of the process cannot be completely automatized.
The processing chain is ready to be used to trace the patterns of different years, and furthermore, it is not limited to S2 and TSX but is ready to ingest images acquired by other sensors because it is based on algorithms that are well consolidated for different sensors.
In addition, the processing chain is strictly dependent on the use of the Copernicus Land Cover Monitoring Service (CLCMS), integrating several products (i.e., Imperviousness map, Riparian Zones product, EU-Hydro database) as input data. Thanks to the continuous acquisition of the Sentinel 2 constellation, the optical component of the pre-automatized processing chain is ready to be operational, but the SAR-based component currently relies on non-operational data (i.e., TSX). The development of the IRIDE constellation, which will provide data with higher spatial and temporal resolution, also offers perspective on the systematic acquisition for the SAR component. The future developments of this work might considerably improve the existing Copernicus portfolio by responding to the user requirements with seasonal and annual crop thematic products. Thus, the obtained annual crop classifications, analyzing the seasonal crop development and status, could be useful to monitor the diversification degree of agri-environmental measures and to assist national paying agencies in monitoring the use of CAP funding.

5. Conclusions

The European Commission (EC) requires each Member State to monitor the proper use of funds and is considering the use of satellite imagery and Earth Observation (EO) techniques to monitor the condition of crops [44]. Thus, using EO data for the requirements related to the emission of direct payments represents a consistent improvement in cost effectiveness [7].
In this framework, the present work develops a Proof of Concept for an automated processing chain based on TerraSAR-X and Sentinel 2 data to retrieve vegetation biomass and vertical structure and to monitor agricultural practices.
Two classification maps were produced and compared, the first based only on S2 data and the second based on the integration between optical and radar sensors (TSX + S2). Integrating two sensors for crop identification provided evidence of improved accuracy. The novelty consists of investigating the multi-temporal pattern of radiometric indices and radar backscatter to detect the different phenological stages of each crop, identifying the Day of the Year (DoY) in which the classes showed greater separability.
The processing chain prototype should aim at demonstrating the technical operational feasibility toward a more automated approach in the provision of EO data processing and pre-operational and operational products, so it can complement and broaden the currently available panoply of existing products and be a benchmark for the development of future operational service for agricultural management. Paying agencies and national institutions can benefit from this operational service to easily monitor the use of Common Agricultural Policy (CAP) funds.
To develop an automatized processing chain, it is necessary to continuously acquire the required input data over time. For these reasons, this work also contributes to defining mission requirements for upcoming and future constellations (e.g., IRIDE) focusing precisely on the synergistic use of different sensors, which is an added value to obtain a more consistent automatized downstream service.

Author Contributions

Conceptualization, A.T., E.V. and E.S.; methodology, S.S. and E.V.; formal analysis, S.S. and B.M.; investigation, E.V. and S.S.; resources, A.T. and E.S.; writing—original draft preparation, S.S.; writing—review and editing, A.T., E.V., M.R. and E.S.; supervision, A.T., E.V. and E.S.; funding acquisition, A.T. All authors have read and agreed to the published version of the manuscript.

Funding

This study has been partially funded by the project NOCTUA “Landscape Monitoring. For Everyone. From Space,” funded under the Regional Operative Program (POR) of Lombardy Region 2014–2020, call hub innovation and research (Project ID: 1179775) and by Italian National Resilience and Recovery Plan (PNRR) M1C2-25—Investment 4.1.2: Satellite Technology and Spatial Economy, Financed by the European Union—NextGenerationEU, Purchase Order N. 5001035211 (Call for Views on the EO PNRR System And/Or Element Architecture).

Data Availability Statement

Sentinel-2 data were obtained from ESA and are freely available from the Copernicus Open Access Hub portal (https://scihub.copernicus.eu, accessed on 21 September 2021), and from Theia catalogue (https://www.theia-land.fr/, accessed on 20 October 2021. TerraSAR-X data were obtained from DLR and are available from the EOWEB GeoPortal (https://eoweb.dlr.de/egp/, accessed on 7 October 2021) with the permission of DLR. Copernicus products were obtained from the Copernicus Land Monitoring Service (https://land.copernicus.eu/en/products, accessed on 20 October 2021). Metereological data and declarations from farmers were obtained from Regione Lombardia and are freely available from the Regione Lombardia—Open Data portal (https://www.dati.lombardia.it/, accessed on 10 December 2021).

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Tao, S.; Xu, Y.; Liu, K.; Pan, J.; Gou, S. Research Progress in Agricultural Vulnerability to Climate Change. Adv. Clim. Change Res. 2011, 2, 203–210. [Google Scholar] [CrossRef]
  2. Monteleone, B.; Borzí, I.; Bonaccorso, B.; Martina, M. Quantifying Crop Vulnerability to Weather-Related Extreme Events and Climate Change through Vulnerability Curves; Springer: Dordrecht, The Netherlands, 2023; Volume 116, ISBN 0123456789. [Google Scholar]
  3. European Commission Commissions Staff Working Document: Analysis of Links between CAP Reform and Green Deal. 2020, pp. 1–22. Available online: https://agriculture.ec.europa.eu/news/cap-reforms-compatibility-green-deals-ambition-2020-05-20_en (accessed on 7 October 2021).
  4. European Court of Auditors Special Report 04/2020: Using New Imaging Technologies to Monitor the Common Agricultural Policy: Steady Progress Overall, but Slower for Climate and Environment Monitoring; Publications Office of the European Union: Luxembourg, 2020; Volume 60.
  5. Fontanelli, G.; Crema, A.; Azar, R.; Stroppiana, D.; Villa, P.; Boschetti, M. Agricultural Crop Mapping Using Optical and SAR Multi-Temporal Seasonal Data: A Case Study in Lombardy Region, Italy. In Proceedings of the International Geoscience and Remote Sensing Symposium (IGARSS), Québec City, QC, Canada, 13–18 July 2014; pp. 1489–1492. [Google Scholar]
  6. Villa, P.; Stroppiana, D.; Fontanelli, G.; Azar, R.; Brivio, P.A. In-Season Mapping of Crop Type with Optical and X-Band SAR Data: A Classification Tree Approach Using Synoptic Seasonal Features. Remote Sens. 2015, 7, 12859–12886. [Google Scholar] [CrossRef]
  7. Schiavon, E.; Taramelli, A.; Tornato, A.; Lee, C.M.; Luvall, J.C.; Schollaert Uz, S.; Townsend, P.A.; Cima, V.; Geraldini, S.; Nguyen Xuan, A.; et al. Maximizing Societal Benefit Across Multiple Hyperspectral Earth Observation Missions: A User Needs Approach. J. Geophys. Res. Biogeosci. 2023, 128, e2023JG007569. [Google Scholar] [CrossRef]
  8. Meroni, M.; D’Andrimont, R.; Vrieling, A.; Fasbender, D.; Lemoine, G.; Rembold, F.; Seguini, L.; Verhegghen, A. Comparing Land Surface Phenology of Major European Crops as Derived from SAR and Multispectral Data of Sentinel-1 and -2. Remote Sens. Environ. 2021, 253, 112232. [Google Scholar] [CrossRef] [PubMed]
  9. Azar, R.; Villa, P.; Stroppiana, D.; Crema, A.; Boschetti, M.; Brivio, P.A. Assessing In-Season Crop Classification Performance Using Satellite Data: A Test Case in Northern Italy. Eur. J. Remote Sens. 2016, 49, 361–380. [Google Scholar] [CrossRef]
  10. Tatsumi, K.; Yamashiki, Y.; Canales Torres, M.A.; Taipe, C.L.R. Crop Classification of Upland Fields Using Random Forest of Time-Series Landsat 7 ETM+ Data. Comput. Electron. Agric. 2015, 115, 171–179. [Google Scholar] [CrossRef]
  11. Hajj, M.; Baghdadi, N.; Belaud, G.; Zribi, M.; Cheviron, B.; Courault, D.; Hagolle, O.; Charron, F. Irrigated Grassland Monitoring Using a Time Series of TerraSAR-X and COSMO-SkyMed X-Band SAR Data. Remote Sens. 2014, 6, 10002–10032. [Google Scholar] [CrossRef]
  12. Thiong’o, K.; Pasternak, R.; Kleusberg, A.; Thonfeld, F.; Menz, G. Separability of Dominant Crop Cultures in Southern Germany Using TerraSAR-X Data. Adv. Remote Sens. 2015, 4, 97–107. [Google Scholar] [CrossRef]
  13. Gella, G.W.; Bijker, W.; Belgiu, M. Mapping Crop Types in Complex Farming Areas Using SAR Imagery with Dynamic Time Warping. ISPRS J. Photogramm. Remote Sens. 2021, 175, 171–183. [Google Scholar] [CrossRef]
  14. Bazzi, H.; Baghdadi, N.; El Hajj, M.; Zribi, M.; Minh, D.H.T.; Ndikumana, E.; Courault, D.; Belhouchette, H. Mapping Paddy Rice Using Sentinel-1 SAR Time Series in Camargue, France. Remote Sens. 2019, 11, 887. [Google Scholar] [CrossRef]
  15. Lopez-Sanchez, J.M.; Ballester-Berman, J.D.; Hajnsek, I. Rice Monitoring in Spain by Means of Time Series of TerraSAR-X Dual-Pol Images. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2011, 4, 412–422. [Google Scholar] [CrossRef]
  16. Yi, Z.; Jia, L.; Chen, Q. Crop Classification Using Multi-Temporal Sentinel-2 Data in the Shiyang River Basin of China. Remote Sens. 2020, 12, 4052. [Google Scholar] [CrossRef]
  17. Zalite, K.; Antropov, O.; Praks, J.; Voormansik, K.; Noorma, M. Monitoring of Agricultural Grasslands with Time Series of X-Band Repeat-Pass Interferometric SAR. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2016, 9, 3687–3697. [Google Scholar] [CrossRef]
  18. Ranghetti, M.; Boschetti, M.; Ranghetti, L.; Tagliabue, G.; Panigada, C.; Gianinetto, M.; Verrelst, J.; Candiani, G. Assessment of Maize Nitrogen Uptake from PRISMA Hyperspectral Data through Hybrid Modelling. Eur. J. Remote Sens. 2022, 56, 2117650. [Google Scholar] [CrossRef]
  19. Bousbih, S.; Zribi, M.; Lili-Chabaane, Z.; Baghdadi, N.; El Hajj, M.; Gao, Q.; Mougenot, B. Potential of Sentinel-1 Radar Data for the Assessment of Soil and Cereal Cover Parameters. Sensors 2017, 17, 2617. [Google Scholar] [CrossRef] [PubMed]
  20. Setiyono, T.D.; Quicho, E.D.; Gatti, L.; Campos-Taberner, M.; Busetto, L.; Collivignarelli, F.; García-Haro, F.J.; Boschetti, M.; Khan, N.I.; Holecz, F. Spatial Rice Yield Estimation Based on MODIS and Sentinel-1 SAR Data and ORYZA Crop Growth Model. Remote Sens. 2018, 10, 293. [Google Scholar] [CrossRef]
  21. Wu, B.; Zhang, M.; Zeng, H.; Tian, F.; Potgieter, A.B.; Qin, X.; Yan, N.; Chang, S.; Zhao, Y.; Dong, Q.; et al. Challenges and Opportunities in Remote Sensing-Based Crop Monitoring: A Review. Natl. Sci. Rev. 2023, 10, nwac290. [Google Scholar] [CrossRef]
  22. Valentini, E.; Taramelli, A.; Marinelli, C.; Piedelobo, L.; Fassari, M.; Troffa, S.; Casa, R.; Pignatti, S. Hyperspectral Mixture Models in the CHIME Mission Implementation for Topsoil Texture Retrieval. JGR Biosci. 2023, 128, e2022JG007272. [Google Scholar] [CrossRef]
  23. Schiavon, E.; Taramelli, A.; Tornato, A.; Pierangeli, F. Monitoring Environmental and Climate Goals for European Agriculture: User Perspectives on the Optimization of the Copernicus Evolution Offer. J. Environ. Manag. 2021, 296, 113121. [Google Scholar] [CrossRef]
  24. ISTAT Popolazione Residente. Available online: https://demo.istat.it/app/?i=POS&l=it (accessed on 7 October 2021).
  25. Musolino, D.; De Carli, A.; Massarutto, A. Evaluation of the Socioeconomic Impacts of the Drought Events: The Case of the Po River Basin. Eur. Countrys. 2017, 9, 163–176. [Google Scholar] [CrossRef]
  26. Zapata-Sierra, A.J.; Zapata-Castillo, L.; Manzano-Agugliaro, F. Water Resources Availability in Southern Europe at the Basin Scale in Response to Climate Change Scenarios. Environ. Sci. Eur. 2022, 34, 75. [Google Scholar] [CrossRef]
  27. Ceppi, A.; Ravazzani, G.; Corbari, C.; Salerno, R.; Meucci, S.; Mancini, M. Real-Time Drought Forecasting System for Irrigation Management. Hydrol. Earth Syst. Sci. 2014, 18, 3353–3366. [Google Scholar] [CrossRef]
  28. Piedelobo, L.; Taramelli, A.; Schiavon, E.; Valentini, E.; Molina, J.; Nguyen Xuan, A.; Gonz-Aguilera, D. Assessment of Green Infrastructure in Riparian Zones Using Copernicus Programme. Remote Sens. 2019, 11, 2967. [Google Scholar] [CrossRef]
  29. Canone, D.; Previati, M.; Bevilacqua, I.; Salvai, L.; Ferraris, S. Field Measurements Based Model for Surface Irrigation Efficiency Assessment. Agric. Water Manag. 2015, 156, 30–42. [Google Scholar] [CrossRef]
  30. ARPA Lombardia. Sintesi Meteo-Climatica 2016. 2016. Available online: https://www.arpalombardia.it/temi-ambientali/meteo-e-clima/clima/il-clima-in-lombardia/ (accessed on 7 October 2021).
  31. Regione Lombardia—SIARL Particelle Agricole Provincia Di Lodi (Dati Generali). Available online: https://www.dati.lombardia.it/Agricoltura/Particelle-agricole-Provincia-di-Lodi-dati-general/hvwv-fgj3 (accessed on 14 November 2021).
  32. Regione Lombardia—SIARL Particelle Agricole Provincia Di Pavia. Available online: https://www.dati.lombardia.it/browse?q=Particelle%20agricole%20Provincia%20di%20pavia&sortBy=relevance (accessed on 14 November 2021).
  33. Regione Lombardia. Linee Guida per La Produzione Integrata. 2021. Available online: https://fitogest.imagelinenetwork.com/it/disciplinari/?anno=2020&idArea=3 (accessed on 14 November 2021).
  34. Regione Lombardia Uso e Copertura Del Suolo 2015 (DUSAF 5.0). Available online: https://www.geoportale.regione.lombardia.it/metadati?p_p_id=detailSheetMetadata_WAR_gptmetadataportlet&p_p_lifecycle=0&p_p_state=normal&p_p_mode=view&_detailSheetMetadata_WAR_gptmetadataportlet_identifier=r_lombar%3Acbc7ef8a-2012-4463-a197-7ffd04760bbb&_j (accessed on 20 October 2021).
  35. European Union Copernicus Land Monitoring Service. Available online: https://land.copernicus.eu/ (accessed on 1 March 2022).
  36. Donadieu, J.; L’Helguen, C. Technical Note: SENTINEL-2A L2A Products Description. 2016. Available online: https://www.theia-land.fr/wp-content/uploads/2018/12/SENTINEL-2A_L2A_Products_Description.pdf (accessed on 1 March 2022).
  37. Roberts, D.A.; Gardner, M.; Church, R.; Ustin, S.; Scheer, G.; Green, R.O. Mapping Chaparral in the Santa Monica Mountains Using Multiple Endmember Spectral Mixture Models. Remote Sens. Environ. 1998, 65, 267–279. [Google Scholar] [CrossRef]
  38. Thorp, K.R.; French, A.N.; Rango, A. Effect of Image Spatial and Spectral Characteristics on Mapping Semi-Arid Rangeland Vegetation Using Multiple Endmember Spectral Mixture Analysis (MESMA). Remote Sens. Environ. 2013, 132, 120–130. [Google Scholar] [CrossRef]
  39. Quintano, C.; Fernández-Manso, A.; Roberts, D.A. Multiple Endmember Spectral Mixture Analysis (MESMA) to Map Burn Severity Levels from Landsat Images in Mediterranean Countries. Remote Sens. Environ. 2013, 136, 76–88. [Google Scholar] [CrossRef]
  40. NV5 Geospatial Software. Available online: https://www.nv5geospatialsoftware.com/docs/MajorityMinorityAnalysis.html (accessed on 2 February 2022).
  41. Rossi, M.; Candiani, G.; Nutini, F.; Gianinetto, M.; Boschetti, M. Sentinel-2 Estimation of CNC and LAI in Rice Cropping System through Hybrid Approach Modelling. Eur. J. Remote Sens. 2022, 56, 2117651. [Google Scholar] [CrossRef]
  42. Mzid, N.; Pignatti, S.; Huang, W.; Casa, R. An Analysis of Bare Soil Occurrence in Arable Croplands for Remote Sensing Topsoil Applications. Remote Sens. 2021, 13, 474. [Google Scholar] [CrossRef]
  43. Sarmap. Running SARscape Analytics Toolbox in ENVI. 2020. Available online: https://www.sarmap.ch/index.php/sarscape-analytics-toolbox/ (accessed on 1 March 2022).
  44. European Union Regulation (EU)—1306/2013 of the European Parliament and of the Council on the Financing, Management and Monitoring of the Common Agricultural Policy. Off. J. Eur. Union 2013, L347/549, 549–607.
Figure 1. The study area is located in the Lombardy Region (Italy) within the Po plain. The red dot represents the location of the Area of Interest.
Figure 1. The study area is located in the Lombardy Region (Italy) within the Po plain. The red dot represents the location of the Area of Interest.
Land 13 00091 g001
Figure 2. Phenological calendar for the main crops in the study area.
Figure 2. Phenological calendar for the main crops in the study area.
Land 13 00091 g002
Figure 3. Satellite dataset. Study area and spatial extent of satellite data.
Figure 3. Satellite dataset. Study area and spatial extent of satellite data.
Land 13 00091 g003
Figure 4. Data pre-processing and processing chain.
Figure 4. Data pre-processing and processing chain.
Land 13 00091 g004
Figure 5. Crop class distribution in the AOI. Pie chart with the % distribution of each class and crop relative areas (in hectares).
Figure 5. Crop class distribution in the AOI. Pie chart with the % distribution of each class and crop relative areas (in hectares).
Land 13 00091 g005
Figure 6. (a) S2-based decision tree; (b) decision tree S2 + TSX based.
Figure 6. (a) S2-based decision tree; (b) decision tree S2 + TSX based.
Land 13 00091 g006
Figure 7. Intermediate crop classification map, resulting from MESMA.
Figure 7. Intermediate crop classification map, resulting from MESMA.
Land 13 00091 g007
Figure 8. Temporal profiles of (a) NDVI and (b) LAI for each class of crop.
Figure 8. Temporal profiles of (a) NDVI and (b) LAI for each class of crop.
Land 13 00091 g008
Figure 9. NDVI temporal behavior of winter cereals. The trend can be explained by the crop rotation between winter cereals (i.e., ray and durum wheat) and forage.
Figure 9. NDVI temporal behavior of winter cereals. The trend can be explained by the crop rotation between winter cereals (i.e., ray and durum wheat) and forage.
Land 13 00091 g009
Figure 10. Temporal profiles of BSI values for each class of crop. The trend lines show the behavior of tree crops (green line) compared to maize (yellow line).
Figure 10. Temporal profiles of BSI values for each class of crop. The trend lines show the behavior of tree crops (green line) compared to maize (yellow line).
Land 13 00091 g010
Figure 11. Average and standard deviation of the coefficient of variation (CoV) for each class of crop.
Figure 11. Average and standard deviation of the coefficient of variation (CoV) for each class of crop.
Land 13 00091 g011
Figure 12. Classification maps and confusion matrices. (a) classification resulting from the use of S2; (b) classification resulting from the integration of S2 and TSX.
Figure 12. Classification maps and confusion matrices. (a) classification resulting from the use of S2; (b) classification resulting from the integration of S2 and TSX.
Land 13 00091 g012
Table 1. Technical specification of sensors.
Table 1. Technical specification of sensors.
SensorSensor
Type
Acquisition ModeProcessing LevelSpatial
Resolution
Revisit TimeSpectral Ranges
Land 13 00091 i001Optical
Multispectral
-Level 2A10–60 m5 daysVIS
NIR
SWIR
Land 13 00091 i002Synthetic Aperture Radar (SAR)Stripmap mode
Single Look Complex
-3 m11 daysX band
Table 2. Satellite dataset: Sentinel 2 (left) and TerraSAR-X (right). The acquisition dates are associated with the Day of the Year (DoY). For S2, the cloud cover for each tile (T32TMR and T32TNR) is reported, together with the average value.
Table 2. Satellite dataset: Sentinel 2 (left) and TerraSAR-X (right). The acquisition dates are associated with the Day of the Year (DoY). For S2, the cloud cover for each tile (T32TMR and T32TNR) is reported, together with the average value.
SensorDataDoYCloud
Cover
T32TMR
Cloud
Cover
T32TNR
Average Cloud
Cover
Sentinel 213/01/2016136%0%3%
TerraSAR-X19/01/201619
TerraSAR-X10/02/201641
TerraSAR-X21/02/201652
TerraSAR-X03/03/201663
Sentinel 223/03/20168315%36%26%
TerraSAR-X25/03/201685
TerraSAR-X05/04/201696
Sentinel 222/04/201611325%20%23%
TerraSAR-X27/04/2016118
TerraSAR-X08/05/2016129
Sentinel 222/05/201614342%8%25%
TerraSAR-X30/05/2016151
TerraSAR-X10/06/2016162
TerraSAR-X21/06/2016173
Sentinel 201/07/201618334%37%36%
TerraSAR-X02/07/2016184
Sentinel 211/07/201619345%36%41%
Sentinel 221/07/201620350%37%44%
TerraSAR-X24/07/2016206
TerraSAR-X04/08/2016217
Sentinel 210/08/20162239%24%17%
TerraSAR-X26/08/2016239
TerraSAR-X06/09/2016250
Sentinel 209/09/201625313%7%10%
Sentinel 219/09/201626318%16%17%
TerraSAR-X28/09/2016272
Sentinel 229/09/201627338%29%34%
TerraSAR-X09/10/2016283
TerraSAR-X20/10/2016294
Sentinel 208/11/20163131%26%14%
TerraSAR-X22/11/2016327
TerraSAR-X14/12/2016349
Sentinel 228/12/20163631%0%1%
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Valentini, E.; Sapio, S.; Schiavon, E.; Righini, M.; Monteleone, B.; Taramelli, A. Development of a Pre-Automatized Processing Chain for Agricultural Monitoring Using a Multi-Sensor and Multi-Temporal Approach. Land 2024, 13, 91. https://0-doi-org.brum.beds.ac.uk/10.3390/land13010091

AMA Style

Valentini E, Sapio S, Schiavon E, Righini M, Monteleone B, Taramelli A. Development of a Pre-Automatized Processing Chain for Agricultural Monitoring Using a Multi-Sensor and Multi-Temporal Approach. Land. 2024; 13(1):91. https://0-doi-org.brum.beds.ac.uk/10.3390/land13010091

Chicago/Turabian Style

Valentini, Emiliana, Serena Sapio, Emma Schiavon, Margherita Righini, Beatrice Monteleone, and Andrea Taramelli. 2024. "Development of a Pre-Automatized Processing Chain for Agricultural Monitoring Using a Multi-Sensor and Multi-Temporal Approach" Land 13, no. 1: 91. https://0-doi-org.brum.beds.ac.uk/10.3390/land13010091

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop