Skip to content

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

@alexyshr

Description

@alexyshr

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

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions