datacarpentry / r-raster-vector-geospatial

Introduction to Geospatial Raster and Vector Data with R
https://datacarpentry.org/r-raster-vector-geospatial
Other
113 stars 111 forks source link

Incorrect use of coord_quickmap() #372

Open djhunter opened 2 years ago

djhunter commented 2 years ago

All occurrences of coord_quickmap() should be replaced with coord_fixed(). In addition, the explanation in Episode 1 should be deleted or revised.

Explanation: It is only appropriate to use coord_quickmap() if your data is in degrees, as the source code for coord_quickmap() makes clear:

# compute distance corresponding to 1 degree in either direction
    # from the center
    x.dist <- dist_central_angle(x.center + c(-0.5, 0.5), rep(y.center, 2))
    y.dist <- dist_central_angle(rep(x.center, 2), y.center + c(-0.5, 0.5))
    # NB: this makes the projection correct in the center of the plot and
    #     increasingly less correct towards the edges. For regions of reasonnable
    #     size, this seems to give better results than computing this ratio from
    #     the total lat and lon span.

In fact, coord_quickmap() calls the helper function dist_central_angle(), which can be found in coord-munch.r.

dist_central_angle <- function(lon, lat) {
  # Convert to radians
  lat <- lat * pi / 180
  lon <- lon * pi / 180

  hav <- function(x) sin(x / 2) ^ 2
  ahav <- function(x) 2 * asin(x)

  n <- length(lat)
  ahav(sqrt(hav(diff(lat)) + cos(lat[-n]) * cos(lat[-1]) * hav(diff(lon))))
}

Sending meter measurements to this function produces nonsensical results.

For example, the first plot in Episode 1 is supposed to be showing square, 1m by 1m pixels: one meter on the horizontal axis should appear the same length as 1 meter on the vertical axis. But since coord_quickmap() chooses an incorrect aspect ratio, the scales on the axes are different and the resulting image is too narrow. The appropriate choice here is to use a 1:1 aspect ratio, which is the default ratio of coord_fixed(). (Equivalently, you could use coord_equal(), which is an alias for coord_fixed(), though this doesn't appear to be stated explicitly in the ggplot documentation.)

History: This issue appears to have been fixed in #242, but then it was subsequently broken again in 9d5314578673e6cfefe5fe777bf653e5796ab560.