Heatmap and contour plots transposed and broken for non-square grids

The plot functions heatmap and contourf give incorrect results when plotting GaussianRandomField objects defined on non-square grids.

For example

using GaussianRandomFields, Plots, Random
Random.seed!(87234632804762304832)
cov = CovarianceFunction(2, Matern(1/2, 5/4))
grf = GaussianRandomField(cov, GaussianRandomFields.Cholesky(), range(0, 2, 50), range(0, 1, 25))
heatmap(grf)

gives the output

while

using GaussianRandomFields, Plots, Random
Random.seed!(87234632804762304832)
cov = CovarianceFunction(2, Matern(1/2, 5/4))
grf = GaussianRandomField(cov, GaussianRandomFields.Cholesky(), range(0, 2, 25), range(0, 1, 50))
heatmap(grf)

and

using GaussianRandomFields, Plots, Random
Random.seed!(87234632804762304832)
cov = CovarianceFunction(2, Matern(1/2, 5/4))
grf = GaussianRandomField(cov, GaussianRandomFields.Cholesky(), range(0, 2, 25), range(0, 1, 50))
contourf(grf)

gives a blank axes as output and the message Arrays have incorrect length or dimension..

I believe the issue is that the Plot.jl heatmap and contourf methods expect the z argument corresponding to the two-dimensional field to be plotted, to be an array of shape (dim_y, dim_x) where dim_y is the number of points in the spatial grid in the vertical (y) dimension and dim_x is the number of points in the spatial grid in the horizontal (x) dimension, such that z[i, j] corresponds to the field value at the jth point along the x-axis and ith point along the y-axis, which I believe is the opposite of the dimension ordering convention in GaussianRandomFields for two-dimensional fields.

I think simply transposing the passed sampled field should give the correct result. For example, calling the plot methods manually

using GaussianRandomFields, Plots, Random
Random.seed!(87234632804762304832)
cov = CovarianceFunction(2, Matern(1/2, 5/4))
grf = GaussianRandomField(cov, GaussianRandomFields.Cholesky(), range(0, 2, 50), range(0, 1, 25))
z = sample(grf)
contourf(grf.pts..., transpose(z))

gives a more reasonable appearing output

I think this can probably just be fixed by changing sample(grf) to tranpose(sample(grf)) in the recipe below

https://github.com/PieterjanRobbe/GaussianRandomFields.jl/blob/5d9257cb452e2b91387e206a1d074da63a81156f/src/plot.jl#L137-L140

If this seems reasonable I am happy to submit a PR making this change.

PieterjanRobbe / GaussianRandomFields.jl

Heatmap and contour plots transposed and broken for non-square grids #41