Generate a syntetic log-normal two dimensional random field I would like to test some functions that I wrote related to the kriging applied to rain data. In order to do that, I would like to generate a synthetic log-normal 2D random field.
The idea is to extract from it some points with their two dimensional coordinates and (of course) their values. These should represent the rain stations data and locations. After that I would like to apply the functions in order to compute the semi-variogram and the kriging. Finally, I will compare the results with the syntetic 2-D random field.
I do not think that is enough to compute a bi-variate log-normal distribution. Am I right?
I have tried to search some packages for "synthetic rain data generator" but what I have found deal with time series and have a link to rain spatial distribution.
I am a newbie but I hope to have myself clear.

An easy option is to assume that there is a smooth varying trend that can be model as a two-dimensional Gaussian Process (GP) $\mathcal{GP}(s,t)$ along two continua $s$ and $t$ and sample directly from that GP. Based on the covariance function used (e.g. squared exponential, Matérn, etc.) we can control the roughness of the resulting surface as well as (mostly importantly) the relation between different points in our covariance surface.
Do note that while we might assume a zero-centred GP. If we have measurements $y_i$ at $(s_i,t_i)$ we would like to use, we could first estimate the mean two-dimensional function ${\mu}(s,t)$ using a separate function, evaluate $\hat{\mu}(s,t)$ along our whole grid and then use those estimates our GP mean. In any case, exponentiation would be done after sampling from our GP.
As always the definite readable reference for GPs is Rasmussen & Williams (2006) "Gaussian Processes for Machine Learning"; Chapt. 4 "Covariance Functions" is the one that is most relevant.
Enough talk, show me the (R) code (I tried to use as "vanilla" code as possible for readability):
x=c(-1.0, 1.0) # minimum and maximum of our continuum 
p=50 # number of grid points along each direction
N=3 # number of "random" realisations 
s=t=seq(x[1], x[2], length.out=p)
s_t_grid=expand.grid(s, t) 

mus=rep(0, p^2) # Assuming a zero-centred GP

dist_matrix=dist(gd, method="eucl")
dist_matrix_full=as.matrix(dist_matrix)

get_K_matrix=function(r, cls=1){
  # r : distance between points
  # cls : characteristic length scale
  K_st=exp(-(r^2)/(2*cls^2)) # squared exponential kernel 
  return(K_st)
}

set.seed(43)
library(MASS)
# we sample a *high*-dimensional covariance, give it a moment. 
ys_0p05=mvrnorm(N, mu=mus, Sigma=get_K_matrix(dist_matrix_full, cls=0.05))  
ys_0p50=mvrnorm(N, mu=mus, Sigma=get_K_matrix(dist_matrix_full, cls=0.50))  
ys_1p50=mvrnorm(N, mu=mus, Sigma=get_K_matrix(dist_matrix_full, cls=1.50))
# we could exponentiate these if needed.

where our ys_* variables hold N realisations of our 2-D GPs.
par(mfrow=c(3,3), pty='s',mar=c(3,0,1.5,0))
cols=hcl.colors(10, "YlOrRd")
for(i in 1:3){
  contour(s, t, matrix(ys_0p05[i,], nrow=p), col=cols, 
          main=paste0('', 'GP sample with CLS: 0.05'))
  contour(s, t, matrix(ys_0p50[i,], nrow=p), col=cols, 
          main=paste0('', 'GP sample with CLS: 0.50'))
  contour(s, t, matrix(ys_1p50[i,], nrow=p), col=cols, 
          main=paste0('', 'GP sample with CLS: 1.50'))
}


