spatial_rotation.R

Rotation of Arbitrary Spatial Data

library(data.table)
library(sp)
library(ggplot2)
library(broom)
library(gridExtra)

set.seed(1)

# Set-up Dummy Data (for example)
n <- 
data3D <- data.table(merge(data.frame(x = runif(n,min = 504007.5,max=504874.9),y = runif(n,min=6251663.6,max=6252643.0)),data.frame(z = seq(0,1,by=1))))
data3D[,avg_phie:=rnorm(1,mean = 0.30,sd = 0.02),by=c("x","y")]
data3D[,avg_swe:=rnorm(1,mean = 0.27,sd = 0.04),by=c("x","y")]
data3D[,avg_h:=runif(1,min=0,max=50)*(1-0.2)+0.2,by=c("x","y")]
data3D[,phie:=rnorm(.N,mean=avg_phie,sd=0.02),by=c("x","y")]
data3D[,swe:=rnorm(.N,mean=avg_swe,sd=0.04),by=c("x","y")]
data3D[z>avg_h,phie:=0,]
data3D <- data.frame(subset(data3D,select=c(x,y,z,phie,swe)))
coordinates(data3D) = ~x+y+z
# Set Projection so we can rotate from here!
#proj4string(data3D) <- "+proj=utm +zone=12 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"

rotate <- function(spobj, angle=0) {
  tmp <- spobj
  suppressWarnings(proj4string(tmp) <-  "+proj=aeqd +ellps=sphere +lat_0=90 +lon_0=0")
  tmp <- spTransform(tmp, CRS(paste0("+proj=aeqd +ellps=sphere +lat_0=90 +lon_0=", -angle)))
  return(tmp)
}

localize <- function(spobj){
  tmp <- spobj
  tmp@coords[,1] <- tmp@coords[,1]-min(tmp@coords[,1])
  tmp@coords[,2] <- tmp@coords[,2]-min(tmp@coords[,2])
  return(tmp)
}

# Compute tests:
dt_prj_0 <- localize(rotate(data3D,0))
dt_prj_45 <- localize(rotate(data3D,45))
dt_prj_90 <- localize(rotate(data3D,90))
dt_prj_135 <- localize(rotate(data3D,135))
dt_prj_180 <- localize(rotate(data3D,180))


# Display Results
grid.arrange(
  ggplot(tidy(data3D),aes(x=x,y=y,colour=swe))+geom_point()+ggtitle("Original CRS")+coord_fixed(),
  ggplot(tidy(dt_prj_0),aes(x=x,y=y,colour=swe))+geom_point()+ggtitle("Transform to Local CRS")+coord_fixed(),
  ggplot(tidy(dt_prj_45),aes(x=x,y=y,colour=swe))+geom_point()+ggtitle("Transform 45deg Local CRS")+coord_fixed(),
  ggplot(tidy(dt_prj_90),aes(x=x,y=y,colour=swe))+geom_point()+ggtitle("Transform 90deg Local CRS")+coord_fixed(),
  ggplot(tidy(dt_prj_135),aes(x=x,y=y,colour=swe))+geom_point()+ggtitle("Transform 135deg Local CRS")+coord_fixed(),
  ggplot(tidy(dt_prj_180),aes(x=x,y=y,colour=swe))+geom_point()+ggtitle("Transform 180deg Local CRS")+coord_fixed(),  
  ncol=3,top="Transformed Data Coordinates")

Sample Output
Allows for arbitrary rotation & zeroing of coordinate reference system:

Hypothetically, If this was done before 3D kriging, could use this to output directly to grid cells for numerical simulation...