如何在R中正确投影和绘制栅格数据。

4

我有一个等面积Behrmann投影的光栅图像,想将其投影到Mollweide投影并绘制。

然而,当我使用以下代码进行操作时,绘图似乎不正确,因为地图延伸到了两侧,并且有各种陆地的轮廓,这些轮廓是我不希望出现的。此外,地图超出了绘图窗口。

有谁能帮我把它漂亮地绘制出来吗?

谢谢!

可以从此链接下载使用的数据文件。

这是我目前拥有的代码:

require(rgdal)
require(maptools)
require(raster)

data(wrld_simpl)

mollCRS <- CRS('+proj=moll')
behrmannCRS <- CRS('+proj=cea +lat_ts=30')

sst <- raster("~/Dropbox/Public/sst.tif", crs=behrmannCRS)

sst_moll <- projectRaster(sst, crs=mollCRS)
wrld <- spTransform(wrld_simpl, mollCRS)

plot(sst_moll)
plot(wrld, add=TRUE)

enter image description here


请查看帖子,我认为它涉及到同样的问题。 - koekenbakker
谢谢!虽然所描述的问题听起来相同,但是建议的修复方法在我的情况下并没有改变任何东西。 - Pascal
1个回答

6

好的,由于在页面的示例似乎有效,我尽可能地模仿它。 我认为问题出在栅格图像的最左侧和最右侧重叠的部分。 裁剪并进行中间的Lat-Lon投影,就像示例中一样,似乎可以解决您的问题。

也许这个解决方法可以成为一个更优雅的解决方案的基础,直接解决问题,因为两次重新投影栅格并不利于性能优化。

# packages
library(rgdal)
library(maptools)
library(raster)

# define projections
mollCRS <- CRS('+proj=moll')
behrmannCRS <- CRS('+proj=cea +lat_ts=30')

# read data
data(wrld_simpl)
sst <- raster("~/Downloads/sst.tif", crs=behrmannCRS)

# crop sst to extent of world to avoid overlap on the seam
world_ext = projectExtent(wrld_simpl, crs = behrmannCRS)
sst_crop = crop(x = sst, y=world_ext, snap='in')

# convert sst to longlat (similar to test file)
# somehow this gets rid of the unwanted pixels outside the ellipse
sst_longlat = projectRaster(sst_crop, crs = ('+proj=longlat'))

# then convert to mollweide
sst_moll <- projectRaster(sst_longlat, crs=mollCRS, over=T)
wrld <- spTransform(wrld_simpl, mollCRS)

# plot results
plot(sst_moll)
plot(wrld, add=TRUE)

enter image description here


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接