3

我想重现小插图中的最后一个示例marmap：太平洋地区的“marmap-DataAnalysis”。该示例显示了以 lon = 50 为中心的世界的正交投影。示例如下：

library(marmap)
library(raster)
# Get data for the whole world. Careful: ca. 21 Mo!
world <- getNOAA.bathy(-180, 180, -90, 90, res = 15, keep = TRUE)

# Switch to raster
world.ras <- marmap::as.raster(world)

# Set the projection and project
my.proj <-   "+proj=ortho +lat_0=0 +lon_0=50 +x_0=0 +y_0=0"
world.ras.proj <- projectRaster(world.ras,crs = my.proj)

# Switch back to a bathy object
world.proj <- as.bathy(world.ras.proj)

# Set colors for oceans and land masses
blues <- c("lightsteelblue4", "lightsteelblue3",
           "lightsteelblue2", "lightsteelblue1")
greys <- c(grey(0.6), grey(0.93), grey(0.99))

# And plot!
plot(world.proj, image = TRUE, land = TRUE, lwd = 0.05,
     bpal = list(c(0, max(world.proj, na.rm = T), greys),
                 c(min(world.proj, na.rm = T), 0, blues)),
     axes = FALSE, xlab = "", ylab = "")
plot(world.proj, n = 1, lwd = 0.4, add = TRUE)

但是，我想将中央更改为太平洋子午线，例如 lon = 155.5。我尝试通过将投影参数更改为，

my.proj <- "+proj=ortho +lat_0=20 +lon_0=155.5 +x_0=0 +y_0=0"

但是之后，

world.ras.proj <- projectRaster(world.ras,crs = my.proj)

结果是：

Error in if (nr != x@nrows | nc != x@ncols) { : 
  missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1],     xy[,  :
  259 projected point(s) not finite
2: In rgdal::rawTransform(projection(raster), crs, nrow(xy), xy[,     1],  :
  4 projected point(s) not finite

我怎么能在太平洋地区绘制“测深世界”？

4

2 回答 2

2

这可以通过包marmap的当前/以前版本来解决raster。您必须使用函数的antimeridian=TRUE参数getNOAA.bathy()和一些技巧来允许通过 raster 包计算投影。

第一个技巧是下载数据，lon1 = lon2 = 0因为反子午线下载 2 个不同的数据集：从反子午线到 lon1 和从 lon2 到反子午线。将 lon1 和 lon2 设置为 0 会下载整个世界。

然后，您必须手动切换回 -180 到 180 之间的经度值（而不是 0 到 360 的animeridian参数产生的getNOAA.bathy()），因此该rownames(world2) <- ...行。

最后，您必须应用相同的-180校正来指定投影。这是代码：

library(marmap)
library(raster)
# Get data for the whole world. Careful: ca. 21 Mo!
world2 <- getNOAA.bathy(0, 0, -90, 90, res = 15, keep = TRUE, antimeridian=TRUE)
rownames(world2) <- as.numeric(rownames(world2))-180

# Switch to raster
world.ras <- marmap::as.raster(world2)

# Set the projection and project
my.proj <-   "+proj=ortho +lat_0=20 +lon_0=155-180 +x_0=0 +y_0=0"
world.ras.proj <- projectRaster(world.ras,crs = my.proj)

# Switch back to a bathy object
world.proj <- as.bathy(world.ras.proj)

# Set colors for oceans and land masses
blues <- c("lightsteelblue4", "lightsteelblue3",
           "lightsteelblue2", "lightsteelblue1")
greys <- c(grey(0.6), grey(0.93), grey(0.99))

# And plot!
plot(world.proj, image = TRUE, land = TRUE, lwd = 0.05,
     bpal = list(c(0, max(world.proj, na.rm = T), greys),
                 c(min(world.proj, na.rm = T), 0, blues)),
     axes = FALSE, xlab = "", ylab = "")
plot(world.proj, n = 1, lwd = 0.4, add = TRUE)

结果如下：

在此处输入图像描述

(2.6-7)的新版本raster解决了跨越日期线的投影问题。但是，由于从 NOAA 服务器下载测深数据时的舍入误差，一些缺失的单元格可能会出现在图中。这是您在原始问题中发布的代码示例：

在此处输入图像描述

这是summary()数据：

summary(world)
# Bathymetric data of class 'bathy', with 1440 rows and 720 columns
# Latitudinal range: -89.88 to 89.88 (89.88 S to 89.88 N)
# Longitudinal range: -179.88 to 179.88 (179.88 W to 179.88 E)
# Cell size: 15 minute(s)

# Depth statistics:
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#  -10635   -4286   -2455   -1892     214    6798 
# 
# First 5 columns and rows of the bathymetric matrix:
#          -89.875 -89.625 -89.375 -89.125 -88.875
# -179.875    2746    2836    2893    2959    3016
# -179.625    2746    2835    2892    2958    3015
# -179.375    2746    2835    2891    2957    3014
# -179.125    2746    2834    2890    2956    3013
# -178.875    2746    2834    2889    2955    3012

因此，使用antimeridian=TRUE上面详述的解决方案应该是最好的。

于 2017-11-21T12:49:38.510 回答
2

我已经简化了你的问题（总是很好，对我来说数据下载不起作用）。在本质上：

library(raster); library(rgdal)
prj1 <- "+proj=ortho +lat_0=0 +lon_0=0 +x_0=0 +y_0=0"
prj2 <- "+proj=ortho +lat_0=20 +lon_0=155.5 +x_0=0 +y_0=0" 
r <- raster()
r <- init(r, 'col')

# works
x1 <- projectRaster(r, crs = prj1)
# fails
x2 <- projectRaster(r, crs = prj2)

这是一个错误。我已经在光栅版本 2.6-2 中修复了它（正在开发中，应该在下周左右可用）

于 2017-11-02T01:37:00.843 回答