R – Correcting Plotted Linestrings Created in R

rsf

I need to create multiple linestrings (transects) based on multiple coordinates, and then plot them (for a study area map) and later on, I will use these lines for joining attributes from another data set and in other spatial analyses.

The first part –creating the linestrings– seems to work fine for me, based on all (great) answers in R – Create multiple linestrings from multiple coordinates & Guzmán's answer in R: sf package points to multiple lines with st_cast. If I use mapview::mapview(), they all return the transects. But when I try to plot it through ggplot2 it returns me:

Error in st_transform.sfc(X[[i]], …) :
cannot transform sfc object with missing crs

When I set a crs and try to plot it, it plots in a completelly different location than expected (pretty much out of the globe, hehe). I've tried setting the crs in different parts of the code, check the ## in the code below, e.g. inside the function in st_linestring, and outside the function in st_sfc, st_sf, or even using st_crs(lines) <- 4326.

Please find below the code I'm using.

library(tidyverse); library(rnaturalearth); library(sf); library(mapview)

df_tbl <- tibble::tribble(
  ~lat,   ~lon,   ~lat2,   ~lon2,
  -34.305, 173.520, -34.461, 173.258,
  -34.461, 173.258, -34.214, 173.378,
  -34.214, 173.378, -34.362, 173.113)

rows <- split(df_tbl, seq(nrow(df_tbl)))

lines <- lapply(rows, function(row) {
  lmat <- matrix(unlist(row[1:4]), ncol = 2, byrow = TRUE)
  st_linestring(lmat) # I've tried crs = 4326 here
})

lines <- st_sfc(lines) # I've tried crs = 4326 here

lines_sf <- st_sf('ID' = 1:3, 'geometry' = lines) # I've tried crs = 4326 here

# st_crs(lines) <- 4326 # And I also tried like this... 

mapview::mapview(lines_sf)

I'm trying to plot it along with the tip of New Zealand's North Island:

nz_sf <- 
  rnaturalearth::ne_countries(scale = "medium", 
                              country = "new zealand", 
                              returnclass = "sf")

## After setting a 'crs' for 'lines_sf'
p <- 
  ggplot() + 
  geom_sf(data = nz_sf) + 
  coord_sf(xlim = c(172.5,173.8), ylim = c(-34.2,-34.9))+ # This just work if the line below is commented 
  geom_sf(data = lines_sf) +
  theme_bw()

Checking attributes(lines_sf$geometry) shows me that, apparently, it's all good with projection.

It is the first time I'm "building" spatial features inside R from scratch, and I might be doing something wrong. What am I missing?

Best Answer

You've got your latitude and longitude swapped. If I plot and label the axes I see you have latitudes greater than 90, which is a big clue:

> plot(lines_sf$geom)
> axis(1)
> axis(2)

enter image description here

st_linestring takes coords in X,Y order, but you are passing lat,long, which is Y,X. Flip those columns and you'll be sorted.

Something like:

lines <- lapply(rows, function(row) {
  lmat <- matrix(unlist(row[c(2,1,4,3)]), ncol = 2, byrow = TRUE)
  st_linestring(lmat) # I've tried crs = 4326 here
})

should get you there.

Related Question