在R语言中如何基于经纬度求解椭圆的笛卡尔中心及长短轴？

Hey there! Let's work through how to calculate the Cartesian center, major axis, and minor axis lengths for each grouped ellipse from your geographic point data in R. Here's a step-by-step solution tailored to your needs:

1. Understand the Approach #

Since you're dealing with geographic points grouped into ellipses (each with ≥3 points), we can leverage covariance matrix eigenvalue decomposition to derive the ellipse parameters. This method works well for fitting ellipses to point clusters, as it captures the spread and orientation of the data. Note: If your data is in raw latitude/longitude (spherical coordinates), we'll first project it to a planar Cartesian system (like UTM) to avoid distortion.

2. Prepare Your Data #

First, let's structure your sample data into a usable format, and add a grouping column (adjust this to match your actual group assignments):

# Sample data with groups (adjust groups to match your actual clustering)
point_data <- data.frame(
  longitude = c(-118.8267, -115.9665, -117.2978, -117.2962, -117.1625),
  latitude = c(33.73430, 33.25514, 34.18589, 34.18449, 34.00642),
  location = paste0("location ", 1:5),
  group = c("cluster_A", "cluster_A", "cluster_B", "cluster_B", "cluster_B") # Ensure your real data has ≥3 points per group!
)

3. Project Geographic Coordinates to Cartesian (If Needed) #

Raw latitude/longitude are spherical, so we'll convert to UTM (a planar coordinate system) using the sf package. Skip this step if your data is already in Cartesian coordinates:

library(sf)
library(dplyr)

# Convert to sf object (WGS84 coordinate system)
sf_points <- st_as_sf(point_data, coords = c("longitude", "latitude"), crs = 4326)

# Auto-detect and convert to UTM (avoids distortion)
utm_crs <- st_crs(paste0("+proj=utm +zone=", 
                         cut(st_coordinates(sf_points)[,1], breaks = seq(-180, 180, 6)) %>% as.integer() + 30,
                         " +datum=WGS84"))
sf_utm <- st_transform(sf_points, crs = utm_crs)

# Extract Cartesian x/y coordinates
planar_data <- sf_utm %>%
  mutate(x = st_coordinates(.)[,1],
         y = st_coordinates(.)[,2]) %>%
  st_drop_geometry()

4. Calculate Ellipse Parameters per Group #

We'll create a custom function to compute the center, major axis, and minor axis for each valid group (≥3 points):

# Function to compute ellipse parameters from a point cluster
get_ellipse_params <- function(cluster_data) {
  # Calculate center coordinates
  center_x <- mean(cluster_data$x)
  center_y <- mean(cluster_data$y)
  
  # Compute covariance matrix of the cluster
  cov_matrix <- cov(cluster_data[, c("x", "y")])
  
  # Eigen decomposition to get axis lengths and orientation
  eigen_vals <- eigen(cov_matrix)$values
  eigen_vecs <- eigen(cov_matrix)$vectors
  
  # Axis lengths (2*sqrt(eigenvalue) gives full axis length for 1σ ellipse)
  major_axis <- 2 * sqrt(eigen_vals[1])
  minor_axis <- 2 * sqrt(eigen_vals[2])
  
  # Optional: Major axis angle (in degrees, from x-axis)
  major_axis_angle <- atan2(eigen_vecs[2, 1], eigen_vecs[1, 1]) * (180 / pi)
  
  # Return results as a tibble
  tibble(
    group = unique(cluster_data$group),
    center_x = center_x,
    center_y = center_y,
    major_axis_length = major_axis,
    minor_axis_length = minor_axis,
    major_axis_angle_deg = major_axis_angle
  )
}

# Apply function to each group (filter out groups with <3 points)
ellipse_results <- planar_data %>%
  group_by(group) %>%
  filter(n() >= 3) %>%
  do(get_ellipse_params(.)) %>%
  ungroup()

# View the results
print(ellipse_results)

5. Notes on Customization #

• Confidence Ellipses: If you want a 95% confidence ellipse instead of a 1σ ellipse, multiply the axis lengths by sqrt(qchisq(0.95, 2)) (~2.4477).
• Least-Squares Ellipse Fit: For a strict geometric least-squares fit (instead of distribution-based), use the conicfit package's EllipseFit() function. Install it with install.packages("conicfit") if needed.

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