Developing Transformation Functions for VENμS and Sentinel-2 Surface Reflectance over Israel

Vegetation and Environmental New micro Spacecraft (VENμS) and Sentinel-2 are both ongoing earth observation missions that provide high-resolution multispectral imagery at 10 m (VENμS) and 10–20 m (Sentinel-2), at relatively high revisit frequencies (two days for VENμS and five days for Sentinel-2). Sentinel-2 provides global coverage, whereas VENμS covers selected regions, including parts of Israel. To facilitate the combination of these sensors into a unified time-series, a transformation model between them was developed using imagery from the region of interest. For this purpose, same-day acquisitions from both sensor types covering the surface reflectance over Israel, between April 2018 and November 2018, were used in this study. Transformation coefficients from VENμS to Sentinel-2 surface reflectance were produced for their overlapping spectral bands (i.e., visible, red-edge and near-infrared). The performance of these spectral transformation functions was assessed using several methods, including orthogonal distance regression (ODR), the mean absolute difference (MAD), and spectral angle mapper (SAM). Post-transformation, the value of the ODR slopes were close to unity for the transformed VENμS reflectance with Sentinel-2 reflectance, which indicates near-identity of the two datasets following the removal of systemic bias. In addition, the transformation outputs showed better spectral similarity compared to the original images, as indicated by the decrease in SAM from 0.093 to 0.071. Similarly, the MAD was reduced post-transformation in all bands (e.g., the blue band MAD decreased from 0.0238 to 0.0186, and in the NIR it decreased from 0.0491 to 0.0386). Thus, the model helps to combine the images from Sentinel-2 and VENμS into one time-series that facilitates continuous, temporally dense vegetation monitoring.


Introduction
We are entering a golden age of high-quality public domain earth observation (EO) data. The increase in the number of EO satellite sensors provides an opportunity to combine the images from the various sensors into a temporally dense time-series. This combination of data significantly improves the temporal resolution of EO and enhances our ability to monitor land surface changes [1,2]. Many studies have previously noted the advantageous use of information at a high spatial and temporal resolution for land cover change [3][4][5], agricultural management [6][7][8], and forest monitoring [9,10]. The availability of public-domain imagery archives on the one hand, and the development of high-performance computing systems on the other, have allowed scientists to work with large volumes of EO data for continuous monitoring and analysis of earth surface phenomena [11][12][13][14][15]. Nevertheless, this progress creates new challenges: To analyze a time-series of images acquired by different sensors, the images must undergo radiometric harmonization-i.e., the spectral differences between their corresponding bands must be minimized [16].
The most common way of integrating data from two sensors is by developing empirical transformation models [17][18][19]. Harmonizing datasets, from two different sensors, requires geometric and radiometric corrections [16]. A BDRF normalization is an important step in the radiometric correction process performed in many studies [16,[19][20][21][22]. Additionally, co-registration of both datasets is important in order to minimize misregistration issues during the comparison of the image pairs [16,18,19]. These studies also emphasize the importance of using a large and representative dataset in the creation of the transformation functions [16,19,21]. Sentinel-2 and the Landsat series, for example, are public-domain optical spaceborne sensors, predominantly used for land cover monitoring. In order to combine imagery from these sensors into a unified time series, previous studies have developed transformation functions between the different Landsat sensors [17,23], and also, between Landsat ETM+ and OLI with Sentinel-2 [16,18,19,[24][25][26]. Nevertheless, many studies have highlighted that the difference in reflectance is not only a function of the band response but also of the target pixels. When there is a time-lag between the acquisitions of the images, during which the land-cover changes, the resulting differences in spectral reflectance between the images are not only caused by the difference in the response of the sensors. Accordingly, this time-lag should be as small as possible, since a latent assumption in the development of the transformation models is that the bias is related to the sensor differences rather than land-cover change. In addition, the prevalent land-cover types in the images used in the model development will determine the model's generality. Thus, the empirical models developed to integrate data from two sensors are often region-specific and less applicable to other locations [18,19,26]. Hence, regional transformation coefficients should be derived in order to combine datasets from different sensors. This paper describes the development of band transformation functions between the VENµS and Sentinel-2 satellite sensors over Israel.
VENµS is a joint satellite of the Israeli and French space agencies (ISA and CNES). Launched in August 2017, it is a near-polar sun-synchronous microsatellite at a 98 • inclination, orbiting at an altitude of 720 km. The satellite has a two-day revisit time and the sensor covers a swath area of 27 km with a constant view angle. The VENµS sensor is a multispectral camera with 12 narrow spectral bands in the range of 415-910 nm. The surface reflectance product is provided at a spatial resolution of 10 m for all bands [27,28]. The major focus of the VENµS mission is vegetation monitoring (with an emphasis on precision agriculture applications that are expected to benefit from the red-edge bands [29]) and the measurement of atmospheric water vapor content and aerosol optical depth [30].
Sentinel-2 is an EO mission from The European Space Agency (ESA) Copernicus program. It includes two satellites, each equipped with a Multi-Spectral Instrument (MSI), namely Sentinel-2A (launched June 2015) and Sentinel-2B (launched March 2017). Both sun-synchronous satellites are orbiting the earth at an altitude of 786 km [31]. Sentinel-2A and Sentinel-2B have a combined revisit time of five days. The push broom MSI sensor has a 20.6 • field of view covering a swath width of 290 km. MSI has 13 spectral bands with varying spatial resolution, 10 m for visible (red, blue, green) and broad near-infrared (NIR); 20 m for red-edge, narrow NIR and short-wave infrared (SWIR); and 60 m for water vapor and cirrus cloud bands. VENµS and Sentinel-2 sensors produce 10-, and 12-bit, radiometric data, respectively.
VENµS has a relatively narrow view angle compared to Sentinel-2, the latter acquiring images at nadir with a wider view angle of ±10.3 • from nadir. This wider field of view may cause bidirectional reflectance effects because most of the land surface consists of non-Lambertian surfaces. For better cross-sensor calibration, bidirectional reflectance distribution effects need to be minimized. Roujean et al. [32] explained an observed the reflectance variation across the swath as the bidirectional reflectance distribution function (BRDF). Roy et al. [21] examined the directional effects on Sentinel-2 surface reflectance in overlapping regions of adjacent image tiles and concluded that the difference in reflectance due to BDRF effects may introduce significant noise for monitoring applications if the BRDF effects are not treated. Studies by Claverie et al. [33] and Roy et al. [20] reported that a single set of global BRDF coefficients has shown satisfying BRDF normalization. These global coefficients have been derived for the visible, NIR and SWIR bands [20] and the red-edge bands [22]. Claverie et al. [16] Remote Sens. 2019, 11, 1710 3 of 18 and Roy et al. [21] reported that the use of these global coefficients for the BDRF correction resulted in a stable and operationally efficient correction for Sentinel-2 data.
The above literature review suggests that a BDRF correction would be required to produce a harmonized product of VENµS and Sentinel-2 surface reflectance. Accordingly, the aim of this study was to develop a transformation model, based on near-simultaneously acquired imagery, from these sensors over Israel. The specific objectives of this study were (1) to create a harmonized (both geometrically and radiometrically corrected) surface reflectance product of VENµS and Sentinel-2 imagery by adapting protocols previously established for other sensors, (2) to derive the transformation model coefficients for the overlapping spectral bands, and (3) to assess the model performance.

Description of the Sentinel-2 and VENµS Dataset
The state of Israel is covered by seven Sentinel-2 tiles and 27 VENµS tiles ( Figure 1). As a first step, same-day acquisitions from VENµS and Sentinel-2 were inventoried. Eleven dates of near-synchronous acquisitions were found for the period from April 2018 to November 2018. In total, 77 Sentinel-2 and 230 VENµS images were used to derive the band transformation model (Table A1). Atmospherically corrected reflectance products from both sensors were used in this analysis. VENµS level-2 products were obtained from the Israel VENµS website maintained by Ben-Gurion University of the Negev (https://venus.bgu.ac.il/venus/). Sentinel-2 level-2A data were obtained from the European Space Agency's Copernicus Open Access Hub website (https://scihub.copernicus.eu/dhus/#/home). Table 1 lists the overlapping spectral bands of the VENµS and Sentinel-2 sensors and their attributes. Figure 2 illustrates the spectral response functions of the VENµS and Sentinel-2 bands. The above literature review suggests that a BDRF correction would be required to produce a harmonized product of VENμS and Sentinel-2 surface reflectance. Accordingly, the aim of this study was to develop a transformation model, based on near-simultaneously acquired imagery, from these sensors over Israel. The specific objectives of this study were (1) to create a harmonized (both geometrically and radiometrically corrected) surface reflectance product of VENμS and Sentinel-2 imagery by adapting protocols previously established for other sensors, (2) to derive the transformation model coefficients for the overlapping spectral bands, and (3) to assess the model performance.

Description of the Sentinel-2 and VENµS Dataset
The state of Israel is covered by seven Sentinel-2 tiles and 27 VENμS tiles ( Figure 1). As a first step, same-day acquisitions from VENμS and Sentinel-2 were inventoried. Eleven dates of nearsynchronous acquisitions were found for the period from April 2018 to November 2018. In total, 77 Sentinel-2 and 230 VENμS images were used to derive the band transformation model (Table A1). Atmospherically corrected reflectance products from both sensors were used in this analysis. VENμS level-2 products were obtained from the Israel VENμS website maintained by Ben-Gurion University of the Negev (https://venus.bgu.ac.il/venus/). Sentinel-2 level-2A data were obtained from the European Space Agency's Copernicus Open Access Hub website (https://scihub.copernicus.eu/dhus/#/home). Table 1 lists the overlapping spectral bands of the VENμS and Sentinel-2 sensors and their attributes. Figure 2 illustrates the spectral response functions of the VENμS and Sentinel-2 bands.  The grey-shaded regions in the overlap between Sentinel-2 tiles were used in the NBAR correction assessment (shown in Figure 4).   VENμS and Sentinel-2 bands were grouped into two categories based on spatial resolution (10 m and 20 m). Most of the bands preserved their native resolution, such that the original reflectance VENµS and Sentinel-2 bands were grouped into two categories based on spatial resolution (10 m and 20 m). Most of the bands preserved their native resolution, such that the original reflectance values were retained. However, the VENµS red-edge and NIR bands were resampled to 20 m to match the Sentinel-2 red-edge and narrow NIR bands.
The following considerations were made during the development of the transformation functions for Sentinel-2 and VENµS reflectance: (1) In order for the regression model input to be representative of the reflectance variance in Israel, large spatial and temporal coverages were considered; (2) the difference in reflectance values between the different sensors over non-lamebrain surfaces was corrected; (3) errors from defective or misregistered pixels were removed [25]. Figure 3 presents the steps in the development of the transformation model. values were retained. However, the VENμS red-edge and NIR bands were resampled to 20 m to match the Sentinel-2 red-edge and narrow NIR bands.
The following considerations were made during the development of the transformation functions for Sentinel-2 and VENμS reflectance: (1) In order for the regression model input to be representative of the reflectance variance in Israel, large spatial and temporal coverages were considered; (2) the difference in reflectance values between the different sensors over non-lamebrain surfaces was corrected; (3) errors from defective or misregistered pixels were removed [25]. Figure 3 presents the steps in the development of the transformation model. Figure 3. The processing chain for the development of transformation models between VENμS and Sentinel-2 reflectance imagery. * The VENμS red-edge bands (8-10) and near infra-red (NIR) band (11) were resampled to 20 m. ** The Sentinel-2 broad NIR and narrow NIR bands (i.e., bands 8 and 8A) were compared with the VENμS NIR band (11).
The Level-2 products used in this study were atmospherically corrected before dissemination by their respective agencies: The level-2A Sentinel-2 images were atmospherically corrected using Sentinel-2 Atmospheric Correction (S2AC) and the VENμS level-2 images were corrected using the MACCS-ATCOR Joint Algorithm (MAJA). Since these atmospherically corrected reflectance products were used to develop the transformation, atmospheric correction is not listed as a step in the process.

BRDF Correction
The BRDF correction was carried out using the c-factor technique that uses global coefficients [20,22]. Table 2 lists global coefficients that have been previously validated for Sentinel-2 and Landsat [16]. Since the VENμS bands are spectrally similar to Sentinel-2, we applied the same coefficients to the VENμS imagery. In the current work, nadir BRDF-adjusted reflectance (NBAR) values were derived for both Sentinel-2 and VENμS.
The NBAR reflectance and c-factor were calculated as:  (11) were resampled to 20 m. ** The Sentinel-2 broad NIR and narrow NIR bands (i.e., bands 8 and 8A) were compared with the VENµS NIR band (11).
The Level-2 products used in this study were atmospherically corrected before dissemination by their respective agencies: The level-2A Sentinel-2 images were atmospherically corrected using Sentinel-2 Atmospheric Correction (S2AC) and the VENµS level-2 images were corrected using the MACCS-ATCOR Joint Algorithm (MAJA). Since these atmospherically corrected reflectance products were used to develop the transformation, atmospheric correction is not listed as a step in the process.

BRDF Correction
The BRDF correction was carried out using the c-factor technique that uses global coefficients [20,22]. Table 2 lists global coefficients that have been previously validated for Sentinel-2 and Landsat [16]. Since the VENµS bands are spectrally similar to Sentinel-2, we applied the same coefficients to the VENµS imagery. In the current work, nadir BRDF-adjusted reflectance (NBAR) values were derived for both Sentinel-2 and VENµS. The NBAR reflectance and c-factor were calculated as: where ρ(λ) is the spectral reflectance for wavelength λ, θ Sensor is the actual sensor's sun-illumination geometry (i.e., angles of view zenith, sun zenith, and view-sun relative azimuth), θ nadir is the sensor's sun-illumination geometry at nadir position (when the view zenith angle equals zero), and c(λ) is the correction factor for wavelength λ. k vol and k geo are the volumetric and geometric kernels and f iso , f vol , and f geo are the constant values of BRDF spectral model parameters ( Table 2). The volumetric and geometric kernels are the functions defined by the view and sun illumination geometry [32]. A detailed explanation of these kernel functions is given in the theoretical document of MODIS BDRF/Albedo product [34]. In the Equation (2) nominator, the view zenith angle is set to nadir (zero) and the average value of the solar zenith angle is applied in order to normalize the VENµS and Sentinel-2 reflectance. This radiometric normalization addresses the difference in reflectance that is the result of BRDF effects.

Spatial Co-Registration and Cloud Masking
The VENµS and Sentinel-2 NBAR products were co-registered to a sub-pixel precision of <0.5 pixels (RMSE) using the AutoSync Workstation tool in ERDAS IMAGINE. A second order polynomial transformation model was used in conjunction with the nearest neighbor resampling method. Tie-points with significant bias were eliminated, while the remaining tie-points were well-distributed throughout the image space to ensure proper geometrical registration. Table 3 shows the number of tie-points used for co-registration with the corresponding RMSE values. Shadow and cloud contaminated pixels were masked out of the analysis by using scene quality information from the VENµS QTL file and the Sentinel-2 scene quality flags.

Transformation Models
In each VENµS image, 3000-6000 random points were generated with a minimum distance of 60m between every two points. The regression model between VENµS and Sentinel-2 reflectance was produced based on 90% of these points, while the remaining 10% were used for validation. Overall, 733,562 pixels from 175 VENµS-Sentinel-2 image pairs were processed using Ordinary Least Square (OLS) regression.
where ρ Sentinel−2 is the Sentinel-2 NBAR; ρ VENµS is the VENµS NBAR; c 0 and c 1 are the OLS regression coefficients.
Since some degree of misregistration is expected, possible outliers in the dataset were removed using Cook's distance method [35]. Cook's distance (D i ) is defined as the sum of all the changes in the regression model when observation i is removed.
whereŶ j is the jth fitted response value,Ŷ j(i) is the jth fitted response value, obtained when excluding i, MSE is the mean squared error, and ρ is the number of coefficients in the regression model. The threshold used to remove the outliers in the training dataset was three times the mean Di. Values above those thresholds were removed. To remove outliers from the validation dataset, the mean Di was set as the threshold. As a result, a higher proportion of the data was removed relative to the training data. This was done to accentuate the differences in model performance-i.e., model performance using the full validation dataset (similar to the real-world data) as compared to model performance when the outliers are removed (similar to the dataset used to train the model, but slightly more refined).
Once the outliers were excluded using Cook's distance, the final coefficients were derived based on a regression model using the remaining training pixels. These values were used to transform the VENµS reflectance in the set of validation pixels. The VENµS pixels were compared with the corresponding Sentinel-2 pixels prior to the transformation, and again post-transformation. The performance of the resulting VENµS to Sentinel-2 transformation model was subsequently assessed in three ways. First, orthogonal distance regression (ODR) was performed to assess the average proportional change between the two reflectance datasets [17]. Unlike the OLS regression, used to derive the model, the ODR slope value does not favor one variable over the other, and is only used to assess the relative divergence between the two datasets. Second, the mean absolute difference (MAD) was used to measure the difference in reflectance before and after transformation. Finally, the similarity index derived from spectral angle mapper (SAM) [36] was used to compare the reflectance values pre-and post-transformation, where smaller angle values denote higher similarity. MAD and SAM values were calculated as: where t i represents the test reflectance (i.e., reflectance values after transformation), r i denotes the reference reflectance (i.e., original reflectance values), and n represents the number of pixels considered in each band.
where t i represents the testing reflectance value of band i (i.e., reflectance values after transformation), r i denotes the reference reflectance (i.e., original reflectance values), and n represents the number of bands. ODR, MAD, and SAM were only calculated for the validation set of pixels. A hypothetical perfect agreement between two sensors is expected to produce a MAD value of zero, an ODR slope of one and a SAM of zero.

Results
A BRDF correction was carried out for the co-acquired VENµS and Sentinel-2 imagery. The NBAR products were evaluated by comparing the Sentinel-2 NBAR reflectance images in the overlapping zones of adjacent tiles marked as grey-shaded regions in Figure 1. The mean absolute difference value for all the pairs of Sentinel-2 NBAR products was derived. To quantify the performance of the BRDF correction, the mean absolute difference found in this study was compared against previously reported values for Sentinel-2 [21,22] and our correction showed better performance ( Figure 4). Accordingly, the uncertainty of this correction for Israel is significantly smaller than for the areas where this correction was originally tested.

Results
A BRDF correction was carried out for the co-acquired VENμS and Sentinel-2 imagery. The NBAR products were evaluated by comparing the Sentinel-2 NBAR reflectance images in the overlapping zones of adjacent tiles marked as grey-shaded regions in Figure 1. The mean absolute difference value for all the pairs of Sentinel-2 NBAR products was derived. To quantify the performance of the BRDF correction, the mean absolute difference found in this study was compared against previously reported values for Sentinel-2 [21,22] and our correction showed better performance ( Figure 4). Accordingly, the uncertainty of this correction for Israel is significantly smaller than for the areas where this correction was originally tested.
. Figure 4. Mean absolute difference between Nadir BRDF Adjusted Reflectance (NBAR) in the overlapping regions of Sentinel-2 imagery (shown in Figure 1) during the months of April and August, which represent the spring and the summer (high vegetation coverage in spring vs. low vegetation coverage in summer). The values found in this study are compared against the values reported in Roy et al. [21,22] for the month of April. RE denotes Red-edge.
The VENμS and Sentinel-2 NBAR products were co-registered with acceptable precision as shown in Table 3. In total, 733,562 training points and 89,198 validation points were randomly distributed over the imagery footprint. An example of outlier removal is shown for the green band reflectance in Figure 5. Cook's distances that were higher than the threshold of three times the average distance were treated as outliers and excluded from the regression model ( Figure 5D). The scatter plots in Figures 5A,C demonstrate the effect of outlier removal on the data distribution in the scatter plot. Table 4 presents the training and validation datasets by the band. The number of outliers for each band is slightly different because the Cook's distance distribution is slightly different for each band.  Figure 1) during the months of April and August, which represent the spring and the summer (high vegetation coverage in spring vs. low vegetation coverage in summer). The values found in this study are compared against the values reported in Roy et al. [21,22] for the month of April. RE denotes Red-edge.
The VENµS and Sentinel-2 NBAR products were co-registered with acceptable precision as shown in Table 3. In total, 733,562 training points and 89,198 validation points were randomly distributed over the imagery footprint. An example of outlier removal is shown for the green band reflectance in Figure 5. Cook's distances that were higher than the threshold of three times the average distance were treated as outliers and excluded from the regression model ( Figure 5D). The scatter plots in Figure 5A,C demonstrate the effect of outlier removal on the data distribution in the scatter plot. Table 4 presents the training and validation datasets by the band. The number of outliers for each band is slightly different because the Cook's distance distribution is slightly different for each band.   OLS regression was performed to derive the transformation function between corresponding VENµS and Sentinel-2 bands in the visible, red-edge and NIR spectral regions. The transformation coefficients are presented in Table 5. A gradual decrease in slope, in conjunction with an increase of the intercept values, can be observed from the blue to the NIR region. Li et al. [37] also observed a similar trend for Landsat 8 to Sentinel-2 transformation. However, this trend does not appear in other Landsat and Sentinel-2 studies [16,18,19]. Similar to the OLS regression, the pre-transformation ODR slopes in Figure 6A were in the range of 0.77 (NIR-865 nm) to 1 (blue). While the plurality of the data prior to the transformation (represented in yellow in Figures 6 and 7) is centered close to the identity line, some scatter is observed in all of the bands. Figure 6B shows the scatter plots of all bands after applying the transformation by using the coefficients in Table 5. The post-transformation values of ODR slopes were all closer to 1, indicating that VENµS reflectance was transformed to become closer to Sentinel-2 reflectance. Thus, this transformation removed part of the systemic bias that is caused by the differences between the sensors.
The model performance was assessed for the full validation dataset following the outlier removal. Figure 7 shows the marginal improvement of the ODR slope values following the removal of outliers from the validation dataset. By lowering the Cook's distance threshold even more, a higher proportion of data points was removed as outliers from the validation dataset compared to the training dataset. This was done to accentuate the differences between the full validation dataset and the remaining data after outlier removal. Despite the emphasis on these differences, the model coefficients did not change significantly. The ODR slopes prior to outlier removal and post-outlier removal, presented in Figures 6B and 7B, respectively, show a minor change of less than 0.05. Accordingly, the coefficients given in Table 5 are expected to perform well for real-world data that contain some outliers. Remote Sens. 2019, 11, x FOR PEER REVIEW 12 of 19  The model performance was assessed for the full validation dataset following the outlier removal. Figure 7 shows the marginal improvement of the ODR slope values following the removal of outliers from the validation dataset. By lowering the Cook's distance threshold even more, a higher proportion of data points was removed as outliers from the validation dataset compared to the training dataset. This was done to accentuate the differences between the full validation dataset and the remaining data after outlier removal. Despite the emphasis on these differences, the model coefficients did not change significantly. The ODR slopes prior to outlier removal and post-outlier removal, presented in Figures 6B and 7B, respectively, show a minor change of less than 0.05. Accordingly, the coefficients given in Table 5 are expected to perform well for real-world data that contain some outliers. The MAD values between the pre-transformed VENµS and Sentinel-2 reflectance show an increasing trend as a function of the wavelength that ranges from 0.024 (Blue) to 0.049 (NIR) for the full dataset, and 0.020 (Blue) to 0.041 (NIR) post-outlier removal. This MAD decreased post-transformation to 0.019 (Blue) and 0.039 (NIR) for the full dataset, and 0.015 (Blue) to 0.029 (NIR) post-outlier removal (Figure 8). It is evident that the transformation reduces the MAD, whether the outliers are removed or not, but that the removal of outliers further reduces the MAD. The SAM angle value for the post-transformation decreased, relative to the pre-transformed VENµS and Sentinel-2 reflectance (Figure 9). This indicates that our model transformation function increased the spectral similarity between the reflectance spectra. Therefore, the transformation developed in this paper seems to decrease systematic bias due to sensor differences, while outlier removal seems to decrease the differences by removing other sources of variation. These include atmospheric conditions that were not completely corrected, residual BRDF effects that remain untreated by the constant coefficients used in the C-factor method, and misregistration between the images.
increasing trend as a function of the wavelength that ranges from 0.024 (Blue) to 0.049 (NIR) for the full dataset, and 0.020 (Blue) to 0.041 (NIR) post-outlier removal. This MAD decreased posttransformation to 0.019 (Blue) and 0.039 (NIR) for the full dataset, and 0.015 (Blue) to 0.029 (NIR) post-outlier removal (Figure 8). It is evident that the transformation reduces the MAD, whether the outliers are removed or not, but that the removal of outliers further reduces the MAD. The SAM angle value for the post-transformation decreased, relative to the pre-transformed VENμS and Sentinel-2 reflectance (Figure 9). This indicates that our model transformation function increased the spectral similarity between the reflectance spectra. Therefore, the transformation developed in this paper seems to decrease systematic bias due to sensor differences, while outlier removal seems to decrease the differences by removing other sources of variation. These include atmospheric conditions that were not completely corrected, residual BRDF effects that remain untreated by the constant coefficients used in the C-factor method, and misregistration between the images.

Discussion
In this study, we developed a transformation model for VENμS and Sentinel-2 sensors over Israel. The new model coefficients provide an opportunity to use observations at high temporal resolutions for land surface monitoring by combining Sentinel-2 and VENμS observations. The availability of the red-edge spectral bands in both sensors is significant for precision agriculture applications like irrigation management [38] and LAI assessment [29]. The same-day VENμS and Sentinel-2 image pairs were acquired during April 2018 to November 2018 and used to calibrate the transformation model. A total of eight spectral bands, namely three bands from the visible region, three red-edge bands and two from the NIR region showed an increased spectral similarity post-

Discussion
In this study, we developed a transformation model for VENµS and Sentinel-2 sensors over Israel. The new model coefficients provide an opportunity to use observations at high temporal resolutions for land surface monitoring by combining Sentinel-2 and VENµS observations. The availability of the red-edge spectral bands in both sensors is significant for precision agriculture applications like irrigation management [38] and LAI assessment [29]. The same-day VENµS and Sentinel-2 image pairs were acquired during April 2018 to November 2018 and used to calibrate the transformation model. A total of eight spectral bands, namely three bands from the visible region, three red-edge bands and two from the NIR region showed an increased spectral similarity post-transformation.
The broad and narrow NIR bands of Sentinel-2 MSI were compared against the VENµS NIR band. Even though the narrow NIR has a better spectral overlap with the VENµS NIR band (Figure 2), it is not very different from the Sentinel-2 broad NIR (Figures 6-8). However, Claverie et al. [16] and Li et al. [37] highlighted that the Sentinel-2 narrow NIR band has shown better performance than has the broad NIR band. Flood [18] highlighted that the difference in reflectance values of Landsat 7-ETM+ and 8-OLI for the Australian landscape was smaller than for the entire globe. In our study, more than 60 percent of the data were non-vegetated surfaces (Table 4). Thus, the transformation function developed in this study is expected to perform better for barren surfaces than for vegetation. This model can be applied for landscapes that are similar to those of Israel, i.e., Mediterranean regions. Its applicability for different environments warrants further examination.
One strength of this study is the use of near-synchronously acquired VENµS and Sentinel-2 image pairs. This minimizes changes to the land-surface, the atmosphere and sun position and, therefore, reduces any bias between the image pairs that may be caused by the temporal delay between the acquisitions of the pair of images. Accordingly, the differences between the images can largely be attributed to systematic sensor bias rather than actual changes to the ground leaving reflectance, which is expected to be minimal [18]. In this respect, this study presents an advantage over studies conducted using image pairs taken a few days apart [16,23,37].
The dissemination of a 5 m native resolution VENµS Level-2 product is expected in the near future. Mandanici and Bitelli [26] pointed out that a difference in spatial resolution can also affect the transformation function. Hence, the transformation function developed here would need to be tested to see if the change in the products' spatial resolution has an effect, especially for the red-edge bands (VENµS-5 m to Sentinel-20 m). In addition, hyperspectral high spatial resolution imagery can be used to assess the influence of spectral resolution differences, as suggested by Claverie et al. [16] and Zhang et al. [19].

Conclusions
A first-of-its-kind cross-sensor calibration study for VENµS and Sentinel-2 surface reflectance data for Israel is presented. An effective processing chain that considers radiometric and geometric corrections was employed to derive the cross-sensor surface reflectance transformation model. Post-transformation, the ODR slopes were close to unity, the spectral similarity has increased as demonstrated by a reduction of the SAM value from 0.093 to 0.071, and the MAD between VENµS and Sentinel-2 reflectance was substantially decreased in all bands. This indicates that the models presented here can successfully be used to create a dense time-series of VENµS and Sentinel-2 imagery. The combined dataset of VENµS and Sentinel-2 provides high-frequency multispectral imagery that can support crop and vegetation monitoring studies, with the added advantage of red-edge bands that are absent from veteran sensors such as the Landsat series. Funding: This study received support from the Ministry of Science, Technology, and Space, Israel, under grant number 3-14559.