如何在R中基于WorldClim数据计算热指数并绘制全球温暖指数？

Hey there! Let's work through your questions about calculating and mapping the Yim & Kira Warmth Index using WorldClim data in R. I'll break this down into clear, actionable steps.

First, let's adapt your existing function to work with raster data. Your original function was built for tibble rows, but we need a version that can handle the 12 monthly temperature values for each raster pixel.

Step 1: Rewrite the Warmth Index Function for Rasters #

Instead of passing separate warm and temp arguments, we'll create a function that takes a vector of 12 monthly temperatures (one for each month) and returns the Warmth Index for that location:

warmth_index_raster <- function(temp_vector) {
  # Filter months where temperature > 5°C
  warm_temps <- temp_vector[temp_vector > 5]
  
  # Handle cases with no warm months (return 0)
  if (length(warm_temps) == 0) {
    return(0)
  }
  
  # Apply Yim & Kira formula: sum(warm temps) - (5 * number of warm months)
  sum(warm_temps) - 5 * length(warm_temps)
}

Step 2: Load WorldClim Monthly Temperature Data #

We'll use the terra package (the modern replacement for raster) to work with WorldClim's TIFF files. Assuming you've downloaded the 12 monthly temperature layers:

library(terra)

# Load all 12 monthly temperature TIFFs into a SpatRaster stack
temp_stack <- rast(list.files("path/to/worldclim/temp", pattern = "\\.tif$", full.names = TRUE))

Make sure the files are ordered from January to December!

Step 3: Compute Global Warmth Index #

Use terra::app() (equivalent to raster::calc()) to apply our function to every pixel in the stack. This will generate a single-layer raster with the Warmth Index for every global location:

wi_global <- app(temp_stack, warmth_index_raster)

Now let's create a graded global map of the Warmth Index. We'll use tidyterra to integrate terra with ggplot2 for clean, customizable plots.

Step 1: Basic Continuous Map #

For a smooth color gradient that shows fine-scale variation:

library(tidyterra)
library(ggplot2)

ggplot() +
  geom_spatraster(data = wi_global) +
  scale_fill_viridis_c(
    name = "Warmth Index",
    option = "plasma", # Colorblind-friendly palette
    na.value = "white" # Mask ocean/NA values
  ) +
  theme_minimal() +
  labs(
    title = "Global Warmth Index (Yim & Kira)",
    x = "Longitude",
    y = "Latitude"
  )

Step 2: Graded (Classified) Map #

To create discrete ecological classes (e.g., for comparing habitat zones), use terra::classify() to bin the Warmth Index values:

# Define class breaks and labels (adjust these to fit your research needs)
wi_classes <- classify(
  wi_global,
  breaks = c(-Inf, 0, 50, 100, 200, 300, Inf),
  labels = c("No Growing Season", "Low", "Moderate", "High", "Very High", "Extreme")
)

# Plot classified raster
ggplot() +
  geom_spatraster(data = wi_classes) +
  scale_fill_viridis_d(
    name = "Warmth Index Class",
    option = "plasma",
    na.value = "white"
  ) +
  theme_minimal() +
  labs(
    title = "Global Warmth Index Classes",
    x = "Longitude",
    y = "Latitude"
  )

Bonus: Apply to Your Plant Population Tibble #

If you want to calculate the Warmth Index for your existing plant location tibble, you can reuse the same function with dplyr::rowwise():

library(dplyr)

plant_data <- raw_data %>%
  rowwise() %>%
  mutate(
    warmth_index = warmth_index_raster(
      c(temp_1, temp_2, temp_3, temp_4, temp_5, temp_6,
        temp_7, temp_8, temp_9, temp_10, temp_11, temp_12)
    )
  ) %>%
  ungroup()

内容的提问来源于stack exchange，提问作者Harry