1.

INTRODUCTION

Although there are many definitions of phenology, from a scientific perspective, phenology is the science that studies the phases of the life cycle of living organisms (plants or animals)¹ and how the seasonal and interannual variations of climate affect them^2,3. Traditionally, farmers or naturalists and some scientists have recorded phenological observations in paper for decades. This practice is part of a more generic movement of involvements of citizens in data capture that is currently referred as citizen science (CS). So far, phenology in-situ observations were reduced to small traditional elements closely located to an observer’s home, providing detailed information at species scale and the volunteers’ efforts were a bit biased only representing a small demarcated area in a very precise and systematic way, as they were usually more likely to collect data from accessible locations (close to the roads or urban areas)⁴.

Indeed, this was a modest localized approximation to cover the strong variability of the vegetation phenology across biomes and to improve the knowledge of vegetation changes as climate change impacts. However, the relevance of phenology as a key indicator and regulator of the impacts of climate change is indisputable and widely recognized through multiple research approaches ^5–7.

Given that the phenological changes correspond to specific vegetation physiological processes⁸, changes in historical seasonality of vegetation phenology (leaf unfolding and flowering of plants in spring, fruit ripening, color changing, and leaf fall in autumn, among others) are considered a very sensitive indicator to demonstrate the local effects of global climate change in terms of ecosystem response⁹. This is why the collection of these observations has become a topic of major relevance. However, this in-situ approach is facing some limitations. On one hand, the strong variability of the vegetation phenology across biomes⁶ (Richardson et al., 2013) requires having more data to improve the knowledge about these changes. Additionally, the current characteristics of the CS approach (traditional species and reduced coverage) are limiting the volunteers’ aims which currently need a better guidance to identify phenological active areas of interest. This situation opens a new challenge for scientists which should develop an approach for guiding the observers in achieving their objective mainly by expanding the phenology report sites. Indeed, remote sensing (RS) has become the optimal technology to overcome limitations in spatial coverage thanks to the several available satellites using different sensors^10–14. Try-outs at large scale have been made by using data from the Advanced Very High Resolution Radiometer (AVHRR) sensor since 1981¹⁵ at 1 and 8 km spatial resolution. Unfortunately, at this low spatial resolution the behaviour of individual species of plants and their particular phenological events are almost undetectable. When moving to medium spatial resolution sensors, one can find sensors, such as MODIS, with more appropriated spectral and spatial resolutions (13 bands, and pixel size from 250 to 1000m) to detect spatial variation during the phenological life cycle. MODIS distributes a Global Land Cover dynamics product⁹ that provides 8-day information related to land surface phenology that has been used in many phenology studies over large geographic regions^8,16–18. Nevertheless, in many areas of the Earth, MODIS spatial resolution (250 m at best) is too coarse to account for small heterogeneous landscape¹³. Alternatively, Landsat legacy sensors present an optimal spectral and spatial resolution of 8 bands and 30 m respectively. In fact, Landsat was the first sensor used to characterize vegetation seasonality at regional scales¹⁹, but its nominal 16-day temporal resolution hampered by clouds might be too low to: a) generate a regular and dense enough time series once cloud cover has been masked and; b) detect some phenological changes such as blossoming that occur quickly over a time period lower than 16 days. Indeed, temporal resolution of remote sensing data will determine the detail of monitored seasonal variation in vegetation dynamics. Fortunately, the Sentinel-2 constellation, consisting in Sentinal-2A and -2B, presents some outstanding characteristics, such as a high spatial (up to 10 m), a high spectral (13 bands from VNIR to SWIR) and a high temporal (almost only 5 days considering Sentinel-2A and Sentinel-2B) resolutions allowing monitoring vegetation phenology at an unprecedented detail as recognized by the European Commission²⁰. Therefore, RS has the potentiality to become a great ally for phenology report in-situ volunteer guidance by making interoperable what is both measured from space and from human visual observations.

It is also worth mentioning, that the result from RS monitoring may greatly vary depending on the remote sensing data extraction method used on the analysis, although some evaluation approaches have been already published²¹. The ingestion of Sentinel-2A and –2B data in an environment that facilitates time series analysis is essential in the context of phenology monitoring of large areas. Indeed, the use of high resolution satellite data implies an increase in the number of files and/or its size, something that may be an issue when the study area is particularly large. Moreover, the analysis of this type of data requires a considerable amount of time spent on data preparation and organization, including data download, radiometric correction²², computation of vegetation derived products²³, as well as a methodical organization of large time series of imagery to facilitate data management and optimization for further data analysis. Since March 2018, Sentinel-2 data provider, ESA, has made a remarkable effort to facilitate the use of its data by distributing it at Level 2A (Bottom of atmosphere reflectance - BOA) and can be downloaded from the Copernicus Open Access Hub²⁴. The downloading step is facilitated through an API that facilitates a regular image update²³. Visualization of Sentinel-2 data is possible through the Sentinel Hub Earth Observation browser but data analysis options are only available for a fee. Regarding the analytical part, phenology analysis functions based on the extraction of long time series, requires visiting a subset of the long list of files for retrieving results such as a time profile of a single position in the space²³, a process that demands to have data well organized and in a way that facilitates spatio-temporal queries. Therefore more is needed.

This study presents a methodological design for harmonizing phenological vegetation monitoring and observation protocols from RS and CS carried out in Catalunya (32 000 km² region located at North-East of Spain) by 1) developing a proper system that ingests Sentinel-2 data in an Open Data Cube environment, the Catalan Data Cube, and; 2) by performing a detailed analysis of the study area to identify and map potential phenological and homogeneous active areas. Maps are then delivered to volunteers through a WMS/WFS service guiding them to targeted areas. Volunteers use an interface properly linked to iNaturalist²⁴, one of the world’s most popular nature biodiversity observation apps, so facilitating recording and sharing their own observations and also allowing users to report feedback via NiMMbus²⁵, the OGC Geospatial User Feedback (GUF). By creating this interoperable approach, in the context of PhenoTandem project (CSEOL initiative funded by ESA), the scientific community also ensures the availability of high-quality research data to be correlated with RS data.

2.

STUDY AREA

This study design was carried out in Catalunya, a region covering 32 000 km², located on the North-East of the Iberian Peninsula (SW Europe). Both, wide altitudinal range (from 0 to more than 3 000 m) and heterogeneous biodiversity landscape, composed of agriculture land (herbaceous and orchards) and mixed with areas of natural and semi-natural vegetation, including coniferous, sclerophyllous and deciduous forest makes it a unique spot for testing an interoperable approach between CS and RS phenology monitoring. Additionally, in Catalunya a wide and strong community of citizen science observers exist, in the FenoCat network (Phenologic Network of Catalonia)²⁶ and the RitmeNatura Citizen Science Observatory²⁷, thus facilitating the engagement of volunteers to participate in using our protocol.

3.

THE CATALAN DATA CUBE

In order to facilitate the access and the organization of Sentinel-2A and -2B atmospheric corrected product (S2L2A dataset) in an optimal way for performing time series analysis, an Open Data Cube (ODC)^28,29 was built. The ODC is an approach that helps in handling multidimensional cube of data (covering x, y and t dimensions). The ODC provides an optimal organization of data by indexing it enabling easy actualization of the cube with new scenes corresponding to new dates. Moreover, the ODC facilitates retrieval of a given subset of data through an API written in Python. The Catalan Data Cube is a local installation of the software that allows full control of data sources and the processes and algorithms applied to them²³ and is freely available at http://datacube.uab.cat/cdc/. Among its features, it automatically downloads Sentinel -2A and -2B L2A data from the ESA Sentinel Data hub covering the whole study area. Later, the data is ingested creating a tile structure of mosaicked granules for a single date in a single time slice, thus, reducing fragmentation of scenes. The tiles series started in March 2018. Subsequently, this infrastructure was used for the analysis of time series and the extraction of phenological profiles of some areas of interest. A Map browser on top of the ODC that used a WMS service that responds binary arrays have been also developed. This is used as a data source for the visualization portal that is described in Section 5 of the present document.

4.

IDENTIFICATION OF PHENOLOGICAL AREAS OF INTEREST FROM REMOTE SENSING

With the aim to guide the volunteers to phenologically active areas of interest in Catalunya, an initial set of habitats with a dominant species presenting different physiognomic structure (forest, shrublands or grassland), including natural habitats as well as agricultural crops, was identified. This first selection was done taking into account the phenologic interest of the species according to botanical expert-knowledge and the requirement that exists sufficient homogeneous patches to observe these species, both at the citizen science and remote sensing level.

The identification process consisted of two steps. The first step was based on the combination of existing vegetation and land cover cartography in order to identify locations with strong evidence of a given species or habitat. Different approaches depending on the physiognomic structure have been followed according to the databases availability.

4.2

Prioritization of areas based on Senintel-2 time series analysis

The second part consists of the filtering of these areas by using time series analysis of remote sensing vegetation indices with the objective to exclude areas with low indices variability, so facilitating identification of phenophases by observers.

4.2.1

Extraction of time-series of Sentinel-2 reflectance data

Vegetation indices have been traditionally used as indicators of the state and conditions of vegetation and have been established to be as sensitive enough for phenology monitoring^35,36. Usually, most of the phenology studies are based on proxies related to photosynthetic activity, easily associated with green-up, maturity or senescence phenophases. However, flowering is an essential and very remarkable phenophase for some species. For instance, fruits and nuts trees are flowering before greening. Therefore, in this study, two different indices were used to capture the green-up and the flowering depending on the species. The Normalized difference Vegetation Index (NDVI), widely recognized as the most commonly used remote sensing index for phenology monitoring was selected as a proxy of phenophases related to green-up mainly in forests and shrub species. Moreover, the Enhanced Bloom Index (EBI)³⁷ was used as a proxy for flowering phenology monitoring mainly in crops and orchards.

To extract the spatio-temporal data from a Sentinel-2 dataset from the Catalan Datacube a script in Python was coded using the Open Data Cube libraries, and selected areas were converted to a GeoJson file. Once reflectance data was extracted for 2018 and 2019 it was aggregated in 5-day composites accounting for the quality control flags and selecting a maximum value for the time interval, thus trying to minimize cloudy pixels. This was meant to optimize the satellite coverage for each time interval, reducing noise from cloud contamination and regularizing the time series into a single synthetic year^38,39.

4.2.2

Computation of vegetation indices and time-series analysis

Vegetation indices have been traditionally used as indicators of the state and conditions of vegetation and have been established to be as sensitive enough for phenology monitoring^35,36. Usually, most of the phenology studies are based on proxies related to photosynthetic activity, easily associated with green-up, maturity or senescence phenophases. However, flowering is an essential and very remarkable phenophase for some species. For instance, fruits and nuts trees are flowering before greening. Therefore, in this study, two different indices were used to capture the green-up and the flowering depending on the species. The Normalized difference Vegetation Index (NDVI), widely recognized as the most commonly used remote sensing index for phenology monitoring was selected as a proxy of phenophases related to green-up mainly in forests and shrub species. Moreover, the Enhanced Bloom Index (EBI)³⁷ was used as a proxy for flowering phenology monitoring mainly in crops and orchards.

Once the time series homogenized, NDVI and EBI temporal profiles for each area were computed and analysed using TIMESAT^40,41 version 3.2 fitting the time series to a logistic function. Additionally, weights were applied to the time series analysis to ensure remaining clouds were screen out. A value of 0 was assigned when there was a cloud high probability and a value of 0.75 when the probability was low or intermediate. For each time series, the beginning, the end and the length of season, the base value, the time of the middle of the season, the maximum value, the seasonal amplitude, the small integrated value, and the large integrated value were computed using TIMESAT. To screen out the priority areas, the amplitude of the area used was defined as the difference between the maximum value and the base level, applying the following criteria: a) only those areas presenting an average flowering amplitude signal for the selected period > 10% or, b) an averaged greening-up amplitude signal for the selected period > 20%. Only areas that met one of these criteria were selected. These criteria ensure the selection of those areas that would be easily spotted by volunteers as in there, the green-up or flowering signal is maximized. Figure 1 shows a time series of Sentinel-2 data fitted to a logistic distribution for two areas: beech (Fagus sylvatica – panels A and C) and almond tree orchards (Prunus dulcis – panels C and D). Panels A and B show maximum NDVI and EBI seasonal amplitudes for beech and almond tree orchards, respectively, being optimal sites for volunteers to gather phenology data. On the other hand, panels C and D show lower variability being discarded areas for volunteers’ phenology monitoring. In particular, panel C is a self-cast shadow beech area in which correlating volunteers’ data and remote sensing data might lead to misinterpretation due to the high amount of unsuitable reflectance values caused by the lack of a topographic correction. Moreover, panel D seems to be a less active almond tree orchard area compared to panel B area. Finally, as a result of the criteria imposed, 50% and 10% of the areas considered to be flowering-up or greening-up were finally discarded, yielding around 32 000 areas along the study area which the volunteers may visit for phenophase reporting (Figure 2).

Figure 1.

Time series for beech and almond tree orchard sites presenting: high amplitude (panels A and B) and low amplitude (panels C and D). Red small and big circles as the applied weights, red line is the applied logistic fit and blue circle is the vegetation index from Sentinel-2 data

Figure 2.

Final selection of the greening and flowering sites (in red) based on the time series amplitude for the study area.

5.

INTEROPERABLE VISUALIZATION PORTAL FOR CITIZEN SCIENCE

The Catalan Data Cube is supplemented with a web map browser interface that provides an easy a friendly way to visualize the structure and the content of the data. The portal is based on MiraMon technology and is using international standards and is composed by a client site viewer and analytical tool written in JavaScript and a server side service written in C. Both elements communicate using Open Geospatial Consortium Standards such as Web Map Server (ISO 19128) and Web Coverage Service (WCS). Even if based on an existing standard, the approach is unique because the WMS server provides binary arrays of values instead of static pictures. Those arrays are transmitted and processed in the client to generate personalized visualization and some analysis on demand. Any other WMS compatible client is able to access the same data in a more classical approach.

A topographic map and areas of interest selected from the combination of existing WMS of the Institut Cartogràfic i Geològic de Catalunya (ICGC) cartography is used as a background of the map browser. The map also show the distribution of the areas selected that are considered homogeneous to be detected by remote sensing as a single habitat or species and volunteers can use the tool to find the areas that are closer to their location (Figure 3). The interface has also an integrated compass function that supports the orientation and positioning of user, when GPS of the device used is activated. This way, citizen volunteers are guided to phenological active and interesting areas other than the typical spots there are used to observe, and their observations acquire a new dimension. By clicking in one of the spots, users see the coordinates and the species dominant in the spot and also a link redirecting them to the data capture interface. By activating the animation functionality, they are able to get a view of the Sentinel-2 time series data as well as the temporal profile of, for instance, the NDVI time series that they can correlate with their past in-situ observations.

Figure 3.

Interoperable visualization portal for guiding citizens to areas of phenological interest.

Another of the elements needed by the users is a Citizen Science capture application that provides them with the capability to create and register observations in the field. In this study, we created a project in iNaturalist²⁴ and included the necessary species-specific phenological phases to observe as an addition field. iNaturalist generic interface might be too generic and we developed a simpler capture interface based on 4 steps: 1) upload a picture, 2) select a species, 3) select the corresponding phenophase and 4) report the position of the observation. The capture app is built on top of iNaturalist API that allows the capture app to use the iNaturalist user account system and to submit observations to the iNaturalist database. This way, during the field campaigns user can register its observation in the platform directly from the web map browser (in this case, the PhenoTandem website) giving the user a fluid experience in the visualization portal.

Another component of the data capture system is the integration with the NiMMbus²⁵ service elaborated in the H2020 NextGEOSS project. NiMMbus is a widget to provide user feedback items specified in the OGC Geospatial User Feedback (GUF) standard such as, comments, ratings, discovered issues, quality reports etc. The NiMMbus system is already integrated if the map browser and can be also integrated in any website using a widget. Internally, the widget makes use of a specific WPS service. Volunteers participating in the data capture campaigns are able to comment on an area of interest, remote sensing products, the observation protocol, and even suggest new relevant areas that have not been initially detected by the remote sensing algorithms. All feedback is declared public by default allowing volunteers to see users’ feedback from other users and complement about it.

6.

RESULTS AND CONCLUSIONS

Phenology monitoring has become one of the main indicators of climate change impacts on vegetation. Traditional observations registered by citizens are limited due to its modest localized approximation and lack of guidance to interesting phenological active areas. Satellite remote sensing observations, especially those presenting high spatial and temporal resolutions, have the potential to spatio-temporally characterize vegetation dynamics, so phenology. In this study, a remote sensing based protocol to identify phenological areas of interest where to direct volunteers’ was developed. Existing cartography combined with Sentinel-2 time series analysis has been demonstrated as an effective methodology to support phenology monitoring by citizen science. Temporal profiles of vegetation indices related to greening or blossoming and its amplitude signal are good variables to filter potential phenological active areas of interest. With this method, around new 5 000 areas presenting phenological changes have been identified. This new approach increases significantly the phenological active areas of interest to be visited by volunteers’ and also the list of species to observe. Therefore, existing observation protocols have been actualized to include new target species. The achievements in the use of this methodology are twofold: 1) citizen science volunteers’ aiming at registering phenological changes are now guided to areas potentially active; 2) the quality of in-situ data useful for remote sensing studies validation increases significantly while the quantity may also increase. Access to data was prepared by implementing an Open Data Cube. The ODC offers a solution for storing big data products in an efficient way with a modest hardware and provides an easy access to data by indexing it. As a result of this study, this technology is recommended for organizing and processing remote sensing products acquired by the main satellite platforms. Moreover, the web map browser implemented allows access data to regular users other than expert users with knowledge on ODC API by using an outstanding technology that works with multidimensional binary arrays of data allowing user to generate time series analysis or animations without any server intervention. The open source code of the MiraMon web map browser presented here can be accessed from: https://github.com/joanma747/MiraMonMapBrowser. Finally, the protocol described in this study is adaptable in any region and is of general use.