0

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?

1
  • Tip: avoid library(tidyverse) - it pulls in so many packages you don't need here. You only need ggplot2 and tribble, and you don't really need tribble for making examples. The lighter you can keep your examples the easier they are for people to run.
    – Spacedman
    Commented Nov 15, 2021 at 6:28

1 Answer 1

0

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.

1
  • Ouch - that was a noobie error! Many thanks, @Spacedman - it worked perfectly now. Thanks also for the comment above. I'll pay more attention to the packages needed next time :)
    – nwd
    Commented Nov 15, 2021 at 9:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.