Converting quartic kernel heatmap into large polygon with R

You can perform your estimation by slightly modifying the kde2d function from the MASS package. To my knowledge there is currently no package in R that implements bivariate KDE estimation with a quartic (biweight) kernel for a bivariate case.

The univariate biweight kernel can be extended to a multivariate kernel in several ways and the simplest is to just use a product kernel, where you use the univariate kernel for each of your dimensions and then multiply the result. You can find the mathematical expression for the biweight product kernel here. When you incorporate this kernel into the kde2d density estimator from the MASS package it looks as follows

kde_biweight_kernel <- function(x,y, bw_x, bw_y, xrange, yrange){
  # This function is based on the kde2d function from 
  # the MASS package. The only difference is that the Gaussian
  # kernel is substituted with a biweight product kernel
  
  # product kernel:
  biweight_kernel <- function(u){
    mask = abs(u) > 1
    kernel_val = (15/16)*((1-u^2)^2)
    kernel_val[mask] = 0
    return(kernel_val)
  }
  
  lims = c(xrange, yrange)
  n = 500
  nx <- length(x)
  n <- rep(n, length.out = 2L)
  # get grid on which we want to estimate the density
  gx <- seq.int(lims[1L], lims[2L], length.out = n[1L])
  gy <- seq.int(lims[3L], lims[4L], length.out = n[2L])
  
  # inputs to kernel
  ax <- outer(gx, x, "-" )/bw_x
  ay <- outer(gy, y, "-" )/bw_y
  
  # evaluate and multiply kernel results along both axes
  res = tcrossprod(biweight_kernel(ax), biweight_kernel(ay))/(nx * bw_x * bw_y)
  return(list(x = gx, y = gy, z = res))
}

Using the kde_biweight_kernelfunction you can compute your desired density as follows

library(MASS)
library(birk)
library(kedd)
library(sf)
library(ks)


# load data
my.data <- read.csv("fake_points_oahu.csv", header = TRUE, stringsAsFactors = FALSE, na.strings = "NA")
# Import shapefile
st_layers("fake_points_oahu_utm.shp")
points_utm <- st_read(dsn = "fake_points_oahu_utm.shp", layer = 'fake_points_oahu_utm')

x = my.data$longitude
y = my.data$latitude

# determine bandwidth for biweight kernel along both axes
bw_x = h.amise(x, deriv.order = 0, kernel = "biweight")$h
bw_y = h.amise(y, deriv.order = 0, kernel = "biweight")$h

# get ranges in which you want to estimate density
xrange = c(min(my.data$longitude), max(my.data$longitude))
yrange = c(min(my.data$latitude),  max(my.data$latitude))

# get 2d density estimate with quartic (biweight) kernel
result = kde_biweight_kernel(x,y, bw_x, bw_y, xrange, yrange)

Note that the bandwidth is especially computed for the biweight kernel case. The resulting density object is a bit different from the ks::kde object. It does, for example, not have contour levels yet. We can get the contour levels by computing the quantiles with a slightly modified version of the kde2dQuantile function from the rmngbpackage

# get quantiles of interest:
kde2dQuantile <- function(d, X, Y, probs = .05) {
  xInd <- sapply(X, function(x) which.closest(d$x, x))
  yInd <- sapply(Y, function(x) which.closest(d$y, x))
  zValues <- d$z[cbind(xInd, yInd)]
  quantile(zValues, probs=probs)
}
# get quantiles
quantiles = kde2dQuantile(result, x, y, seq(0,1,by=0.001))

From your question I am not sure which quantile you are interested in so I just chose the 1% quantile. To be able to plot the data in the same way as in the question, we have to format the density result in the same way as objects from the kde class:

# to make the kde estimate compatible with the other density estimates
# from the ks package, the result can be converted to a named list.
# -> create ks::KDE object:
axes = matrix(c(result$x,result$y), ncol = 2)
colnames(axes) = c('longitude', 'latitude')

my.fhat_biweight = list('x' = axes,
                        'eval.points' = list(result$x, result$y),
                        'estimate' = result['z']$z,
                        'gridtype' = 'linear', 'gridded' = TRUE,
                        'binned' = TRUE, 'names' = c("longitude","latitude" ))

# add quantile to ks::KDE object
my.fhat_biweight$cont = quantiles

# change class (make sure ks package is loaded for this)
class(my.fhat_biweight) <- "kde"

Finally plot the biweight kernel density over the data

plot(my.data$longitude, my.data$latitude)
plot(my.fhat_biweight, lwd = 3, display = "filled.contour", cont = cont=c(96.5), add = TRUE)

this outputs: