2

I am trying to create a map with spatial data plotted on top. I am following the example using ggmap and then adding the New Zealand coastline using geom_sf. I have two problems:

  1. The projection of map_data(nz) doesn't match the projection of ggmap and so the coastline is not aligned to the map, even though both are in WGS84 lon, lat.
  2. When I try to apply coord_sf(xlim, ylim) to the plot I get a st_cast.POINT error:

The reprex requires a google key I'm afraid. The method in this post doesn't work for me

library(sf)
library(dplyr)
library(ggmap)

nz <- map_data("nz") %>%
  st_as_sf(coords = c("long", "lat"), remove = FALSE, crs = st_crs(4326)) %>%
  group_by(group) %>%
  summarise(
    region = region[1],
    do_union = FALSE
  ) %>%
  st_cast("LINESTRING") %>%
  ungroup()

gkey <- readLines("sjrw_google_key.dat")
register_google(key = gkey)
basemap <- get_map(location = c(lon = 175.5, lat = -38),
                   zoom = 8,
                   maptype = 'terrain-background',
                   source = 'stamen')

attr(basemap, "bb")
#>           ll.lat   ll.lon  ur.lat   ur.lon
#> bottom -39.37417 173.7449 -36.604 177.2606

ggmap(basemap) +
  geom_sf(data = nz, inherit.aes = FALSE) +
  coord_sf(crs = st_crs(4326)) +
  coord_sf(xlim = c(174.5, 176.5), ylim = c(-39.2, -36.6))
#> Coordinate system already present. Adding new coordinate system, which will replace the existing one.
# Error in st_cast.POINT(x[[1]], to, ...) : 
#  cannot create MULTILINESTRING from POINT

Created on 2021-10-21 by the reprex package (v2.0.1)

hugh-allan
  • 1,170
  • 5
  • 16
Simon Woodward
  • 1,946
  • 1
  • 16
  • 24
  • 3
    It looks like you are performing a coordinate conversion with `st_as_sf(coords = c("long", "lat"), remove = FALSE, crs = st_crs(4326))` are you absolutely positive with the CRS 4326 coordinate system? A slightly different coordinate system may be the issue here? – Dave2e Oct 20 '21 at 21:52
  • @Dave2e as far as I know the data from map_data and from ggmap are both in WGS84 (4326). – Simon Woodward Oct 20 '21 at 23:07
  • 2
    Are you 100% sure that the two layers do align? Obviously in theory the coastline of NZ should match the outline of the land area, but it looks like you are retrieving the coastline and land area from different sources. Without running your example, I wonder if they are simply different resolution/quality data, and in fact do not align perfectly anyway (the black line looks pretty coarse in your map). Maybe you could share a link to the example you are using? – hugh-allan Oct 21 '21 at 07:11
  • 1
    Thanks. It's all there in the reprex. – Simon Woodward Oct 21 '21 at 08:01
  • 1
    Based on the warning message, I'd try and combine the two `coord_sf()` lines as a start - sounds like the first one is replaced by the second anyway. What I was thinking by "sharing your example" - does the example you're following show the basemap and the coastline matching each other? Where are you following the example from? – hugh-allan Oct 29 '21 at 08:04
  • @hugh-allan this is a simplified version of my actual plot. I'm not sure why the nz coastline is misaligned, but I don't actually need it in my plot. My main problem is the inability to crop the figure. Yes it's an issue of inconsistent crs, ggmap does not seem to correctly set the crs and maybe this causes the error later. I'm also unclear whether my several geom_sf calls all need to use the same crs or not. Thanks for your help. – Simon Woodward Oct 29 '21 at 19:02
  • 1
    Ok, so does the code produce the image in the example or return an error? If it produces the image, then that's great. If you want to crop the image _exactly_ to the limits you supply, add `expand = FALSE` to `coord_sf()`, otherwise it will add a small buffer. See https://ggplot2.tidyverse.org/reference/ggsf.html . I don't know if the data all needs the same crs, but you could try setting the crs of `nz` outside `ggplot()`, before plotting maybe? – hugh-allan Oct 29 '21 at 20:38
  • 1
    And going in a slightly different direction, if don't _have to_ use `ggmap`, I would suggest you investigate `tmap` for plotting spatial data. It's pretty easy to use and is similar to `ggplot` in many ways. It might be able to do what you require https://geocompr.robinlovelace.net/adv-map.html – hugh-allan Oct 29 '21 at 20:42
  • Good suggestion thanks! The last line of the code returns an error. Without the last line it returns the figure. But I want to crop the figure. As you can see I do set the crs of nz. – Simon Woodward Oct 30 '21 at 03:16

1 Answers1

1

This is a bit tricky because I can't recreate the example without an API. I don't use ggmap myself so I'm not sure, but it's sounding like coord_sf() (or it's limits) doesn't like ggmap?

How about you retrieve exactly the area you want to plot by supplying the limits when you get the basemap, like in this example R: cropping/zooming a map

Or perhaps this example which describes downloading a specific area as the basemap https://gis.stackexchange.com/questions/155334/ggmap-clip-a-map

Additionally, you should not have two separate lines for coord_sf(); the second will override the first and you may aswell not have the first. Maybe because the crs is being overridden/ignored, the limits don't match the data coordinates? Combine these to one line and see what happens. Alternatively, don't set the crs inside coord_sf(); try using data = nz %>% st_transform(4326) in geom_sf() instead.

hugh-allan
  • 1,170
  • 5
  • 16