R：如何在相等的 area/compromise 世界地图上投影全球二维场？

Question

我愿意

采用具有 lon/lat 个值（例如地表气温）的 global data 矩阵，缺失值给出为 NAs。
在等面积或折衷（例如罗宾逊）投影中绘制它（热带和极地地区都应表示）
在顶部添加大陆轮廓
在侧面添加网格和标签

编辑：我正在寻找的投影看起来像这样：

使用 lat/lon 值获得顶部有地图的矩形图像的基本图很简单，但是一旦我尝试投影矩阵数据，事情就会变得复杂得多。它不可能真的这么复杂——我一定是遗漏了什么。试图理解 openproj 、 spTransform 和 SpatialPoints 给我的印象是我必须转换我的矩阵坐标，以便我有一种网格对象。

下面的例子是基于对 R - Plotting netcdf climate data and uses the NOAA air surface temperature data from http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.html 的回答。

编辑：似乎（至少）有三种方法可以解决这个问题，一种使用 mapproj（并将矩阵转换为多边形，如下面的回复中指出），另一种使用 openproj、CRS 和 spTransform', and maybe a third one using 光栅`。

library(ncdf)
library(fields)
temp.nc <- open.ncdf("~/Downloads/air.1999.nc")
temp <- get.var.ncdf(temp.nc,"air")
temp[temp=="32767"] <- NA # set missing values to NA
temp.nc$dim$lon$vals -> lon
temp.nc$dim$lat$vals -> lat
temp.nc$dim$time$vals -> time
temp.nc$dim$level$vals -> lev
lat <- rev(lat)
temp <- temp[, ncol(temp):1, , ] #lat being our dimension number 2
lon <- lon -180
temp11 <- temp[ , , 1, 1] #Level is the third dimension and time the fourth.

所以数据是lon x lat x temp11并且可以绘制成图像

image.plot(lon,lat,temp11-273.15,horizontal=TRUE) 
library(maptools)
map("world",add=TRUE,wrap=TRUE)

到目前为止一切都很好并且可以正常工作。但我希望将地图和图像都投影到更相关的网格上。所以接下来是我尝试找到一种方法来做到这一点。

尝试使用 mapproj 投影地图（失败）

library("mapproj")
eg<-expand.grid(lon,lat)
xyProj<-mapproject(eg[,1],eg[,2],projection="vandergrinten")
im<-as.image(temp11,x=xyProj$x,y=xyProj$y) # fails
image.plot(im)

而这个项目至少勾勒出了大陆的轮廓（但有些线条被扭曲了）

map("world",proj="vandergrinten",fill=FALSE)

两者如何结合？谢谢！

Answer 1

这个问题我遇到过很多次了。我最终编写了一个有助于从矩阵创建多边形的函数。然后使用相同的投影设置将它们投影到地图上。这里唯一没有很好完成的是网格线的标记——我已经删除了它们，因为它们变得很乱：

library(ncdf4)
library(fields)
temp.nc <- nc_open("air.1999.nc")
lon <- ncvar_get(temp.nc, "lon")
lat <- ncvar_get(temp.nc, "lat")
time <- ncvar_get(temp.nc, "time")
lev <- ncvar_get(temp.nc, "level")
temp <- ncvar_get(temp.nc, "air")
nc_close(temp.nc)

temp[temp=="32767"] <- NA # set missing values to NA
lat <- rev(lat)
temp <- temp[, ncol(temp):1, , ] #lat being our dimension number 2
lon <- lon -180
temp11 <- temp[ , , 1, 1] #Level is the third dimension and time the fourth.

# plot
library("mapproj")
image(lon, lat, temp11); map("world",add=TRUE,wrap=TRUE)

library(sinkr) # https://github.com/marchtaylor/sinkr
polys <- matrixPoly(lon, lat, temp11)

png("vandergrinten.png", width=6, height=6, res=400, units="in")
op <- par(mar=c(1,1,1,1))
map("world", 
  proj="vandergrinten", parameters = NULL, orient=c(90,0,0),
  xlim=c(-180,180), ylim=c(-90,90),
  fill=FALSE, wrap=TRUE
)
COLS <- val2col(temp11, col = jetPal(100))
for(i in seq(polys)){
  tmp <- mapproject(polys[[i]],
    proj="vandergrinten", parameters = NULL, orient=c(90,0,0)
  )
  polygon(tmp$x, tmp$y, col=COLS[i], border=COLS[i], lwd=0.1)
}
map("world",
  proj="vandergrinten", parameters = NULL, orient=c(90,0,0),
  add=TRUE
)
map.grid(nx = 18, ny=9, col=8, labels = FALSE)
par(op)
dev.off()

更新：

定义其他限制有点棘手，但这里有一个示例（您会看到我必须调整设备的尺寸以便绘制所有内容。）。重要提示：asp=1 应设置为 plot 以防止失真：

xlim=c(-180,180)
ylim=c(-80,80)
LIM <- mapproject(x=xlim, y=ylim, proj="vandergrinten", orient=c(90,0,0))
 
png("vandergrinten2.png", width=6, height=4.5, res=400, units="in")
op <- par(mar=c(1,1,1,1))
plot(1, t="n", asp=1 ,axes=FALSE, xlab="", ylab="", xlim=range(LIM$x), ylim=LIM$y)
map("world", 
  proj="vandergrinten", orient=c(90,0,0),
  fill=FALSE, wrap=TRUE, add=TRUE
)
COLS <- val2col(temp11, col = jetPal(100))
for(i in seq(polys)){
  tmp <- mapproject(polys[[i]],
    proj="vandergrinten", parameters = NULL, orient=c(90,0,0)
  )
  polygon(tmp, col=COLS[i], border=COLS[i], lwd=0.1)
}
map("world", proj="vandergrinten", orient=c(90,0,0), add=TRUE)
map.grid(nx = 18, ny=9, col=8, labels = FALSE)
par(op)
dev.off()

R：如何在相等的 area/compromise 世界地图上投影全球二维场？

R: How to project a global 2D field on an equal area/compromise world map?

maps

plot

r

geospatial

netcdf

更新：