I am trying to plot some home ranges, and I'd like to "zoom" to a smaller area. I can plot the range polygon without issue, but the problem occurs when I try to look closer at a specific area.
When I set the X and Y limits to the area of interest (using coord_sf()), the polygon is removed. I think this is because ggplot automatically removes data outside of the plotting area. For regular data, I know coord_cartesian() solves this problem, but that seems not to be an option when using polygons and geom_sf() ("Error: geom_sf() must be used with coord_sf()").
Is there some equivalent for this type of data, or perhaps a simple way to crop the polygon based on X and Y limits?
Example below.
library(sf)
library(ggmap)
library(adehabitatHR)
library(rnaturalearth)
library(lwgeom)
data <- data.frame(x = c(-50.3, -49.9, -50.0, -50.6, -55.3, -55.4, -55.5, -55.3, -54.9, -54.4, -51.5, -51.2, -50.8, -50.3),
y = c(50.3, 48.8, 48.1, 47.4, 48.2, 47.4, 50.1, 48.1, 47.5, 50.7,50.4, 50.7, 50.5, 48.3))
data$id <- as.factor("a")
#create a SpatialPointsDataFrame by defining coordinates
coordinates(data) <- c("x", "y")
# currently there is no CRS for this data. Since it is lat/long, set it as such:
proj4string(data) = CRS("+proj=longlat +datum=WGS84 +no_defs")
# now it needs to be converted to UTM
data <- spTransform(data, CRS("+init=epsg:32621 +proj=utm +zone=21 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 "))
# NOTE: I know proj4strings are no longer current but this is old code that still works for the purpose of this example
# run code to determine minimum convex polygon
mcp <- mcp(data, percent = 95, unin = c("m"),
unout = c("km2"))
# view resulting polygons by plotting them
plot(data, col = as.factor(data@data$id), pch = 16, axes = TRUE)
plot(mcp, col = alpha(1:7, 0.5),add = TRUE)
# now plot polygons on map
# import map and check projection
canada <- ne_countries(country = "canada", scale = "large", returnclass = "sf") # in rnaturalearth package
st_crs(canada)
# convert to matching utm
nl <- st_transform_proj(canada, "+init=epsg:32621 +proj=utm +zone=21 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
# now plot larger area
ggplot(data = nl) +
geom_sf(data = nl, fill = "gray80") +
coord_sf(xlim = c(100000, 1200000), ylim = c(5000000, 6100000)) + # map coordinates in utm
labs(x = "Longitude", y = "Latitude", colour = "") + # labels
geom_polygon(data = mcp, fill = "goldenrod", alpha=0.7, aes(x = long, y = lat, group = group))
# zoom to area of interest
ggplot(data = nl) +
geom_sf(data = nl, fill = "gray80") +
coord_sf(xlim = c(600000, 920000), ylim = c(5100000, 5400000)) + # map coordinates in utm
labs(x = "Longitude", y = "Latitude", colour = "") + #labels
geom_polygon(data = mcp, fill = "goldenrod", alpha = 0.7, aes(x = long, y = lat, group = group))
# polygon removed
