projectRaster：太平洋のマーマップを持つ太平洋 - ワールドマップの地形データ（ "NOAA.nc"）のラスタ投影

Question

太平洋地域の場合、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に中心を変更したいIによりこれを試みました。投影パラメータをに変更する。

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 

がどのように私は太平洋地域における '海底地形の世界' をプロットすることができことができますか？

Answer 1

私はあなたの質問を簡略化しました（いつもうまくいって、私にとってはデータのダウンロードがうまくいかなかった）。本質的に：

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で修正しました（開発中です。来週利用可能になるはずです）

Answer 2

これは、rasterパッケージの現在/以前のバージョンでmarmapに解決できます。 getNOAA.bathy()関数の引数antimeridian=TRUEを使用し、ラスタパッケージによる投影の計算を許可するにはトリッキーを使用する必要があります。

アンチメリディアンがantimeridianからlon1に、そしてlon2からantimeridianに2つの異なるデータセットをダウンロードするので、最初のトリックはlon1 = lon2 = 0でデータをダウンロードすることです。 lon1とlon2を0に設定すると、世界全体がダウンロードされます。

次に、-180〜180の間の値の値に手動で戻す必要があります（animeridianの引数である0〜360ではなく、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) 

そして、ここではその結果である：ここで はコードである

がrasterの新バージョン（2.6-7）日付ライン全体に投影する問題を解決します。ただし、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を使用したソリューションは、最高でなければなりません。