r - 如何在 ggplot2 中正确绘制投影网格数据？

Question

我一直在使用 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)))

如果你想完全理解问题，也许你想尝试大域和小域。代码是相同的，只是点大小和 map 限制发生了变化。下面的值适用于大的完整域。
好的，现在让我们开始吧!
从瓷砖开始
最明显的解决方案是使用瓷砖。咱们试试吧。

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...

好的，但那些不是真正的模型方 block ，这是一个黑客。此外，模型框在域的顶部发散，并且都以相同的方式定向。不太好。让我们自己投影正方形，即使我们已经知道这不是正确的做法……也许它看起来不错。

#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

好吧，让我们尝试用 lat-lons 代替...

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 无法用于出版物。

Answer 1

在深入挖掘之后，您的模型似乎基于“兰伯特圆锥”投影中的 50 公里规则网格。但是，您在 netcdf 中的坐标是“单元格”中心的 lat-lon WGS84 坐标。

鉴于此，一种更简单的方法是重建原始投影中的单元格，然后在转换为 sf 后绘制多边形。对象，最终经过重投影。像这样的东西应该可以工作(请注意，您需要从 github 安装 devel 版本的 ggplot2 才能工作):

load(url('https://files.fm/down.php?i=kew5pxw7&n=loadme.Rdata'))
library(raster)
library(sf)
library(tidyverse)
library(maps)
devtools::install_github("hadley/ggplot2")

#   ____________________________________________________________________________
#   Transform original data to a SpatialPointsDataFrame in 4326 proj        ####

coords = data.frame(lat = values(s[[2]]), lon = values(s[[3]]))
spPoints <- SpatialPointsDataFrame(coords, 
                                   data = data.frame(data = values(s[[1]])), 
                                   proj4string = CRS("+init=epsg:4326"))

#   ____________________________________________________________________________
#   Convert back the lat-lon coordinates of the points to the original      ###
#   projection of the model (lcc), then convert the points to polygons in lcc
#   projection and convert to an `sf` object to facilitate plotting

orig_grid = spTransform(spPoints, projection(s))
polys = as(SpatialPixelsDataFrame(orig_grid, orig_grid@data, tolerance = 0.149842),"SpatialPolygonsDataFrame")
polys_sf = as(polys, "sf")
points_sf = as(orig_grid, "sf")

#   ____________________________________________________________________________
#   Plot using ggplot - note that now you can reproject on the fly to any    ###
#   projection using `coord_sf`

# Plot in original  projection (note that in this case the cells are squared): 
my_theme <- theme_bw() + theme(panel.ontop=TRUE, panel.background=element_blank())

ggplot(polys_sf) +
  geom_sf(aes(fill = data)) + 
  scale_fill_distiller(palette='Spectral') +
  ggtitle("Precipitations") +
  coord_sf() + 
  my_theme 

# Now Plot in WGS84 latlon projection and add borders: 

ggplot(polys_sf) +
  geom_sf(aes(fill = data)) + 
  scale_fill_distiller(palette='Spectral') +
  ggtitle("Precipitations")  +
  borders('world', colour='black')+
  coord_sf(crs = st_crs(4326), xlim = c(-60, 80), ylim = c(15, 75))+
  my_theme 

但是，要在原始投影中添加边界，您必须将 loygon 边界提供为 sf。目的。从这里借:

Converting a "map" object to a "SpatialPolygon" object

像这样的东西会起作用:

library(maptools)
borders  <- map("world", fill = T, plot = F)
IDs      <- seq(1,1627,1)
borders  <- map2SpatialPolygons(borders, IDs=borders$names, 
                               proj4string=CRS("+proj=longlat +datum=WGS84")) %>% 
            as("sf")

ggplot(polys_sf) +
  geom_sf(aes(fill = data), color = "transparent") + 
  geom_sf(data = borders, fill = "transparent", color = "black") +
  scale_fill_distiller(palette='Spectral') +
  ggtitle("Precipitations") +
  coord_sf(crs = st_crs(projection(s)), 
           xlim = st_bbox(polys_sf)[c(1,3)],
           ylim = st_bbox(polys_sf)[c(2,4)]) +
  my_theme

作为旁注，既然我们“恢复”了正确的空间引用，也可以构建正确的 raster数据集。例如:

r <- s[[1]]
extent(r) <- extent(orig_grid) + 50000

会给你一个合适的 raster在 r :

r
class       : RasterLayer 
band        : 1  (of  36  bands)
dimensions  : 125, 125, 15625  (nrow, ncol, ncell)
resolution  : 50000, 50000  (x, y)
extent      : -3150000, 3100000, -3150000, 3100000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_1=30. +lat_2=65. +lat_0=48. +lon_0=9.75 +x_0=-25000. +y_0=-25000. +ellps=sphere +a=6371229. +b=6371229. +units=m +no_defs 
data source : in memory
names       : Total.precipitation.flux 
values      : 0, 0.0002373317  (min, max)
z-value     : 1998-01-16 10:30:00 
zvar        : pr 

看到现在分辨率是50Km，范围是公制坐标。因此，您可以使用 r 进行绘图/工作使用 raster 的函数数据，例如:

library(rasterVis)
gplot(r) + geom_tile(aes(fill = value)) + 
  scale_fill_distiller(palette="Spectral", na.value = "transparent") +
  my_theme  

library(mapview)
mapview(r, legend = TRUE)