library(sp)
> library(gstat)
> library(rgdal)

>xy<-as.data.frame(cbind(nutrients$Easting,nutrients$Northing))
>xy_1<-na.omit(xy)
>names(xy_1)=c("Easting","Northing")
>
>coordinates(xy_1)=~Easting+Northing


>class(xy_1)
[1] "SpatialPoints"
attr(,"package")
[1] "sp"


>Nmin<-na.omit(nutrients$Nmin)
> 
>Nmin<-as.data.frame(Nmin)
> class(Nmin)
[1] "data.frame"

locationsDD<-SpatialPointsDataFrame(xy_1,Nmin)

str(locationsDD)
Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 85 obs. of  1 variable:
  .. ..$ Nmin: num [1:85] 9.04 10.38 9.05 13.07 11.73 ...
  ..@ coords.nrs : num(0) 
  ..@ coords     : num [1:85, 1:2] 2989189 2989108 2989027 2988946 2988864 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:85] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "Easting" "Northing"
  ..@ bbox       : num [1:2, 1:2] 2988020 338723 2989189 339360
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "Easting" "Northing"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr NA


grid<-expand.grid(x=seq(from=338723.43,to=339359.52,by=10),y=seq(from=2988020.21,to=2989188.87,by=10))


vgm=variogram(Nmin~1,locationsDD)
> vgm
    np      dist    gamma dir.hor dir.ver   id
1   45  27.02613 3.098949       0       0 var1
2   30  54.12824 2.726203       0       0 var1
3  149  82.89347 3.773270       0       0 var1
4  123 100.54355 4.661074       0       0 var1
5   29 127.05775 4.861359       0       0 var1
6  216 166.31922 3.529234       0       0 var1
7   65 191.29360 4.583603       0       0 var1
8   41 220.80042 4.527452       0       0 var1
9  260 248.58705 3.750910       0       0 var1
10  76 277.89292 4.447972       0       0 var1
11  91 312.90805 5.292077       0       0 var1
12 233 333.81945 3.641161       0       0 var1
13  64 369.49771 6.194402       0       0 var1
14 201 405.77421 4.126315       0       0 var1
15 107 424.84415 4.474050       0       0 var1

fit=fit.variogram(vgm,model=vgm(1,"Exp",150,1))
> plot(vgm,fit)

> kriged=krige(Nmin~1,locationsDD,grid,model=fit)
Error in bbox(newdata) : object not a >= 2-column array