Relationship between MRPV Model Parameters from MISRL2 Land Surface Product and Land Covers: A Case Study within Mainland Spain

Abstract

In this study, we showed that the multi-angle satellite remote sensing product, MISR L2 Land Surface (MIL2ASLS), which has a scale of 1.1 km, could be suitable for improving land-cover studies. Using seven images from this product, captured by the multi-angle imaging spectroradiometer sensor (MISR), we explored the values reached by the three parameters (ρ0, Θ, and k) of the Rahman–Pinty–Verstraete model, which was modified by Martonchick (MRPV). Thereafter, we compared the values and behaviors shown in seven Co-ordination of Information on the Environment (CORINE) land cover categories, in the red and near infrared (NIR) bands, over the seven MISR orbits captured in 2006 for Mainland Spain. Furthermore, we used Normalized Difference Vegetation Index (NDVI), Leaf Area Index (LAI), and Fraction of Photosynthetically Active Radiation (FPAR) ancillary data and the illumination angles from the same pixels, which made up the images. These ancillary data were also provided by the MISR products. An inferential statistic test was performed to evaluate the relationship between each parameter–band combination, and the land cover in every MISR orbit used. The results suggested that the ρ0 parameters of this product seemed to be the most related to photosynthetic activity, and it should be comparable with the widely-used NDVI. On the other hand, the k and Θ parameter values were not related, or at least not entirely related, to the phenology of land coverage. These seemed to be more influenced by the anisotropy behavior of the studied land cover pixels. Additionally, we observed, by constructing analysis of variance, how the mean of each MRPV parameter–band differed statistically (p < 0.01) by land covers and orbits. This study suggested that the MISR MRPV model parameter data product has great potential to be used to improve land cover applications.

1. Introduction

For many years now, there has been a need to retrieve both spatial and vertical structural information from remote sensing in order to provide critical inputs for ecosystem and process models [1]. Descriptions of canopy structure are vital to obtain accurate radiative transfer modeling joined to water and energy balances, which are usually a poorly-modeled element of global ecosystem models [2]. Currently, remote sensing is the only realistic way to quickly obtain land surface information across wide regions and with a broad range of temporal scenes. Traditionally, sources of remote sensing that have been used the most have been multispectral signatures. However, non-photosynthetic canopy elements represent major biomass components, and the vertical distribution and structural density can vary extensively under canopies with similar multispectral signatures [3].

Quantitative description of vegetation structure has been proposed as one of the key requirements for a significant improvement in modeling the carbon cycle and global land characterizations [4]. The use of multi-spectral data, acquired during the growing season, is the most commonly-used approach for mapping vegetation and land cover. However, this method is usually limited, as it is only capable of indirectly capturing canopy height [4]. Therefore, active remote sensing technology, such is airborne Light Detection and Ranging (LiDAR) sensors , or the Geoscience Laser Altimeter System (GLAS) in satellite platforms are nowadays popular techniques for examining the structure of vegetation canopies. In this sense, there are several studies based on these techniques, and on different areas around the Earth, and within various fields of work. For instance, there exist studies that prove the utility of LiDAR metrics for modeling biomass and improving forest inventories [5,6,7], or their use to estimate urban tree damage in an individual tree estimation approach [8], or even to study plant phenotyping in the agricultural field [9]. There are also examples where researchers have fused LiDAR acquisitions with radar or spectral images. One of these is described in Reference [10], where the authors mapped and estimated forest area and aboveground biomass in the “Miombo Woodlands” in Tanzania using data from airborne laser scanning (ALS), interferometric synthetic aperture radar (InSAR) derived from TanDEM-X, RapidEye optical imagery, and global forest map products derived from Landsat and ALOS PALSAR L-band radar imagery, which yielded improvements in terms of precision in the estimations. However, the high cost of airborne LiDAR data collection over large areas ensures that much of the terrestrial Earth surface will not be covered by airborne acquisitions in the near future, and also makes it difficult for easy regular updates. Although LiDAR in a satellite platform could be a solution, it still cannot provide comprehensive spatially-contiguous LiDAR datasets. Some studies have demonstrated the utility of multiangular reflectance data for mapping vegetation structural characteristics [4,11,12,13]. The multiangular approach could provide potential in forest height mapping over large regions at a lower cost than LiDAR technologies, thereby allowing for more regular updates.

From the first data obtained by a satellite in the 1970s, it was possible to evaluate the fact that the reflectance of different land covers was dependent on observation and illumination angles. That dependence implied that land covers were not isotropic; that is, they did not reflect the same amount of energy in all directions. Therefore, the source of noise must be corrected [14]. In those years, the Bidirectional Reflectance Function (BRDF) was defined as a variable that explained that reflectance was dependent on both the observation and illumination angles, as well as the wavelength of reflected radiation [5]. However, it was also possible to verify how magnitude and contour of this function expressed land cover composition, shape, and geometric structure where the incident energy was reflected [5].

Next, the anisotropy that the different terrestrial surfaces presented were considered as the noise to correct, or as an opportunity to explore land surface 3D characteristics. For the study of the land covers anisotropy, the launch of satellites with onboard sensors was promoted. These sensors allowed the study of different directional behaviors when a target reflected solar radiation. Out of these sensors, it is important to name the POlarization and Directionality of the Earth’s Reflectances instrument (POL-DER), the Multi-angle Imaging SpectroRadiometer (MISR), and the Compact High-Resolution Imaging Spectrometer (CHRIS Proba). The latter two have the ability to capture images of almost the same location (pixel) from several different angles during the same flight path [15]. In this sense, it was easy to find several studies where MISR products were used. Recently, different authors (e.g., References [15,16,17,18,19,20,21]) have proven the MISR utility for characterizing vegetal canopy, and for mapping vegetation and land cover at medium- and large-spatial scales.

The effects of the anisotropic behavior of the different targets can be observed in all remote sensing data, but the magnitude of these effects changes with object reflection properties. Different radiative transfer models (RTM) have been developed to try and model BRDF effects of surfaces, depending on their features. Vegetation RTM describes the interaction of light and vegetal canopy by combining different view and illumination conditions with their structural and biochemical properties. Today, there are several RTM models, including the Rahman–Pinty–Verstraete model (RPV) [22]. MISR level 2 land product (MI2ASLS) provides the parameters for the modified RPV model (MRPV). It is a semi-empirical, multiplicative, three-parameter model that returns Bidirectional Reflectance Factor (BRF) values for an arbitrary geophysical medium [23,24]. Furthermore, the MRPV model parameters have been directly related to the density of trees, proportions of shadow and bright surfaces, anisotropic scattering, and solar zenith angles from the scenes taken [19,22,25,26,27,28,29,30].

Against this background, our objective was to study the values taken by three MRPV MISR L2 land surface parameters in the most representative seven agrarian and forest Co-ordination of Information on the Environment (CORINE) land cover categories in Mainland Spain (non-irrigated arable land, permanently irrigated band, agroforestry areas, broad-leaved forest, coniferous forest natural grasslands, and transitional woodland shrub) in the red and near infrared (NIR) spectral bands, and in various MISR schemes, captured in 2006 from Mainland Spain. Thereby, we sought to assess the relationship between the parameter values and land covers. This aim is considered crucial for future research, based on land cover studies in Mainland Spain through the use of the not-widely-studied product, MIL2ASLS, and its MRPV parameters.

2. Materials and Methods

2.1. Study Site

The study area for this research was the intersection of MISR path 200, 201, and Mainland Spain. This area starts in Block 56 and ends in Block 61 of MISR (see Figure 1). Mainland Spain has a wide range of different natural surfaces and is an example of the Mediterranean zone.

2.2. Co-Ordination of Information on the Environment Land Cover 2006

The CORINE (Co-ordination of Information on the Environment) Land Cover (CLC) consists of a set of maps of the European environmental landscape, based on photo-interpretation of satellite images (SPOT, Landsat TM, and The Multispectral Scanner System or better known MSS). Ancillary data (aerial photographs, topographic or vegetation maps, statistics, and local knowledge) were used to refine the interpretation and the assignment of the territory into the categories of CORINE Land Cover nomenclature. The smallest surfaces mapped were 25 hectares. Linear features less than 100 m in width were not considered [31].

The CORINE program provides comparable digital maps of land cover (and partly land use) for each country, and much of Europe. These maps shows 44 different land cover categories organized hierarchically in three levels [22].

We used a map from Spain that consisted of a 1:100,000 vectorial layer, from 2006. On that layer, seven CORINE land covers were studied (see

Table 1. The Co-ordination of Information on the Environment (CORINE) land cover categories used.

CLC_2006 Code	Category Name	Abbreviations
211	Non-irrigated arable land	Non Irri
212	Permanently irrigated land	Per Irri
244	Agroforestry areas	Agro For
311	Broad-leaved forest	Br For
312	Coniferous forest	Con For
321	Natural Grasslands	Grass
324	Transitional Woodland-shrub	Wood Sh

). We selected these natural land cover categories, as we considered them, to be the most representative of our study area. Furthermore, they also represented different 3D structures and phenological characteristics, which are very important when considering if the three MRPV parameters from the MISR L2 land surface were somehow related to these land cover properties, as some authors have noted [10,13,16,17,18,19,20,21].

2.3. MISR and MISR Data Used

The Multi-angle Imaging Spectroradiometer (MISR) was launched on 18 December 1999, and has a sensor with a 705 km sun-synchronous Earth orbit, which is onboard the Earth Observing System (EOS) Terra spacecraft. The MISR instrument has nine push cameras that view the Earth’s surface from different angles (nominal and design values); 0° (An), 26.1° (Af, Aa), 45.6° (Bf, Ba), 60.0° (Cf, Ca), and 70.5° (Df, Da), relative to nadir. In both, (f) corresponds to forward and (a) is relative to afterward, along the direction of the flight. This sensor gathers Earth information in four spectral bands, positioned at 446, 558, 672, and 866 nm, with a temporal resolution from 2 to 9 days, and a spatial resolution of 1100 or 275 m [32].

The MISR team provides several elaborate data products, including MIL2ASLS. This product provides bi-hemispherical reflectance (BHR), hemispherical directional reflectance factor (HDRF), directional-hemispherical reflectance (DHR), bi-hemispherical reflectance (BHR), and the bidirectional reflectance factor (BRF). MIL2ASLS product also includes parameters for reflectance factors based on a semi-empirical modified MRPV model (ρ0, Θ, and k), Leaf Area Index (LAI), fraction of photosynthetically-active radiation (FPAR), and the normalized difference vegetation index (NDVI) [32]. MIL2ASLS shows a 1.1-km spatial resolution, and is provided without data in some locations. These no-data zones are a consequence of glitter contamination, topographical accidents, or the presence of clouds [33].

As mentioned in the Introduction, we worked with the three parameters (ρ0, Θ, and k) of the RPV model, which were modified by Martonchik et al. [29] (MRPV). In addition to these parameters, we used some ancillary data, also provided by the MIL2ASLS product. These auxiliary data consisted of NDVI, FPAR, LAI, and zenith and azimuth illumination angle measurements from the same MISR orbits as the MRPV parameters. The MISR dataset was from 2006, as the CORINE land cover map used here was for that year. We used seven MISR images, as we wanted to maximize the number of pixels without errors and clouds; moreover, these pixels must be the same along the dates. We also needed to ensure that at least one MISR image from every season was included in the dataset (see the MIL2ASLS MISR orbits used in

Table 2. MISR schemes used. All of these are from the F07-0022 version and captured in 2006.

Date	DOY	Orbit and Path
12-Feb	43	O32738 p201
26-Apr	117	O33801 p200
4-Jun	156	O34369 p201
22-Jul	204	O35068 p201
7-Aug	220	O35301 p201
11-Nov	316	O36699 p201
13-Dec	348	O37165 p201

and Figure 1).

2.4. Rahman–Pinty–Verstraete Modified (MRPV) Model

Rahman et al. [22] introduced a nonlinear, empirical, three-parameter model for vegetation canopy, which is abbreviated as the RPV (Rahman–Pinty–Verstraete) model. The RPV model splits a BRF field into an amplitude component (ρ0). The BRF field is represented by the product of three separate functions, accounting for both the viewing and illumination directions. The RPV model formulates the BRF of a surface as:

ρ_sfc (z₀, Ω₀, → Ω; ρ0, ρ_c, Θ, k) = ρ0 M_I (θ₀, θ; k) F_HG (g; Θ) H (ρ_c; G),

(1)

where

$${\mathrm{M}}_{\mathrm{I}}\left({\mathsf{\theta }}_0,\mathsf{\theta };\mathrm{k}\right)=\frac{{\cos }^{\mathrm{k}-1}{\mathsf{\theta }}_0{\cos }^{\mathrm{k}-1}\mathsf{\theta }}{{({\cos \mathsf{\theta }}_0+\cos \mathsf{\theta })}^{1-\mathrm{k}}},$$

(2)

$${\mathrm{F}}_{\mathrm{HG}}\left(\mathrm{g}:\mathsf{\Theta }\right)=\frac{1-{\mathsf{\Theta }}^2}{{(1+2\mathsf{\Theta }\cos \mathrm{g}+{\mathsf{\Theta }}^2)}^{3/2}},$$

(3)

$$\mathrm{H}\left({\mathsf{\rho }}_{\mathrm{c}};\mathrm{G}\right)=1+\frac{1-{\mathsf{\rho }}_{\mathrm{c}}}{1+\mathrm{G}},$$

(4)

cos g = cos θ cos θ₀ + sin θ sin θ₀ cos ɸ,

(5)

G = [tan² θ₀ + tan² θ₀ – 2 tan θ₀ tan θ cos ɸ]^1/2

(6)

where ρ0 yields the overall reflectance level; and θ and θ₀ are the observation and illumination zenith angles, respectively. The relative azimuth angle (ɸ) is zero when the source of illumination is behind the sensor. M_I is the modified Minnaert function, which permits the mathematical representation of all shapes of the angular field through the k parameter. The k parameter represents bowl- (indicated by k < 1) or bell- (indicated by k > 1) shaped anisotropy. The function F_HG is based on the Henyey–Greenstein function, and the parameter Θ establishes the degree of forward (positive Θ) or backward (negative Θ) scattering. The H function allows the modeling of a hotspot with the ρ0 parameter. In some cases, the parameter accounting for the hot spot effect, ρ_c, can be set at a fixed value or equaled to the scalar amplitude parameter ρ0, being three parameters instead of four. The MIL2ASLS product only includes three parameters (ρ0, k and Θ).

Martonchik et al. [20] modified the RPV model by replacing its Henyey–Greenstein function with an exponential function after applying a logarithmic transformation, which resulted in the MRPV quasi-linear model being obtained. The fact that it is a quasi-linear model makes it easier to use in an operational context [23,34].

2.5. Data Processing and Data Analysis

The MIL2ASLS data from the seven MISR schemes in the red and NIR band were downloaded from the NASA Langley Distributed Active Archive Center (DAAC) and re-projected from the Space Oblique Mercator (SOM) to the Universal Transverse Mercator (UTM) system. We worked with red and NIR bands as these have the wavelengths most related to natural land covers in the available MISR bands, and, for that reason, they are the most studied in the literature [35].

We observed cloud contamination in some pixels from the MIL2ASLS product; therefore, an additional cloud mask was applied on each of the seven MISR schemes. This cloud mask was created through RGB compositions, which, in turn, were constructed by the BRF values from the red channel of Da, An, and Df cameras of each of the study area pixels. These BRF values were obtained from the MI1B2T product, which can be downloaded from the same site as MIL2ASLS. In this manner, it was easy to detect cloud contamination pixels by eye, and were therefore easy to remove [36].

We also carried out a buffer of 500 m on the CORINE land cover map with the seven CLC categories shown in Table 1. The buffer operation was carried out to avoid mixing the two or more types (at least 78% homogeneity in the land cover category of every pixel) that could have different behaviors in their phenology or 3D structural properties, and may translate into noise in the parameter values of every MISR image.

Next, the center of the pixels of the three parameters (ρ0, Θ, and k) and the ancillary data (NDVI, FPAR, LAI, and zenith angles) obtained from the MIL2ASLS images were spatially joined with the CORINE land cover map that resulted after applying the buffer and selecting the categories. As a consequence of this operation, we obtained a composition map made up of points with the parameter, as well as the ancillary variable values, CLC code, MISR band, and orbit (date of capture) information from the study area. Then, we assessed the reliability pixel value included in MIL2ASLS. Finally, pixels were selected at random to obtain a sample made up of 30 different pixels from each MISR orbit and CLC category. The latter was created to have a balanced analysis to perform other inferential analyses with robust results using this kind of sample. It is important to note that some of the FPAR and LAI showed no data, as it is difficult to get a random sample with 30 data pixels for these two variables.

With the final sample of pixels, we built a box plot by grouping the variables by CLC 2006 codes and MISR orbits. Then, the Spearman correlation coefficients between each combination parameter–band were calculated, as well as the NDVI, FPAR, and LAI, in order to evaluate the relationship between the MRPV parameters of this MISR product and the variables that express the phenological character of land covers.

Finally, we constructed linear multi-regression models and the variance analysis (ANOVA) and Kruskal–Wallis test. The linear regression models were made using as dependent variable (Y) each parameter–band combination, and using the CLC category and orbit as explanatory variables (X_n). A two-factor test of variance (ANOVA and Kruskal-Wallis) was conducted, using each combination parameter band as a numerical variable, and in which the categorical variables were the CLC category and the MISR orbit. These tests were made first to assess the relationship between land cover and orbit and the values achieved by the parameters. On the other hand, variance tests were performed to evaluate how the parameter–band combination differed between land covers and/or between MISR orbits captured on different days and belonging to different seasons. We decided to apply ANOVA and the Kruskal Wallis test (the non-parametrical test of variance) as some analyses of variance did not meet the homoscedasticity assumption of the parametrical test. Finally we used a post hoc pairwise test (Tukey) to know what pairs of CLC categories or orbits were statistically different.

3. Results and Discussion

3.1. Rahman–Pinty–Verstraete Modified (MRPV) Model Parameters Behavior in Different Land Uses and Captures in Mainland Spain

Figure 2 shows the ρ0, Θ and k values in the red and NIR spectral bands of MISR orbits from 2006 in different classes of land use CLC for Peninsular Spain.

The ρ0 red values from the MISR scenes ranged between 0.01 and 0.14, while those related to ρ0 NIR were between 0.05 and 0.35. Meanwhile, the k red average values were from 0.35 to 1, and, for the k NIR, oscillated between 0.5 and 0.9. Finally, the average values of Θ in the spectral band of red were between –0.1 and –1.4, and those found in the NIR band were from 0 and –0.9 (see Figure 2).

After reviewing the values presented by other authors (who have also worked with the parameters of the MRPV model), we verified that the values reported in our study coincided with those examined in the literature [25,37,38].

The graphs in Figure 2 show that the k values were never greater than one; furthermore, the Θ parameter values were always negative. The latter indicated a predominance of backscattering [26,37,38]. Pinty et al. [17] were the first to report how the anisotropic character of terrestrial targets is influenced by the structure and density of vegetation on a sub-pixel scale. The authors of the papers cited before demonstrated how, for a light soil with dark and vertical vegetation, the presence of trees at medium-high densities created a bell shape (values on the BRF curve along the viewing angles). This pattern was observed in the red spectrum, as in NIR, this spectral contrast between the bright background and the trees occurred on the contrary [17]. The k values less than one indicated how the BRF curve, along the zenith observation angles, presented a bowl shape, the most common in natural land covers overall, at a 1.1-km spatial resolution [26]. This parameter was highly dependent on spatial resolution, and when we have a medium resolution (as the MIL2ASLS product is), it appeared harder to find the bell form BRF curve [26], which occurred in our results.

On the one hand, Figure 2 shows how parameter ρ0, in the Br For, Con For, and Wood Sh CLC categories, presented more similar median values in all the MISR orbits than the other land cover classes. These results were consistent with those exposed by Alcaraz-Segura [39], who used NDVI values in the Iberian Ecosystems, instead of ρ0 ones. Figure 2 also shows how the average values from the Con For CLC were still more similar in the study of MISR orbits than values belonging to Br For, or even more so than those belonging to Wood Sh. This is because Con For on the CLC 2006 map is more homogenous than Br For CLC, which does not show separate deciduous and evergreen plants. The last point explains the fact that, in Figure 2a,b, we can observe ρ0 red and NIR values from Br For CLC with more dispersion within the same orbit (interquartile range found in the boxplots) than in Con For CLC. Additionally, this also explains the fact that we could see more steady behavior in the values of ρ0 red and NIR pixels belonging to the Br For and Wood Sh CLC categories.

In Figure 2a,b, it can be seen how the crops (Non Irri, Per Irri), pasture (Grass), and agroforestry (Agro For) CLC ρ0 red, as well as ρ0 at NIR values, presented more variability through the MISR orbits, taken on seven different dates, than the Br For, Con For and Wood Sh classes did. In general, the trend of the ρ0 red and NIR values seemed to be related to the land covers phenological pattern. To further investigate this idea, we focused on Figure 3, which shows four graphs with the NDVI, FPAR, NDVI calculated using the ρ0 red and NIR values (in Figure 2a,b) and the LAI variables. From Figure 3a,c, the median and the interquartile ranges were extracted, and, in general, the behavior of the ρ0 and NDVI values along the seven MISR orbits were quite similar. The FPAR (Figure 3b) and the LAI value (Figure 3d) dispersion within one CLC and orbit were higher than the dispersion taken using NDVI (Figure 3a) and NDVI from ρ0 (Figure 3c). However, the median values that demonstrated the behaviors of these four variables (Figure 3) were similar. Thereby, in general, it can be noted that the mean values may support the idea that the ρ0 parameter has a relationship with spectral properties and land cover phenological features, as other authors have already described [19,25].

The k parameter in the NIR and the red band seemed to present lower median values in the orbits belonging to winter and autumn than in those of summer and spring (this fact was slightly more marked in the red wavelength than in the NIR one). That is to say, the covers present, in general terms and from the k parameter, a more isotropic behavior (k values close to the unity) from the dates in summer and spring than in those of autumn and winter. Furthermore, the zenith angles from the MISR scenes used here were higher in the winter images than in those belonging to summer and spring, as seen in Figure 4.

The parameter k shows the difference between observing BRF from the nadir camera in comparison to the BRF from other positions, further away from this nadir position [18,40]. The fact that finding average values of k expressed greater anisotropy in the winter and autumn schemes than in summer and spring indicated the influence of the zenith illumination angle on these values. That is, in winter, the sun position is farther from the zenith than in spring and summer. This fact also implied that, in winter and autumn, the sensor captured much more shade from the nadir in comparison with other positions more distant from it, than in summer and winter. Thereby, it could be found greater anisotropy expressed by k in winter and autumn than in summer and spring. The influence of the zenith angle of the data collection on the values reached by k has already been reported in Pinty et al. [17] or in Reference [37]; furthermore, their results coincided with those presented in our current work.

From Figure 3c,d, it can be seen that the k parameter values did not follow the same behavior as ρ0, NDVI, FPAR, and LAI, i.e., the values reached by k did not seem to be dependent on the phenology.

Figure 2e,f shows the Θ parameter pixel values in the different orbits and CLC categories used in this study. The Θ values given from the land use classes in the seven MISR orbits were more anisotropic in the red band than in the NIR, i.e., the values reached were lower in red than in NIR. This fact was due to the effect of multiple dispersion on the NIR channel, which caused the anisotropic behavior to be less than in the red band. This physical effect was observed and explained by Sandmeier et al. [32]. Furthermore, these results were consistent with those obtained by Armston et al. [19], where the authors studied different natural covers and dates in Queensland, Australia, and found this effect on bare soil, as well as in soil in addition to vegetation.

A comparison of these graphs (Figure 2e,f) to the NDVI, FPAR, and LAI ones (Figure 3) showed that the trends that followed the median and the interquartile ranges of each variable were not very similar. For instance, the Θ red and NIR values of the Grass CLC category appeared to be steadier through the orbits than the others. The Θ red values of the Br For CLC class were slightly more stable through the orbits than the values belonging to the Con For CLC category; however, the NDVI plot (Figure 3a) occurred on the contrary.

With regards to the anisotropy described by the Θ parameter in our study, it could be appreciated how the Br For and Con For CLC categories were those that presented Θ values more different than zero, i.e., they were the most anisotropic from the perspective of the Θ parameter. Camacho de Coca et al. [41] explained that the natural land covers made up of vegetation with a predominance of vertical structure 3D were the most anisotropic, and even more so when the background was bright. For that reason, it also explained why we found that the values of the Θ parameter were even more different than zero (more anisotropic) in the orbits used that belonged to summer or near summer. It is in those orbits where the background as the driest of all seven dates used for this study (see Figure 2e,f).

3.2. MRPV Model Parameters Related to NDVI, LAI, FPAR

Table 3. Spearman correlation coefficient from the NDVI and each of the MRPV parameters. In the second column, the Rho’s from all seven MISR orbit pixels is shown. Missing cells indicate non-significant Spearman correlation. The Non Irri-Wood Sh are the abbreviations to the CLC categories.

NDVI	All	Non Irri	Per Irri	Agro For	Br For	Con For	Grass	Wood Sh
ρ0 red	–0.8346	–0.8825	–0.9641	–0.9403	–0.7485	–0.815	–0.8269	–0.805
ρ0 NIR		0.5805	0.8276	0.72			0.646	0.1805
k red	–0.3129		–0.4634	–0.6789			–0.3762
k NIR	–0.277			–0.3853
Θ red	–0.1991	–0.5469	–0.5201
Θ NIR	0.1505	0.6451	0.6085	0.6451			0.514

,

Table 4. Spearman correlation coefficient from the FPAR and each of the MRPV parameters. In the second column, the Rho’s from all seven MISR orbit pixels is shown. Missing cells indicate non-significant Spearman correlation. The Non Irri-Wood Sh are the abbreviations to the CLC categories.

FPAR	All	Non Irri	Per Irri	Agro For	Br For	Con For	Grass	Wood Sh
ρ0 red	–0.4899	–0.7665	–0.91	–0.8412			–0.6926	–0.2941
ρ0 NIR	0.2628	0.6179	0.8669	0.759			0.6697	0.3957
k red	–0.3198		–0.3925	–0.6586
k NIR				–0.3414
Θ red		–0.5661	–0.5037
Θ NIR	0.2667	0.5139	0.6129	0.6746			0.5002

and

Table 5. Spearman correlation coefficient from the LAI and each of the MRPV parameters. In the second column, the Rho’s from all seven MISR orbit pixels is shown. Missing cells indicate non-significant Spearman correlation. The Non Irri-Wood Sh are the abbreviations to the CLC categories.

LAI	All	Non Irri	Per Irri	Agro For	Br For	Con For	Grass	Wood Sh
ρ0 red	–0.5864	–0.7888	–0.9146	–0.8404	0.3965	–0.4148	–0.7268	0.5205
ρ0 NIR	0.1502	0.616	0.877	0.7071	–0.1565		0.5743	0.2713
k red	–0.1895		–0.4093	–0.5231
k NIR				–0.3177
Θ red	–0.2438	–0.6513	–0.5009
Θ NIR	0.149	0.4759	0.6483	0.5578			0.3612

show the Spearman Rho’s coefficient between the MRPV parameters from the MIL2ASLS and the NDVI, FPAR, and LAI. In the tables, no Rho’s coefficient value meant that the relationship was not significant. The results suggested that the correlation was stronger between the ρ0 and NDVI than between the other three parameters and the NDVI, as well as occurring cases in the correlation between these MRPV parameters and the FPAR and LAI variables. The correlation between ρ0 red and NDVI, FPAR and LAI was slightly stronger than that between the ρ0 NIR and NDVI, FPAR, and LAI. Finally, it appeared that the relationship between Θ and the last variables mentioned were higher than the relationship between those and the k parameter. These results proved (as seen in Figure 2 and Figure 3) the relationship between the ρ0 parameter and the NDVI, FPAR, and LAI, and matched those described by other authors using different sites than our study area [19,25].

3.3. Inferential Statistics

Table 6. Results from the linear multi-regression models and the ANOVA test (Df is the Degrees of freedom, F is the Fisher statistical value, lm refers to the linear model, * implying significance at p < 0.01 level).

Response	Factor	Df	F	Sig	R2 Adj. (lm)	F Statistic (lm)
NDVI (ρ0)	CLC	6	138.482	*	0.67	57.13
	Orbit	6	146.738	*
	CLC & Orbit	36	28.837	*
ρ0 red	CLC	6	351.08	*	0.73	75.3
	Orbit	6	176.01	*
	CLC & Orbit	36	13.04	*
ρ0 NIR	CLC	6	270.506	*	0.74	70.27
	Orbit	6	119.823	*
	CLC & Orbit	36	28.855	*
k red	CLC	6	19.517	*	0.49	28.12
	Orbit	6	184.4	*
	CLC & Orbit	36	3.859	*
k NIR	CLC	6	65.27	*	0.57	37.16
	Orbit	6	191.4225	*
	CLC & Orbit	36	7.09	*
Θ red	CLC	6	62.68	*	0.37	17.53
	Orbit	6	43.16	*
	CLC & Orbit	36	5.7	*
Θ NIR	CLC	6	62.69	*	0.41	20.08
	Orbit	6	43.17	*
	CLC & Orbit	36	5.725	*

shows the results of the linear model regressions and the parametric analyses of variance performed using the MRPV parameters (Y) and the categorical variables, named MISR orbit and CLC category (X_n).

Table 7. Results from the Kruskal–Wallis test (Df is the degrees of freedom, * implying significance at p < 0.01 level).

Response	Factor	Df	K-W Chi-Square	Sig
NDVI (ρ0)	CLC	6	284.5169	*
	Orbit	6	279.8697	*
	CLC & Orbit	48	897.9475	*
ρ0 red	CLC	6	696.3339	*
	Orbit	6	225.8915	*
	CLC & Orbit	48	1007.171	*
ρ0 NIR	CLC	6	543.24	*
	Orbit	6	175.5939	*
	CLC & Orbit	48	937.2772	*
k red	CLC	6	54.94	*
	Orbit	6	596.8069	*
	CLC & Orbit	48	702.8155	*
k NIR	CLC	6	158.9761	*
	Orbit	6	504.9536	*
	CLC & Orbit	48	786.1684	*
Θ red	CLC	6	242.0182	*
	Orbit	6	150.5548	*
	CLC & Orbit	48	515.93	*
Θ NIR	CLC	6	328.6615	*
	Orbit	6	54.9	*
	CLC & Orbit	48	596.662	*

presents the results from the non-parametric test of variance. At a glance, the results showed how the two categorical variables explained the parameters in a significant way (p < 0.01) (see Table 6). The models yielded adjusted R-square values, which rated between 0.37 and 0.74. Taking both the adjusted R-square and F statistics values into account, it appeared that the relationship between the parameters and the MISR dates and CLC category was stronger in the case of ρ0 than in the k, and was weakest in the models explaining the Θ parameter.

All analysis of variance performed showed significantly different mean and median using Orbit, CLC, and Orbit and CLC (see Table 6 and Table 7). The resulting F from the parametrically analyses of variance, as well as the Chi-square results from the non-parametrically analyses of variance (see Table 6 and Table 7) were higher for the CLC than for the MISR orbit factor in tests related to ρ0 and Θ, and was the opposite in the case of the k parameter test. The NDVI F and Chi-square values were similar between the CLC and Orbit cases.

For the post hoc test performed using the CLC categories shown in Figure 5, the results illustrated that the majority of these groups were significantly different, using NDVI (ρ0) mean values. This was a good result, since the aim was to use this parameter, for instance, to separate between land covers, but was not a striking finding, as NDVI has already been utilized in some studies [35] to classify or to detect changes in the land uses. For this reason, we focused on the pair of images that did not significantly differ from the mean NDVI (ρ0) values. These pairs were built using the Wood Sh and the Grass (transitional woodland-shrub and natural grassland), and the pair formed by the Non Irri and Per Irri (non-irrigated arable land and permanently irrigated arable land) CLC categories. Table 7 shows, on the one hand, that the pair conformed by Wood Sh and Grass were significantly different from the mean values of the k red and NIR and by the Θ red and Θ NIR. On the other hand, the pair built by Non Irri and Per Irri differed statistically from the mean values of the k red. Thus, these results may indicate that the k and Θ parameters from the MIL2ASLS product, in addition to ρ0, which matched well with classical NDVI (like the results shown in Section 3.2), could help to differentiate between those land covers, which are difficult to differentiate with NDVI even, including data for different seasons of the year.

For the post hoc test performed using the MISR orbit (Figure 6a) it can be seen that the NDVI (ρ0) values were grouped into three categories. The first category was made up of the MISR orbit from April, which presented the highest mean value, the other was made up of the orbits from December, November, and June. Finally, the last category was confirmed by the orbits of February, July, and August which showed the lowest value of NDVI (ρ0). Effectively, the vegetation in those three categories, in general, presented very marked different phenological behaviors within Mainland Spain. However, for the k and the Θ parameter, the classes resulting from the post hoc test were not conformed in the same way as in the NDVI (ρ0).

In general, the parameter-band CLC categories in the different studied orbits had a high degree of separability, at least given that the mean values belonged to its CLC or orbit. This is why other authors have included multi-angular information in their work on the classification of images, among other utilities. For example, Brown et al. [42] observed how the inclusion of multi-angular information, in addition to only multi-spectral from the POLDER sensor, produced an improvement in the overall accuracy of close to five percent, achieving a maximum improvement between the classes of wooded savanna and grassland. On the other hand, Mahtab et al. [43] also obtained increases in the degree of separability of classes, introducing multi angular information in addition to multispectral. Su et al. [44] increased the overall accuracy by 10–15% using multi-angle information, which was also used to classify between different types of coverage. MIL2ASLS with its MRPV parameters provide us with a source of multi spectral and multi angular information that can be successfully used for great- and medium-scale land cover studies. This is because this product provides a multispectral view of classical remote sensing, thanks to the ρ0 parameter, and, on the other hand, a multi-angular view with the k and Θ parameters.

4. Conclusions

According to the results presented in the sections of this study, it can be said that the parameter ρ0 (at the red and NIR channel) seemed to be the most related to the NDVI, FPAR, and LAI, i.e., it appeared to be the parameter that could better express the photosynthetic activity or the phenological behavior in the land covers studied in Mainland Spain. We noted this first at a glance when comparing the boxplots of the ρ0 by CLC category in the different seven MISR orbits also captured on seven dates, with NDVI, FPAR, and LAI; and, second, by studying the Spearman correlation coefficient between the ρ0 and the other three variables. The k and the Θ parameter values (also at the red and NIR channel) from the pixels explored showed different patterns than those followed by NDVI, FPAR, and LAI. We compared the k and Θ with the illumination angles and noted how these angles could control, or at least somehow influence the values reached by the k and Θ parameters, and, therefore, the grade of anisotropy expressed by them.

From the inferential test, linear models were constructed to explain each combination of parameter and band, and their analyses of variance were executed using the CLC category and the MISR orbit as a categorical variable to allow us to assess how these two explicative variables could explain each parameter variance in a significant way. The CLC and the Orbit explained more variance of ρ0 than k or Θ, and post hoc tests, performed by pair groups, of these two variables, showed how ρ0, in general, alone or used as the classic NDVI form, showed many differences from its mean values. However, the key was to use ρ0 in addition to k or Θ, to assure the separability of the CLC types that could only not be differenced with ρ0.

This study suggested that, although the number of studies up to the present, used for the MIL2ASLS MRPV model parameter data product, is not very large, it has great potential to be used for improving land cover research. For this reason, amplifying the number of future studies and applications, which involve the parameters of MIL2ASLS, is worthwhile.