---
output: github_document
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)
```

# VirtualPals 🐨: Virtual Species and Niche Modeling in R

[![R-CMD-check](https://github.com/your-username/VirtualPals/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/your-username/VirtualPals/actions/workflows/R-CMD-check.yaml)
[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)

**VirtualPals** is an R package for creating virtual species by defining their ecological niches in a multi-dimensional environmental space. It allows you to project these niches onto real-world geographic landscapes to simulate species distributions. This tool is ideal for researchers and students who need to test Species Distribution Model (SDM) algorithms, explore ecological theory, or generate plausible data for educational purposes.

The core of the package is its ability to define a species' fundamental niche as a fully-rotatable, N-dimensional ellipsoid using real environmental data.

---

## 🚀 Key Features

* **Flexible Niche Definition**: Define species niches as 2D or 3D ellipsoids with full control over the center (optimum), axes (tolerances), and rotation.
* **Geographic Projection**: Generate geographic suitability maps from environmental raster data based on your defined niche.
* **Multiple Output Types**: Create maps of binary (presence/absence), probabilistic, or distance-based suitability.
* **Data Filtering**: Easily extract points from any dataset that fall within a defined niche.
* **Rich Visualization**: Integrates with `plotly` for interactive 3D visualizations of niches and with base R plotting for clear geographic maps.
* **Bundled Data**: Comes with sample WorldClim bioclimatic data so you can run examples out of the box.

---

## 📦 Installation

You can install the development version of VirtualPals from GitHub:

```{r}
# install.packages("devtools")
devtools::install_github("paanwaris/virtualPals")
```

---

## 🗺️ Core Workflow: From Niche to Map

This example demonstrates the full workflow using the bioclimatic data bundled with the package. We will define a niche for a hypothetical species and map its potential global distribution.

### Step 1: Load Libraries

```{r}
library(virtualPals)
library(terra)
library(plotly)
library(rnaturalearth) # For plotting world boundaries
```

### Step 2: Load Packaged Environmental Data

The package includes sample bioclimatic raster data. We load them using `system.file()` which finds the files inside the installed package directory.

```{r}
# Find the path to the bundled data files
path_bio1 <- system.file("extdata", "Bio1.tif", package = "virtualPals")
path_bio12 <- system.file("extdata", "Bio12.tif", package = "virtualPals")
path_bio4 <- system.file("extdata", "Bio4.tif", package = "virtualPals")

# Load the rasters
bio1 <- rast(path_bio1)   # Annual Mean Temperature
bio12 <- rast(path_bio12) # Annual Precipitation
bio4 <- rast(path_bio4)   # Temperature Seasonality
```

### Step 3: Prepare Environmental Layers

We'll stack the rasters, name them, and mask them to only include terrestrial areas.

```{r}
# Combine layers into a single SpatRaster
bio_stack <- c(bio1, bio12, bio4)
names(bio_stack) <- c("AnnualTemp", "AnnualPrecip", "TempSeasonality")

# Get world boundaries and mask the environmental data
world_sf <- ne_countries(scale = "medium", returnclass = "sf")
bio_stack_terrestrial <- mask(bio_stack, vect(world_sf))
plot(bio_stack_terrestrial)
```

### Step 4: Define the Species Niche

Let's define a niche for a species that prefers temperate, moderately wet environments with noticeable seasons. We use a **named vector** for `center` so the function can reliably match niche dimensions to the raster layers.

```{r}
# Define the 3D niche using a rotated ellipsoid
species_niche <- build_ellipsoid(
  center = c(AnnualTemp = 300, AnnualPrecip = 1500, TempSeasonality = 2000),
  axes = c(100, 500, 2500), # Tolerances for Temp, Precip, and Seasonality
  angles = c(0, 0, pi/12)   # We'll add a slight rotation
)

bio_df <- as.data.frame(bio_stack_terrestrial, xy = TRUE, na.rm = TRUE)
df <- bio_df[sample(1:nrow(bio_df), size = 100000, replace = FALSE), ]

plot_ly(df, x = ~AnnualTemp, y = ~AnnualPrecip, z = ~TempSeasonality,
                     type = 'scatter3d', mode = 'markers',
                     marker = list(color = "lightgrey",
                                   size = 3),
                     name = "Points") %>%
  add_trace(data = species_niche$surface, x=~x, y=~y, z=~z,
            type="scatter3d", mode="lines",
            line=list(color="blue"),
            name="Ellipsoid surface", inherit = FALSE) %>%
  add_markers(x = species_niche$center[1], y = species_niche$center[2], z = species_niche$center[3],
              marker = list(color = 'red', size = 5),
              name = "Center") %>%
  layout(
    title = "Ellipsoid in 3D with Axes and Center",
    scene = list(
      xaxis = list(title = "X"),
      yaxis = list(title = "Y"),
      zaxis = list(title = "Z")
    ),
    legend = list(x = 0.05, y = 0.95)
  )
```

### Step 5: Generate the Species Distribution Map

Now, we project this abstract niche onto our global environmental map to find all suitable habitats.

```{r}
# Generate a binary presence/absence map
suitability_map <- generate_suitability(
  env_rasters = bio_stack_terrestrial,
  niche = species_niche,
  output_type = "binary"
)
```


### Step 6: Visualize the Geographic Distribution

Let's plot our final map and add country borders for context.

```{r}
plot(suitability_map,
     col = c("grey90", "darkgreen"),
     legend = FALSE,
     main = "Potential Global Distribution of Virtualis palsicus")

# Add world boundaries
lines(vect(world_sf), col = "black", lwd = 0.5)
```

### Step 7: Visualize the Niche in 3D Environmental Space

Finally, let's create an interactive 3D plot to see how the species' niche (the ellipsoid) fits within the world's available environments.

```{r}
# Convert rasters to a data frame, subsampling for performance
env_df <- as.data.frame(bio_stack_terrestrial, na.rm = TRUE)
env_df_sample <- env_df[sample(nrow(env_df), 50000), ]

# Find which of the sampled points fall inside our niche
points_inside <- points_in_ellipsoid(species_niche, env_df_sample)
points_inside <- as.data.frame(points_inside)

# Define your colors
my_colors <- c("grey", "darkgreen", "blue")

plot_ly(colors = my_colors) %>%
  # Add all available environmental points (the background)
  add_markers(data = env_df_sample,
              x = ~AnnualTemp, y = ~AnnualPrecip, z = ~TempSeasonality,
              marker = list(color = "grey", size = 2, opacity = 0.4),
              name = "Global Environments") %>%
  # Highlight the suitable points that fall inside the niche
  add_markers(data = points_inside,
              x = ~AnnualTemp, y = ~AnnualPrecip, z = ~TempSeasonality,
              marker = list(color = "darkgreen", size = 2.5),
              name = "Suitable Environments") %>%
  # Draw the niche boundary as a 3D mesh
  add_trace(data = species_niche$surface,
            x = ~x, y = ~y, z = ~z,
            type = "scatter3d", 
            mode = "lines",
            opacity = 0.15, 
            line = list(color = "blue", width = 2),
            name = "Niche Boundary") %>%
  layout(
    title = "Species Niche in 3D Environmental Space",
    scene = list(
      xaxis = list(title = "Annual Temp (°C)"),
      yaxis = list(title = "Annual Precip. (mm)"),
      zaxis = list(title = "Temp. Seasonality")
    )
  )

```


---

## 🛠️ Function Overview

* `build_ellipsoid()`: Creates a 2D or 3D ellipsoid object with specified center, axes, and rotation.
* `generate_suitability()`: Projects an `ellipsoid` niche onto `SpatRaster` environmental layers to create a suitability map.
* `points_in_ellipsoid()`: Filters a `data.frame` or `matrix` to return only the points that fall inside an `ellipsoid`.