R:如何在相等的 area/compromise 世界地图上投影全球二维场?
R: How to project a global 2D field on an equal area/compromise world map?
我愿意
- 采用具有
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)
两者如何结合?谢谢!
这个问题我遇到过很多次了。我最终编写了一个有助于从矩阵创建多边形的函数。然后使用相同的投影设置将它们投影到地图上。这里唯一没有很好完成的是网格线的标记——我已经删除了它们,因为它们变得很乱:
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()
我愿意
- 采用具有
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)
两者如何结合?谢谢!
这个问题我遇到过很多次了。我最终编写了一个有助于从矩阵创建多边形的函数。然后使用相同的投影设置将它们投影到地图上。这里唯一没有很好完成的是网格线的标记——我已经删除了它们,因为它们变得很乱:
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()