scallops <- read.table("scallops.txt", header=T)

myscallops <- scallops
myscallops[,"lgcatch"] <- log(scallops$tcatch+1)

library(geoR)

R <- 6371
lat <- myscallops$lat*pi/180
lon <- myscallops$long*pi/180

u <- cbind(R*cos(lat)*cos(lon), R*cos(lat)*sin(lon), R*sin(lat))

obj <- cbind(myscallops$long,myscallops$lat,myscallops$lgcatch)

scallops.geo <- as.geodata(obj,coords.col=1:2,data.col=3)

scallops.var <- variog(scallops.geo, estimator.type="classical")
scallops.var.robust <- variog(scallops.geo, estimator.type="modulus")

par(mfrow=c(2,1))
plot(scallops.var)
plot(scallops.var.robust)

scallops.var.fit <- variofit(scallops.var.robust,ini.cov.pars=c(1.0,2.0),cov.model="exponential", fix.nugget=FALSE,nugget=1.0)

scallops.lik.fit <- likfit(scallops.geo,ini.cov.pars=c(1.0,2.0),cov.model="exponential",trend="cte", fix.nugget=FALSE,nugget=1.0, nospatial=FALSE)

detach(package:geoR)