projectRaster : 태평양 지역의 marmap을 사용하여 태평양 세계지도에서 수심 측정 데이터 ("NOAA.nc")의 래스터 투영

Question

태평양 지역의 경우 marmap 짤막음에서 마지막 예제를 재현하고 싶습니다 : 'marmap-DataAnalysis'. 예는 경도 여기 = (50)을 중심으로 전 세계의 직교 투영의 예이다 보여준다. I는 태평양 자오선, 즉 경도 = 155.5에 중심을 변경하려는 그러나

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) 

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 인수와 래스터 패키지에 의한 투영 계산을 허용하는 몇 가지 속임수를 사용해야합니다.

첫 번째 트릭은 안티 메리 디안이 2 개 데이터 세트를 다운로드하기 때문에 lon1 = lon2 = 0으로 데이터를 다운로드하는 것입니다. 안티 메리 디안에서 lon1로, 그리고 lon2에서 antimeridian으로 다운로드합니다. lon1과 lon2를 0으로 설정하면 전 세계를 다운로드합니다.

그런 다음, (getNOAA.bathy()animeridian의 인수에 의해 제조 된 (360)가 아닌 0) -180 수동 경도 180 사이의 값으로 전환하는 따라서 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를 사용하여 솔루션을 가장해야한다 .