Satellite remote sensing dataset of Sentinel-2 for phenology metrics extraction from sites in Bulgaria and France

Site Description:

In this dataset, there are seventeen production crop fields in Bulgaria where winter rapeseed and wheat were grown and two research fields in France where winter wheat – rapeseed – barley – sunflower and winter wheat – irrigated maize crop rotation is used. The full description of those fields is in the database "In-situ crop phenology dataset from sites in Bulgaria and France" (doi.org/10.5281/zenodo.7875440).

Methodology and Data Description:

Remote sensing data is extracted from Sentinel-2 tiles 35TNJ for Bulgarian sites and 31TCJ for French sites on the day of the overpass since September 2015 for Sentinel-2 derived vegetation indices and since October 2016 for HR-VPP products. To suppress spectral mixing effects at the parcel boundaries, as highlighted by Meier et al., 2020, the values from all datasets were subgrouped per field and then aggregated to a single median value for further analysis.

Sentinel-2 data was downloaded for all test sites from CREODIAS (https://creodias.eu/) in L2A processing level using a maximum scene-wide cloudy cover threshold of 75%. Scenes before 2017 were available in L1C processing level only. Scenes in L1C processing level were corrected for atmospheric effects after downloading using Sen2Cor (v2.9) with default settings. This was the same version used for the L2A scenes obtained intermediately from CREODIAS. 

Next, the data was extracted from the Sentinel-2 scenes for each field parcel where only SCL classes 4 (vegetation) and 5 (bare soil) pixels were kept. We resampled the 20m band B8A to match the spatial resolution of the green and red band (10m) using nearest neighbor interpolation. The entire image processing chain was carried out using the open-source Python Earth Observation Data Analysis Library (EOdal) (Graf et al., 2022).

Apart from the widely used Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI), we included two recently proposed indices that were reported to have a higher correlation with photosynthesis and drought response of vegetation: These were the Near-Infrared Reflection of Vegetation (NIRv) (Badgley et al., 2017)  and Kernel NDVI (kNDVI) (Camps-Valls et al., 2021). We calculated the vegetation indices in two different ways: 

First, we used B08 as near-infrared (NIR) band which comes in a native spatial resolution of 10 m. B08 (central wavelength 833 nm) has a relatively coarse spectral resolution with a bandwidth of 106 nm.

Second, we used B8A which is available at 20 m spatial resolution. B8A differs from B08 in its central wavelength (864 nm) and has a narrower bandwidth (21 nm or 22 nm in the case of Sentinel-2A and 2B, respectively) compared to B08.

The High Resolution Vegetation Phenology and Productivity (HR-VPP) dataset from Copernicus Land Monitoring Service (CLMS) has three 10-m set products of Sentinel-2: vegetation indices, vegetation phenology and productivity parameters and seasonal trajectories (Tian et al., 2021). Both vegetation indices, Normalized Vegetation Index (NDVI) and Plant Phenology (PPI) and plant parameters, Fraction of Absorbed Photosynthetic Active Radiation (FAPAR) and Leaf Area Index (LAI) were computed for the time of Sentinel-2 overpass by the data provider. 

NDVI is computed directly from B04 and B08 and PPI is computed using Difference Vegetation Index (DVI = B08 - B04) and its seasonal maximum value per pixel. FAPAR and LAI are retrieved from B03 and B04 and B08 with neural network training on PROSAIL model simulations. The dataset has a quality flag product (QFLAG2) which is a 16-bit that extends the scene classification band (SCL) of the Sentinel-2 Level-2 products. A “medium” filter was used to mask out QFLAG2 values from 2 to 1022, leaving land pixels (bit 1) within or outside cloud proximity (bits 11 and 13) or cloud shadow proximity (bits 12 and 14). 

The HR-VPP daily raw vegetation indices products are described in detail in the user manual (Smets et al., 2022) and the computations details of PPI are given by Jin and Eklundh (2014). Seasonal trajectories refer to the 10-daily smoothed time-series of PPI used for vegetation phenology and productivity parameters retrieval with TIMESAT (Jönsson and Eklundh 2002, 2004).

HR-VPP data was downloaded through the WEkEO Copernicus Data and Information Access Services (DIAS) system with a Python 3.8.10 harmonized data access (HDA) API 0.2.1. Zonal statistics [’min’, ’max’, ’mean’, ’median’, ’count’, ’std’, ’majority’] were computed on non-masked pixel values within field boundaries with rasterstats Python package 0.17.00.

The Start of season date (SOSD), end of season date (EOSD) and length of seasons (LENGTH) were extracted from the annual Vegetation Phenology and Productivity Parameters (VPP) dataset as an additional source for comparison. These data are a product of the Vegetation Phenology and Productivity Parameters, see (https://land.copernicus.eu/pan-european/biophysical-parameters/high-resolution-vegetation-phenology-and-productivity/vegetation-phenology-and-productivity) for detailed information.

File Description:

4 datasets:

1_senseco_data_S2_B08_Bulgaria_France; 1_senseco_data_S2_B8A_Bulgaria_France; 1_senseco_data_HR_VPP_Bulgaria_France; 1_senseco_data_phenology_VPP_Bulgaria_France

3 metadata:

2_senseco_metadata_S2_B08_B8A_Bulgaria_France; 2_senseco_metadata_HR_VPP_Bulgaria_France; 2_senseco_metadata_phenology_VPP_Bulgaria_France

The dataset files “1_senseco_data_S2_B8_Bulgaria_France” and “1_senseco_data_S2_B8A_Bulgaria_France” concerns all vegetation indices (EVI, NDVI, kNDVI, NIRv) data values and related information, and metadata file “2_senseco_metadata_S2_B08_B8A_Bulgaria_France” describes all the existing variables. Both “1_senseco_data_S2_B8_Bulgaria_France” and “1_senseco_data_S2_B8A_Bulgaria_France” have the same column variable names and for that reason, they share the same metadata file “2_senseco_metadata_S2_B08_B8A_Bulgaria_France”.

The dataset file “1_senseco_data_HR_VPP_Bulgaria_France” concerns vegetation indices (NDVI, PPI) and plant parameters (LAI, FAPAR) data values and related information, and metadata file “2_senseco_metadata_HRVPP_Bulgaria_France” describes all the existing variables. 

The dataset file “1_senseco_data_phenology_VPP_Bulgaria_France” concerns the vegetation phenology and productivity parameters (LENGTH, SOSD, EOSD) values and related information, and metadata file “2_senseco_metadata_VPP_Bulgaria_France” describes all the existing variables.

Bibliography

G. Badgley, C.B. Field, J.A. Berry, Canopy near-infrared reflectance and terrestrial photosynthesis, Sci. Adv. 3 (2017) e1602244. https://doi.org/10.1126/sciadv.1602244.

G. Camps-Valls, M. Campos-Taberner, Á. Moreno-Martínez, S. Walther, G. Duveiller, A. Cescatti, M.D. Mahecha, J. Muñoz-Marí, F.J. García-Haro, L. Guanter, M. Jung, J.A. Gamon, M. Reichstein, S.W. Running, A unified vegetation index for quantifying the terrestrial biosphere, Sci. Adv. 7 (2021) eabc7447. https://doi.org/10.1126/sciadv.abc7447.

L.V. Graf, G. Perich, H. Aasen, EOdal: An open-source Python package for large-scale agroecological research using Earth Observation and gridded environmental data, Comput. Electron. Agric. 203 (2022) 107487. https://doi.org/10.1016/j.compag.2022.107487.

H. Jin, L. Eklundh, A physically based vegetation index for improved monitoring of plant phenology, Remote Sens. Environ. 152 (2014) 512–525. https://doi.org/10.1016/j.rse.2014.07.010.

P. Jonsson, L. Eklundh, Seasonality extraction by function fitting to time-series of satellite sensor data, IEEE Trans. Geosci. Remote Sens. 40 (2002) 1824–1832. https://doi.org/10.1109/TGRS.2002.802519.

P. Jönsson, L. Eklundh, TIMESAT—a program for analyzing time-series of satellite sensor data, Comput. Geosci. 30 (2004) 833–845. https://doi.org/10.1016/j.cageo.2004.05.006.

J. Meier, W. Mauser, T. Hank, H. Bach, Assessments on the impact of high-resolution-sensor pixel sizes for common agricultural policy and smart farming services in European regions, Comput. Electron. Agric. 169 (2020) 105205. https://doi.org/10.1016/j.compag.2019.105205.

B. Smets, Z. Cai, L. Eklund, F. Tian, K. Bonte, R. Van Hoost, R. Van De Kerchove, S. Adriaensen, B. De Roo, T. Jacobs, F. Camacho, J. Sánchez-Zapero, S. Else, H. Scheifinger, K. Hufkens, P. Jönsson, HR-VPP Product User Manual Vegetation Indices, 2022.

F. Tian, Z. Cai, H. Jin, K. Hufkens, H. Scheifinger, T. Tagesson, B. Smets, R. Van Hoolst, K. Bonte, E. Ivits, X. Tong, J. Ardö, L. Eklundh, Calibrating vegetation phenology from Sentinel-2 using eddy covariance, PhenoCam, and PEP725 networks across Europe, Remote Sens. Environ. 260 (2021) 112456. https://doi.org/10.1016/j.rse.2021.112456.

站点描述：

本数据集包含保加利亚的17个生产型农田（种植冬油菜与冬小麦），以及法国的2个研究型农田，其轮作制度分别为冬小麦-油菜-大麦-向日葵，以及冬小麦-灌溉玉米。上述农田的完整说明收录于数据集《保加利亚与法国站点原位作物物候数据集（In-situ crop phenology dataset from sites in Bulgaria and France）》（doi.org/10.5281/zenodo.7875440）。

方法与数据描述：

遥感数据提取自哨兵二号（Sentinel-2）瓦片：保加利亚站点对应瓦片35TNJ，法国站点对应瓦片31TCJ，数据获取时间为卫星过境当日。其中，哨兵二号衍生植被指数的数据时段为2015年9月起，高分辨率植被物候与生产力（High Resolution Vegetation Phenology and Productivity, HR-VPP）产品的数据时段为2016年10月起。

为抑制地块边界处的光谱混合效应（Meier等，2020），所有数据集的值按农田分组，随后聚合为单一中位数值用于后续分析。

哨兵二号数据从CREODIAS（https://creodias.eu/）下载，测试站点的下载版本为L2A处理级别，设置最大场景云量覆盖阈值为75%。2017年之前的影像仅支持L1C处理级别，下载后的L1C级别影像使用Sen2Cor（v2.9）默认设置进行大气校正，该版本与从CREODIAS中间获取的L2A影像所用版本一致。

随后，从每个农田地块的哨兵二号影像中提取数据，仅保留场景分类（Scene Classification, SCL）类别为4（植被）与5（裸土）的像素。将20m分辨率的B8A波段重采样至与绿光、红光波段（10m分辨率）一致的空间分辨率，重采样方法为最近邻插值。整套影像处理流程通过开源Python地球观测数据分析库（EOdal）完成（Graf等，2022）。

除广泛使用的归一化差分植被指数（Normalized Difference Vegetation Index, NDVI）与增强型植被指数（Enhanced Vegetation Index, EVI）外，本数据集纳入两项新近提出的植被指数——据报道二者与植被光合作用及干旱响应的相关性更高：分别为植被近红外反射率（Near-Infrared Reflection of Vegetation, NIRv）（Badgley等，2017）与核化归一化差分植被指数（Kernel NDVI, kNDVI）（Camps-Valls等，2021）。植被指数的计算采用两种方式：
1. 使用原生空间分辨率为10m的B08波段作为近红外（Near-Infrared, NIR）波段。B08波段中心波长为833nm，光谱分辨率相对较粗，带宽为106nm。
2. 使用空间分辨率为20m的B8A波段。与B08相比，B8A的中心波长为864nm，带宽更窄（哨兵二号A、B版本的带宽分别为21nm与22nm）。

哥白尼陆地监测服务（Copernicus Land Monitoring Service, CLMS）提供的高分辨率植被物候与生产力（HR-VPP）数据集包含三款10m分辨率的哨兵二号产品：植被指数、植被物候与生产力参数，以及季节轨迹（Tian等，2021）。数据提供商已针对哨兵二号过境时刻，计算了植被指数（归一化植被指数（Normalized Vegetation Index, NDVI）与植物物候指数（Plant Phenology Index, PPI））以及植物参数（光合有效辐射吸收比（Fraction of Absorbed Photosynthetic Active Radiation, FAPAR）与叶面积指数（Leaf Area Index, LAI））。

NDVI直接通过B04与B08波段计算得到；PPI通过差分植被指数（Difference Vegetation Index, DVI = B08 - B04）及其逐像素季节最大值计算得到。FAPAR与LAI通过B03、B04与B08波段反演得到，反演方法为基于PROSAIL模型模拟的神经网络训练。该数据集附带质量标志产品（QFLAG2），为16位数据，用于扩展哨兵二号L2级产品的场景分类波段（SCL）。使用“中值”滤波器掩膜去除QFLAG2值为2至1022的像素，仅保留陆地区域像素（位1），无论其处于云邻近区域（位11与13）或云影邻近区域（位12与14）。

HR-VPP每日原始植被指数产品的详细说明见用户手册（Smets等，2022），PPI的计算细节见Jin与Eklundh（2014）的研究。季节轨迹指用于植被物候与生产力参数反演的PPI每10日平滑时间序列，反演工具为TIMESAT（Jönsson与Eklundh，2002，2004）。

HR-VPP数据通过WEkEO哥白尼数据与信息访问服务（Data and Information Access Services, DIAS）系统下载，使用Python 3.8.10与统一数据访问（Harmonized Data Access, HDA）API 0.2.1。使用rasterstats Python包0.17.00，针对农田边界内未被掩膜的像素值计算分区统计指标，包括['min', 'max', 'mean', 'median', 'count', 'std', 'majority']。

物候始期（Start of Season Date, SOSD）、物候终期（End of Season Date, EOSD）与季长（Length of Seasons, LENGTH）从年度植被物候与生产力参数（Vegetation Phenology and Productivity Parameters, VPP）数据集中提取，作为额外的对比数据源。上述数据为植被物候与生产力参数产品的产物，详细信息见：https://land.copernicus.eu/pan-european/biophysical-parameters/high-resolution-vegetation-phenology-and-productivity/vegetation-phenology-and-productivity。

文件描述：

本数据集包含4个数据文件与3个元数据文件：
- 数据文件：
1. "1_senseco_data_S2_B08_Bulgaria_France"
2. "1_senseco_data_S2_B8A_Bulgaria_France"
3. "1_senseco_data_HR_VPP_Bulgaria_France"
4. "1_senseco_data_phenology_VPP_Bulgaria_France"
- 元数据文件：
1. "2_senseco_metadata_S2_B08_B8A_Bulgaria_France"
2. "2_senseco_metadata_HR_VPP_Bulgaria_France"
3. "2_senseco_metadata_phenology_VPP_Bulgaria_France"

数据集文件"1_senseco_data_S2_B08_Bulgaria_France"与"1_senseco_data_S2_B8A_Bulgaria_France"包含所有植被指数（EVI、NDVI、kNDVI、NIRv）的数据值及相关信息，元数据文件"2_senseco_metadata_S2_B08_B8A_Bulgaria_France"描述了所有现有变量。由于二者的列变量名称完全一致，因此共享同一份元数据文件"2_senseco_metadata_S2_B08_B8A_Bulgaria_France"。

数据集文件"1_senseco_data_HR_VPP_Bulgaria_France"包含植被指数（NDVI、PPI）与植物参数（LAI、FAPAR）的数据值及相关信息，元数据文件"2_senseco_metadata_HR_VPP_Bulgaria_France"描述了所有现有变量。

数据集文件"1_senseco_data_phenology_VPP_Bulgaria_France"包含植被物候与生产力参数（LENGTH、SOSD、EOSD）的数据值及相关信息，元数据文件"2_senseco_metadata_phenology_VPP_Bulgaria_France"描述了所有现有变量。

参考文献：

1. Badgley G, Field C B, Berry J A. 冠层近红外反射率与陆地光合作用[J]. Science Advances, 2017, 3(3): e1602244. https://doi.org/10.1126/sciadv.1602244.
2. Camps-Valls G, Campos-Taberner M, Moreno-Martínez Á, et al. 用于量化陆地生物圈的统一植被指数[J]. Science Advances, 2021, 7(7): eabc7447. https://doi.org/10.1126/sciadv.abc7447.
3. Graf L V, Perich G, Aasen H, et al. EOdal: 一款用于基于地球观测与网格化环境数据开展大规模农业生态研究的开源Python包[J]. Computers and Electronics in Agriculture, 2022, 203: 107487. https://doi.org/10.1016/j.compag.2022.107487.
4. Jin H, Eklundh L. 一种用于改进植物物候监测的基于物理原理的植被指数[J]. Remote Sensing of Environment, 2014, 152: 512-525. https://doi.org/10.1016/j.rse.2014.07.010.
5. Jönsson P, Eklundh L. 基于函数拟合的卫星传感器时间序列季节提取[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(8): 1824-1832. https://doi.org/10.1109/TGRS.2002.802519.
6. Jönsson P, Eklundh L. TIMESAT——一款用于分析卫星传感器时间序列的程序[J]. Computers & Geosciences, 2004, 30(8): 833-845. https://doi.org/10.1016/j.cageo.2004.05.006.
7. Meier J, Mauser W, Hank T, et al. 欧洲区域高分辨率传感器像素尺寸对共同农业政策与智慧农业服务的影响评估[J]. Computers and Electronics in Agriculture, 2020, 169: 105205. https://doi.org/10.1016/j.compag.2019.105205.
8. Smets B, Cai Z, Eklund L, et al. HR-VPP Product User Manual Vegetation Indices[R]. 2022.
9. Tian F, Cai Z, Jin H, et al. 基于涡度协方差、物候相机与欧洲PEP725网络的哨兵二号植被物候校准[J]. Remote Sensing of Environment, 2021, 260: 112456. https://doi.org/10.1016/j.rse.2021.112456.