prepare a data from a data.table ready for printing as a 3d grid

dt2xyz.R

# Prepare wide data for 3D printing
# Expects to be passed 3 columns: x,y,z
# x and y will be binned
# a function of z will be used to prepare the z coordinate

library(data.table)
library(ggplot2)
library(r2stl)
library(Hmisc)
library(assertthat)


# Prepare bins
into_bins <- function(x, bins=30) {
    xmin <- min(x, na.rm=TRUE)
    xmax <- max(x, na.rm=TRUE)
    cuts <- seq(xmin, xmax, length.out=bins+1)
    res <- Hmisc::cut2(x, cuts)
    return(as.integer(res))
}

# Prepare matrix based on bins (assumes square)
bins2matrix <- function(dd, bins, fun=median, show.plot=FALSE){

    # now create grid
    grid <- setDT(expand.grid(x=seq(bins), y=seq(bins)))
    # collapse by median
    dd <- dd[, .(z=fun(z, na.rm=TRUE)), by=,.(x,y)]
    grid <- dd[grid, on=c(x='x', y='y')]
    # replace missing with zero
    grid[is.na(z), z := 0]
    if (show.plot) {
        # subtract one else plots with a rounding error?
        gg <- ggplot(grid, aes(x=x-1, y=y-1, z=z)) + stat_summary_2d(fun = function(x) x)
        print(gg)
    }
    grid <- as.matrix(dcast.data.table(grid, x ~ y))
    # get rid of rownames
    grid <- grid[ , 2:bins+1]
    return(grid)
}


# via flowingdata.com
# https://flowingdata.com/2018/05/07/3-d-printing-how-to-prepare-the-data-in-r/
flattenSurface <- function(z, n_per_pt = 4) {
    # Generate squares for each cell value.
    newz <- matrix(0, nrow=dim(z)[1]*n_per_pt, ncol=dim(z)[2]*n_per_pt)
    for (i in 1:dim(z)[1]) {
        for (j in 1:dim(z)[2]) {
            curr_val <- z[i, j]

            curr_x <- ( (i-1) * n_per_pt + 1 ) : ( (i-1) * n_per_pt + n_per_pt )
            curr_y <- ( (j-1) * n_per_pt + 1 ) : ( (j-1) * n_per_pt + n_per_pt )

            coords <- expand.grid(x=curr_x, y=curr_y)
            newz[coords$x, coords$y] <- curr_val
        }
    }

    return(newz)
}

getClosedSurface <- function(z, border=1) {
    # Close surface with zeros around border.
    closedz <- matrix(nrow=dim(z)[1]+border*2, ncol=dim(z)[2]+border*2)
    closedz[1:border,] <- 0
    closedz[,1:border] <- 0
    closedz[(dim(closedz)[1]-border+1):dim(closedz)[1],] <- 0
    closedz[,(dim(closedz)[2]-border+1):dim(closedz)[2]] <- 0
    closedz[(border+1):(dim(closedz)[1]-border), (border+1):(dim(closedz)[2]-border)] <- z

    return(closedz)
}


# Wrapper function
dt2xyz <- function(dt,
                    x='x', y='y', z='z',
                    bins=30, fun=median, show.plot=FALSE, grid=4, border = NA) {
    # convert data.table of arbitrary length into binned matrix suitable for r2stl
    if (is.na(border)) {
        border = grid
    }
    # check data.table suitable vars and suitable size
    assert_that(is.data.table(dt))
    assert_that(all(c(x, y, z) %in% names(dt)))
    assert_that(bins > 1)
    assert_that(nrow(dt) > bins)

    dd <- dt[,.(xx=get(x),yy=get(y),zz=get(z)), with=TRUE]
    dd <- dd[, .(x=into_bins(xx, bins), y=into_bins(yy, bins), z=zz)]
    dd <- bins2matrix(dd, bins, show.plot=show.plot)
    dd <- flattenSurface(dd, n_per_pt = grid )
    dd <- getClosedSurface(dd, border = border)

    return(dd)

}


# Vignette / example
# Prepare dummy data
dt <- data.table(x=rnorm(1000), y=rnorm(1000), z=abs(rnorm(1000)))
grid <- dt2xyz(dt, show.plot=TRUE)
persp(grid, theta=30)