2

对于一个 Shiny uni 项目，我面临着折线穿越太平洋日期变更线的问题。从一个中心（武汉）我想创建 72 条线，显示它们连接不同区域。我已经创建了一个循环，以确保经度 > 0 已更改为大于 0 的经度。有没有人有解决方案使折线正确穿过日期线？

在照片中，您可以看到我当前的情节没有正确越界。

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

输出闪亮

4

2 回答 2

2

您可以尝试几件事。一种是使用以下breakAtDateLine = TRUE选项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
  )

如您所见，它打破了这条线，但在屏幕左侧继续它，这并不理想。

在此处输入图像描述

我们可以尝试的另一件事是运行 gcIntermediate 函数，并在运行该函数后breakAtDateLine = FALSE手动调整纬度和经度。如果我们设置，这将是最简单的。sp=FALSE

请注意，我们只需要更正从我们所在位置向东的线 - 这些是唯一穿过日期线的线。每个位置的情况都不相同，但逻辑应该相似。

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)

在此处输入图像描述

当它到达地图顶部时仍然有点奇怪，但希望这是你可以忍受的东西

于 2020-02-25T12:11:46.817 回答
0

我最近在仪表板上工作并且遇到了同样的问题。我得到的解决方案与上面克里斯的评论有点混合。差异：

  1. 不是创建 p1 和 p2，而是直接在 gcIntermediate 中调用 df 列

  2. 保持sp = T和breakAtDateLine = T

  3. 对生成的 sp 对象执行 for 循环，不修改 df 并使用 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)

输出（类似于 Chris 的） 循环利用了这样一个事实，即gcIntermediate()通常只为不跨越 ante-meridiem 的特征创建一条线。但是当线确实交叉时，它将为该要素制作两条线，在负坐标中设置的线始终在列表中排在第二位。因此，循环查找具有 2 个特征的行，取第二个特征的长并添加 +360，将第二行坐标 rbind 到第一行，并将第二行归零。

这个解决方案做的事情几乎相同，而且速度并不快（我针对 Chris 的第二个答案 1000 倍做了一个基准测试，这个解决方案平均只快 0.3%）。但是，它确实减少了代码、对象混乱和通过转换处理数据。

于 2020-09-08T21:41:40.073 回答