1

我正在尝试使用 mapView 在 R 中绘制跨越经度 180 度的太平洋网格的数据。原生 crs 是“+proj=merc +lon_0=150 +lat_ts=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0”。等效数据的示例是：

test.coords<-as.data.frame(cbind(c(runif(15,-180,-130),runif(5,160,180)),runif(20,40,60)))
test.sp <- SpatialPointsDataFrame(coords = cbind(test.coords$V1,test.coords$V2), data = test.coords,
                                 proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
test.sp <- spTransform(test.sp, 
          "+proj=merc +lon_0=150 +lat_ts=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
plot(test.sp)
mapview(test.sp)

使用绘图功能时，点居中，而使用 mapView 时，它们会在地图的两侧（链接图像）之间分割。使用 mapView 的视图可以以不同的经度为中心吗？

对此地图使用本机 crs 无济于事。我得到一个错误。

mapview(plot, native.crs=TRUE)

警告信息：sf 层不是 long-lat 数据

谢谢

地图查看图片 太平洋

4

2 回答 2

1

这是多边形的解决方法-它会产生一些警告，但总体上可以。

# transform grid to standard latitude and longitude
# original grid in crs "+proj=merc +lon_0=150 +lat_ts=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
grid<-spTransform(grid,"+init=epsg:4326")
summary(coordinates(grid))

#       V1               V2       
# Min.   :-179.8   Min.   :37.17  
# 1st Qu.:-168.8   1st Qu.:51.70  
# Median :-154.3   Median :55.02  
# Mean   :-118.0   Mean   :54.65  
# 3rd Qu.:-131.4   3rd Qu.:58.31  
# Max.   : 179.8   Max.   :65.56  

out <- lapply(1:length(grid), function(i) grid[i, ])

for (i in 1:length(grid)) {

  cds <- slot(slot(slot(out[[i]], "polygons")[[1]], "Polygons")[[1]], "coords")

  cds_polygon <- Polygons(list(Polygon(cds)), ID="out.x")
  out.x <- SpatialPolygons(list(cds_polygon))
  proj4string(out.x) <- CRS("+init=epsg:4326")

  if (cds[1,1]>0) { # depending to which side one want to move polygons that are on both sides of International Date Line
    out.y <- elide(out.x, shift=c(-360,0))
  } else {
    out.y<-out.x}

  if(i==1){grid.shifted<-out.y} else {
    grid.shifted<-raster::union(grid.shifted,out.y)
  }
}

# add data in the original grid
grid.shifted <- SpatialPolygonsDataFrame(grid.shifted, grid@data, match.ID = F)

mapview(grid.shifted)
于 2020-03-13T23:42:09.810 回答
1

好的，因此这是您问题的半最佳解决方案，原因如下：

  1. 这仅适用于点数据
  2. 对于大型点集合，它可能无法很好地扩展
  3. 投影到另一个 CRS 可能不再起作用（尽管这并不真正相关，因为手头的问题是关于可视化）

对于一组线或多边形，我们需要另一种方法，但我目前不知道 sp/sf 原生解决方案来解决对象的往返坐标。

library(sp)
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(mapview)

test.coords<-as.data.frame(cbind(c(runif(15,-180,-130),runif(5,160,180)),runif(20,40,60)))
test.sp <- SpatialPointsDataFrame(coords = cbind(test.coords$V1,test.coords$V2), data = test.coords,
                                  proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

test_sf = st_as_sf(test.sp)

shift = function(x) {
    geom = st_geometry(x)
    st_geometry(x) = st_sfc(
        lapply(seq_along(geom), function(i) {
            geom[[i]][1] = ifelse(geom[[i]][1] < 0, geom[[i]][1] + 360, geom[[i]][1])
            return(geom[[i]])
        })
        , crs = st_crs(geom)
    )
    return(x)
}

mapview(shift(test_sf)) +
    mapview(test_sf, col.regions = "orange")

由reprex 包（v0.3.0）于 2019 年 12 月 11 日创建

我选择使用 sf 而不是 sp ，因为mapview在内部使用 sf 并且希望找到一个可以集成的通用方法来解决这个问题。

于 2019-12-11T18:10:06.103 回答