geom_tile squares aren't square in geom_sf plate carrée

[Kodiologist]

When a I draw a square with geom_tile in a plot with coord_sf(crs = 4326), the result appears taller than it is wide. Maybe there's a problem with geom_tile, or maybe coord_sf(crs = 4326) isn't actually using plate carrée (in which one degree of longitude is the same length as one degree of latitude). I'm not sure which.

#3654 is possibly related.

library(ggplot2)

crs = 4326 # https://epsg.io/4326

borders = rnaturalearth::ne_states(iso_a2 = "US", returnclass = "sf")
borders = borders[tolower(borders$postal) == "ny",]

ggplot() +
    geom_tile(aes(x, y, width = 5, height = 5),
        data = data.frame(y = 40.7831, x = -73.9712)) +
    geom_sf(data = borders, fill = NA) +
    coord_sf(datum = NA, crs = crs) +
    theme_void()

ggsave("/tmp/test.png", width = 1, height = 1)

[clauswilke]

I cannot run your example. The rnaturalearth::ne_states() call fails. It also doesn't help to hide the graticule. Below is a reprex that avoids rnaturalearth and doesn't hide the graticule. It shows that geom_tile() does exactly the right thing, but the aspect ratio of the entire plot is indeed distorted.

library(ggplot2)

crs = 4326 # https://epsg.io/4326

nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)

ggplot() +
  geom_sf(data = nc) +
  geom_tile(
    aes(x, y, width = 2, height = 2),
    data = data.frame(x = -83, y = 34)) +
  coord_sf(crs = crs) 

^{Created on 2020-01-07 by the reprex package (v0.3.0)}

Not sure whether this is a problem with ggplot2 or with sf. @edzer?

[edzer]

geom_sf() doesn't use plate carree, but uses equidistant rectangular. The aspect ratio of long and lat is chosen such that in the middle of the plot area, 1 m north equals 1 m east. Only when the middle of your plot area coincides with the equator this corresponds to plate carree.

[Kodiologist]

Huh. Is it possible to get plate carrée by some other means, then, or is that just not implemented for ggplot2?

[edzer]

You can convert your data to plate carree.

x = st_transform(nc, "+proj=eqc")
ggplot()+ geom_sf(data = x[1])

[Kodiologist]

Thanks. Oddly, I can't get the transformed geom_tile to display. Perhaps I'm getting something wrong.

library(ggplot2)
library(sf)

crs.lonlat = 4326 # https://epsg.io/4326

nc = sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
nc = st_transform(nc, "+proj=eqc")

d = data.frame(x = -83, y = 34)
d = as.data.frame(st_coordinates(st_transform(
  st_as_sf(d, coords = c("x", "y"), crs = crs.lonlat),
  "+proj=eqc")))

ggplot() +
  geom_sf(data = nc) +
  geom_tile(
    aes(X, Y, width = 2, height = 2),
    data = d) +
  coord_sf(crs = "+proj=eqc")

ggsave("/tmp/test.png", width = 3, height = 2)

[edzer]

Units are now meter, and you're drawing a tile of 2 x 2 m. Try values of 200000 or so.

[Kodiologist]

That's confusing. It looks like the correspondence between meters and degrees is set at the equator, so each degree is equal to about 40075017 / 360 of the projection's nominal meters, where the numerator is the circumference of the equator in meters.

At any rate, I guess my question is answered, so thank you, but you may want to clarify these issues in the documentation of geom_sf or coord_sf.

[clauswilke]

With PR #3659 that I'm currently finalizing, this becomes very simple. By default, non-sf geoms now work in long-lat coordinates.

library(ggplot2)
library(sf)
#> Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)

ggplot() +
  geom_sf(data = nc) +
  geom_tile(
    aes(x, y, width = 2, height = 2),
    data = data.frame(x = -83, y = 34)
  ) +
  coord_sf(crs = "+proj=eqc")

^{Created on 2020-02-01 by the reprex package (v0.3.0)}

cc @yutannihilation

[clauswilke]

The same using a different projection.

library(ggplot2)
library(sf)
#> Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)

ggplot() +
  geom_sf(data = nc) +
  geom_tile(
    aes(x, y, width = 2, height = 2),
    data = data.frame(x = -83, y = 34)
  ) +
  coord_sf(crs = 3347)

^{Created on 2020-02-01 by the reprex package (v0.3.0)}

[yutannihilation]

Now that #3659 gets merged, I thought we can close this issue, but it seems crs = "+proj=eqc") no longer works because of the change in PROJ, iiuc. I found EPSG:4087 is Equidistant Cylindrical but the unit is in metre, so I have no idea.

I'm closing this anyway, not mainly because the issue is resolved, but because there seems nothing we can do on ggplot2's side. Please correct me if I'm wrong!

[edzer]

FYI - I don't see any problems using crs = "+proj=eqc" - which error do you see?

[yutannihilation]

Oh, thanks. I see no error, but a different result with the same code as #3698 (comment). I thought this was somehow related to the deprecation of proj4string-notation, but I have no idea.

library(ggplot2)
library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUE

nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)

ggplot() +
  geom_sf(data = nc) +
  geom_tile(
    aes(x, y, width = 2, height = 2),
    data = data.frame(x = -83, y = 34)
  ) +
  coord_sf(crs = "+proj=eqc")

^{Created on 2022-08-21 with reprex v2.0.2}

[edzer]

I saw that too, but this doesn't seem a PROJ problem: if we use + coord_sf(crs = st_crs(nc)) or + coord_sf(crs = st_crs(4326)) things look good, if we for instance use + coord_sf(crs = st_crs(3857)) (web Mercator) things look bad.

My guess is that the x = -83, y = 34 are interpreted in the CRS of st_coord(), meaning they practically land at geodetic coordinates long/lat 0, 0, which might actually be the right thing to do:

> st_sfc(st_point(c(-83, 34)), crs = st_crs(3857)) |> st_transform(st_crs(nc))
Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -0.0007456017 ymin: 0.0003054272 xmax: -0.0007456017 ymax: 0.0003054272
Geodetic CRS:  NAD27
POINT (-0.0007456017 0.0003054272)
> st_sfc(st_point(c(-83, 34)), crs = st_crs("+proj=eqc")) |> st_transform(st_crs(nc))
Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -0.0007456017 ymin: 0.0003054272 xmax: -0.0007456017 ymax: 0.0003054272
Geodetic CRS:  NAD27
POINT (-0.0007456017 0.0003054272)

[edzer]

I think this is all in agreement with the ?coord_sf documentation, and, as the documentation says, the plot above is obtained by:

library(ggplot2)
library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUE

nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)

ggplot() +
  geom_sf(data = nc) +
  geom_tile(
    aes(x, y, width = 2, height = 2),
    data = data.frame(x = -83, y = 34)
  ) +
  coord_sf(crs = "+proj=eqc", default_crs = st_crs(4326))

[yutannihilation]

Yes, you are right. Sorry, it seems I was a bit confused. I tired some combination of crs and default_crs, but couldn't get that result, probably because I specified them oppositely... Anyway, glad that this issue is confirmed to be resolved. Thanks so much for your help!

library(ggplot2)
library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUE

nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)

ggplot() +
  geom_sf(data = nc) +
  geom_tile(
    aes(x, y, width = 2, height = 2),
    data = data.frame(x = -83, y = 34)
  ) +
  coord_sf(crs = st_crs(4087), default_crs = st_crs(4326))

^{Created on 2022-08-21 with reprex v2.0.2}