Bdtjtk

library(rgdal)

library(maptools)

library(gstat)

library(sp)



data <- read.table("meuse.txt", sep="", header=TRUE) # read txt file



# transform the data frame into a spatial data frame

coordinates(data) <- ~ x + y



## Set the coordinate system

proj4string(data) <- CRS("+init=epsg:4326")



## the epsg numbers can be found here: http://spatialreference.org/ref/



# import the border shp file

border <- readOGR("meuse_area.shp", "meuse_area")

proj4string(border) <- CRS("+init=epsg:4326")



# import a raster from a ArcInfo ASCII format

zinc <- read.asciigrid("zinc.asc")

proj4string(zinc) <- CRS("+init=epsg:4326")



# Let's first create a prediction grid for the interpolation, starting from 

# the shape file

vals <- border@bbox

deltaLong <- as.integer((vals[1, 2] - vals[1, 1]) + 1.5)

deltaLat <- as.integer((vals[2, 2] - vals[2, 1]) + 1.5)

gridRes <- 0.5   # change this value to change the grid size (in metres)

gridSizeX <- deltaLong / gridRes

gridSizeY <- deltaLat / gridRes

grd <- GridTopology(vals[, 1], c(gridRes, gridRes), c(gridSizeX, gridSizeY))

pts <- SpatialPoints(coordinates(grd))

pts1 <- SpatialPointsDataFrame(as.data.frame(pts), 

                               data=as.data.frame(rep(1, nrow(as.data.frame(pts)))))



Overlay <- overlay(pts1, border)

pts1$border <- Overlay

nona <- na.exclude(as.data.frame(pts1))

coordinates(nona) <- ~ x + y

gridded(nona) <- TRUE

proj4string(nona) <- CRS("+init=epsg:4326")  # remember to set the coordinate 

                                             # system also for the prediction grid

writeAsciiGrid(nona, "prediction_grid.asc")





# For the Co-kriging we need to obtain the value of the covariate for each 

# observation

over <- overlay(zinc, data)

data$zinc <- over$zinc.asc

str(as.data.frame(data))



# also the prediction grid need to be overlayed with the covariate

over <- overlay(zinc, nona)

nona$zinc <- over$zinc.asc



# for the cokriging, the first thing to do is create an object with the

# function gstat() that contains both the variable and the covariate

str(data)

complete.cases("data")

str(zinc)

complete.cases("zinc")

g <- gstat(id="lead", formula=lead ~ 1, data=data)

g <- gstat(g, id="zinc", formula=zinc ~ 1, data=data)





# Fitting the variogram

# first, plot the residual variogram

vario <- variogram(g)

Error in na.fail.default(list(zinc = c(NA, NA, NA, NA, NA, NA, NA, NA,
: missing values in object

I know that there is no missing in zinc when I edit the object in Notepad. What did I miss? There is no NA in zinc.asc.These are my data.

I want to perform cokriging and I am stuck with variogram.

One cause of the confusion here may be that sp used to include it's own overlay function.

https://www.rdocumentation.org/packages/sp/versions/0.9-7/topics/overlay

However the most recent version of the sp package replaced the "overlay" function with "over". I do no believe raster::overlay behaves the same as the overlay function from sp 0.9-7