Mesh (ndgrid) to points list

There’s also Iterators.product, with or without collect (usually without, if you can help it), which is essentially the same as the ((xg, yg) for xg in x, yg in y) solution above but easier to write in an abitrary number of dimensions.

julia> collect(Iterators.product(x, y))
50×50 Matrix{Tuple{Float64, Float64}}:
 (-2.0, -2.0)      (-2.0, -1.91837)      …  (-2.0, 1.91837)      (-2.0, 2.0)
 (-1.91837, -2.0)  (-1.91837, -1.91837)     (-1.91837, 1.91837)  (-1.91837, 2.0)
 (-1.83673, -2.0)  (-1.83673, -1.91837)     (-1.83673, 1.91837)  (-1.83673, 2.0)
 (-1.7551, -2.0)   (-1.7551, -1.91837)      (-1.7551, 1.91837)   (-1.7551, 2.0)
 (-1.67347, -2.0)  (-1.67347, -1.91837)     (-1.67347, 1.91837)  (-1.67347, 2.0)
 (-1.59184, -2.0)  (-1.59184, -1.91837)  …  (-1.59184, 1.91837)  (-1.59184, 2.0)
 ⋮                                       ⋱
 (1.67347, -2.0)   (1.67347, -1.91837)   …  (1.67347, 1.91837)   (1.67347, 2.0)
 (1.7551, -2.0)    (1.7551, -1.91837)       (1.7551, 1.91837)    (1.7551, 2.0)
 (1.83673, -2.0)   (1.83673, -1.91837)      (1.83673, 1.91837)   (1.83673, 2.0)
 (1.91837, -2.0)   (1.91837, -1.91837)      (1.91837, 1.91837)   (1.91837, 2.0)
 (2.0, -2.0)       (2.0, -1.91837)          (2.0, 1.91837)       (2.0, 2.0)

Although it doesn’t sound like quite what you’re after here, sometimes it can be useful to use CartesianIndices for things like this:

gridaxes = (x,y)
for ci in CartesianIndices(eachindex.(gridaxes))
    point = getindex.(gridaxes, Tuple(ci))
    # do stuff with point
end