LD-SDM: Language-Driven Hierarchical Species Distribution Modeling

## Abstract

Abstract We focus on the problem of species distribution modeling using global-scale presence-only data. Most previous studies have mapped the range of a given species using geographical and environmental features alone. To capture a stronger implicit relationship between species, we encode the taxonomic hierarchy of species using a large language model. This enables range mapping for any taxonomic rank and unseen species without additional supervision. Further, we propose a novel proximity-aware evaluation metric that enables evaluating species distribution models using any pixel-level representation of ground-truth species range map. The proposed metric penalizes the predictions of a model based on its proximity to the ground truth.
We describe the effectiveness of our model by systematically evaluating on the task of species range prediction, zero-shot prediction and geo-feature regression against the state-of-the-art.
Results show our model outperforms the strong baselines when trained with a variety of multi-label learning losses.

## 1 Introduction

The ever-increasing availability of crowdsourced and remote sensing data has opened opportunities to solving challenging tasks that were not possible before. One such task is global-scale joint species distribution modeling (SDM). Understanding species ranges is not only essential for biodiversity conservation but also for understanding the effects of climate change. However, complex biophysical and geographic processes make SDM a challenging task.

Traditionally, species ranges were mapped using geographical and environmental features. While the current approaches have their merits, we contend they fall short in accurately representing the inter-species relationships inherent in the taxonomic hierarchy. For example, *species* belonging to the same *genera* tend to be found in similar locations. They often share many common characteristics, including their general habitat preferences []. Nevertheless, this pattern diminishes as one moves up in the taxonomic hierarchy. To capture these complex hierarchical relationships between the species, we use large language models (LLMs) to encode their taxonomic information in text (as shown in Figure [1](#S1.F1)). Intuitively, this should enhance the identification of rare taxa []. To the best of our knowledge, no previous work has incorporated large language models for SDM.

Figure: Figure 1: Species Distribution Modeling using LLM Guidance. We use text prompts in the form of key-value pairs for training our species distribution model. This kind of prompt captures the taxonomy hierarchy of species. We use only the exhaustive prompt (d) for training while the rest are used for zero-shot predictions.
Refer to caption: extracted/5293161/imgs/front_page.jpg

In order to incorporate species-specific embeddings into SDM, it is necessary to reformulate the SDM problem. Previous approaches to SDM [, , ] tried to answer the question, “which species are likely to be found at a given location?”. In contrast, our approach to SDM models, “which locations are likely for a given species to be found?”. Assuming uniform prior over the species distribution [, ], this is an equivalent formulation as before. But by doing so, we can now incorporate a variety of information, including species-specific embeddings or any other metadata, which opens up new opportunities for more accurate modeling. Alternatively, this formulation makes it easier to incorporate unseen species that were not present in the training data.

Previously, SDMs were trained using environmental features while ignoring the geographic location aspect. Recently, several works showed the benefit of incorporating geographic locations in various applications including SDM [, , ] and fine-grained image classification [, , ]. A series of works [, ] proposed spherical harmonics representation of geographic location.
Inspired by these works, we propose to use the Spherical Fourier Neural Operators (SFNO’s) [] as a geographical and environmental feature extractor for SDM.

Figure: Figure 2: Proposed Architecture. During training, features extracted from the environmental covariates are added to a learnable positional embedding and passed to the spherical Fourier neural network (SFNO). The species range map is obtained by computing cross-attention between the representations from SFNO and embeddings obtained from LLaMA-2 for a text description of a given species. Inference requires only cross-attention computation if embeddings are precomputed.
Refer to caption: extracted/5293161/imgs/arch_ihsdm.png

When it comes to predicting species range, the evaluation process typically involves using range maps that are defined by experts or curated carefully [, ]. However, these maps can often be sparse and difficult to obtain. Alternatively, using any crowdsourced species observation data can also be challenging as there is no precise definition of species absence. Intuitively, models with false predictions close to their true locations are better than those with distant false predictions. This is because models in the former case learn *spatially localized* species distributions. Previous studies have used the area under the curve (AUC) [, ] and the mean average precision (MAP) [] to evaluate SDMs. However, these metrics ignore the spatial distribution of predictions and weigh them equally, making them highly sensitive to spatial shifts in the predictions. Even a single pixel-level shift in predictions could add significant bias to the scores.

To this end, we propose a proximity-adjusted evaluation metric that penalizes predictions based on their distance to the closest true locations. We construct a One-Way Probability-Weighted Chamfer Distance (PWCD) metric, that measures the dissimilarity between the true and predicted species range maps, which is robust to small shifts. This enables an effective evaluation of SDMs using pixel-level range maps and crowdsourced data (where absence is not confirmed).

The contribution of our work is threefold:

- 1.
We use large language models to encode taxonomic information of species, which allows mapping over any taxonomic level.
- 2.
We reformulate the SDM problem which enables the use of species-specific embeddings. This allows inference on unseen species.
- 3.
We propose a novel proximity-aware metric for assessing the performance of SDMs using pixel-level range maps.

## 2 Related Works

Species distribution modeling (SDM) is a set of techniques that uses environmental, geographic, and species observation data to predict species distribution across geographic space and time [, ]. In this context, we restrict ourselves to the task of species range mapping over geographic space only. In the past, SDM relied on hand-crafted geographic and environmental features as input to capture the implicit spatial correlations in data [, , ]. These approaches were limited to one species per model due to their reliance on traditional supervised learning algorithms. With advancements in deep neural network architectures, end-to-end frameworks are now able to model joint species distributions and complex implicit relations in data [, , ].

Another important aspect of SDM is the type of species occurrence data. Species occurrence data may exist in one of the two forms: *presence-only* and *presence-absence*. Collecting presence-absence observations is a challenging task, requiring expert surveys especially to confirm the absence of a given species in some region []. However, once collected, this kind of data can be easily fed into supervised learning frameworks, without requiring sophisticated modifications [, , , , ]. With no expert knowledge, presence-only data can be recorded when any species of interest is encountered. This way of collecting data does not require the data collector to confirm absences. However, since there are no negative samples, it is difficult to use standard supervised learning frameworks [, , ].

Most of the recent works formulate the problem of SDM using presence-only data as a multi-label classification task [, ]. Here, the goal of SDM is to predict the likelihood over species for a given location and environmental features. A variety of multi-label learning losses have been proposed in the literature that can easily be plugged in for this task. For instance, the AN-full (full assume negative) loss [] randomly samples locations that are assumed to be negative during training. The ME-full (full maximum entropy) [] tries to maximize the entropy of predictions of negative samples and at randomly sampled locations during training. Other losses such as ASL (asymmetric loss) [], Hill loss [], and RAL (robust asymmetric loss) [] weigh the negative samples in a way that easy negative and potential false negative samples are down-weighted. The authors of [] presented an extensive comparison of the performance of SDMs trained with losses involving pseudo-negative sampling. Here, we extend their work by additionally comparing their SDMs trained with the ASL and RAL losses.

Several studies have attempted to assess the performance of species distribution models trained on presence-only data [, , , , ]. They measure performance using Kappa index, Kullback–Leibler divergence, area under the curve (AUC), or mean average precision (MAP). However, these metrics are inadequate in quantifying the performance of models on pixel-level data. It is worth noting that none of these metrics consider the spatial proximity of predictions to their ground-truth observations, hence making them sensitive to spatial shifts in predictions.

Finally, another line of work uses high-resolution satellite imagery at each species occurrence location to guide SDM methods [, , , , , ]. For example, [] uses satellite imagery and geolocation to directly predict the likelihood over species while [] uses a contrastive learning framework to establish similarity between ground-level images of species and corresponding satellite imagery. These methods are limited in geographic scale due to the high-resolution nature of raw data []. Except [], all the works view this problem as an information retrieval task rather than a species prediction task.

## 3 Method

In this section, we describe our approach to hierarchical species distribution modeling with LLMs.

### 3.1 Problem Definition

We consider presence-only data in the form of occurrence records: $\mathcal{D}=\{(l_{i},y_{i})|\>i=1,...,N\}$, where $l_{i}\in\mathbb{R}^{2}$ is the geolocation of each observation and $y_{i}\in[1,...,S]$ is the species of the observed organism. Additionally, a set of environmental features $e_{i}\in\mathbb{R}^{n}$ is included for each location $l_{i}$. Many previous studies [, , ] have tackled the task of SDM by training a neural network to model the conditional probability of $p(y|\>l,e)$. This translates to the problem of “how likely a *species* can be found at a given *location*”. However, our approach involves modeling the probability of $p(r|\>y,l_{r},e_{r})$, which answers the question “Is it likely for a given species to be found at given location?”.
Here, $r\in[1,...,L]$ is a set of discrete locations that are of interest when predicting the range of species belonging to some class. In this framework, we define the set of locations $r$ as a contiguous grid that is laid over the surface of the globe. The variables $l_{r}$ and $e_{r}$ represent fixed parameters in the model, namely the location vectors and environmental features defined at each location in $r$ respectively. It is easy to see that the input $y$ and the parameters $(l_{r},e_{r})$ can have any arbitrary representations. So, this new formulation allows for flexible encoding of both the species class and geographic location, compared to the previous formulation which only allowed flexibility in representing the geographic location.

### 3.2 Architecture Overview

Our objective is to train a neural network of the form: $f_{\theta}(y;l_{r},e_{r}):\mathbb{R}^{d}\rightarrow[0,1]^{L}$, where $d$ is the dimensionality of the species feature vector.
$f_{\theta}$ is composed of three modules: 1) a location encoder network $f_{\alpha}$, 2) a species encoder network $g_{\beta}$, and 3) a multi-label classifier $h_{\gamma}$. See Figure [2](#S1.F2) for a high-level architecture diagram.
The network $f_{\alpha}$ consists of convolutional feature extractors and the Spherical Fourier Neural Operator (SFNO) [] blocks that take as input the fixed parameters: $l_{r}$ and $e_{r}$. The location vectors are used as learnable positional encodings for the SFNO blocks. The species encoder $g_{\beta}$ takes as input the text descriptions corresponding to a particular species of class $y$. The multi-label classifier $h_{\gamma}$ consists of a cross-attention network that fuses features from the SFNO and text encoder and performs multi-label classification using a linear layer. We call our model Language-Driven Species Distribution Model (LD-SDM).

### 3.3 Species-specific embeddings

Many popular crowdsourcing science platforms, such as iNaturalist [], contain not just occurrence records, but also valuable species-specific information. This information may include the complete taxonomic hierarchy of each species (such as class, order, family, etc.) or the IUCN red list category (such as endangered, near threatened, etc.). While often ignored, this fine-grained information could be incorporated to enhance the performance and interpretability of SDMs. In this work, we incorporate the taxonomic hierarchy of species into our SDM in the form of text. Each species is uniquely identified by its taxonomic classification ranks, which are represented as key-value pairs. For instance, a key-value pair that represents the *wild duck* is shown below:

t = {class: *Aves*, order: *Anseriformes*, family: *Anatidae*, genus: *Anas*, species: *Anas platyrhynchos*}

Recently, Large Language Models (LLMs) have demonstrated excellent capabilities in representing implicit knowledge, linguistic structure, and hierarchical concepts when trained on massive amounts of data []. With no finetuning, pretrained LLMs have shown impressive zero-shot performance with zero in-domain data. In this work, we use the recently released, LLaMA-2 [] to encode the species-specific key-value pairs as text. We experiment with three variants of LLaMA-2: 7B, 13B, and 70B.

### 3.4 Loss Functions

Our formulation of the SDM falls in the domain of multi-label long-tailed learning, where the model learns to predict the likelihood probability over *each* location for a given species. Since we are modeling the likelihood over locations, our problem differs from the Single Positive Multi Label (SPML) task []. Here, we are given multiple positive locations of occurrence for each species. We explore different types of losses for our models that do not require pseudo-negative sampling, such as ASL and RAL. However, for the state-of-the-art models, we use all types of losses, including AN-full and ME-full, thereby extending their work. We describe the losses in detail below.

Full Assume-Negative (AN-Full). This loss is a kind of pseudo-negative sampling multi-label learning loss based on the argument that the majority of the species are absent at any given location []. This is typically used in SPML problems. The loss is given by:

$$ $\displaystyle\mathcal{L}_{\text{AN-full}}(\hat{y},y)=-\frac{1}{S}\sum_{s=1}^{S }[\mathbbm{1}_{[y_{s}=1]}\alpha log(\hat{y}_{s})$ (1) $\displaystyle+\mathbbm{1}_{[y_{s}=0]}log(1-\hat{y}_{s})+log(1-\hat{r}_{s})]$ $$

where $\hat{r}$ is a uniformly randomly sampled location from the spatial domain of the globe.

Full Maximum-Entropy (ME-Full).
This loss was introduced in [] to solve SPML problems. The idea of this loss is to maximize the entropy of predictions for negative samples and for randomly sampled pseudo-negative samples. The loss is given by:

$$ $\displaystyle\mathcal{L}_{\text{ME-full}}(\hat{y},y)=-\frac{1}{S}\sum_{s=1}^{S }[\mathbbm{1}_{[y_{s}=1]}\alpha log(\hat{y}_{s})$ (2) $\displaystyle+\mathbbm{1}_{[y_{s}=0]}H(\hat{y}_{s})+H(\hat{r}_{s})]$ $\displaystyle H(\hat{y})=-[\hat{y}log(\hat{y})+(1-\hat{y})log(1-\hat{y})]$ $$

where again $\hat{r}$ is a uniformly randomly sampled location.

Asymmetric Loss (ASL). This loss, proposed by [], weights positive and negative samples differently. The positive and hard negative samples are given more weight. It is a fairly general loss for multi-label learning which can also be used for SPML, without the need for psuedo-negative sampling. This loss is given by:

$$ $\displaystyle\mathcal{L}_{\text{ASL}}(\hat{y},y)=-\frac{1}{S}\sum_{s=1}^{S}[ \mathbbm{1}_{[y_{s}=1]}(1-\hat{y_{s}})^{\gamma^{+}}log(\hat{y_{s}})$ (3) $\displaystyle+\mathbbm{1}_{[y_{s}=0]}(p_{s}^{\tau})^{\gamma^{-}}log(1-p_{s}^{ \tau})]$ $$

where $p^{\tau}=\text{max}(\hat{y}-\tau)$ is used for hard thresholding easy negatives.

Robust Asymmetric Loss (RAL). Recently, [] proposed an improved version of the ASL loss by combining the Asymmetric Polynomial Loss (APL) [] and the Hill loss [], which makes it less sensitive to hyperparameters. This loss is given by:

$$ $\displaystyle\mathcal{L}_{\text{RAL}}(\hat{y_{s}},y_{s})=[\sum_{m=1}^{M} \mathbbm{1}_{[y_{s}=1]}\alpha_{m}(1-\hat{y_{s}})^{m+\gamma^{+}}$ (4) $\displaystyle+\sum_{n=1}^{N}\mathbbm{1}_{[y_{s}=0]}\psi(\hat{y_{s}})\beta_{n}( p_{s}^{\tau})^{n+\gamma^{-}}]$ $$

where $\psi(\hat{y}_{s})=\lambda-\hat{y}_{s}$ is the Hill loss term. $M$, $N$, $\alpha_{n}$ and $\beta_{n}$ are balance parameters for positive and negative samples.

## 4 Experiments

In this section, we describe the implementation details and performance of our models compared to the state-of-the-art models on the task of species range prediction and geo-feature regression.

### 4.1 Implementation Details

As described in section [3.2](#S3.SS2), we use SFNO to extract geographic and environmental-specific features. The architecture is borrowed from []. We use two spectral and encoder layers in the SFNO. The encoder dimension used is 128 while the number of SFNO blocks are two. For the coordinate-only models, we drop the environmental features during the training. For the text encoder, we use *frozen* LLaMA-2 [] (70B variant) and extract its encoder outputs for each species. Before training, we precompute and save the embeddings for each species in the disk. For the multi-label classifier, we use a single-layer multi-headed cross-attention network followed by a fully-connected layer with sigmoid activations. We compare LD-SDM with two state-of-art-models: SINR [] and SIREN(SH) []. SINR learns an end-to-end neural representation of geographic locations that is used to predict the likelihood over species. SIREN(SH) uses a spherical harmonic representation of geographic locations which is then used to predict the likelihood over species. The rest of the details about our models and the implementation of state-of-the-art models are described in Appendix [C](#A3).

### 4.2 Proximity-Aware Evaluation

Evaluating SDMs is a challenging task since it requires the comparison of spatially distributed probability values. Pixel-wise comparison methods are not well-suited for this task since they are susceptible to spatial prediction shifts. Intuitively, false positive predictions that are geographically distant from true positive predictions are more significant than those that are close.
This is not inherent in the evaluation metrics used by previous studies such as the Kullback–Leibler divergence, AUC, or MAP, which are invariant of the spatial arrangement of probability mass.

A well-known spatially explicit distance used for comparing any two distributions is the Earth-Mover (Wasserstein) Distance. Due to its computational requirements and infeasibility on large-scale data points, Chamfer Distance is a popularly used proxy []. However, it has only been developed for point clouds and there doesn’t exist a chamfer distance measure to compare two spatial probability distributions. Inspired by [], we make a simple modification to their Density-Aware Chamfer Distance (DCD) metric that works for comparing a *predicted* probability distribution map ($P_{1}$) with a *binary* probability map ($P_{2}$) representing the ground truth. We add a likelihood probability term which weights the distance term. We call this metric One-Way Probability Weighted Chamfer Distance (PWCD), which measures the false positive rate in a proximity-aware manner. This metric is between [0,1] and could also be used to compute the number of false positives instead of the rate. It is given by:

$$ $d_{\text{PWCD}}(P_{1},P_{2})=\frac{1}{|N|}\sum_{x\in P_{1}}\left(1-e^{-\alpha p (x)||x-y||_{2}}\right)$ (5) $$

$$ $FP_{\text{PWCD}}(P_{1},P_{2})=\sum_{x\in P_{1}}\left(1-e^{-\alpha p(x)||x-y||_ {2}}\right)$ (6) $$

where, $y=\underset{z\in P_{2}\cap p(z)=1}{argmin}||x-z||_{2}$. This represents the nearest neighbor ground-truth observation of a given pixel $x$ in the predicted probability distribution map. $N=\{z\in P_{2}|\>p(z)=0\}$ is the set of all negative observations in the ground truth. $\alpha$ is a temperature parameter that controls the exponential weighting of the distance. See Appendix [B](#A2) for more details. For our case, $|P_{1}|=|P_{2}|=L$. Notice that we have omitted the density term in the metric as our probability maps are a contiguous grid of pixels and equally dense everywhere. Since it is a proxy for false positives, it can be plugged into any other metric of choice such as the MAP and AUC. In this work, we report proximity-adjusted MAP scores denoted as MAP(pa) by replacing false positives with our metric. All the scores reported are an average over scores computed for each species category.

**Table 1: Comparison of performance of SDMs over a variety of multi-label learning losses on two tasks: species range prediction (and zero-shot range prediction) and geo-feature regression. The proximity-aware metrics proposed are defined in Section [4.2](#S4.SS2).**
|  |  |  | GBIF’21 | GBIF’22 | Zero Shot | Geo Feature |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Loss | Env. | Method | MAP(pa) $\uparrow$ | $d_{\text{PWCD}}$ $\downarrow$ | MAP(pa) $\uparrow$ | $d_{\text{PWCD}}$ $\downarrow$ | MAP(pa) $\uparrow$ | $d_{\text{PWCD}}$ $\downarrow$ | (Mean $R^{2}$) |
| $\mathcal{L}_{\text{AN-full}}$ [] | ✗ | SINR [] | 62.64 | 0.177 | 60.28 | 0.178 | - | - | 0.446 |
| $\mathcal{L}_{\text{ME-full}}$ [] | ✗ | SINR | 61.88 | 0.189 | 58.81 | 0.187 | - | - | 0.518 |
| $\mathcal{L}_{\text{ASL}}$ [] | ✗ | SINR | 67.11 | 0.152 | 64.55 | 0.155 | - | - | 0.619 |
| $\mathcal{L}_{\text{RAL}}$ [] | ✗ | SINR | 62.02 | 0.291 | 60.42 | 0.297 | - | - | 0.536 |
| $\mathcal{L}_{\text{AN-full}}$ | ✗ | SIREN(SH) [] | 61.99 | 0.203 | 60.67 | 0.215 | - | - | 0.562 |
| $\mathcal{L}_{\text{ME-full}}$ | ✗ | SIREN(SH) | 64.44 | 0.197 | 64.16 | 0.212 | - | - | 0.630 |
| $\mathcal{L}_{\text{ASL}}$ | ✗ | SIREN(SH) | 62.56 | 0.189 | 62.43 | 0.206 | - | - | 0.481 |
| $\mathcal{L}_{\text{RAL}}$ | ✗ | SIREN(SH) | 63.13 | 0.192 | 62.98 | 0.209 | - | - | 0.488 |
| $\mathcal{L}_{\text{ASL}}$ | ✗ | LD-SDM (ours) | 73.88 | 0.174 | 71.45 | 0.182 | 59.27 | 0.648 | 0.583 |
| $\mathcal{L}_{\text{RAL}}$ | ✗ | LD-SDM (ours) | 73.46 | 0.176 | 71.52 | 0.184 | 55.22 | 0.694 | 0.587 |
| $\mathcal{L}_{\text{AN-full}}$ | ✓ | SINR | 64.25 | 0.184 | 62.43 | 0.181 | - | - | - |
| $\mathcal{L}_{\text{ME-full}}$ | ✓ | SINR | 64.13 | 0.200 | 60.67 | 0.222 | - | - | - |
| $\mathcal{L}_{\text{ASL}}$ | ✓ | SINR | 67.23 | 0.172 | 69.49 | 0.179 | - | - | - |
| $\mathcal{L}_{\text{RAL}}$ | ✓ | SINR | 64.98 | 0.224 | 63.77 | 0.252 | - | - | - |
| $\mathcal{L}_{\text{ASL}}$ | ✓ | LD-SDM (ours) | 74.16 | 0.222 | 72.54 | 0.231 | 61.14 | 0.702 | - |
| $\mathcal{L}_{\text{RAL}}$ | ✓ | LD-SDM (ours) | 75.26 | 0.219 | 74.13 | 0.224 | 59.65 | 0.715 | - |

**Table 2: Comparison of performance of our models using three different variants of LLaMA-2. The 70B variant achieves the best performance followed by the 7B and 13B variants.**
|  |  | GBIF’21 | GBIF’22 | Zero Shot | Geo Feature | FLOPS |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Text Encoder | encoding dim | MAP(pa) $\uparrow$ | $d_{PWCD}$ $\downarrow$ | MAP(pa) $\uparrow$ | $d_{PWCD}$ $\downarrow$ | MAP(pa) $\uparrow$ | $d_{PWCD}$ $\downarrow$ | (Mean $R^{2}$) | TFLOPS |
| LLaMA-2-7B | 4096 | 71.24 | 0.177 | 70.25 | 0.178 | 56.27 | 0.719 | 0.579 | 0.174 |
| LLaMA-2-13B | 5120 | 70.11 | 0.184 | 69.09 | 0.193 | 55.70 | 0.732 | 0.563 | 0.342 |
| LLaMA-2-70B | 8192 | 73.88 | 0.174 | 71.45 | 0.182 | 59.27 | 0.648 | 0.583 | 1.848 |

### 4.3 Training and Evaluation

We focus on bird species occurrence data which is abundantly available yet challenging to model as compared to other species due to complex migration patterns. We compiled 131.7 million observations of 4141 species of birds using the GBIF ([gbif.org](gbif.org)) platform which includes data from various community science platforms. This data is collected for the years 2012-2020. For evaluation, we use curated research-grade observations from the iNaturalist [] and eBird [] for 2021 and 2022. These contain 17.9 million and 17.2 million observations respectively. Using the same platform, we created a dataset of 84 species that are unseen (not present in training data) and rare ($<$ 1000 observations) for zero-shot performance evaluation. This contains a total of 26k observations. These platforms are useful since errors in misclassification are highly mitigated due to community consensus [], hence being a strong representative of the ground truth. Further, we use the same procedure as described in [] for the geo-feature regression task.

All the training and evaluations are done at a spatial resolution of $0.2^{o}$. This implies $L=1,620,000$. During training, we create ground-truth pixel-level range maps for each species at $0.2^{o}$ resolution using the observations present in the occurrence data. For this, we use the histogram binning technique and clip all the values greater than one to one []. These maps are then used to calculate the various multi-label learning losses described previously. Finally, we report MAP(pa) and $d_{PWCD}$ using $\alpha=0.1$ for the species range prediction task and $R^{2}$ for the geo-feature regression task. See Appendix [C](#A3) for more details on the training settings and hyperparameters used.

Figure: Figure 3: Predicted Range Maps (Probability Values). The predictions are generated at a spatial resolution of 0.2 degrees. The best performing models are used in each case. For the SIREN(SH) and SINR models, predictions are generated at the center of each raster cell. The visualizations show that our model performs better than other models at resolving higher frequencies over the sphere and localizing the species distribution.
Refer to caption: extracted/5293161/imgs/Results.jpeg

### 4.4 Results

Quantitative Evaluation. In Table [1](#S4.T1), we report the performance of all the models considered in this work. LD-SDM has shown better performance than the other models on the task of species range prediction. The zero-shot performance of LD-SDM on unseen species is reasonably good. Although it has an increased false-positive rate, it can be attributed to the model using higher-level taxonomic relationships. $\mathcal{L}_{\text{RAL}}$ tends to produce overconfident predictions as noticeable from the high MAP(pa) scores corresponding to high $d_{PWCD}$ scores (especially for SINR and SIREN(SH)). On average, $\mathcal{L}_{\text{ASL}}$ loss is likely to be better performing than the other losses over all the models. For SINR, $\mathcal{L}_{\text{ASL}}$ performs the best as compared to the other losses by a large margin. This is also evident in the geo-feature regression task. Our model has also produced competitive scores on the task of geo-feature regression. The coordinate-only models perform nearly as well as models incorporating environmental features. So, dropping environmental features could result in faster training times.

**Table 3: Comparison of taxonomic mapping performance of our models using the three different variants of LLaMA-2 for taxonomic encoding. These performances are reported on the combined presence-only dataset from GBIF for the years 2012-2022. It is worth noting that LLaMA-2-13b has a lower performance than LLaMA-2-7b, possibly due to differences in their training data.**
|  | LLaMA-2-7B | LLaMA-2-13B | LLaMA-2-70B |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Taxonomic Level | MAP(pa) $\uparrow$ | $d_{PWCD}$ $\downarrow$ | MAP(pa) $\uparrow$ | $d_{PWCD}$ $\downarrow$ | MAP(pa) $\uparrow$ | $d_{PWCD}$ $\downarrow$ |
| Class | 97.10 | 0.117 | 95.2 | 0.204 | 97.52 | 0.163 |
| Order | 56.15 | 0.450 | 55.94 | 0.472 | 58.31 | 0.414 |
| Family | 38.18 | 0.659 | 36.47 | 0.688 | 40.13 | 0.632 |
| Genus | 24.95 | 0.710 | 24.97 | 0.711 | 27.86 | 0.692 |

**Table 4: Comparison of performance of three different backbones used in the parameterization network. All the models are trained with the same loss with environmental features included in the inputs.**
|  | Params | GBIF’21 | GBIF’22 |  |  |
| --- | --- | --- | --- | --- | --- |
| Backbone | M | MAP(pa) $\uparrow$ | $d_{PWCD}$ $\downarrow$ | MAP(pa) $\uparrow$ | $d_{PWCD}$ $\downarrow$ |
| Segformer [] | 20.81 | 71.64 | 0.247 | 70.54 | 0.255 |
| DINOv2 [] | 41.34 | 72.19 | 0.232 | 72.01 | 0.240 |
| SFNO [] | 19.19 | 74.16 | 0.222 | 72.54 | 0.231 |

**Table 5: SFNO ablation in LD-SDM. All experiments run with the $\mathcal{L}_{\text{ASL}}$ loss and the LLaMA-2-70B variant.**
|  | GBIF’21 | GBIF’22 |  |  |
| --- | --- | --- | --- | --- |
| encoder dim | MAP(pa) $\uparrow$ | $d_{PWCD}$ $\downarrow$ | MAP(pa) $\uparrow$ | $d_{PWCD}$ $\downarrow$ |
| 64 | 68.55 | 0.147 | 67.28 | 0.149 |
| 128 | 73.88 | 0.174 | 71.45 | 0.182 |
| 256 | 73.67 | 0.178 | 71.82 | 0.188 |

Qualitative Examples. In Figure [3](#S4.F3), we show global-scale species distribution maps generated for five bird species. To make fair comparisons, we used the best-performing coordinate-only models to ensure no bias is added into models from the environmental features. Although SINR and SIREN(SH) could generate species distribution maps at any spatial resolution, they tend to learn a very low-frequency function over the sphere as evident from the blob-like maps. SIREN(SH) produces a lot of false positive predictions possibly due to the limitation of the number of Legendre polynomials it can fit. Our model on the other hand learns to predict species ranges at a more fine-grained level. This is evident from the visualization of the geographic location features learned by the location encoder of each model (shown in Appendix [D](#A4)). Fine-grained geographic features are important for better geolocalization of species. SINR and SIREN(SH) tend to learn similar features over large geographic regions while LD-SDM produces higher variability in the features at a finer scale.

Effect of LLaMA-2 variants on LD-SDM. We assess the impact of using different variants of LLaMA-2 for LD-SDM in Table [2](#S4.T2). In particular, we compare the 7B, 13B, and 70B variants. Surprisingly, the 13B variant performs the worst. This is possibly due to the differences in the training data of the variants. Empirically, this can be explained by examining the T-SNE plots of text embeddings. In Appendix [A](#A1), we present T-SNE plots along with mean intra-cluster and inter-cluster distances of text embeddings. Ideally, both intra-cluster and inter-cluster distances should be high to ensure separability within a particular taxonomic class and between taxonomic classes respectively. 13B variant has the poorest separability at every taxonomic rank. The 7B and 70B variants have almost identical separatabilities. The impressive performance of the 70B variant can be attributed to its high dimensionality of the species-specific encoding. Although 70B has noticeably higher computational overheads, it is compensated by precomputing the text embeddings before training. Overall, all the variants are able to perform better than the state-of-the-art.

Taxonomic range mapping using LLaMA-2 variants. We perform range prediction over the taxonomic ranks: *class*, *order*, *family*, and *genus*. The procedure is the same as for generating species range maps. We precompute the embeddings for each rank as shown in Figure [1](#S1.F1). The range map for the rank *class* represents the unconditional distribution of finding any bird *species* in the world. As one moves down in the hierarchy, the range maps become more and more fine-grained. As expected, the performance of LD-SDM drops as one moves down in the hierarchy (Table  [3](#S4.T3)). This can be explained by examining the increasing $d_{\text{PWCD}}$ scores. The T-SNE plots of the embeddings of various ranks in Appendix [A](#A1) help in better understanding the increase in false positive rates as one moves down in the hierarchy. As it can be seen, the intra-cluster distances decrease substantially when moving down in the hierarchy. As a result, it becomes challenging for the model to distinguish between birds for example of *genera* belonging to the same *family*. On the other hand, when training in a fully supervised manner at the *species* rank, the model is able to achieve impressive performance. This shows the importance of supervision in SDM.

Spherical Harmonics are important for SDM. We compare the performance of LD-SDM when using different backbones in the location encoder network (Table [4](#S4.T4)). All the models are trained from scratch using the LLaMA-2-70B variant. We use the same training inputs including the environmental features. All the models are trained using the $\mathcal{L}_{\text{ASL}}$ loss. Notice that SFNO gives the best performance over Segformer and DINOv2. With SFNO, we observe a significant gain in performance and a reduction in false positives. Overall, a robust geographic feature extractor is beneficial for SDM and SFNO proves to be one such architecture.

Lightweight SFNO performs surprisingly well. In Table [5](#S4.T5), we show the effect on the performance of species range prediction due to different components of the SFNO. It is noteworthy that an encoder dimension of 128 for SFNO yields results on par with an encoder dimension of 256 while reducing computational complexity. On the other hand, having fewer encoder layers in SFNO is beneficial since performance gain from adding additional encoder layers is marginal. In conclusion, a lightweight SFNO is a strong geographic feature extractor that can be used in various downstream global-scale geospatial tasks.

## 5 Conclusions

We introduced a novel approach for species distribution modeling that uses a large-language model to generate a representation of a species. This provides flexibility to generate distributions at different levels of the taxonomic hierarchy and for unseen species.
In addition, we introduced a new evaluation metric, Probability Weighted Chamfer Distance. This metric computes the false positive rate in species range mapping tasks in a proximity-aware manner.
While our focus was on species distribution modeling, the integration of large-language models points to some intriguing potential future use cases. It might be possible, for example, to generate distribution models from arbitrary text, such as “a large white bird” or “a flock of red-winged blackbirds”. This is potentially more useful for non-expert users. Finally, alternative text prompts could be explored as future work which may better reason about a species and its dependence on a particular geographic location.

## Appendix A Species-Specific Embeddings

As described in the paper, we encode the taxonomic information of species in the form of text. Text is a dictionary consisting of key-value pairs representing each taxonomic rank. We use hugging-face’s library transformers to obtain all the LLaMA-2 variants used in the paper. We separately compute and save the embeddings of the key-value pairs for each variant. In Table [6](#A1.T6), we report the number of tokens generated by the LLaMA-2 variants for each taxonomic rank.

**Table 6: Total tokens generated by LLaMA-2 for each taxonomic rank. The dimensionality of each of these tokens depends on the type of LLaMA-2 variant used.**
| Taxonomic Rank | #tokens | #classes | #total tokens |
| --- | --- | --- | --- |
| order | 13 | 40 | 520 |
| family | 19 | 211 | 4009 |
| genus | 25 | 1393 | 34,825 |
| species | 34 | 4141 | 140,794 |

We create T-SNE plots in 2-D for visualizing the embeddings at each rank resulting from the variants (Figure [4](#A0.F4)). To create the plot, we do an average pool across the number of tokens’ dimension for each text embedding. We then use scikit-learn library to compute the T-SNE embeddings for each species and at every taxonomic rank. We color code the embeddings based on the taxonomic rank one level higher than the current rank. For example, the T-SNE embeddings for the rank family are color-coded by their order. Doing this provides a simple intuition for interpreting the separatability of species and taxonomic ranks.

We compute average inter-cluster and intra-cluster distances based on the color coding for each taxonomic rank. This is also done using the scikit-learn package. The inter-cluster distance represents the average dissimilarity between species (or any other rank) of different classes. The intra-cluster distance represents the average dissimilarity between species (or any other rank) of the same class. Intuitively, both distances should be high enough to ensure good separatability for SDM. As we see, the distances decrease as we move down in the hierarchy. This is expected since the distribution of each species becomes more and more fine-grained. The 13B variant has abruptly low distances compared to the other variants suggesting poor performance in SDM. This can possibly be due to differences in training data and the training procedure between the variants.

## Appendix B Probability Weighted Chamfer Distance

We proposed the Probability Weighted Chamfer Distance (PWCD) which is a spatially-aware proxy for false-positive rates when comparing two probability maps. This distance could also be used in other metrics such as precision, MAP, or AUC to make them spatially aware. In Listing LABEL:pseudocode, we describe NumPy-style pseudocode for computing this distance metric.

Figure: Listing 1: NumPy-style pseudocode for computing $d_{\text{PWCD}}$

We use the xrspatial’s proximity function for computing the nearest-neighbor distance map given the binary ground-truth range map. All the distances are measured in pixel units using Euclidean distance. The minimum possible distance is 0 which is for the positive pixels having species presence. The maximum possible distance depends on the spatial resolution of the range map. In our case, the spatial resolution of $0.2^{o}$ gives a maximum distance of 2012.46. It is noteworthy that this distance function could be replaced with any other distance function for example the haversine distance.

We visualized the pixel-level $d_{\text{PWCD}}$ map in Figure [5](#A0.F5). Notice that the predictions spatially close to the true presence of species are weighted less than the predictions that are distant. This weighting is due to the nearest-neighbor distances of the pixels making the metric spatially-aware. The chamfer distance map describes the weighting of predictions as a function of distance from the nearest presence pixel. The weight gradually increases as the distance from the nearest presence pixel increases. Note that the nearest-neighbor distance map differs from the chamfer distance map, which is in the range of [0,1]. The chamfer distance map is computed using the nearest-neighbor distance map, using the formula for $d_{\text{PWCD}}$ and a constant likelihood term equal to one.

In Figure [6](#A0.F6), we show the effect of temperature parameter $\alpha$ on the $d_{\text{PWCD}}$ metric. The figure shows chamfer distance as a function of nearest-neighbor distance at a fixed likelihood probability of one and varying values of $\alpha$. The parameter $\alpha$ controls the exponential weighting of the nearest-neighbor distance in the metric. Increasing $\alpha$ gradually increases the weight of nearby pixels. The metric starts penalizing more of the nearby pixels in the computation of false positive rates. A large value of $\alpha$ will make the metric more sensitive to spatial shifts in the predictions while a lower value of $\alpha$ will be more robust towards spatial shifts. However, an extremely low value of $\alpha$ will make the metric ignore distant false positives. As a result, an optimal value of $\alpha$ is essential to evaluate SDMs.

Figure: Figure 9: Features Learned by Location Encoder. We show the visualization of location embeddings learned by the location encoder of various models. For visualization, the embeddings are projected to a 3-dimensional space using Independent Component Analysis. SINR learns a low-frequency representation of location while SIREN(SH) and LD-SDM learn to localize at a fine-grained level.
Refer to caption: extracted/5293161/imgs/pos_embeds.jpg

In order to evaluate the localization of species distributions, we have plotted $d_{\text{PWCD}}$ scores with varying $\alpha$ in Figure [7](#A0.F7). The first plot shows the relationship between $d_{\text{PWCD}}$ and the nearest neighbor distance for varying $\alpha$ values. By analyzing this graph, one can determine the optimal value of $\alpha$ for their specific application. For our purposes, we selected $\alpha=0.1$, as it is robust to spatial shifts and accurately detects the distribution of species.

In the second plot, we have shown the $d_{\text{PWCD}}$ scores as a function of the likelihood probability, at a constant nearest neighbor distance of 100 pixels. As the value of $\alpha$ increases, the metric starts penalizing low likelihood values. This gives us more control over the trade-off between the distance and likelihood scores, allowing us to determine the optimal $\alpha$ value for our application.

## Appendix C Training

Below we describe the details of the datasets used. We also describe the implementation details of the models and the losses. Each training run was done on a single NVIDIA A40 GPU.

### C.1 Data

Our training and evaluation data was collected using various citizen science platforms such as the iNaturalist ([inaturalist.org](https://www.inaturalist.org/)) [] and eBird ([ebird.org](https://ebird.org/home)) []. The maps in Figure [8](#A0.F8) display the geographical locations of the observations for each dataset across the globe. They were compiled using the GBIF ([gbif.org](https://gbif.org/)) [] platform, which hosts occurrence records of species from several citizen science platforms. For compilation, we used the following conditions:

- 1.
All observations must have valid labels for each taxonomic rank and geographical location values.
- 2.
All observations must have the *research grade* status label.

To avoid bias, we excluded species with less than 100 observations. The evaluation datasets only contain observations from iNaturalist and eBird. For zero-shot evaluation, we selected the species that were not present in the training dataset. Further, we filtered out species with more than 1000 observations to focus on rare species. Table [7](#A3.T7) shows the details of each dataset.

**Table 7: Number of observations in each dataset compiled in this work.**
| Dataset | Year | #Observations | #species |
| --- | --- | --- | --- |
| Train | 2012-2020 | 131.7M | 4141 |
| Test | 2021 | 17.9M | 3960 |
| Test | 2022 | 17.2M | 4010 |
| Zero-Shot | 2021-2022 | 26k | 84 |

We use bioclimatic variables and elevation data as environmental covariates in the SDM. The variables are derived from the WorldClim dataset ([worldclim.org](https://test2.biogeo.ucdavis.edu/wc2/#/)) [], which are at a spatial resolution of $0.01^{o}$. The dataset is resampled to $0.2^{o}$ resolution using bilinear interpolation. In total, the dataset contains 20 different variables (channels). We use nine different geographic features in the geo-feature regression task as done by []. These include above-ground carbon, elevation, etc. All the features are resampled to $0.2^{o}$ resolution.

**Table 8: Training settings for SINR**
| Config | Value |
| --- | --- |
| backbone | ResNet-style MLP |
| encoding dim | 256 |
| depth | 4 |
| input encoding | sin-cos |
| optimizer | AdamW |
| base learning rate | 5e-4 |
| learning rate scheduler | exponential decay |
| batch size | 2048 |
| epochs | 10 |

**Table 9: Training settings for SIREN(SH)**
| Config | Value |
| --- | --- |
| backbone | SIREN |
| encoding dim | 256 |
| depth | 4 |
| input encoding | spherical-harmonics |
| L | 20 |
| optimizer | AdamW |
| base learning rate | 1e-4 |
| learning rate scheduler | exponential decay |
| batch size | 2048 |
| epochs | 20 |

**Table 10: Training settings for LD-SDM**
| Config | Value |
| --- | --- |
| backbone | SFNO |
| encoding dim | 128 |
| SFNO blocks | 2 |
| multi-label classifier | hidden_dim: 64 |
|  | num_heads: 16 |
| optimizer | AdamW |
| base learning rate | 1e-5 |
| batch size | 1 |
| accumulate grad batches | 64 |
| epochs | 10 |

### C.2 Models

We experimented with three different SDMs: 1) SINR []; 2) SIREN(SH) []; LD-SDM (ours). For the baseline models, we used the same hyperparameters as reported by authors. Tables [8](#A3.T8), [9](#A3.T9) and [10](#A3.T10), show the hyperparameters used for each model. We use $L=20$ for SIREN(SH) since experiments with larger $L$ were computationally infeasible.

### C.3 Losses

We evaluated four different multi-label learning losses in our experiments, namely: 1) AN-full []; 2) ME-full []; 3) ASL []; 4) RAL []. For AN-full and ME-full, we used a positive weight of 2048. The pseudo-negative sampling is done according to [], by using a spherical buffer around each observation. In the case of ASL and RAL, we used 0 and 4 for $\gamma^{+}$ and $\gamma^{-}$ respectively. In RAL, we set $\lambda$ to 1.5.

## Appendix D Location Embeddings

We present the location-specific features learned by the SDMs in Figure [9](#A2.F9). The location embeddings displayed are the output of the location encoder network for each SDM. We use Independent Component Analysis (ICA) to project the high-dimensional embeddings to a 3-dimensional space. All the maps are generated at a spatial resolution of $0.2^{o}$. SINR produces a smooth map indicating that it has failed to learn high-frequency functions. SIREN(SH) produces noisy location features spread across the globe. LD-SDM produces a map that resembles the distribution of training data. We see that in the locations where observations are absent, LD-SDM has learned similar features.

## Appendix E Limitations

As pointed out in [], the data obtained from the process of crowdsourcing is biased towards regions with heavy human traffic. Most often, rich biodiverse regions are underrepresented limiting the performance of SDMs. In this work, we do not explicitly deal with the observation bias in presence-only SDM. However, our goal was to incorporate taxonomic information into SDM to enhance unseen species discovery.

Previous studies have evaluated species SDMs using expert range maps, which are only available for a limited number of species. Our proposed evaluation metric allows for the evaluation of SDMs using crowdsourced data, which provides a more comprehensive representation of numerous species. However, we note that none of these methods are completely free of errors. Yet, with our proposed evaluation metric, innovative methods could be developed for semi-automatic evaluation of SDMs in the future, to strike a balance between accuracy and efficiency.

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: What exactly is $\hat{r}$ in lines 223-224 and 234-235? A location embedding? 

Have the authors tried llama 3 for the embeddings? Or maybe CLIP for better visual alignment?
soundness: 2
strengths: **(S1)**: The problem tackled by this paper is underexplored in most ML works. The proposed solution appears to be straightforward to implement.

**(S2)**: The method can handle zero-shot prediction for species not seen in training data. This is a valuable advantage as a new classifier doesn’t need to be re-trained for each new query species. 

**(S3)**: The new metric introduced accounts for spatial proximity which makes a more fine-grained analysis of performance possible.
summary: The authors introduce a reformulation of species-distribution modeling to allow species embeddings via LLMs and introduce a new metric that takes spatial proximity into account.
weaknesses: **(W1)**: No theoretical guarantees are given for the proposed metric. Since the authors are introducing a new metric that is posited to be more robust to variations in location, it would be helpful to establish the bounds of tolerance for variability as well as formal guarantees. Nothing as presented about the metric makes it specific to species-distribution modeling (SDM), and so the authors must be mathematically specific about where and how this metric should be used.

**(W2)**: Lack of clarity in problem formulation. Why is the region variable $r$ discretized (section 3.1)? It would be more useful if $r$ was allowed to continuously vary which would more realistically model geospatial coordinates (i.e. latitude, longitude). One can still model the problem as $p(r|y,l_r, e_r)$ with continuously varying $r$. By discretizing $r$ as a uniform grid across the globe, there might be lots of values of $r$ for which the probability mass would be 0, and in regions of high probability, the discretization might hurt the granularity of the approximation. Modeling the problem with continuously varying $r$ seems more natural and yet is not considered in this work.

**(W3)**: Lack of presentation clarity. Portions of the paper are not well explained. The loss section particularly (3.4) is not explicitly tied to the variables and terms introduced in the problem definition section. As far as I can see, the task description for “species range prediction” is not explicitly defined (what is the input, the output etc.).

**(W4)**: Unclear experimental evaluation. It seems that the authors have compiled their own training and testing datasets. Why weren’t comparisons made directly on the datasets used for SINR (or other baselines) for example? Also, it looks like SINR has better geo-feature regression which suggests better geospatial representation learning, which is not discussed in detail in the experimental evaluation.

Overall, I think that the weaknesses of the paper outweigh its benefits and so I lean towards rejection. The proposed method is not especially novel (eg: the species embeddings are extracted from a pre-trained Llama 2 via natural language queries) and while a new metric is introduced, it is not thoroughly justified. I believe the authors could improve their work by being more precise in their presentation.

---

confidence: 4
contribution: 4
initial_rating: 3
presentation: 4
questions: Since my questions closely parallel weaknesses, they have been stated in the "Weaknesses" section.
soundness: 2
strengths: * The paper is generally clearly written and the figures are well-designed. 
* The incorporation of language into species distribution modeling is novel to the best of my knowledge. 
* The proposed method significantly outperforms other methods (at least on the new metric introduced in this paper - see comments on this below).
* Comprehensive evaluation of different combinations of prior methods and losses (including the $ASL$ and $RAL$ multi-label losses, which had not been previously studied for this problem). 
* Separating out evaluation on "unseen" and "rare" species is and interesting and valuable idea. 
* The paper introduces a new proximity-aware metric which seems interesting and intuitively reasonable. 
* The hyperparameters are clearly specified and there is enough information to reproduce the results of the paper. 
* Species distribution modeling is an interesting and important problem.
summary: This paper considers the important problem of large-scale species distribution modeling. Different from prior works, this paper uses an LLM to encode taxonomic information for each species. The proposed model is also capable of zero-shot prediction, which is a novel capability. This work also introduces a novel evaluation metric that is meant to be more sensitive to spatial proximity when considering prediction errors. The paper also curates new training and evaluation datasets. Generally, the paper is interesting and valuable, but there are some experimental methodology issues that should be resolved.
weaknesses: Missing pieces:
* Computational resources seem to be a major confounder (in terms of model sizes and potentially training times). Can you explicitly compare the sizes  and # training steps for different models to contextualize the performance differences we observe between the proposed model and prior works? 
* We know prompt engineering is very important for language models. Were there any attempts at prompt engineering? What if we just use the species name without any other taxonomic information? There doesn't seem to be any evidence that the proposed text encoding is a better choice than other simple alternatives. 
* A new metric was proposed, but the problem it is meant to solve was not demonstrated to exist. (Though I agree that intuitively it seems reasonable that the problem might exist.) Finding a few nice examples would fix this. It was also not demonstrated that the new metric is providing different information than existing ones in practice. For instance, does the new metric rank models differently than previous commonly used metrics? (Showing pre-existing metrics e.g. mAP alongside the new metrics would help.)
* This work only makes one direct comparison to other methods using existing metrics: the Geo Feature results in Table 2. In that case, existing models outperform the proposed model. This raises the concern that the proposed model might underperform prior models on other pre-existing tasks as well. Since the newly proposed metric is not yet accepted as the gold standard, it is important to demonstrate how the proposed model performs on prior metrics as well. 
* How does the proposed metric handle the case of "vagrants" - species that are spotted far from their normal ranges? It seems like this is a problematic case (especially for birds!) that is not discussed. 
(Apologies if any of these are not missing and I simply missed it, please point it out if so!)

Slightly misleading or unclear bits:
* L42-50 seems to indicate that the framework of Mac Aodha et al. cannot use "species-specific rich information". Is this true? What would stop their model from using fixed species embeddings (e.g. those based on taxonomic string embeddings) and learning the location encoder? Or using some sort of regularization to encourage learned species embeddings to have relative distances similar to those based on taxonomic string embeddings? It's not clear to me that the "reformulation" of the species distribution problem is necessary to achieve these goals, as is suggested by the paper. 
* L184: I'm not quite sure I'd characterize the approach in this paper as "natural language" - the paper focuses on a highly constrained and artificial sort of language (taxonomic strings). 
* "To avoid bias, we excluded species with less than 100 observations." What kind of bias are we talking about, and how does this avoid it? 

Minor comments (no need to respond):
* The citation for iNaturalist (the platform) should probably be www.inaturalist.org instead of Van Horn et al. (which is a dataset derived from iNaturalist data) . 
* L102: Strange capitalization for MacKenzie et al. 
* L105: "With no expert knowledge" - I'd clarify that no expert knowledge is required to capture the raw data, but it is required to identify the species. 
* Why not include common names? Might be closer to the distribution the language model was trained on.

---

confidence: 3
contribution: 1
initial_rating: 3
presentation: 2
questions: As follows from the weaknesses:
1) What is the justification for assuming that taxonomically similar species tend to have similar geographic ranges, given that speciation is often promoted by geographic isolation?
2) How does using full taxonomic names facilitate zero-shot predictions for unseen species, especially considering that species names may not inherently convey information about geographic range?
3) Are there benefits of using the full taxonomic name over the binary name, or of using the binary name compared to a one-hot encoding? Additionally, how does the conclusion that "Spherical Harmonics are important for SDM" follow from the results in Table 5?
soundness: 2
strengths: 1) The paper is written and structured.
2) The usage of an adapted version of the chamfer distance is interesting, since eveluation metrics are typically a bottleneck in developing SDMs.
3) The results show that their method is competitive against two competing methods.
summary: The authors tackle the problem of predicting the geographic range of a biological species given a set of opportunistic observations. They do so by leveraging a model the uses spherical harmonics to better fit functions defined on the Earth’s surface. Their main novelty is to encode the species with their full taxonomic name, thus allowing species that share part of their life tree to be encoded closer together.
weaknesses: 1) I was a bit surprised by their choice of main contribution: the use of full taxonomic names based on language models. This seems to assume that taxonomically similar species tend to have similar ranges, but I could not find a justification for this. I’m not an ecologist, but this seems to go counter the idea that speciation is promoted by the geographic isolation of populations.
2) Although I do understand that such an approach may be useful for obtaining range maps of higher-level taxa based on training only with species-level data, I could not grasp the rationale of why this would help performing zero-shot with unseen species. Why would the species name itself carry any information about its range?
3) The experimental setup does not manage to show the value of the different elements in the method: is there an advantage to using the full taxonomic name vs the binary name? Does the binary name provide better results that a one-hot encoding of the species? Even the last conclusion, “Spherical Harmonics are important for SDM”, is not really based on the provided results, since it is not clear why the models compared in Table 5 are actually comparable, or that they differ only in the use of spherical harmonics.

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 3
questions: [Q1] Line 35 states that “species belonging to the same genera tend to be found in similar locations”. This seems like a big assumption, as I would have assumed that the most common reason historically that new species were formed is because populations became geographically isolated from one another. It makes sense that they would have an affinity for the same habitats, but it is not clear if they would be found in the same location.  Is there evidence to support this claim e.g., from the literature or your training data? 

[Q2] How is the zero-shot evaluation performed? Is the full taxonomic name for the unseen species just provided to the model at inference time? 

[Q3] Can you provide statistics of how many of the species in the test sets are from the same genus as those in the train set. Also, report on average how many species are in each genus in the train set.
soundness: 2
strengths: The main strengths are:  
[S1] This is the first SDM approach, that this reviewer is aware of, that incorporates language. However, only a restrictive subset of all possible input text is used, i.e., a string encoding the taxonomic hierarchy for the species of interest      
[S2] Unlike existing methods that typically predict for a single location of interest, the proposed model makes a prediction for all locations simultaneously for one species at a time     
[S3] In addition to the standard losses used in this problem space, the authors also evaluate different training losses from the single positive multi label literature    
[S4] The quantitative results in Table 2 demonstrate a non-trivial improvement over existing state of the art methods
summary: This work addresses the problem of species distribution modelling, i.e., estimating the spatial distribution of a species from a limited set of point observations. A new model, LD-SDM, is proposed that directly predicts a probability for a discretized set of locations indicating if a species is present there or not. The model is conditioned on a set of environmental features and uses a pre-trained language model to encode the identity of the species of interest. Results are presented on a new dataset adapted from public crowdsourced data where the proposed model outperforms existing work.
weaknesses: The main weaknesses are:  
[W1] The central assumption that “species belonging to the same genera tend to be found in similar locations” is a big assumption and not backed up with evidence. See Q1 below for a request for more information about this.   
[W2] It is not clear if it is necessary to use an LLM to encode the taxonomy/hierarchy. A simpler approach could be to train a per-species embedding as in SINR. For unseen species at inference time, one could just use the average embedding of all species from the train set that belong to the same genus. At the very least, an ablation akin to this would allow the readers to understand if the LLM adds value.   
[W3] A new training and evaluation dataset is proposed. However, it is not clear why the data from Cole et al. ICML 2023 is not used. Line 72 states that hand crafted range maps are “not widely available”. However, this is exactly the type of evaluation data used in Cole et al.   
[W4] The proposed evaluation metric would appear to be very sensitive to missing presences in the evaluation data. As indicated in the paper, species observation data is highly spatially biassed. Thus, only evaluating models on locations where data has been observed runs the risk of incorrectly penalising correct predictions from a model when there may in fact actually be presences for the species in that location, but they are just not in the evaluation dataset. The advantage of the evaluation protocol used in Cole et al. is that the evaluation data is spatially dense. 

Addition minor comments, that do not require a response in the rebuttal:   
L25 “remote-sensing task” -> SDM is not a remote sensing task    
L27 goal is to “produce large-scale range maps” -> the SDM problem represents a wider set of tasks than just range estimation, e.g., abundance     
L51 It is not clear what this sentence means, i.e., there are models that have been trained with location only, location and environmental features, and location and remote sensing images. It should be rephrased    
L71 “presence-only” data is introduced here for the first time and not defined  
L142 “n” is used for the dimensionality of the environmental features and “N” is used for the number of training observations  
L344 should make it clear in the text that most of the data in GBIF is actually from iNaturalist or eBird  
L345  Van Horn et al. 2018 is the incorrect reference, i.e., it is the iNaturalist 2017 image classification dataset, not the iNaturalist website    
L361 this spatial binning will not result in an even distribution around the globe, it would be better to use H3 cells so the poles are not overrepresented. Perhaps that is what is being done, but it is not clear?   
L516 should note in the limitations the big assumption that taxonomic is a good proxy for range similarities  
Fig 9 the spatial distribution of unseen species is very different from the others. This should be commented on

Meta-review: This submission received four reviews, with reviewers expressing concerns about the validity of several assumptions and the absence of key comparisons. As the authors did not provide a rebuttal to address these points, the reviewers maintained their initial negative scores. We encourage the authors to carefully incorporate the feedback received to strengthen their work for future submissions.

Decision: Reject