Generating Data

This code generates data for flight paths at 3 separate “levels” or heights, with random “relative humidity”" data sampled at each level.

suppressPackageStartupMessages(library(tidyverse)) # dplyr and other tools

# make fake data
t <- seq(0, 10, by=0.01)
x <- t*cos(t)
y <- t*sin(t)
z1 <- rnorm(n = 1001, mean = 5, sd = .5)
rh1 <- rnorm(n=1001, mean=34, sd=1)
z2 <- rnorm(n = 1001, mean = 10, sd = .5)
rh2 <- rnorm(n=1001, mean=33.9, sd=1)
z3 <- rnorm(n = 1001, mean = 15, sd = .5)
rh3 <- rnorm(n=1001, mean=33.8, sd=1)
dfw1 <- data.frame(x, y, z=z1, rh=rh1, lev=1)
dfw2 <- data.frame(x, y, z=z2, rh=rh2, lev=2)
dfw3 <- data.frame(x, y, z=z3, rh=rh3, lev=3)
dfw <- bind_rows(dfw1, dfw2, dfw3)

Plotting in 3d

Turns out it’s not very difficult…a quick google search turned up a help page, which then helped me create the following. I show what the 2d plot looks like, then the 3D plot.

library(plotly)
library(viridis)

# first this is what it looks like in 2d

# 2d plot
ggplot(data=dfw, aes(x=x,y=y, group=lev))+
  geom_path() + theme_bw() +
  labs(title = "A Spiral Flight Path in 2D")

# 3d plot
plot_ly(dfw, x=~x, y=~y, z=~z) %>%
  add_markers(marker=list(color=~rh, 
                          colorscale=viridis(n = 12),
                          colorbar=list(title='Relative. \n Humidity'))) %>% 
  layout(title="Simulated Spiral Flight Path at 3 Heights",
         scene = list(zaxis=list(range=c(0,20))))

Linear Interpolation

The following uses akima and fields for the interpolation.

library(akima) # interpolation
library(fields) # plotting/interp

# pull out variables for one level
x1 <- dfw1$x 
y1 <- dfw1$y
z1 <- dfw1$rh

# using the interp linear method:
interp1 <- interp(x1, y1, z=dfw1$rh, linear = T, nx = 90, ny = 90) # output grid 90x90

glimpse(interp1) # quick check of data

## List of 3
##  $ x: num [1:90] -9.48 -9.3 -9.12 -8.94 -8.77 ...
##  $ y: num [1:90] -5.44 -5.29 -5.14 -4.99 -4.84 ...
##  $ z: num [1:90, 1:90] NA NA NA NA NA NA NA NA NA NA ...

We can then view our newly interpolated layer using the raster package, and convert the data to a raster, and extract the XY points back out from the interpolated data.

library(raster) # for converting to raster

# Raster plot
image.plot(interp1, nlevel = 20, col = viridis(n=20))
points(x1, y1, bg="white", col="gray80",pch=21, cex=0.5) #add the original points

# Contour plot: Nice, but can't overlay points/layers
filled.contour(interp1$x,interp1$y,interp1$z, color.palette = viridis)

# Convert Data to a Raster (to extract XY)
interp1_raster <- raster(interp1) # make a raster

# Extract Data to Dataframe from Raster
rh.linear <- as.data.frame(rasterToPoints(interp1_raster)) %>% 
  rename(rh=layer) %>% # rename to the RH variable
  mutate(height=mean(dfw1$z))

Finally, we can use plotly again to visualize!

# Plotly Plot
plot_ly(rh.linear, x=~x, y=~y, z=~height, 
        marker=list(color=~rh, colorscale=viridis(n = 12), show.scale=TRUE)) %>%
  add_markers() 

Since this is based on a partial spiral, linear interpolation is pretty messy, and it looks kind of weird, but it works. Just for fun, zoom way in on the layer until you can see the individual points that make up the grid. You can make a grid more or less dense depending on what you need, or the surface density you are interested in using.

Let’s see how the spline interpolation compares to this.

Spline Interpolation

Ok, so this uses a different method, but same idea, same data. Notice there are specific arguments for spatial data (lon.lat, miles). Could and probably should convert these XY points to lon/lat, or spatial data, but for now, let’s move forward (and ignore the warnings).

# Here we spline each layer or height separately
spline1 <- Tps(data.frame(x1,y1), z1, miles = F, lon.lat = F)

## Warning: 
## Grid searches over lambda (nugget and sill variances) with  minima at the endpoints: 
##   (GCV) Generalized Cross-Validation 
##    minimum at  right endpoint  lambda  =  5059.179 (eff. df= 3.00103 )

spline2 <- Tps(data.frame(dfw2$x,dfw2$y), dfw2$rh, miles = F, lon.lat = F)

## Warning: 
## Grid searches over lambda (nugget and sill variances) with  minima at the endpoints: 
##   (GCV) Generalized Cross-Validation 
##    minimum at  right endpoint  lambda  =  5059.179 (eff. df= 3.00103 )

spline3 <- Tps(data.frame(dfw3$x,dfw3$y), dfw3$rh, miles = F, lon.lat = F)

## Warning: 
## Grid searches over lambda (nugget and sill variances) with  minima at the endpoints: 
##   (GCV) Generalized Cross-Validation 
##    minimum at  right endpoint  lambda  =  5059.179 (eff. df= 3.00103 )

# get some warnings, but that's ok.

# make grid 90 x 90 (can make smaller if you want fewer points)
fullRH.grid1 <- predictSurface(spline1, nx = 90, ny = 90) 
fullRH.grid2 <- predictSurface(spline2, nx = 90, ny = 90) 
fullRH.grid3 <- predictSurface(spline3, nx = 90, ny = 90) 

Now we can view our grid and convert to a raster.

# raster plot
image.plot(fullRH.grid1, nlevel = 20, col = viridis(n=20), legend.args = list(text="RH", line=0.5), main=paste0("Relative Humidity at: ", round(mean(dfw1$z), digits = 2), "m"))
points(x1, y1, cex=.5, pch=16) # add points from orig dataset

# Contour plot: notice grid and contours
filled.contour(fullRH.grid1, color.palette = viridis) 

# Convert to Raster (to extract XY)
rh1.sp <- raster(fullRH.grid1)
rh1.spline <- as.data.frame(rasterToPoints(rh1.sp)) %>% 
  rename(rh=layer) %>% # rename to the RH variable
  mutate(height=mean(dfw1$z)) # take mean of the height for plotting purposes

rh2.sp <- raster(fullRH.grid2)
rh2.spline <- as.data.frame(rasterToPoints(rh2.sp)) %>% 
  rename(rh=layer) %>% # rename to the RH variable
  mutate(height=mean(dfw2$z))

rh3.sp <- raster(fullRH.grid3)
rh3.spline <- as.data.frame(rasterToPoints(rh3.sp)) %>% 
  rename(rh=layer) %>% # rename to the RH variable
  mutate(height=mean(dfw3$z))

Now let’s view this spline interpolation using plotly.

# set axis 
axx <- list(
  #nticks = 5,
  title="X",
  autorange = "reversed"
  #range = c(-121.525, -121.508)
)

axy <- list(
  #nticks = 5,
  title="Y",
  autorange = "reversed"
  #range = c(38.191, 38.194)
)

axz <- list(
  nticks = 5,
  title="Height"
  #range = c(5.5,6.5)
)

# put it together
plot_ly(rh1.spline, x=~x, y=~y, z=~height, 
        marker=list(color=~rh, colorscale=viridis(n = 12), 
                    showscale=TRUE)) %>% 
  layout(scene = list(xaxis=axx,
                      yaxis=axy,
                      zaxis=axz)) %>%
  add_markers()

This definitely appears to be a smoother interpolation of our simulated data.