3D Surface Generator

This code will create a 3D surface based on the input of only 3 points.

library(gstat) # for the IDW function (idw)
library(spatial)
library(fields)
library(raster)
library(marmap)
library(lattice)

The main functions we’ll need for interpolating points are below:

triangleFill2D <- function(df, num.splits){
  
  # First set the plotting window to display multiple plots
  par(mfrow = c(2, 3))
  # Loop through, keeping the points in diff dfs so that points aren't compared twice
  # This loop works
  first.df <- as.data.frame(df[, 1:2]) 
  second.df <- as.data.frame(df[, 1:2])
  
  for(j in 1:num.splits){
    
    num1 <- nrow(first.df) 
    num2 <- nrow(second.df)
    
    # Create a blank df
    third.df <- matrix(0, nrow = 0, ncol = 3, byrow = TRUE)
    third.df <- as.data.frame(third.df)
    colnames(third.df) <- c("x", "y") # t is the triangle number
    
    # Loop
    for(n in 1:num1){
      for(i in 1:num2){
        newx <- mean(c(first.df$x[n], second.df$x[i]))
        newy <- mean(c(first.df$y[n], second.df$y[i]))
       
        # make new entry to paste on end
        new.row <- c(newx, newy) 
        third.df <- rbind(third.df, new.row)
      }
    }
    
    # Remove duplicates from the third list
    colnames(third.df) <- c("x", "y")
    third.df <- unique(third.df[ , 1:2])
    third.df <- as.data.frame(third.df)
    
    # Create the new first list
    first.df <- rbind(first.df, second.df)
    first.df <- unique(first.df[ , 1:2])
    first.df <- as.data.frame(first.df)
    plot(first.df[ , 1], first.df[ , 2], xlim = c(0,5), ylim = c(0, 5), 
         pch = ".", main = paste0("Split #", (j-1)))
    
    # Create the new second list
    second.df <- third.df
    
    # Now make sure the second.df doesn't have any of the same points at the first.df
    # So combine first and second, then remove duplicates
    combo.df <- as.data.frame(rbind(first.df, second.df))
    combo.df <- as.data.frame(unique(combo.df[ , 1:2]))
    combo.df2 <- as.data.frame((rbind(first.df, combo.df)))
    combo.df2 <- combo.df2[!(duplicated(combo.df2)|duplicated(combo.df2, fromLast=TRUE)), ] # eliminate duplicates
    second.df <- combo.df2
    
  }
  
  final.df <- rbind(first.df, second.df)
  final.df <- as.data.frame(unique(final.df[, 1:2]))
  print(nrow(final.df))
  
  # Add the original 3 z coordinates
  final.df$z <- 0
  final.df$z[1:nrow(df)] <- df[,3]
  
  return(final.df)
}

surfaceCreator <- function(df, num_neighbors, stdev, runs){
  for(j in 1:runs){
    for(i in 1:nrow(df)){
      # Get the row
      kk <- df[i, ]
      
      # Calculate the distance between the points and add as a dist column
      df$diffx <- df$x - kk$x
      df$diffy <- df$y - kk$y
      df$dist <- sqrt((df$diffx)^2 + (df$diffy)^2)
      
      # Get the 6 points that are closest to it
      # order based on x, then order that based on y, then take top 6. That should be the points?
      ordered_by_dist <- df[order(df$dist), ]
      nearest6 <- ordered_by_dist[2:(num_neighbors+1), ] # the nearest point is the point itself
      
      # plot(x = kk$x, y = kk$y)
      # points(x = nearest6$x, y = nearest6$y)
      
      
      # Then take mean of the 6 points and the current point and use that as the mean and then have some stdev
      avg <- mean(c(nearest6$z, kk$z))
      newz <- new.row <- abs(rnorm(n = 1, mean = avg, sd = stdev)) # Here, I'm allowing the Z value to vary
      
      # Assign that value back to the df
      df$z[i] <- newz
    }
  }
  return(df)
}

Below, I create the three points that will be used to generate the surface. Each point has an X, Y and Z coordinate.

values <- c(1,1,3,
            2,4,7,
            3,2,1)

tri <- matrix(values, nrow = 3, ncol = 3, byrow = TRUE)
tr <- as.data.frame(tri)
colnames(tr) <- c("x", "y", "z") 

Take a look at the points:

Here, I decide how many times I want to split the points and then run the function to do just that.

# How many splits?
num.splits <- 5

# Run the function to carry out all the splits 
filled.triangle <- triangleFill2D(df = tr, num.splits = num.splits) 

## [1] 561

plot(filled.triangle[ , 1], filled.triangle[ , 2], xlim = c(0,5), ylim = c(0, 5), pch = ".",
     main = paste0("Split #", num.splits))

The elevation of all the points except the first three are set at zero. The first three act to seed the surface generator. The generator creates random changes from the seed elevation to assign values to the points we just created above. The roughness of the surface is controlled by setting the standard deviation. Higher standard deviations create rougher surfaces.

By running the surface generator function several times on the same data, we grow the landscape. Here, I run the function 50 times by inputting that in the “runs” variable. The more runs, the more dramatic the landscape.

# The plot window has been changed, so let's change it back
 par(mfrow = c(1, 1))

# Now run the surface generator
trial <- surfaceCreator(df = filled.triangle, num_neighbors = 10, stdev = 7, runs = 50)
plot(trial[ , 1], trial[ , 2], xlim = c(0,5), ylim = c(0, 5), pch = ".", col = trial[ , 3])

Now, we’ll start making this data spatial by creating coordinates and assigning a coordinate reference system

# CRS for the US in WGS 1984 (world geographic 1984)
wgs1984.proj <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

# Designate x and y columns as coordinates
coordinates(trial) = ~x + y

From there, I create a grid that will have values assigned from the triangle. The grid can then be made into a raster and that raster will eventually be our 3D surface

# Get min and max coord values
minx <- min(trial$x)
maxx <- max(trial$x)
miny <- min(trial$y)
maxy <- max(trial$y)

# Create a data.frame from all combinations of the supplied vectors or factors.
grd <- expand.grid(x = seq(from = minx, to = maxx, by = 0.005), y = seq(from = miny, to = maxy, by = 0.005))  # expand points to grid
coordinates(grd) <- ~x + y
gridded(grd) <- TRUE

Now, let’s look at how the triangle’s points overlap with the grid

# Plot the points and grid
plot(grd, cex = 1.5, col = "grey")
points(trial, pch = 1, col = "red", cex = 1)

Now, we’ll use inverse distance weighting to fill the grid with elevation (z) values interpolated from the triangle’s data.

# Apply the IDW
idw <- idw(formula = z ~ 1, locations = trial, newdata = grd)  # apply idw model for the data

## [inverse distance weighted interpolation]

idw.ras <- raster(idw)

# Rasterize and crop the IDW (might not be needed)
proj4string(idw.ras) <- wgs1984.proj
plot(idw.ras)

Might as well try a thin plate spline to see if things look much different.

m <- Tps(coordinates(trial), trial$z)
tps <- interpolate(idw.ras, m)
plot(tps)

The IDW doesn’t look so good. We could work to make it look a little better, but why waste our time? Let’s run with the TPS, it looks great!

Now, we’ll make a cool plot to really see the surface we’ve generated.

# Create hillshade
slope <- terrain(tps, opt = 'slope')
aspect <- terrain(tps, opt = 'aspect')
hill <- hillShade(slope, aspect, angle = 45, direction = 0)

plot(hill, col=grey(0:100/100), legend=FALSE)
plot(tps, col=rev(topo.colors(50, alpha=0.35)), add=TRUE)

Look at a contour map:

contour(tps,col = "blue")

Now, look at it from the side. By playing with the aspect, you can rotate the image to get the view you like.

bmap <- as.bathy(tps)
wireframe(unclass(bmap), shade = TRUE, aspect = c(2/1, 0.7), ylab = "y", xlab = "x", zlab = "z")

That’s it! You can keep running the code and each time you’ll end up with different surfaces. Enjoy!