【发布时间】：2017-09-22 14:23:33
【问题描述】：
多年来，我一直在使用ggplot2 绘制气候网格数据。这些通常是投影的 NetCDF 文件。单元格在模型坐标中是方形的，但根据模型使用的投影，在现实世界中可能并非如此。
我通常的做法是先将数据重新映射到合适的规则网格上，然后再进行绘图。这会对数据进行小幅修改，通常这是可以接受的。
但是，我认为这已经不够好了：我想直接绘制投影数据，而不需要重新映射，因为如果我没记错的话，其他程序（例如 ncl）可以这样做，而无需触摸模型输出值。
但是，我遇到了一些问题。我将在下面逐步详细说明可能的解决方案，从最简单到最复杂，以及它们的问题。我们能克服它们吗？
编辑：感谢@lbusett 的回答，我得到了包含解决方案的this nice function。如果喜欢请点赞@lbusett's answer！
初始设置
#Load packages
library(raster)
library(ggplot2)
#This gives you the starting data, 's'
load(url('https://files.fm/down.php?i=kew5pxw7&n=loadme.Rdata'))
#If you cannot download the data, maybe you can try to manually download it from http://s000.tinyupload.com/index.php?file_id=04134338934836605121
#Check the data projection, it's Lambert Conformal Conic
projection(s)
#The data (precipitation) has a 'model' grid (125x125, units are integers from 1 to 125)
#for each point a lat-lon value is also assigned
pr <- s[[1]]
lon <- s[[2]]
lat <- s[[3]]
#Lets get the data into data.frames
#Gridded in model units:
pr_df_basic <- as.data.frame(pr, xy=TRUE)
colnames(pr_df_basic) <- c('lon', 'lat', 'pr')
#Projected points:
pr_df <- data.frame(lat=lat[], lon=lon[], pr=pr[])
我们为每个模型单元创建了两个数据框，一个带有模型坐标，一个带有真实的经纬交叉点（中心）。
可选：使用较小的域
如果您想更清楚地看到单元格的形状，您可以对数据进行子集化并仅提取少量模型单元格。请注意，您可能需要调整点大小、绘图限制和其他便利设施。您可以像这样子集，然后重做上面的代码部分（减去load()）：
s <- crop(s, extent(c(100,120,30,50)))
如果您想完全理解问题，也许您想同时尝试大域和小域。代码是相同的，只是点大小和地图限制发生了变化。下面的值适用于大的完整域。 好的，现在让我们开始吧！
从瓷砖开始
最明显的解决方案是使用瓷砖。让我们试试吧。
my_theme <- theme_bw() + theme(panel.ontop=TRUE, panel.background=element_blank())
my_cols <- scale_color_distiller(palette='Spectral')
my_fill <- scale_fill_distiller(palette='Spectral')
#Really unprojected square plot:
ggplot(pr_df_basic, aes(y=lat, x=lon, fill=pr)) + geom_tile() + my_theme + my_fill
好的，现在更高级：我们使用真正的 LAT-LON，使用方形瓷砖
ggplot(pr_df, aes(y=lat, x=lon, fill=pr)) + geom_tile(width=1.2, height=1.2) +
borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_theme + my_fill +
coord_quickmap(xlim=range(pr_df$lon), ylim=range(pr_df$lat)) #the result is weird boxes...
好的，但那些不是真正的模型方块，这是一个 hack。此外，模型框在域的顶部发散，并且都以相同的方式定向。不是很好。让我们自己投影正方形，即使我们已经知道这不是正确的做法……也许看起来不错。
#This takes a while, maybe you can trust me with the result
ggplot(pr_df, aes(y=lat, x=lon, fill=pr)) + geom_tile(width=1.5, height=1.5) +
borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_theme + my_fill +
coord_map('lambert', lat0=30, lat1=65, xlim=c(-20, 39), ylim=c(19, 75))
首先，这需要很多时间。不能接受的。另外，再次重申：这些不是正确的模型单元。
尝试使用点，而不是图块
也许我们可以用圆点或方点代替瓦片，也可以投影它们！
#Basic 'unprojected' point plot
ggplot(pr_df, aes(y=lat, x=lon, color=pr)) + geom_point(size=2) +
borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_cols + my_theme +
coord_quickmap(xlim=range(pr_df$lon), ylim=range(pr_df$lat))
我们可以使用方点...和投影！我们越来越近了，尽管我们知道它仍然不正确。
#In the following plot pointsize, xlim and ylim were manually set. Setting the wrong values leads to bad results.
#Also the lambert projection values were tired and guessed from the model CRS
ggplot(pr_df, aes(y=lat, x=lon, color=pr)) +
geom_point(size=2, shape=15) +
borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_theme + my_cols +
coord_map('lambert', lat0=30, lat1=65, xlim=c(-20, 39), ylim=c(19, 75))
不错的结果，但不是完全自动的，而且绘图点不够好。我想要真实的模型细胞，它们的形状被投影改变了！
多边形，也许吧？
如您所见，我正在寻求一种正确绘制模型框的方法，以正确的形状和位置投影。当然，模型中为正方形的模型框，一旦投影出来，就不再是规则的形状了。所以也许我可以使用多边形并投影它们？
我尝试使用rasterToPolygons 和fortify 并关注this 帖子，但没有这样做。我试过这个：
pr2poly <- rasterToPolygons(pr)
#http://mazamascience.com/WorkingWithData/?p=1494
pr2poly@data$id <- rownames(pr2poly@data)
tmp <- fortify(pr2poly, region = "id")
tmp2 <- merge(tmp, pr2poly@data, by = "id")
ggplot(tmp2, aes(x=long, y=lat, group = group, fill=Total.precipitation.flux)) + geom_polygon() + my_fill
好的，让我们尝试用经纬度替换...
tmp2$long <- lon[]
tmp2$lat <- lat[]
#Mh, does not work! See below:
ggplot(tmp2, aes(x=long, y=lat, group = group, fill=Total.precipitation.flux)) + geom_polygon() + my_fill
（对不起，我改变了图中的色标）
嗯，甚至不值得尝试投影。也许我应该尝试计算模型单元角的纬度，并为此创建多边形，然后重新投影？
结论
- 我想在其原生网格上绘制投影模型数据，但我无法这样做。使用瓦片是不正确的，使用点是不正确的，并且使用多边形似乎由于未知原因而不起作用。
- 通过
coord_map()投影时，网格线和轴标签错误。这使得投影的 ggplots 无法用于出版物。
【问题讨论】：
-
你试过
ggmap吗？ -
@ulfelder 我记得过去使用过它，但最近没有。我忘记了它的存在。我会调查的。有什么特别要指出的吗？
-
我无法加载您的数据；你确定网址正确吗？
-
@SymbolixAU 我更改了下载提供程序，请重试。这应该更可靠地工作。
-
您的数据是从磁盘调用的：“/home/adriano/EURO-CORDEX_59_tas_pr_prc_monmean_1998-2000.nc” 所以您将一个没有真实数据但只有一个 url 的对象保存到您的机器中。您可以将 RData 保存在栅格图层中吗？
标签： r ggplot2 netcdf map-projections proj