data("nuvel1")
# North America relative to Pacific plate:
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
data("san_andreas")
san_andreas$azi.PoR <- PoR_shmax(san_andreas, PoR)
# convert back to geo CRS
PoR2Geo_azimuth(san_andreas, PoR) |> head()
Run the code above in your browser using DataLab