"CRS object has comment, which is lost in output" in krige

Question

I am trying to krige for water quality dataset with Latitude, longitude which using CRS("+init=epsg:4326").

GGT <- read.csv("C:/Users/user/Data/newdata2019.csv")
coordinates(GGT) = ~Lon+Lat
GGT <- st_as_sf(GGT)
st_crs(GGT) <- 4326

GGTgrid <- readOGR('C:/Users/user/Desktop/FisheryScience/Data/Maps/GGTgrid.shp')
GGTgrid1 <- st_as_stars(GGTgrid, crs = 4326)
st_crs(GGTgrid1) <- 4326

-This is what GGT dataset looks like

vario <- variogram(log(DO_S)~1, GGT)
model_GGT <- fit.variogram(vario, model=vgm(psill = 1, model= 'Sph', range= 200, nugget =1))
plot(vario, model = model_GGT)

-It seems to work fine by here

Then it shows error when I run the code block below

krige_result <- krige(formula = log(DO_S)~1, GGT, GGTgrid1, model = model_GGT)

with long lines of error

Warning message in proj4string(obj):
"CRS object has comment, which is lost in output"Warning message in proj4string(obj):
"CRS object has comment, which is lost in output"Warning message in proj4string(obj):

Answer 1

The warnings you get are indicating that your work may be affected by the change introduced with PROJ 6 (and GDAL 3), adopted by R-spatial and rspatial. You can get all the details using this two links:

https://rgdal.r-forge.r-project.org/articles/CRS_projections_transformations.html

https://r-spatial.org/r/2020/03/17/wkt.html

To make these warning messages disappear, you just have to use objects of type sf and stars (by installing/loading the packages of the same names) which take into account these recent changes. So, I suggest you use the following few lines of code at the beginning of your script to replace your first six lines of code. This will give you two objects (i.e. GGT of type sf and GGTgrid1 of type stars) :

GGT <- read.csv("C:/Users/user/Data/newdata2019.csv")
coordinates(GGT) = ~Lon+Lat
GGT <- st_as_sf(GGT)
st_crs(GGT) <- 4326

GGTgrid <- readOGR('C:/Users/user/Desktop/maps/GGTgrid.shp')
coordinates(GGTgrid) <- ~x+y
GGTgrid1 <- st_as_stars(GGTgrid1, crs = 4326)
st_crs(GGTgrid1) <- 4326

It is easier for me to work on real data than to work "virtually" with the name of your objects as I don't have your original files. So I prefer to show you how to proceed for your analysis with the "meuse" data contained in the sp package.

By analogy with the reprex I give you, I think you should be able to manage with your own files. And you will see, no more warning message will appear :-)

Please, find below my reprex.

Reprex

• Loading the library and the data

library(sp)
library(sf)
library(stars)
library(gstat)

data(meuse)                # loading the data (equivalent of your csv file)
coordinates(meuse) = ~x+y  # you already know this step ;-)

# Just a look to the class of original data 
class(meuse)
#> [1] "SpatialPointsDataFrame"
#> attr(,"package")
#> [1] "sp"                # "meuse" is an object of class 'sp'


data(meuse.grid)           # loading the data (equivalent of your shp file)
gridded(meuse.grid) = ~x+y

# Just a look to the class of original data 
class(meuse.grid)
#> [1] "SpatialPixelsDataFrame"
#> attr(,"package")
#> [1] "sp"                # "meuse" is an object of class 'sp'

• Converting the meuse data into sf object and the meuse.grid data into stars object

# Convert 'sp' object 'meuse' (i.e. SpatialPointsDataFrame) into 'sf' object
meuse <-  st_as_sf(meuse)

class(meuse) 
#> [1] "sf"         "data.frame"    # meuse is indeed of class 'sf'

# Convert 'sp' object 'meuse.grid' (i.e. SpatialPixelDataFrame) into 'stars' object
meuse.grid <-  st_as_stars(meuse.grid)

class(meuse.grid) 
#> [1] "stars"                       # meuse.grid is indeed of class 'stars'

• Compute and plot the variogram

vario <-  variogram(log(zinc)~1, meuse)
model_meuse <-  fit.variogram(vario, model = vgm(psill = 1, model = "Sph", range = 200, nugget = 1))
plot(vario, model = model_meuse)

• Krige and plot predictions and variances

krige_result <- krige(formula = log(zinc)~1, meuse, meuse.grid, model = model_meuse)
#> [using ordinary kriging]
class(krige_result)
#> [1] "stars"
krige_result
#> stars object with 2 dimensions and 2 attributes
#> attribute(s):
#>                  Min.   1st Qu.    Median      Mean   3rd Qu.      Max. NA's
#> var1.pred  4.77655207 5.2376428 5.5728908 5.7072284 6.1717618 7.4399908 5009
#> var1.var   0.08549102 0.1372838 0.1621815 0.1853301 0.2116141 0.5002793 5009
#> dimension(s):
#>   from  to offset delta refsys point values x/y
#> x    1  78 178440    40     NA    NA   NULL [x]
#> y    1 104 333760   -40     NA    NA   NULL [y]

plot(krige_result[1]) # plot predictions

plot(krige_result[2]) # plot variances

^{Created on 2021-10-19 by the reprex package (v2.0.1)}