我想使用ggplot
和facet_wrap
生成12张Map,这些Map显示了1月至12月苏格兰大陆架的海洋温度。我得到了一个包含25,000个观测值的“.csv”文件,其具有2016-2018年的每月温度数据以及相关联的纬度和经度值。我曾尝试使用geom_raster
Map此数据,但收到以下错误消息
Error: cannot allocate vector of size 39.7 Gb
In addition: Warning messages:
1: In f(...) :
Raster pixels are placed at uneven horizontal intervals and will be shifted. Consider using geom_tile() instead.
2: In f(...) :
Raster pixels are placed at uneven vertical intervals and will be shifted. Consider using geom_tile() instead
当然,我也尝试过使用geom_tile
,但当我运行代码时,我得到了一个没有颜色的空白Map,尽管我已经指定了填充和颜色。下面是一些类似于我的 Dataframe 的示例数据,请注意,我的latitude
和longitude
值的真实的数据不是由seq(...,...,0.1)
均匀隔开的(抱歉,我不知道如何在不进行排序的情况下生成示例数据),因此您不会遇到与我相同的错误,代码也将正常工作
#In my dataframe the lat and long values are not equally spaced!
mapoc_temp = expand.grid(data.frame(Longitude= seq(-64.5,-62.4,0.1),
Latitude= seq(42.7,44.8,0.1),
year = sample(c(2016,2017,2018), 22, replace = T),
month = sample(month.abb, 22, replace = T)))
mapoc_temp$Temp = runif(nrow(mapoc_temp))
下面是geom_raster
的代码,它给出了我提到的错误
library(mapdata)
library(ggplot2)
#map I'm using for ggplot
canada = map_data("worldHires", "Canada")
ggplot(mapoc_temp, aes(x=Longitude, y=Latitude)) +
#try to map my temp on ocean
geom_raster(aes(fill = Temp, x = Longitude), interpolate = TRUE) +
geom_polygon(data = canada, aes(x=long, y=lat, group=group), colour="grey50", fill = 'grey55')+
coord_sf(xlim=c(-64.5,-62.8), ylim=c(42.7,45)) +
#get months
facet_wrap(vars(month))
这是我尝试geom_tile
时的代码,我再次注意到它在我的样本数据中有效,但在我的真实的数据中无效,我相信这与我的坐标之间的距离不相等有关。
ggplot(mapoc_temp, aes(x=Longitude, y=Latitude)) +
geom_tile(aes(fill = Temp, x = Longitude, y = Latitude), colour = mapoc_temp$Temp) +
geom_polygon(data = canada, aes(x=long, y=lat, group=group), colour="grey50", fill = 'grey55')+
coord_sf(xlim=c(-64.5,-62.8), ylim=c(42.7,45)) +
facet_wrap(vars(month))
这是使用geom_tile获得的图片
这是我想要制作的大致内容,但显然是12张Map,因为我希望每个月都有一张分辨率更高的Map(如果可能,但在这一点上,我会选择任何用我的温度数据着色的Map)。
我有一种感觉,我必须沿着这类代码(我不久前发现的代码片段)做更多的事情,但我已经尝试过操纵它,没有成功。有什么建议吗?
#my CRS value!
latlong = "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"
#for transforming
planar ="+proj=utm +zone=20 +datum=NAD83 +units=km +no_defs +ellps=GRS80 +towgs84=0,0,0"
out <- x%>%
st_as_sf(coords=c("Longitude","Latitude"),crs=latlong)%>%
st_transform(planar)%>%
st_coordinates()%>%
raster::rasterize(.,grid,field=xyz[,3],fun=rasterFun)%>%
raster::projectRaster(.,crs=latlong)%>%
raster::rasterToPolygons(.)%>% # this part is slow
st_as_sf()%>%
rename(MAP = layer)
3条答案
按热度按时间fae0ux8s1#
看起来geom_polygon(Map图层)覆盖了geom_tile图层。下面是使用示例数据的示例,geom_tile称为last,alpha设置为0.4以显示下面的Map。
由reprex package(v0.3.0)于2020年2月22日创建
thtygnil2#
如果geom_polygon图层位于geom_raster之上,则还可以在geom_polygon函数中设置fill= NA
wfsdck303#
你已经用你提供的示例数据回答了你自己的问题。
geom_raster
和geom_tile
函数有它们的使用限制。一个很大的限制是,因为它们要求数据在均匀的网格上,对于非常高的分辨率数据,机器将在绘图前尝试创建均匀的网格,因此需要39.7GB内存。通过将lon/lat值舍入到0.1x0.1网格,计算机现在可以绘制本文其他答案中的数据。尝试对数据运行以下代码块。您可以根据自己的喜好更改四舍五入的级别,但除非绝对必要,否则应避免使用太细的比例。一旦您的数据四舍五入到偶数网格中,您就可以使用以下块来绘制它们。
没有必要创建任何sf形状多边形或添加任何CRS信息到绘图中。tidyverse的设计不是为了利用这些信息,在我看来,它只会使事情不必要地复杂化,以带来sf包。我希望这有帮助!