library("TDA")

X = read.table("crater.txt")  # read data
plot(X[,1], X[,2], "p")

# grid
Xlim <- c(min(X[,1]), max(X[,1])); Ylim <- c(min(X[,2]), max(X[,2]))
byx <- (Xlim[2]-Xlim[1])/200; byy <- (Ylim[2]-Ylim[1])/200 
Xseq <- seq(Xlim[1], Xlim[2], by = byx); Yseq <- seq(Ylim[1], Ylim[2], by = byy)
Grid <- expand.grid(Xseq, Yseq)

# kNN density estimator
k=200
kNN <- knnDE(X = X, Grid = Grid, k = k)

# plot
dev.new()
#persp(Xseq, Yseq, matrix(kNN, ncol = length(Yseq), nrow = length(Xseq)))

# [optional] nice plot
persp(Xseq, Yseq, matrix(kNN, ncol = length(Yseq), nrow = length(Xseq)),ylab = "", zlab = "", theta = -70, phi = 65, ltheta = 70,col = 2, border = NA, main = "kNN", expand = 1, shade = 0.9)