时间戳

首页 > 解决方案 > R中的Mapview和太平洋地区的预测/ crs问题

r - R中的Mapview和太平洋地区的预测/ crs问题

问题描述

我正在尝试使用 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 数据

谢谢

地图查看图片 太平洋

标签： rr-mapview

解决方案

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

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

推荐阅读