I'm working on two sets of different data in the same script; eHf and eNd. I want to use universal kriging (UK) from package gstat to interpolate both data and produce 2 different maps. But at the end, I want to produce one final map of Universal Co_kriging for both eHf and eNd.

variogram and UK maps for eHf

#Load data

    data <- read.csv("C:/Users/23030542/Desktop/Universal kriging/2627.csv")
    data.df <- data.frame(data)
    LONG <- data.df$LONG
    LAT <- data.df$LAT
    eHf <- data.df$EHFI_new
    x <- data.df$x
    y <- data.df$y

    coordinates(data.df) <- c("LONG","LAT")
    proj4string(data.df) <- CRS("+init=epsg:4326")
    SPoints <- spTransform(data.df,CRSobj = "+proj=utm +zone=51 +south +datum=WGS84 +units=m
    +no_defs")
    SP <- SpatialPoints(SPoints@coords, CRS("+proj=utm +zone=51 +south +datum=WGS84 +units=m
    +no_defs"))

#Grid

    x.range_1 <- as.numeric(c(-608380.858, 860489.815))
    y.range_1 <- as.numeric(c(6080795, 7383164))
    grd = expand.grid(x = seq(from = x.range_1[1], to = x.range_1[2], by = 10000), 
              y = seq(from = y.range_1[1], to = y.range_1[2], by = 10000))
    coordinates(grd) = ~x + y
    grd <- SpatialPixels(grd, proj4string = "+proj=utm +zone=51 +south +datum=WGS84 +units=m
    +no_defs")

#define the first order polynomial equation

    first_order = as.formula(eHf ~ x+y)

#Variogram Model

    vgm_resd_cloud <- variogram
 (first_order, SPoints, cloud = F)
    vmodel = vgm(8, "Gau", 200000, 2)
    fit.var <- fit.variogram(vgm_resd_cloud, vmodel)

#Universal Kriging

    UK = krige(first_order, SP, grd, model=vmodel)

This code works well. I proceed with the eNd data by using the same code. However, I need to distiguished the formula of the first order polynomial by using object x1 and x2 instead of x and y.

variogram and UK maps for eNd

#Load data

    data_Nd <- read.csv("C:/Users/23030542/Desktop/Universal kriging/2627_Nd.csv")
    data.df_Nd <- data.frame(data_Nd)
    LONG_Nd <- data.df_Nd$LONG
    LAT_Nd <- data.df_Nd$LAT
    eNd <- data.df_Nd$ENDI_new
    x1 <- data.df_Nd$x1
    y1 <- data.df_Nd$y1

    coordinates(data.df_Nd) <- c("LONG","LAT")
    proj4string(data.df_Nd) <- CRS("+init=epsg:4326")
    SPoints_Nd <- spTransform(data.df_Nd,CRSobj = "+proj=utm +zone=51 +south +datum=WGS84 
    +units=m
    +no_defs")
    SP_Nd <- SpatialPoints(SPoints_Nd@coords, CRS("+proj=utm +zone=51 +south +datum=WGS84 
    +units=m
    +no_defs"))

#define the first order polynomial equation

    first_order_Nd = as.formula(eNd ~ x1+y1)

#variogram for eNd

    vgm_resd_Nd <- variogram
 (first_order_Nd, SP_Nd, cloud = F)
    vmodel_Nd = vgm(4, "Gau", 200000, 2)
    fit_var_Nd <- fit.variogram(vgm_resd_Nd, vmodel_Nd)

#Universal kriging

    UK_Nd = krige(first_order_Nd, SP_Nd, grd, model=vmodel_Nd)

when I run the UK_Nd, this error occured.

Error in gstat.formula.predict(d$formula, newdata, na.action = na.action, : NROW(locs) != NROW(X): this should not occur In addition: Warning messages: 1: 'newdata' had 19257 rows but variables found have 413 rows 2: 'newdata' had 19257 rows but variables found have 413 rows

If I use x and y for eNd at the beginning instead of x1 and y1, it works just fine.

    x <- data.df_Nd$x1
    y <- data.df_Nd$y1

But like I said earlier, I need to distiguished those two because I'm going to do Co-Kriging later on. Could Someone explain this error to me. Thanks.