how to make a topoplot in R with coordinates variable distribution

Question

I would like to adapt the example I have found in this page and readapt it to my data in order to get the final topoplot realized with ggplot2

topoplot in ggplot2 – 2D visualisation of e.g. EEG data

Unfortunately it seems I having some problem with the interpolation system. Indeed when I try to plot the coordinates, they have spatial distribution very different from what they are supposed to be (my case is on the right and there is a different scale too)

This creates a problem to make work the interpolation example posted there. I tried to find some workaround and I used a different system with another package, called MBA

dat_interp <- mba.surf(cbind(dat$x, dat$y, dat$signal), no.X = 1000, no.Y = 1000)$xyz.est
mba_result <- MBA::mba.surf(cbind(dat$x, dat$y, dat$signal), 
                            no.X = 1000, no.Y = 1000)

dat_interp_df <- expand.grid(x = dat_interp$x, y = dat_interp$y)
dat_interp_df$signal <- as.vector(dat_interp$z)
ggplot() +
  geom_tile(data = dat_interp_df, aes(x = x, y = y, fill = signal)) +
  geom_point(data = dat, aes(x = x, y = y, color = signal),
             shape = 21, colour = 'black', fill = 'white', size = 2) +
  scale_fill_viridis_c() +
  scale_color_viridis_c() +
  ggtitle("Topography") +
  theme_void()

which is a result far away from what I would liek to recreate. Is there anyone that can help me with this out. Here as follows my data:

dput(dat)
structure(list(X = 1:64, label = c("AF3", "AF4", "AF7", "AF8", 
"AFz", "C1", "C2", "C3", "C4", "C5", "C6", "CP1", "CP2", "CP3", 
"CP4", "CP5", "CP6", "CPz", "Cz", "F1", "F2", "F3", "F4", "F5", 
"F6", "F7", "F8", "FC1", "FC2", "FC3", "FC4", "FC5", "FC6", "FCz", 
"FT7", "FT8", "Fp1", "Fp2", "Fpz", "Fz", "Iz", "O1", "O2", "Oz", 
"P1", "P10", "P2", "P3", "P4", "P5", "P6", "P7", "P8", "P9", 
"PO3", "PO4", "PO7", "PO8", "POz", "Pz", "T7", "T8", "TP7", "TP8"
), x = c(-33.706, 35.738, -54.873, 55.806, 0.23117, -36.269, 
37.675, -65.418, 67.104, -80.273, 83.413, -35.631, 38.16, -64.113, 
66.3, -79.863, 83.134, 0, 0.47908, -27.501, 29.475, -50.221, 
51.912, -64.49, 67.878, -70.233, 73.024, -34.049, 34.793, -60.191, 
62.341, -77.199, 79.561, 0.37692, -80.826, 81.786, -29.391, 29.902, 
0.11233, 0.31257, 0, -30.195, 29.937, 0, -29.247, 74.343, 32.545, 
-53, 55.408, -66.98, 68.381, -72.811, 72.81, -73.635, -36.844, 
36.608, -54.282, 55.984, 0, 0, -84.117, 85.068, -85.242, 85.591
), y = c(76.785, 77.701, 68.493, 69.657, 80.761, -9.7183, -9.3934, 
-11.535, -10.868, -13.721, -12.764, -47.284, -47.123, -46.581, 
-46.424, -46.109, -46.082, -47.183, -9.1413, 56.889, 57.599, 
53.107, 54.323, 48.069, 49.861, 42.535, 44.399, 26.032, 26.41, 
22.744, 23.681, 18.676, 19.984, 27.439, 14.106, 15.454, 83.929, 
84.91, 88.283, 58.525, -118.78, -112.69, -111.73, -115.07, -80.356, 
-74.343, -80.553, -78.576, -79.131, -77.052, -75.945, -72.811, 
-72.81, -73.635, -101.23, -100.58, -97.927, -96.966, -102.2, 
-81.305, -16.351, -15, -45.324, -45.51), signal = c(0.356912525947854, 
0.441154040587387, 0.127495728211593, 0.209844264592143, 0.471423541381149, 
0.920244846231276, 1.20067778812602, 0.901957655651085, 1.26828922845346, 
0.913597538349512, 0.961279110570437, 1.70669913750882, 1.71209354564514, 
1.54924715577636, 1.69030506815152, 1.57151232049429, 1.52222537801321, 
1.7682950526965, 1.00052912202385, 0.258051875687978, 0.473055705981418, 
0.281706344592271, 0.425239318349764, 0.391951492462168, 0.265733842589671, 
0.207140129446265, 0.237923276586573, 0.284985208364728, 0.587702684164806, 
0.337358918978001, 0.722237508242702, 0.532691597382293, 0.613087103723537, 
0.395202773422108, 0.719928099089211, 0.256588947956852, 0.433203966266369, 
0.360419814621179, 0.490884077829908, 0.391080737042044, 0.494551323229936, 
1.33255331306293, 1.1550019409092, 0.943146310863848, 2.13507543110418, 
0.605380856482758, 2.05759696362809, 2.16419480467961, 1.95371325977124, 
1.99675377322271, 1.67495247573419, 1.64695492742097, 1.33445815225302, 
0.772128251322973, 2.07429881383391, 1.60452882592484, 1.53601914533575, 
1.38352896104063, 1.83063042527717, 2.17298717559065, 0.841149491013589, 
0.443649019354451, 1.28666049614812, 0.982024863842714)), row.names = c(NA, 
-64L), class = c("tbl_df", "tbl", "data.frame"))

Answer 1

Base R reducing package dependencies.

v0, eggheadplot

1. Re-scale xy-coords.
2. Set n, x, y, calculate r.
3. Initialise MBA::mba.surf() with it's argument n=r around 1.2, set extend=FALSE, see help file. Not of relevance here, but maybe in the future.
4. Plot with image() + points() + contour() + polygon() (+ chull()).

 dat[c('x', 'y')] = lapply(dat[c('x', 'y')], \(x) (x - min(x))/(max(x) - min(x)))
 n = 1e3; x = y = seq(0, 1, length.out=n)
 r = diff(range(dat$y)) / diff(range(dat$x))
 Z = MBA::mba.surf(dat[c('x', 'y', 'signal')], length(x), length(y), 
                   n=r, extend=FALSE)$xyz.est
 par(xpd = TRUE)
 image(Z, col=hcl.colors(n = n, palette = 'viridis'), asp = 1/r)
 polygon(dat[c('x', 'y')][chull(dat[c('x', 'y')]), ])
 points(dat$x, dat$y, col='white', pch=20)
 contour(x, y, Z$z, add=TRUE, col='white', labcex=1)
 polygon(c(.47, .5, .53), c(.99, 1.03, .99), col='white', border='black') # nose 

Notice that the contour lines are somewhat off. This is due to extend=FALSE.

Notes

(1) signal range

> range(dat$signal)
[1] 0.1274957 2.1729872

(2) hcl.colors(n = n, palette = 'viridis') is not the default colour palette used in image(). A brilliant reference.

(3) We might want to use different spatial interpolation techniques.

Edit

v1, eggheadplot

eggheadplot = \(d, 
                n=100, 
                extend=FALSE, 
                col=NULL, 
                # reversed matlab-like colour palette blue2green2red 
                hull=TRUE, 
                nose=TRUE) {
   
   stopifnot(is.data.frame(d), all(c('x', 'y', 'signal') %in% names(d)))
   d = d[c('x', 'y', 'signal')]
   if (missing(n)) n = 1e3
   x = y = seq.int(0, 1, length.out=n)
   d[c('x', 'y')] = lapply(d[c('x', 'y')], \(x) (x - min(x)) / (max(x) - min(x)))
   r = diff(range(d$y)) / diff(range(d$x))
   Z = MBA::mba.surf(d, n, n, n=r, extend=extend)$xyz.est
   def_par = par(no.readonly = TRUE)
   on.exit(def_par)
   layout(t(1:2), widths=c(10,1))
   par(mar=rep(.5, 4), oma=rep(3, 4))
   p = .05
   if (missing(col)) col = hcl.colors(n = n, palette = 'viridis')
   image(Z, col=col, asp=1/r, ylim=c(0, max(Z$y) + p * diff(range(Z$y))), axes=FALSE)
   if (hull) polygon(d[c('x', 'y')][chull(d[c('x', 'y')]), ], lwd = 2)
   points(d$x, d$y, col='white', pch=20)
   contour(x, y, Z$z, add=TRUE, col='white', labcex=1)
   if (nose) polygon(c(.45, .5, .55), c(.999, 1.05, .999), col='white', border='black')
   axis(1, seq(0, 1, .1))
   axis(2, seq(0, 1, .1), las=1)
   # if (col_legend) {
      v = sort(unlist(Z$z))
      image(1, v, t(v), col=col, axes=FALSE)
      axis(4, seq(0, max(v) + .25, .25))
   #}
}

eggheadplot(dat)

eggheadplot(dat, extend = TRUE)

eggheadplot(dat, col=rev(colorRamps::matlab.like2(100)))

The exporting tool crops the top of the colour legend. You will see the complete colour legend locally.

Answer 2

You can mask your interpolation grid with an ellipse and draw a nose polygon. I also added some white contour lines. Note, that if we use this masked approach, we need to extend the bounds of mba.surf in bbox. I added theme_bw() for the white background and coord_fixed() to retain the axis-ratio provided in data.

library(MBA)
library(ggplot2)

# Nose
hr <- min(max(abs(dat$x)), max(abs(dat$y))) *1.1
nose_df <- data.frame(x = c(-0.1, 0, 0.1) * hr, 
                      y = c(0.98, 1.1, 0.98) * hr) # maybe adjust thes triangle points!

# head ellipse
theta <- seq(0, 2*pi, length.out = 100)
a <- hr * 0.97   # horizontal radius
b <- hr * 1.15   # vertical radius
y_offset <- 15
head_df <- data.frame(x = a * cos(theta), y = b * sin(theta) - y_offset)

dat_interp <- mba.surf(cbind(dat$x, dat$y, dat$signal), 
                       no.X = 1000, 
                       no.Y = 1000,
                       extend = TRUE, 
                       b.box=c(min(dat$x)-20, max(dat$x)+20, min(dat$y)-20, max(dat$y)+20) # I extended the grid a little
                       )$xyz.est
mba_result <- MBA::mba.surf(cbind(dat$x, dat$y, dat$signal), 
                            no.X = 1000, no.Y = 1000,
                            extend = TRUE)
dat_interp_df <- expand.grid(x = dat_interp$x, y = dat_interp$y)
dat_interp_df$signal <- as.vector(dat_interp$z)

# Create mask for ellipse
dat_interp_df <- dat_interp_df[((dat_interp_df$x / a)^2 + ((dat_interp_df$y + y_offset) / b)^2) <= 1, ]

# plot
ggplot() +
  geom_tile(data = dat_interp_df, aes(x = x, y = y, fill = signal)) +
  geom_contour(data = dat_interp_df, aes(x = x, y = y, z = signal), colour = 'white', alpha = 0.5) +
  geom_point(data = dat, aes(x = x, y = y, color = signal),
             shape = 21, colour = 'black', fill = 'white', size = 2) +
  geom_path(data = head_df, aes(x = x, y = y), color = "black", size = 0.5) + # head
  geom_polygon(data = nose_df, aes(x = x, y = y), fill = "white", colour = 'black') + # nose
  ggtitle("Topography") +
  scale_fill_viridis_c() +
  scale_color_viridis_c() +
  theme_bw() +
  coord_fixed()