2

For a Shiny uni project I'm facing problems with polylines crossing the date line in the Pacific Ocean. From one center (Wuhan) I would like to create 72 lines showing their connection different regions. I have already created a loop to make sure that Longitudes > 0 have been changed to a longitude larger than 0. Does anyone have a solution to make the polylines properly cross the date line?

In the photo you can see my current plot is not properly crossing the line.

library(leaflet)
library(geosphere)
library(sp)

df <- read.table(text = "
Longitude.Central Latitude.Central Longitude Latitude
112.2707         30.97564  117.2264 31.82571
112.2707         30.97564  116.4142 40.18238
112.2707         30.97564    4.4699 50.50390
112.2707         30.97564  -71.0589 42.36010
112.2707         30.97564 -123.1210 49.28270
112.2707         30.97564  104.9910 12.56570",
           header = TRUE)

p1 <- as.matrix(df[,c(1,2)])
p2 <- data.frame(df[,c(3,4)])

p3 <- matrix(nrow = length(p2))
for (j in 1:nrow(p2)) {
  if (p2[j,]$Longitude < 0) {
    p3[j] <- p2[j,]$Longitude + 360
  } else {
    p3[j] <- p2[j,]$Longitude
  }
}

p2 <- as.matrix(cbind(p3, p2$Latitude))
df2 <- gcIntermediate(
  p1, p2, 
  n = 72, 
  addStartEnd = TRUE,
  sp = T)

leaflet( ) %>% 
    setView(df$Longitude.Central[[1]], lat = df$Latitude.Central[[1]], 1) %>%
    addTiles(urlTemplate = "https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png") %>% 
    addPolylines(data = df2, weight = 1
    )

# Head of the data
> head(df)
# A tibble: 6 x 4
  Longitude.Central Latitude.Central Longitude Latitude
              <dbl>            <dbl>     <dbl>    <dbl>
1          112.2707         30.97564  117.2264 31.82571
2          112.2707         30.97564  116.4142 40.18238
3          112.2707         30.97564    4.4699 50.50390
4          112.2707         30.97564  -71.0589 42.36010
5          112.2707         30.97564 -123.1210 49.28270
6          112.2707         30.97564  104.9910 12.56570

Output Shiny

Chris
  • 3,836
  • 1
  • 16
  • 34
  • 1
    Hi @Chris, Thanks for responding. I have adjusted the code and added the head of my data, which should clarify the data. Is the current information sufficient for answering my question? – Onno van der Horst Feb 25 '20 at 08:10

2 Answers2

2

There are a couple of things you could try. One is to use the breakAtDateLine = TRUE option in gcIntermediate:

df2 <- gcIntermediate(
  p1, p2, 
  n = 72, 
  addStartEnd = TRUE,
  sp = T,
  breakAtDateLine = T)

leaflet( ) %>% 
  setView(lng = df$Longitude.Central[[1]], lat = df$Latitude.Central[[1]], 1) %>%
  addTiles(urlTemplate = "https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png") %>% 
  addPolylines(data = df2, weight = 1
  )

As you can see, it breaks the line, but continues it on the left side of the screen which isn't ideal.

enter image description here

Another thing we could try is to run the gcIntermediate function with breakAtDateLine = FALSE and manually adjust the latitude and longitude after we run the function. This will be easiest if we set sp=FALSE.

Notice how we only need to correct the lines that head east from our location - these are the only ones that cross the date line. This will not be the same for every location but logic should be similar.

p1 <- as.matrix(df[,c(1,2)])
p2 <- data.frame(df[,c(3,4)])

df2 <- gcIntermediate(
  as.matrix(p1),
  as.matrix(p2), 
  n = 72, 
  addStartEnd = TRUE,
  breakAtDateLine = F,
  sp = F)

# correct the longitudes
res <- lapply(df2, function(x) {
  # if direction is east, correct lon, else don't
  if(x[,'lon'][2] - x[,'lon'][1]  > 0){
    cbind(lon =ifelse(x[,'lon']<0, x[,'lon'] + 360, x[,'lon']), lat = x[,'lat'])
  } else {
    x
  }
})

# and convert back to sp (I just took this code from the gcIntermediate function)
for (i in 1:length(res)) {
  if (!is.list(res[[i]])) {
    res[[i]] <- Lines(list(Line(res[[i]])), ID = as.character(i))
  }
  else {
    res[[i]] <- Lines(list(Line(res[[i]][[1]]), Line(res[[i]][[2]])), 
                      ID = as.character(i))
  }
}

res <- SpatialLines(res, CRS("+proj=longlat +ellps=WGS84"))


leaflet() %>% 
  setView(df$Latitude.Central[[1]], lat = df$Longitude.Central[[1]], 1) %>%
  addTiles(urlTemplate = "https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png") %>% 
  addPolylines(data = res, weight = 1) 

enter image description here

It still goes a little strange when it hits the top of the map but hopefully that's something you can live with

Chris
  • 3,836
  • 1
  • 16
  • 34
  • 1
    Hi Chris, Thanks for the help. I getting an error while running this: Error in x[, "lon"] : no 'dimnames' attribute for array Have you got any idea what I'm doing wrong? – Onno van der Horst Feb 25 '20 at 12:52
  • 1
    Copy and paste the example code into a new session to confirm it works, then see what is different between the example and your actual data, it should definitely work on the example above. – Chris Feb 25 '20 at 13:05
  • 1
    Shit, I left Wuhan itself (and its coordinates) in the data causing the error. Thanks for the help!! – Onno van der Horst Feb 25 '20 at 13:40
0

I was working on a dashboard recently and had the same issue. The solution I came to is a bit of a hybrid to Chris's comment above. Differences:

  1. Instead of creating p1 and p2, the df columns are called directly in gcIntermediate

  2. keep sp = T and breakAtDateLine = T

  3. for loop is performed on the resulting sp object, no modifications to a df and converting with SpatialLines

    df2 <- gcIntermediate(
                       df[,c(1,2)],
                       df[,c(3,4)], 
                       n = 72, 
                       addStartEnd = TRUE,
                       sp = T, 
                       breakAtDateLine = T)
    
    for (l in 1:length(df2@lines)){  
       if (length(df2@lines[[l]]@Lines) > 1){
         df2@lines[[l]]@Lines[[2]]@coords[ ,1] <- df2@lines[[l]]@Lines[[2]]@coords[ ,1] +360
    
         df2@lines[[l]]@Lines[[1]]@coords  <- rbind(df2@lines[[l]]@Lines[[1]]@coords, 
                                            df2@lines[[l]]@Lines[[2]]@coords)
    
         df2@lines[[l]]@Lines[[2]] <- NULL
       }
     }  
    
    leaflet() %>% 
      setView(df$Latitude.Central[[1]], lat = df$Longitude.Central[[1]], 1) %>%
      addTiles(urlTemplate = "https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png") %>% 
      addPolylines(data = df2, weight = 1) 
    

output (similar to Chris's) The loop is taking advantage of the fact that gcIntermediate() normally creates only one line for features that don't cross the ante-meridiem. But when lines do cross, it will make two lines for the feature, with the line set in the negative coordinates always being second in the list. So, the loop looks for lines that have 2 features, takes the long of the second feature and adds +360, rbinds the second lines coords to the first line and Nulls the second line.

This solution does pretty much the same thing and isn't much faster (I did a benchmark running this against Chris's second answer 1000x, with this solution only being 0.3% faster on average). However, it does cut down on code, object clutter, and handing the data through conversions.

dla06c
  • 1