Add `coord_bc()` function to specify ggplot2 bounding box in lat/long

[schckngs]

Hi, This is a great package and I'm starting to use it more often. I'm posting here to get input on mapping sub-regions of BC using ggplot with the bc_bound_hres(). Currently for quickly mapping I convert the lat/lon of my study area to UTM and enter in coord_sf(), e.g.: ggplot() + geom_sf(bc_bound_hres(), fill = "grey", colour = "black") + coord_sf(xlim = c(600000, 800000), ylim = c(745000, 1000000)) If it is possible to set limits in lat/lon it would be faster and more intuitive, but I haven't figured out how best to do that yet. Pardon if there is an obvious way of doing this. Cheers!

[ateucher]

Hi @schckngs, this is a good question!

You can set it using lat/long if you specify crs=st_crs(4326) in the coord_sf() call, but it then reprojects the whole map to lat/long WGS84, which doesn't give a faithful representation of BC (think of how it looks in Google Maps):

library(sf)
#> Linking to GEOS 3.8.1, GDAL 3.1.1, PROJ 6.3.1
library(bcmaps)
library(ggplot2)
ggplot() + 
  geom_sf(data = bc_bound(), fill = "grey", colour = "black") + 
  coord_sf(xlim = c(-132, -129), ylim = c(51.5, 54), crs = st_crs(4326))

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

I think there may be space in bcmaps for a small function (coord_bc()?) that converts it on the fly for you, it could be pretty handy! Here is a prototype. It requires a bit more testing, but let me know if this is what you are thinking of:

library(sf)
library(bcmaps)
library(ggplot2)

coord_bc <- function(xlim, ylim, crs = sf::st_crs(4326), ...) {
  box <- sf::st_bbox(c(xmin = min(xlim), xmax = max(xlim), 
                       ymin = min(ylim), ymax = max(ylim)), 
                     crs = crs)
  box <- sf::st_bbox(
    sf::st_transform(sf::st_as_sfc(box), sf::st_crs(3005))
  )
  
  ggplot2::coord_sf(xlim = box[c("xmin", "xmax")], 
                    ylim = box[c("ymin", "ymax")], 
                    ...)
}

ggplot() + 
  geom_sf(data = bc_bound(), fill = "grey", colour = "black") + 
  coord_bc(xlim = c(-132, -129), ylim = c(51.5, 54))

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

Thoughts @stephhazlitt @boshek @gcperk?

[schckngs]

Thank you! This is great.

[stephhazlitt]

I think this would be a great additional helper function in the bcmaps :package: @ateucher. I added this enhancement issue to the v0.2.0 project board. Great idea, thanks @schckngs.

[boshek]

Using this a little I am seeing that I am really choosing one of the features (or collection of features) to defining this window. Having to find the bounding box of that feature and then manually extract the limits is a little clunky. What if coord_bc directly accepted a bounding features then transformed that:

coord_bc <- function(bounding_feature = NULL, ...) {
  
  if (!inherits(bounding_feature, "sf")) {
    stop("bounding_feature is not an sf object")
  }

  box <- sf::st_bbox(
    sf::st_transform(bounding_feature, sf::st_crs(3005))
  )

  ggplot2::coord_sf(xlim = box[c("xmin", "xmax")],
                    ylim = box[c("ymin", "ymax")],
                    ...)
}

[ateucher]

I think that would be a different (but welcome) function @boshek - I think it's quite common to just want your plot zoomed in on a particular area by specifying lat/long

[boshek]

Would it make sense to have coord_bc accept both lat/long or a bbox?

[ateucher]

Possibly, though a couple of things:

1. I sort of feel that coord_bc should mimic as much as possible coord_sf, i.e., have xlim and ylim args. We could dispatch on xlim (it could be a two element numeric vector or a bbox object), but that feels clunky.

2. If you're clipping based on a specific bbox or another feature, why not do that before plotting with st_intersection? I see using coord_sf/bc as more iterative when you're plotting, which is why I like the easy flexibility of xlim and ylim (I.e., zooming in vs a-priori selecting a specific area to show).

[boshek]

I think you are definitely right that we don't want to change the syntax or behaviour of coord_sf. The attractive thing for me was that the coord_bc performed that intersection on all the layers on the fly.

[stephhazlitt]

Bumping this issue, as I came across another user looking for the feature.

[ateucher]

Just coming back to this - I realized that the original use case defined by @schckngs can be achieved by specifying xlim and ylim in lat/long, and setting the default_crs argument to 4326 (WGS84). As you can see this keeps the map in the CRS of the original data:

library(sf)
library(bcmaps)
library(ggplot2)
ggplot() + 
  geom_sf(data = bc_bound(), fill = "grey", colour = "black") + 
  coord_sf(xlim = c(-132, -129), ylim = c(51.5, 54), default_crs = 4326)

^{Created on 2023-11-20 with reprex v2.0.2}

[ateucher]

So is it worth implementing @boshek's suggestion to specify an sf[c] object or bbox? It's a more general problem than just B.C. though

[schckngs]

Hey @ateucher and all - it has been a while since I looked at this page. Thanks for keeping the conversation going! Reflecting on this - personally, being more familiar now with the sf package and coord_sf(), having a function like this is not as useful. But it sounds like it could be helpful for users new to R or sf/bcmaps (I recall getting going with sf being quite a struggle). Perhaps it's worth adding some of the code above to the documentation instead of making a new function?