paleolimbot / ggspatial

Enhancing spatial visualization in ggplot2
https://paleolimbot.github.io/ggspatial
368 stars 35 forks source link

Limiting extent with coord_sf() and annotation_scale() #67

Closed kstierhoff closed 4 years ago

kstierhoff commented 4 years ago

This may be similar to Issue #62, but I've recently encountered a problem when adding annotation_scale() to a working ggplot2 object with limits specified in coord_sf(). For example, this works:

library(ggplot2)
library(ggspatial)

# Load system file
nc <- st_read(system.file("shape/nc.shp", package="sf"))

# Define bounding box
map.bounds <- st_bbox(nc)

# Make basic plot, setting xlim/ylim
ggplot() + 
  geom_sf(data = nc) + 
  coord_sf(xlim = c(map.bounds["xlim"], map.bounds["xmax"]),
           ylim = c(map.bounds["ylim"], map.bounds["ymax"])) 

but this fails:

# Add annotation scale
ggplot() + 
  geom_sf(data = nc) + 
  coord_sf(xlim = c(map.bounds["xlim"], map.bounds["xmax"]),
           ylim = c(map.bounds["ylim"], map.bounds["ymax"])) + 
  annotation_scale()

with the following error message: Error in st_as_sf.data.frame(point_coords, coords = c("x", "y"), crs = sf_crs) : missing values in coordinates not allowed.

It appears that the problem exists with the named vector elements supplied by the result of st_bbox() passed to the ggplot2 object, because this seems to work just fine:

# Unname xlim/ylim
ggplot() + 
  geom_sf(data = nc) + 
  coord_sf(xlim = unname(c(map.bounds["xlim"], map.bounds["xmax"])),
           ylim = unname(c(map.bounds["ylim"], map.bounds["ymax"]))) + 
  annotation_scale()

, but I don't understand the internals of ggplot2::coord_sf() to know why.

Any ideas?

For now using unname() in the coord_sf() bit is a fine workaround, but I'm curious to know if this is expected behavior or a bug.

Thanks!

paleolimbot commented 4 years ago

Sorry for taking so long here, but I think this typo on your end:

library(sf)
#> Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
# Define bounding box
nc <- st_read(system.file("shape/nc.shp", package="sf"))
#> Reading layer `nc' from data source `/Library/Frameworks/R.framework/Versions/4.0/Resources/library/sf/shape/nc.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 14 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> CRS:            4267

map.bounds <- st_bbox(nc)

# I think you meant:
c(map.bounds["xmin"], map.bounds["xmax"])
#>      xmin      xmax 
#> -84.32385 -75.45698
c(map.bounds["xmax"], map.bounds["ymax"])
#>      xmax      ymax 
#> -75.45698  36.58965

# Not:
c(map.bounds["xlim"], map.bounds["xmax"])
#>      <NA>      xmax 
#>        NA -75.45698
c(map.bounds["ylim"], map.bounds["ymax"])
#>     <NA>     ymax 
#>       NA 36.58965

Created on 2020-07-10 by the reprex package (v0.3.0)

kstierhoff commented 4 years ago

Well that's embarrassing. I'll have to revisit my original code to see if that was the problem outside of my failed attempt at creating a reprex. Thanks for the reply!

paleolimbot commented 4 years ago

No worries! It could easily have been my problem!