solution to reading .raw images into R?

well, your original data is 16-bit integers so the file is size prod(wh) * 16 + 18. When we plot in R we are assuming a colour model - i.e. we are mapping those raw numeric values to some way of encoding colours - as the ?tiff says

TIFF is a meta-format:
the default format written by ‘tiff’ is lossless and stores RGB
values uncompressed-such files are widely accepted, which is their
main virtue over PNG.

So, potentially for every numeric value we have a 3-bytes to encode RGB (0:255 along 3 axes), so we might expect prod(wh) * 3. Even though we used grey colours, those are expanded (redundantly) into RGB (all the axes have the same value, so greyscale - we have no way of knowing if our encoding is sensible for the data you have, we just stretched out white-to-black along the range of values you have in this particular file - you'd want a different data-value min-max to give absolute colours across all files).

There's some extra bytes for the TIFF format itself:

file.info("output.tif")$size
[1] 4285004
prod(wh) * 3
[1] 4284864

I don't know if tiff() can be used to produce greyscale (a single band), but it can compress so

tiff("output.tif", width=wh[1], height=wh[2], compression = "lzw")
par(c(0,0,0,0))
image(matrix(v, wh[1], wh[2])[wh[1]:1,], useRaster = TRUE, col = grey.colors(256))
dev.off()
file.info("output.tif")$size
[1] 1467824

You would also have some axis cruft around the edges I expect, so use image(, axes = FALSE, xlab = "", ylab = "") to get the cleanest output.

Oh, also beware of useRaster - in this case it probably won't do any interpolation/resampling because you are targetting a device with the dimensions of your matrix (and no margin), but to be safe I'd set useRaster = FALSE, though it will be slower.

I don't know how to render more directly to image files, it's something I wish we could do more easily (magick has some ability here and maybe does what I'm talking about).