R:如何在等面积/折衷世界地图上投影全局二维场？

Question

我愿意

• 获取具有lon/lat值(例如地表气温)的全局数据矩阵，缺失值以NA形式给出。
• 在等面积或折衷(例如罗宾逊)投影中绘制它(热带和极地地区都应表示)
• 在顶部添加大陆轮廓
• 在侧面添加网格和标签

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

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

以下示例基于对 R - Plotting netcdf climate data 的回答并使用来自 http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.html 的 NOAA 空气表面温度数据.

编辑:似乎(至少)有三种方法可以解决这个问题，一种使用 mapproj(并将矩阵转换为多边形，如以下响应中所指出)，另一种使用 openproj、CRS 和 spTransform'，也许第三个使用 raster`。

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)

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

尝试使用 mapproj 投影 map (失败)

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

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

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