Why does gstat.predict() function often return NaN values (GSTAT Package)? (R version 3.3.2, Windows 10)

Question

I am trying to simulate a combination of two different random fields (yy1 and yy2 with different mean and correlation length) with an irregular boundary using Gstat package in R. I have attached the picture of my expected outcome. The code is not giving such output consistently and I am frequently getting atleast one of the yy1 and yy2 as NaNs, which results in the Undesired output as shown in image.

The key steps I used are:

1) Created two gstat objects with different means and psill (rf1 and rf2) 2) Created two computational grids (one for each random field) in the form of data frame with two variables “x” and “y” coordinates. 3) Predicted two random fields using unconditional simulation.

Any help in this regard would be highly appreciated.

Attachments: 2 images (link provided) and 1 R code

1) Expected Outcome

2) Undesired Outcome

    library(gstat)

    xy <- expand.grid(1:150, 1:200) # 150 x 200 grid is created in the form of a dataframe with x and y vectors

   names(xy)<-c('x','y') # giving names to the variables

# creating gsat objects

    rf1<-gstat(formula=z~1,locations=~x+y,dummy = T,beta=c(1,0,0),  model=vgm(psill=0.025, range=5, model='Exp'), nmax=20)  # dummy=T treats this as  a unditional simulation
    rf2<-gstat(formula=z~1,locations=~x+y,dummy = T,beta=c(4,0,0),  model=vgm(psill=0.025, range=10, model='Exp'), nmax=20)  # dummy=T treats this as a unditional simulation

# creating two computational grid

    rows<-nrow(xy)

    xy_shift <- expand.grid(60:90, 75:100)
    names(xy_shift)<-c('x','y')

    library(dplyr) # for antijoin

    xy1<-xy[1:(rows/2),]
    xy1<-anti_join(xy1, xy_shift, by = c("x","y")) # creating the irregular boundary

    xy2<-rbind(xy[(rows/2+1):rows,],xy_shift)

    library(sp)

    yy1<- predict(rf1, newdata=xy1, nsim=1) # random field 1
    yy2<- predict(rf2, newdata=xy2, nsim=1) # random field 2

    rf1_label<-gl(1,length(yy1[,1]),labels="field1")
    rf2_label<-gl(1,length(yy2[,1]),labels="field2")

    yy1<-cbind(yy1,field=rf1_label)
    yy2<-cbind(yy2,field=rf2_label)

    yy<-rbind(yy1,yy2)

    yyplot<-yy[,c(1,2,3)]

# plotting the field

    gridded(yyplot) = ~x+y
    spplot(obj=yyplot[1],scales=list(draw = TRUE))