我在等面积贝尔曼投影中有一个栅格,我想将其投影到莫尔韦德投影和绘图上。
然而,当我使用以下代码执行此操作时,绘图似乎不正确,因为地图延伸到两侧,并且存在我意想不到的各种陆地的轮廓。此外,地图延伸到绘图之外窗户。
谁能帮我把这个情节画好吗?
Thanks!
使用的数据文件可以从以下位置下载这个链接 https://dl.dropboxusercontent.com/u/34644229/sst.tif.
这是我到目前为止的代码:
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)
好吧,因为这个例子在this http://r-sig-geo.2731867.n2.nabble.com/Difficulties-with-Eckert-IV-Mollweide-projection-of-RasterLayer-with-global-extent-td7580864.html页面似乎有效,我尝试尽可能地模仿它。我认为出现问题是因为光栅图像的最左侧和最右侧重叠。如示例中所示,裁剪和中间重新投影到经纬度似乎可以解决您的问题。
也许这种解决方法可以成为直接解决问题的更优雅解决方案的基础,因为重新投影栅格两次并没有什么好处。
# 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)
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)