README

This package offers a method of fitting a smooth path to a given set of noisy spherical data observed at known time. It implements a piecewise geodesic curve fitting method on the unit sphere based on a velocity-based penalization scheme. The proposed approach is implemented using the Riemannian block coordinate descent algorithm. To understand the algorithm, you can refer to the paper titled “Intrinsic spherical smoothing method based on generalized Bézier curves and sparsity inducing penalization” by Bak, K. Y., Shin, J. K., & Koo, J. Y. (2023) for the case of order 1. Additionally, this package includes various functions necessary for handling spherical data.

Installation

install.packages("spheresmooth")

Example: APW path fitting

A polar wander dataset presented in Kent and Irving (2010). The 17 Triassic/Jurassic cratonic poles from other major cratons are rotated into North American coordinates and combined with the 14 observations from North America. Our method is applied to these 31 observations ranging in age from 243 to 144 Ma (millions of years ago), which covers the late Triassic and Jurassic periods. The first column represents the time points, and the remaining two columns provides the observed spherical coordinates.

library(spheresmooth)
data(apw_spherical)

apw_cartesian = Spherical_to_Cartesian(apw_spherical[, 2:3])

Generate a knot sequence along the time domain with a specified dimension and a lambda sequence.

t = apw_spherical[, 1]
dimension = 15
initial_knots = knots_quantile(t, dimension = dimension)
lambda_seq = exp(seq(log(1e-07), log(1), length = 40))

Fit the penalized piecewise gedesic curve to the data and obtain the control points for the optimal complexity parameter.

fit = penalized_linear_spherical_spline(t = t, y = apw_cartesian,
                                        dimension = dimension,
                                        initial_knots = initial_knots,
                                        lambdas = lambda_seq)
# choose a curve that minimizes the BIC
best_index = which.min(fit$bic_list)
best_index
#> [1] 33
# obtained control points for the piecewise geodesic curve
fit[[best_index]]$control_points
#>             [,1]       [,2]      [,3]
#> [1,] -0.48039447 -0.2266981 0.8472480
#> [2,] -0.11927705  0.1300343 0.9843089
#> [3,] -0.06869155  0.1967739 0.9780396
#> [4,]  0.01640654  0.2026538 0.9791130
#> [5,]  0.03231127  0.3512434 0.9357265
#> [6,] -0.03372888  0.5085848 0.8603510
#> [7,] -0.08310953  0.5246722 0.8472378
#> [8,] -0.45295783  0.3987197 0.7974032

# install.packages("rworldmap")
# install.packages("ggplot2")
library(rworldmap)
#> Loading required package: sp
#> ### Welcome to rworldmap ###
#> For a short introduction type :   vignette('rworldmap')
library(ggplot2)

worldMap = getMap()
world.points = fortify(worldMap)
#> Warning: `fortify(<SpatialPolygonsDataFrame>)` was deprecated in ggplot2 3.4.4.
#> ℹ Please migrate to sf.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
#> Regions defined for each Polygons
world.points$region = world.points$id
world.df = world.points[, c("long","lat","group", "region")]
colnames(world.df)[1] = "longitude"
colnames(world.df)[2] = "latitude"

cp_best = Cartesian_to_Spherical(fit[[best_index]]$control_points)
cp_long_lat = cp_best * 180 / pi
cp_long_lat_df = data.frame(latitude = 90-cp_long_lat[, 1],
                            longitude = cp_long_lat[,2])

apw_spherical_df = data.frame(apw_spherical)
colnames(apw_spherical_df) = c("time", "latitude", "longitude")
apw_spherical_df$latitude = 90 - apw_spherical_df$latitude * 180 / pi
apw_spherical_df$longitude = apw_spherical_df$longitude * 180 / pi

worldmap = ggplot() +
  geom_polygon(data = world.df, aes(x = longitude, y = latitude, group = group),colour = "grey", fill = "antiquewhite") +
  scale_y_continuous(breaks = (-2:2) * 30) +
  scale_x_continuous(breaks = (-4:4) * 45) +
  geom_point(data = apw_spherical_df, aes(x = longitude, y = latitude), cex = 0.8) +
  geom_point(data = cp_long_lat_df, aes(x = longitude, y = latitude), shape = 23, col = 'blue', cex = 4) +
  geom_line(data = cp_long_lat_df[1:5,], aes(x = longitude, y = latitude), col = 'red', lwd = 1) +
  geom_line(data = cp_long_lat_df[5:8,], aes(x = longitude, y = latitude), col = 'red', lwd = 1) +
  coord_map("ortho", orientation=c(38, 120, 0))
worldmap