R中的Mapview和太平洋地区的投影/crs问题

提问者：小点点

R中的Mapview和太平洋地区的投影/crs问题

我正在尝试使用mapView在R中绘制跨越经度180度的太平洋网格的数据。本机crs是"proj=merclon_0=150lat_ts=0x_0=0y_0=0 data=WGS84单元=mno_defsellps=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层不是经度数据

谢谢你。

mapView图片太平洋

共2个答案

匿名用户

好的，所以这是你问题的半最佳解决方案，原因有几个：

< li >这只适用于点数据 < li >它可能无法很好地适应大型点集合 < li >投影到另一个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，并希望找到一种通用的方法来解决这个问题，这种方法可能会被集成。

匿名用户

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

# 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)