【问题标题】:Plotting country borders with maptools - R用 maptools 绘制国家边界 - R
【发布时间】:2013-05-16 04:45:36
【问题描述】:

根据 Joe Wheatley (http://joewheatley.net/ncep-global-forecast-system/) 的精彩帖子,我设法制作了一个 温度全球地图。但是,我尝试使用 maptools 包来绘制国家边界,而不是只绘制海岸线。 当只绘制东半球国家边界时,问题就来了。我应该错过一些我无法弄清楚的东西, 仍在寻找 stackoverflow 和谷歌。希望你能帮忙。

这是我正在使用的代码(大部分来自 Joe 的帖子)

loc=file.path("ftp://ftp.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/gfs.2013052100/gfs.t00z.sfluxgrbf03.grib2")
download.file(loc,"temp.grb",mode="wb")

system("wgrib2 -s temp.grb | grep :TMP: | wgrib2 -i temp.grb -netcdf TMP.nc",intern=T)
system("wgrib2 -s temp.grb | grep :LAND: | wgrib2 -i temp.grb -netcdf temp.nc",intern=T)

library(ncdf)
landFrac <-open.ncdf("LAND.nc")
lon <- get.var.ncdf(landFrac,"longitude")
lat <- get.var.ncdf(landFrac,"latitude")

temp=open.ncdf("TMP.nc")
t2m.mean <- get.var.ncdf(temp,"TMP_2maboveground")


library("fields")
library("sp", lib.loc="/usr/lib/R/site-library")
library("maptools", lib.loc="/usr/lib/R/site-library")

day="DIA"

png(filename="gfs.png",width=1215,height=607,bg="white")

rgb.palette <- colorRampPalette(c("snow1","snow2","snow3","seagreen","orange","firebrick"), space = "rgb")#colors
image.plot(lon,lat,t2m.mean,col=rgb.palette(200),main=as.expression(paste("GFS 24hr Average 2M Temperature",day,"00 UTC",sep="")),axes=T,legend.lab="o C")
data(wrld_simpl)
plot(wrld_simpl, add = TRUE)

dev.off()

这是生成的图像

这是一张全球地图,我应该在 image.plot 中使用 xlim 和 ylim 来提取区域(即欧洲)

编辑:为 temp.nc 文件添加 url

http://ubuntuone.com/29DKAeRjUCiCzLblgfSLc9

任何帮助将不胜感激,谢谢

【问题讨论】:

    标签: r netcdf maptools grib


    【解决方案1】:

    数据集wrld_simpl 在经度 [-180, 180] 上对齐,但您的网格数据显然是 [0,360]。由于您已加载地图工具,请尝试将修改后的副本添加到您的绘图中:

    plot(elide(wrld_simpl, shift = c(360, 0)), add = TRUE)
    

    还有其他工具可以像这样裁剪和移动/合并和重新定位数据,包括 ?rotate 在包栅格中。这真的取决于你需要什么。

    另一种图形选择是在地图包中使用“world2”

    library(maps)
    map("world2", add = TRUE)
    

    其中任何一个都可以完成上面开始的情节。

    这是另一个与此相关的讨论:Fixing maps library data for Pacific centred (0°-360° longitude) display

    【讨论】:

    • 优秀的答案,像往常一样。
    • 嗨,plot(elide(wrld_simpl, shift = c(360, 0)), add = TRUE) 可以替代 plot(wrld_simpl, add = TRUE) 吗?谢谢。
    • 几乎完美的解决方案。但是现在,边界被绘制在地图之外,就在图例的右侧。欧洲边界不仅出现在地图的左侧,也出现在地图的右侧。
    • 它远非完美,它有几个解决方法。我会在 [0,360] 范围内得到一张不错的地图 - 例如,rgdal 或 maptools 中有用于读取 shapefile 的工具,或者将栅格重新配置为 [-180, 180]。改用 raster 包将是一个好的开始,它简化了很多。如果您可以发布您的 NetCDF 文件,我会详细说明,但我无权访问 wgrib2。
    • 您可以在ubuntuone.com/29DKAeRjUCiCzLblgfSLc9 找到数据文件 temp.nc,还在我的问题中添加了 url。我会看看光栅包,再次感谢。
    【解决方案2】:

    你可以这样做:

    library(raster)
    r <- raster("temp.nc")
    r <- rotate(r)
    plot(r)
    

    【讨论】:

    • 嗨@RobertH,这正是我想要的。现在地图很好。非常感谢。
    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2017-12-19
    • 1970-01-01
    • 2023-01-30
    • 1970-01-01
    • 2019-08-13
    • 2014-05-15
    相关资源
    最近更新 更多