Plant Species Richness Prediction from DESIS Hyperspectral Data: A Comparison Study on Feature Extraction Procedures
and Regression Models

This study focused on two different habitat types in southeast New South Wales, Australia, namely the Southern Tablelands (34°12’26”–34°39’07”S, 150°05’57”–150°40’51”E) and Snowy Mountains (35°43’58”–36°16’30”S, 148°23’16”–148°39’02”E), as shown in Fig. 1. Both regions are located within the climate zone of Cfb (oceanic climates), according to the Köppen–Geiger climate classification system.

The Southern Tablelands region is located to the southwest of Sydney (Fig. 1). It is characterised by high-altitude plains with a rich biodiversity. There are more than 1200 plant species within Southern Tablelands, of which 30 are listed as threatened [Fallding, 2002]. A large part of the landscape has been transformed into suburbs for residential developments and pastures for grazing purposes. Considering the high degree of human interference and habitat alteration, conservation efforts have been undertaken in order to preserve endangered plant species by improving habitat connectivity and condition.

The Snowy Mountains region is located to the southwest of Canberra (Fig. 1). It encompasses the highest mountain ranges of the Australian Alps, serving as an important habitat for alpine-exclusive species. There are 212 species of vascular plants, of which 21 are endemic [Pickering et al., 2008]. Due to its unique status in Australia’s ecosystem, the plant biodiversity in Snowy Mountains has drawn consistent interest from the research community (e.g. Körner [1995], Pickering et al. [2008], and Pickering and Green [2009]).

For on-ground measures of vascular plant species richness, we obtained plant community survey data from the NSW BioNet Vegetation Information System database [Government, 2019]. Field surveys were conducted to collect species richness samples in 2016 and 2017 for the two regions, with a sampling plot area of 400 m² (20 m $\times$ 20 m) at each surveying location. A significant bushfire event occurred during the 2019–2020 summer, with some of the sampling points situated within the affected areas. These bushfire-affected samples were excluded from the data set, based on the National Indicative Aggregated Fire Extent Datasets (NIAFED) provided by the Australian Government Department of Agriculture, Water and the Environment. After the exclusion, a total of 44 and 29 samples were used in this study for analysis for the Southern Tablelands and Snowy Mountains regions, respectively. The locations of these samples are shown in Fig. 1, and their associated information is summarised in Table 1. For each sampling plot, the number of native vascular plant species was calculated and used as the response variable in our analyses.

The histograms of species richness distribution for sampling plots in the two regions are shown in Fig. 2. Generally, the sampling plots in Southern Tablelands show a higher richness of species than Snowy Mountains, with a mean richness value of 44.4 for the former and 23.5 for the latter. The difference in species richness can be mainly attributed to the fact that the Southern Tablelands is located at lower altitudes with a relatively warmer climate than the mountainous region of Snowy Mountains.

The DESIS spectrometer [Krutz et al., 2019] is embedded in the MUSES platform onboard the International Space Station at an altitude of approximately 400 km. It operates in a push-broom imaging mode featuring state-of-the-art radiometric and spectral specifications. It delivers hyperspectral images with 235 spectral bands over the visible and near-infrared regions of 400 $\sim$ 1000 nm, with a spectral resolution of 2.55 nm and a spatial resolution of 30 m. The radiometric resolution for each band is 12 bit with 1 bit gain. The signal-to-noise ratio is 195 at the wavelength of 550 nm.

In our study, the DESIS Level-2A product was used. It consisted of surface reflectance images with atmospheric correction having been applied. The correction was conducted with DLR’s PACO (Python Atmospheric COrrection) software [De los Reyes et al., 2018] where the MODTRAN^® radiative transfer model [Berk, 2016] served as the module for simulating atmospheric effects. As inputs for atmospheric simulation, aerosol optical thickness and water vapour content were retrieved per pixel using reflectance in the red and NIR bands, and bands around the water absorption features of 820 nm, respectively. DESIS spectra intersecting with locations of the species richness samples were queried within CSIRO’s Earth Analytics and Science Innovation (EASI) platform. We selected spectra captured in January 2020 for our analysis. In order to moderate the random noise present in the original spectra, the spectral resolution were down-sampled from 2.55 nm into 10.2 nm bins with the assumption of a Gaussian-shaped spectral response function. The atmospherically affected bands of 759, 769, 933.4, 943.4, and 953.2 nm, and the low quality bands of 402.8, 410.3, and 999.5 nm at the left and right ends of the spectrum were removed. A total of 52 bands were retained after the removal. Pixels flagged as cloud by the DLR Level-1 processing were masked out.

The Bidirectional Reflectance Distribution Function (BRDF) effect in DESIS data is unneglectable, given the $4.1^{\circ}$ field of view in conjunction with the $\pm 15^{\circ}$ along track pointing capability of the DESIS sensor, and the $\pm 25^{\circ}$ along track and $-45^{\circ}$ $\sim$ $+5^{\circ}$ cross track tilting capability of the MUSES platform. Following the approach adopted in Green and Craig [1985] and Ong and Cudahy [2014] for correcting BRDF effect in hyperspectral data, in our study each DESIS spectra was mean-normalised with each spectral band being divided by the mean value over all bands.

For comparison purpose, cloud-free Sentinel-2 multispectral data observed closest to the sensing time of DESIS spectra were also downloaded. These Sentinel spectra were downloaded as Level-2A surface reflectance with atmospheric effects being corrected. Each spectrum consisted of 12 bands covering the visible and near-infrared spectral regions. The Sentinel-2 data were re-sampled into a spatial resolution of 30 m to be consistent with that of the DESIS data. A comparison between the DESIS and Sentinel-2 spectra at one of the ground sampling plots is shown in Fig. 3. Both spectra cover the visible, near-infrared, and short-wave-infrared (SWIR) regions, including the red-edge region that is critical for vegetation mapping. The Sentinel-2 and DESIS spectra show a consistent shape in these spectra regions, with DESIS having a much denser band coverage. The SWIR region is covered by Sentinel-2 only with its two bands. Though information in SWIR is not contained in the DESIS data, it will be provided by the upcoming DLR mission of EnMAP (for which DESIS has been served as a preparation mission), as shown in Fig. 3.

The DESIS hyperspectral data provide rich information in the abundant and spectrally continuous bands, but these bands are highly collinear. We performed dimensionality reduction to address this problem—testing three different approaches, namely the Principal Component Analysis (PCA), Canonical Correlation Analysis (CCA), and Partial Least Squares analysis (PLS). These methods aim to find a linear transformation to project the DESIS spectra from the original space of $n$ spectral bands to a new space of a reduced dimensionality defined by $k$ uncorrelated components, with $k$ being smaller than $n$.

Mathematically, the transformation for dimensionality reduction is written as:

where the input data for the transformation is $\mathbf{X}$, which is a $n\times m$ matrix consisting of the original spectra, with $n$ being the number of observed spectra and $m$ being the number of spectral bands; the output of the transformation is $\mathbf{T}=[{\bm{t}}_{1},{\bm{t}}_{2},\cdots,{\bm{t}}_{k}]$, which is a $n\times k$ matrix consisting of $k$ components of $\mathbf{X}$; $\mathbf{W}=[{\bm{w}}_{1},{\bm{w}}_{2},\cdots,{\bm{w}}_{k}]$ is a $m\times k$ matrix transforming $\mathbf{X}$ from the original $n$-dimensional space into a new space of $k$ components. The weights in $\bm{w}_{i}$ are a measure of relative contributions of the original bands to the transformed component ${\bm{t}}_{i}=\mathbf{X}{\bm{w}}_{i}$. The number of components, $k$, needs to be preset. In this study, $k$ is selected as the one that achieves the highest cross-validation accuracy.

The PCA is an unsupervised algorithm that finds orthogonal components as the ones that explain the maximum variance in the spectral data, disregarding the target values of species richness. In contrast, CCA and PLS are supervised with both spectral data and species richness values being taken into account in the computation of components. The difference between CCA and PLS is that CCA seeks to maximise the correlation between computed components and species richness values, while PLS aims to maximise the covariance between the two.

After reducing the dimensionality of spectral data from the original $n$ bands to $k$ components, regression is conducted to predict species richness from the components, such that the mismatch between model-predicted and ground-truth species richness is minimised. The Kernel Ridge Regression (KRR), Gaussian Process Regression (GPR), and Random Forest Regression (RFR) algorithms are employed and compared, covering respectively the deterministic, Bayesian, and ensemble approaches to statistical regression. Here we focus on formulating our task of estimating species richness within the frameworks of these regression approaches.

The DESIS spectra $\left\{\bm{x}_{i}\right\}_{i=1}^{n}$ can be represented by their components $\left\{\bm{t}_{i}\right\}_{i=1}^{n}$ using one of the dimensionality reduction methods described in Subsection 2.3. The KRR transforms $\left\{\bm{t}_{i}\right\}_{i=1}^{n}$ into a feature space of high dimensionality (potentially infinite dimensionality) via a function $\varphi(\bm{t}_{i})$. A linear model in the high-dimensional feature space is non-linear when it projects back to the original space, thus enabling capturing non-linear relationships in the data. The number of parameters for a linear model increases with space dimensionality, posing a risk of over-fitting. In order to constrain the model complexity, a regularisation term is adopted in KRR to penalise the norm of the coefficient vector $\bm{b}$. The optimisation problem is:

where $\lambda$ is the regularisation parameter controls the relative importance of the regularisation term $\left\|\bm{b}\right\|_{2}^{2}$.

The input data for Eq. 2 are the transformed components $\left\{\bm{t}_{i}\right\}_{i=1}^{n}$. Solving Eq. 2 does not involve calculating $\varphi(\bm{t}_{i})$. Instead, it only requires computation of the inner product $k(\bm{t}_{i},\bm{t}_{j})=\varphi(\bm{t}_{i})\varphi(\bm{t}_{j})$, where $\bm{t}_{i}$ and $\bm{t}_{j}$ are pairs from the input data $\left\{\bm{t}_{i}\right\}_{i=1}^{n}$. In this study, the kernel function $k(\bm{t}_{i},\bm{t}_{j})$ is specified as a combination of a dot-product kernel $k_{d}(\bm{t}_{i},\bm{t}_{j})$, a radial-basis function kernel $k_{r}(\bm{t}_{i},\bm{t}_{j})$, and a white kernel $k_{w}(\bm{t}_{i},\bm{t}_{j})$:

where $\sigma$, $l$, and $\delta$ are hyperparameters that need to be selected with grid search. The dot-product and radial-basis function kernels account for the linearity and non-linearity of the data, respectively, while the white kernel explains the noise in the data. For a new spectrum $\bm{x}$ with its extracted components $\bm{t}$, the predicted species richness is $\hat{y}=\sum_{i=1}^{n}{\alpha}_{i}k(\bm{t}_{i},\bm{t})$, where ${\alpha}_{i}=(\mathbf{K}+\lambda\mathbf{I})^{-1}y_{i}$ with $\mathbf{K}_{i,j}=k(\bm{t}_{i},\bm{t}_{j})$.

In contrast to KRR that is formulated in a deterministic form, GPR is a Bayesian approach for regression. The underlying function correlating DESIS components and species richness is assumed to be distributed probabilistically as a Gaussian Process (GP):

where $m(\bm{t})$ is the mean function which is often set to 0, and $k(\bm{t},\bm{t}^{\prime})$ is the covariance function that can be specified as a kernel function. The input data for Eq. 4 are the transformed components $\left\{\bm{t}_{i}\right\}_{i=1}^{n}$. The covariance function $k(\bm{t},\bm{t}^{\prime})$ for GPR is set in the same way as that for KRR (Eq. 3), i.e., $k(\bm{t},\bm{t}^{\prime})=k(\bm{t}_{i},\bm{t}_{j})$ where $\bm{t}_{i}$ and $\bm{t}_{j}$ are pairs from the input data $\left\{\bm{t}_{i}\right\}_{i=1}^{n}$. In contrast to KRR in which the optimal kernel parameters is found by grid search, the parameters of the covariance function in GPR can be automatically determined based on gradient-ascent on the marginal likelihood function. After the posterior likelihood is determined for the GP, the predicted species richness is $\hat{y}=\sum_{i=1}^{n}{\alpha}_{i}k(\bm{t}_{i},\bm{t})$ for a new spectrum, where ${\alpha}_{i}=(\mathbf{K}+\epsilon^{2}\mathbf{I})^{-1}y_{i}$ with $\mathbf{K}_{i,j}$ being $k(\bm{t}_{i},\bm{t}_{j})$ and $\epsilon$ being a pre-set parameter explaining the noise in the data.

In addition to KRR and GPR that are respectively belong to the deterministic and Bayesian regression categories, we also tested the ensemble method of RFR. A random forest combines a number of decision tree regressors and takes the average regression result over all trees as the final estimate. Due to the ensemble structure, RFR tends to produce robust regression results with high resistance to overfitting and data noise. The training of a random forest regressor minimises the following optimization function:

where $h(\cdot)$ is the decision tree regressor with $\bm{\theta}_{j}$ being the parameters of the $j$th tree; $d$ is the number of decision trees, which was set to 100 in this study.

Hyperparameters that need to be pre-set included the number of components $k$ in the dimensionality reduction methods PCA, CCA, and PLS, and the kernel parameters $\sigma$, $l$, and $\delta$ in the non-linear regression method KRR. The number of components $k$ determines how much information to retain after dimensionality reduction. An optimal selection of $k$ is able to reduce data redundancy without excessively discarding useful information in the original DESIS spectra. In this study, we tested different $k$ values ranging from 1 to 10. The optimal $k$ value was selected based on the accuracy of species richness prediction, and the amount of variance in the spectral data that the retained components could explain. The tunable hyperparameters $\sigma$, $l$, and $\delta$ for the kernel function define the non-linearity structure of the regression model KRR. These kernel parameters were selected based on grid search. A grid of values was tested with each parameter varying from $10^{-5}$ to $10^{5}$ on a logarithmic scale. The combination of $\sigma$, $l$, and $\delta$ that produced the best performance in species richness prediction was selected.

A two-fold validation scheme was used in this study for assessing the modelling accuracy. For each of the Southern Tablelands and Snowy Mountains regions, the whole data set was randomly partitioned into two subsets (Subsets I and II), with Subset I dedicated to training and Subset II for validation (Round I), followed by Subset II for training and Subset I for validation (Round II). This procedure was repeated 100 times with the data set being partitioned differently each time. Correlation diagrams were plotted between the ground-truth species richness and the predicted values from the DESIS spectra. The coefficient of correlation ($r$) and Root-Mean-Square Error (RMSE) were calculated to evaluate the performance of the models. Results with different feature extraction procedures (PCA, CCA, and PLS) and regression models (KRR, GPR, and RFR) were computed and compared.

The DESIS data consist of spectral measurements over the visible and near-infrared bands from 400 nm to 1000 nm. Importance analysis was conducted in order to analyse which spectral bands provide more explanatory power than others in predicting plant species richness. We used the vector length of contribution values of each band to all components used in the regression weighted by the partial correlation coefficients as the importance index of the band:

where $I_{i}$ is the importance index for the $i$th band; $w_{i,j}$ is the $(i,j)$th element of the weight matrix $\mathbf{W}$ in Eq. 1, representing the contribution of the $i$th band to the $j$th component; $p_{j}$ is the partial correlation coefficient of the $j$th component with the target variable of species richness; $k$ is the number of components used in regression. In our study, the importance indices $I_{i}$ are normalised to relative values $\tilde{I}_{i}$ with sum over all bands equal to one:

where $m$ is the number of spectral bands.

Spaceborne multispectral imagery such as Sentinel-2 is more readily available than hyperspectral. Though Sentinel-2 images have less bands than the DESIS hyperspectral data, they are delivered on a more stable and systematic basis. In this study, an analysis is conducted to see if Sentinel-2 multispectral data are able to achieve comparable results. Through this analysis, we hoped to examine how well Sentinel-2 could serve as a substitute for plant biodiversity mapping in instances where hyperspectral data are unavailable.

The Sentinel-2 data set described in Section 2.2 was used for the comparison. In order to conduct a fair comparison between hyperspectral and multispectral data with minimised differences in instrumental specifications and acquisition conditions, a Sentinel-2-like synthetic data set was simulated from the DESIS data. The simulation involved resampling the DESIS spectra using the spectral response functions of Sentinel-2’s visible and near-infrared bands. Both the real and synthetic Sentinel-2 sets were used for species richness prediction, with results being compared to those achieved with DESIS data.

The blue, green, and red curves in Fig. 4 show the DESIS spectra averaged over ground sampling plots with species richness falling into the low, intermediate, and high tertiles, respectively. Results for the Southern Tablelands and Snowy Mountains regions are displayed in Figs 4a and b, respectively. For each region, it was observed that plots of higher richness showed a lower reflectance in the visible range of $400\sim 680$ nm. Considering that major absorption features of chlorophyll are located within the visible region [Zhao et al., 2014b], the lower reflectance in this spectral portion might indicate a higher concentration of chlorophyll. It was also observed that plots of higher richness showed a higher reflectance in the near-infrared plateau of $780\sim 1000$ nm and a steeper red edge between $680\sim 780$ nm, which might suggest a larger Leaf Area Index (LAI) [Delegido et al., 2013] and a greater vegetation vigor [Boochs et al., 1990] for those plots. On the basis of these observations, the spectral shape of high richness plots, as compared with spectra of intermediate and low richness plots, may imply a generally richer vegetation. This is consistent with findings reported in literature that high species richness enhances primary productivity [Wang et al., 2016, Grace et al., 2016] and biomass [Malhi et al., 2020, Tilman et al., 1997].

Fig. 5 shows the $r$ and RMSE values achieved with different numbers of components (ranging from 1 to 10) being selected as features. The results were averaged over the two study regions. It was observed that two components achieved the best performance for all the three feature extraction methods of PCA (Fig. 5a), CCA (Fig. 5b), and PLS (Fig. 5c). When only one component was used with more information in the original spectral data being discarded, lower $r$ values and higher RMSE values were also observed, indicating a poorer performance compared with that achieved by two components. When more than two components were retained with a higher degree of data redundancy presenting, weaker results were also observed. The eigenvalue, percentage of explained variance, and cumulative percentage of explained variance with different components are shown in Table 2. It was seen that with two components, the PCA, CCA, and PLS could explain 93.81$\%$, 89.34$\%$, and 87.38$\%$ of variance in the DESIS data, respectively. Based on these results, the number of components, $k$, was set to two in our experiments for dimensionality reduction with PCA, CCA, and PLS. The first and second components are plotted in Fig. 6. The first component depicted the general shape of the spectral brightness, while the second component highlighted more on local spectral features of the spectrum.

The grid searching result for the optimal KRR hyperparameters $\sigma$, $l$, and $\delta$ is shown in Fig. 7. The best combination of kernel parameters was $\sigma=10^{3}$, $l=10^{3}$, and $\delta=10$. This combination of kernel parameter values was used as default values for KRR in our experiments.

The assessment of different dimensionality reduction and regression methods (Table 3) showed that, for the Southern Tablelands region, the best result was achieved with a combination of PLS for dimensionality reduction and GPR for regression, with $r$ being 0.76 and RMSE being 5.89. Results with PLS as the dimensionality reduction method performed better than those with PCA and CCA. The PLS-based models achieved $0.75\sim 0.76$ for $r$ and $5.89\sim 5.92$ for RMSE, respectively, better than the PCA-based models with $0.70\sim 0.71$ for $r$ and $5.99\sim 6.02$ for RMSE, and the CCA-based models with $0.69\sim 0.73$ for $r$ and $5.98\sim 6.02$ for RMSE.

The assessment results for the Snowy Mountains region in shown in Table 4. The Snowy Mountains region is located at high altitudes with less human interference and a generally lower species richness than the Southern Tablelands (Fig. 2). Compared with the results for the Southern Tablelands in Table 3, it can be seen that generally the $r$ values were lower and RMSE were higher for the Snowy Mountains region (Table 4). The best $r$ and RMSE were 0.68 and 5.95, respectively, achieved with a combination of PLS for dimensionality reduction and RFR for regression.

The correlation diagrams in Fig. 8 show the relationship between the ground-truth species richness and the predicted values from the DESIS spectra, with data samples from the validation set. These results were produced with the best performing models in Tables 3 and 4.

Figs 8a and 8b show results with modelling being conducted for the Southern Tablelands and Snowy Mountains regions separately. When we modelled the pooled data from both regions, we found that accuracy decreased (Fig. 9). The prediction result was $r=0.61$ and $\text{RMSE}=10.1$, which was lower than modelling with data from only one region (Figs 8a and 8b). This indicates that location-specific modelling performs better than using one model to describe multiple regions.

The difficulty of modelling relationships between hyperspectral data and plant species richness for multiple regions, compared modelling for each region separately, is worth further investigation. Different regions may differ in the assemblages of plant species, and their compositional and structural properties, resulting in extra variations in hyperspectral data in addition to those induced by the richness of species. These additional variations may add complexity in exploring useful information in hyperspectral data to predict plant species richness. The location-specific relationship hyperspectral data and plant species richness calls for location-dependant modelling or encoding location information into the input spectra in future studies when mapping plant species richness at continental or global scales is attempted.

The relative importance of DESIS bands in predicting plant species richness is shown in Fig. 10, with the aim to analyse which parts of the spectrum had the most explanatory power. Subplots 10a, b, and c display the results with PCA, CCA, and PLS being used as the feature extraction procedure, respectively. From these subplots, it was observed that bands in the red-edge spectral region of approx. 700 $\sim$ 720 nm showed the highest importance. This may suggest that the slope of the red-edge is an important spectral feature for plant species richness prediction, considering that the low, intermediate, and high species richness plots showed varied red-edge slopes in Fig. 4. This observation might be supported by literature that the red-edge is a critical spectral region for vegetation mapping as it is closely related to biological variables such as leaf area index [Delegido et al., 2013], plant vigour [Boochs et al., 1990], and biochemical contents [Mutanga and Skidmore, 2007]. Followed by the red-edge, the visible range of approx. 400 $\sim$ 700 nm also showed high importance. The importance of visible bands in plant species richness prediction might be justified by the fact that, this portion of the spectrum, especially bands in red and blue, is the major leaf pigment absorption range. It is sensitive to mainly chlorophyll a and b contents, according to the sensitivity analysis result reported in Zhao et al. [2014b, a]. Though less important than the red-edge and visible regions, the near-infrared region with wavelengths longer than 720 nm also showed some explanatory power. This indicated that the near-infrared region also provided contributory information in predicting plant species richness, as this portion of the spectrum, often characterised by high reflectance for vegetation, is sensitive to leaf thickness [Zhao et al., 2014a] and the amount, arrangement, and inclination of leaves in the canopy [Knipling, 1970].

Figs 11a and 11b show correlation diagrams between the ground-truth species richness and values predicted from the real Sentinel-2 multispectral data for the Southern Tablelands and Snowy Mountains regions, respectively. It is seen that, a prediction result of $r=0.66$ and $\text{RMSE}=6.27$ was achieved for the Southern Tablelands region (Fig. 11a), and $r=0.57$ and $\text{RMSE}=6.31$ for the Snowy Mountains region (Fig. 11b). The prediction results with the Sentinel-2-like synthetic data set are shown in Fig. 12. For the Southern Tablelands, the $r$ and RMSE values are 0.65 and 6.19, while for the Snowy Mountains, the $r$ and RMSE values are 0.57 and 6.39.

These multispectral results were slightly worse than the hyperspectral results shown in Fig. 8. This could be explained by the fact that hyperspectral data contain richer information in the spectral domain. However, when taking data availability and stability into consideration, spaceborne multispectral imagery could serve as a reliable alternative for plant biodiversity mapping. Also, data such as Sentinel has much more consistent and regular sampling over time, and harnessing the added information in temporal changes in multispectral signal could provide further explanatory power for biodiversity analyses.

Predicting plant biodiversity from remotely sensed measurements helps with making evidence-informed environmental policies and conducting effective conservation activities. Hyperspectral imagery delivered by the recently launched DESIS instrument opens the potential for plant biodiversity monitoring at a finer spatial scale and with a higher accuracy. In this study, we take the the Southern Tablelands and Snowy Mountains regions in southeast Australia as experimental sites to test the relationship between DESIS data and plant species richness. A two step approach is proposed, where feature extraction techniques are first used to reduce the dimensionality of hyperspectral data, followed by regression models with kernel functions to account for the linearity, non-linearity, and noise in data.

Obtaining informative features from hyperspectral data is important to the success of subsequent interference of plant species richness. In previous studies, features have been primarily selected as a subset of the original bands, or as spectral indices computed from a subset of the original bands (e.g., Malhi et al. [2020], Peng et al. [2018], and Wang et al. [2016]). The bands or indices are often selected based on our a priori knowledge of the spectral properties, such as absorption features of biochemical contents. Though selected features usually offer good explainability, it means we have to discard information in unselected bands. Considering the large number of bands in hyperspectral data, many of them would be discarded as the high collinearity of hyperspectral data often requires a considerable reduction of dimensionality. In contrast, the feature extraction approach, as adopted in this work, makes use of the information in all original bands by transforming them into a new feature space of lower and non-collinear dimensionality.

Linear models have been primarily employed in previous studies to relate features to species richness. For example, multiple linear regression has been adopted by Wang et al. [2016] and Malhi et al. [2020], and stepwise linear regression by Peng et al. [2018]. However, the relationship between features and richness might not necessarily follows simply a linear pattern. In our study, a novel kernel function is proposed (Eq. 3), with the dot-product and radial-basis function kernels account for the linearity and non-linearity of the data, respectively, and the white kernel explains the noise in the data. The ability of our model to explore both linear and non-linear patterns distinguishes our work from aforementioned studies.

Though detailed field surveys have been deemed as the most accurate and reliable way for assessing plant biodiversity, remote sensing data can serve as a proxy for large-scale and cost-effective biodiversity mapping. The mechanism justifying the use of remote sensing has been widely discussed in literature. It is worth noting the exact mechanism is dependant on which method is used and the sensor specifications (e.g., pixel size and spectral range). The richness of species can be estimated either directly from raw spectra, or via the extracted plant functional traits or types [Wang and Gamon, 2019]. Methods based on the spectral variation hypothesis is also intriguing, whereby the spectral variation across spatially adjacent pixels is employed as the proxy [Palmer et al., 2002, Fassnacht et al., 2022]. Here we focus on discussing the underlying mechanism that underpins our study where the DESIS spectra with a pixel size of 30 m and a spectral range of 400–1000 nm are linked to on-ground species richness via feature extraction and statistical regression.

It has been reported in literature that plant communities of high diversity tend to have an enhanced primary productivity [Wang et al., 2016, Grace et al., 2016] and a higher above ground biomass [Malhi et al., 2020, Tilman et al., 1997]. Though the reason to explain the richness–productivity/biomass relationship is a subject of debate, a common theory is that the complementary roles played by different species lead to lower nutrient losses and more sustainable soils [Tilman et al., 1996]. Species complementarity allows plants to capture resources in ways that are complementary in space or/and time, leading to increased biomass production [Cardinale et al., 2007]. For example, a high number of species allows plants to reside in various partitions of niches, resulting in a denser occupation of space and a higher efficiency of water, nutrition, and sunlight usage.

Based on this relationship, many studies have successfully estimated plant species richness from hyperspectral measurements (e.g., Wang et al. [2016], Malhi et al. [2020], and Peng et al. [2018]), given that remotely sensed hyperspectral data is a good proxy of vegetative biomass and primary productivity. A positive and dynamic productivity–diversity relationship is observed in a prairie grassland experiment at Cedar Creek, Minnesota, USA, with NDVI being employed as a proxy of vegetation productivity to estimate species richness [Wang et al., 2016]. The correlation between spectral indices and plant species diversity is also reported in Peng et al. [2018] for a semi-arid sandland ecosystem in Inner Mongolia, China.

It is worth noting that previous studies are primarily focused on ground-measured (e.g., [Peng et al., 2018] and Wang et al. [2016]) or airborne (e.g., Asner [2008]) hyperspectral data, with a limited spatial range. The recently launched DESIS and PRISMA [Pignatti et al., 2013, Verrelst et al., 2021], and the upcoming EnMAP (Environmental Mapping and Analysis Program) [Guanter et al., 2015] missions, enable us to test the potential of spaceborne hyperspectral measurements in plant species richness mapping. Our study shows a promising correlation between the two, and finds that the correlation is location-dependent.

Over the past decades, in-situ samples of plant species richness have been collected via various survey campaigns. In total, more than 188 thousands of samples have been gathered in Australia as of the year of 2018 [Gellie et al., 2018]. However, most of the samples are not able to be matched with a DESIS observation that is temporally close enough to them, as DESIS has not been in operation until 2018. In this study, in order to avoid large temporal discrepancies, we have limited our on-ground samples to those that are less than three years apart with their associated DESIS spectra. The limited spatial coverage of DESIS images, and bushfires during the 2019–2020 summer, have further reduced the number of available samples. As a result, analyses in this study are conducted on a relatively small number of samples. Nevertheless, the ever increasing amount of both satellite images and on-ground samples will enable more comprehensive assessment of the potential of spaceborne hyperspectral remote sensing for plant biodiversity mapping. Moreover, though the data set used in this work does not have information on which strata the observed species come from, it is worth splitting out richness for different strata in future field surveys. With strata information on record, we might be able to model tree richness and understory richness separately, and then combine them once predicted to get total richness. The results reported in this study may serve as a basis for future studies.

It is important to note that the DESIS pixel size ($30\times 30$ m) does not match exactly with the plot area of ground species richness samples ($400\;\text{m}^{2}=20\times 20\;\text{m}$). Though richness measurements could be scaled to a larger or smaller plot area using an assumed power relationship ($S=cA^{z}$) between species richness ($S$) and plot area ($A$) [Rosenzweig, 1995], the location-specific power parameter ($z$) is often hard to be accurately determined without adequate knowledge about the experiment site. It is suggested that, in future field campaigns, it is worthwhile to gather information on the richness-area relationship, in order to facilitate accurate up- and down-scaling of plot areas. Potential inaccuracies in geo-registration of DESIS images and in geo-positioning measurements during field surveys may also result in geo-mismatch between pixels and ground plots. Better geo-registration and geo-positioning accuracies would help reduce uncertainties in predicting species richness in future works.

The DESIS sensor covers the visible and near-infrared (VNIR) portions of the spectrum. As shown in this study, this spectral range is informative in plant biodiversity mapping. However, it is also worthwhile to leverage the potential of the short-wave-infrared (SWIR) bands. The upcoming EnMAP imaging spectroscopy mission [Guanter et al., 2015] is scheduled to be launched in 2022. The EnMAP dual-spectrometer instrument, covering both VNIR and SWIR from 420 nm to 2450 nm (as shown in Fig. 3), will provide an opportunity to integrate information from SWIR for plant biodiversity mapping. In addition to hyperspectral optical imagery, the combination of data from other sensor types, such synthetic-aperture radar (SAR) or LiDAR data, could also be explored to improve our ability in remote mapping of plant biodiversity.