r-spatial/gstat

Behaviour of gstat::variogram when using geographic coordinates: defining CRS and setting CRS to NA

Opened this issue · 0 comments

Greetings Dr. Edzer!

I have learned that I must use projected coordinates to do geostatistics (instead of using geographic coordinates). Nonetheless, I just want to be sure that the behavior of gstat::variogram when using geographic coordinates is right and not a bug.

Why the distances of the semivariogram are between 0 and 150 considering that the maximum distance between points is 4.590449?

If I use geographic coordinates without defining the coordinate system st_crs(st_grd_sample) <- NA, then the semivariogram distances seem fine.

Thank you!

library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(ggplot2)

# Polygon of the state of Ohio, USA
ohio = sf::st_as_sf(maps::map("state", "ohio", plot = FALSE, fill = TRUE))

#Sample 200 points
st_grd_sample = sf::st_sf(sf::st_sample(ohio, 200))
#It has WGS84 geographic coordinates
sf::st_crs(st_grd_sample)$input
#> [1] "EPSG:4326"
#Review max distances between points
the_coords = sf::st_coordinates(st_grd_sample)
the_dist = dist(the_coords)
max(the_dist)
#> [1] 4.590449

#assign random values between 0 and 100
st_grd_sample$value = runif(dim(st_grd_sample)[1], 0, 100)

#
#Geographic Coordinates using CRS
#
y.v2 = gstat::variogram(value~1, loc=st_grd_sample)
#The distances for the semivariogram are between 0 and more than 150 (wrong?)
plot(y.v2$dist, y.v2$gamma, 
     xlab="Distance", ylab=expression(gamma(h)))

#
#Geographic Coordinates without defining CRS
#
st_crs(st_grd_sample) <- NA
y.v2 = gstat::variogram(value~1, loc=st_grd_sample)
#The distances for the semivariogram are ok
plot(y.v2$dist, y.v2$gamma, 
     xlab="Distance", ylab=expression(gamma(h)))

Created on 2023-04-11 with reprex v2.0.2