如何使用R中的非均匀间隔坐标生成高分辨率温度图?

vpfxa7rd  于 2023-02-06  发布在  其他
关注(0)|答案(3)|浏览(173)

我想使用ggplotfacet_wrap生成12张Map,这些Map显示了1月至12月苏格兰大陆架的海洋温度。我得到了一个包含25,000个观测值的“.csv”文件,其具有2016-2018年的每月温度数据以及相关联的纬度和经度值。我曾尝试使用geom_rasterMap此数据,但收到以下错误消息

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 的示例数据,请注意,我的latitudelongitude值的真实的数据不是由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)
fae0ux8s

fae0ux8s1#

看起来geom_polygon(Map图层)覆盖了geom_tile图层。下面是使用示例数据的示例,geom_tile称为last,alpha设置为0.4以显示下面的Map。

#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))

  library(mapdata)
#> Loading required package: maps
  library(ggplot2)

  #map I'm using for ggplot
  canada = map_data("worldHires", "Canada")

  ## Use geom_raster as the last layer in the ggplot2 call,
  ##  otherwise the polygons plot over the tiles.
  ## Below alpha is set on the raster layer to show underlying map.
  ggplot(mapoc_temp, aes(x=Longitude, y=Latitude)) + 
    #try to map my temp on ocean
    geom_polygon(data = canada, aes(x=long, y=lat, group=group), colour="grey50", fill = 'grey55')+
    geom_raster(aes(fill = Temp, x = Longitude),alpha = .4, interpolate = TRUE) +
    coord_sf(xlim=c(-64.5,-62.8), ylim=c(42.7,45)) +
    #get months
    facet_wrap(vars(month))

reprex package(v0.3.0)于2020年2月22日创建

thtygnil

thtygnil2#

如果geom_polygon图层位于geom_raster之上,则还可以在geom_polygon函数中设置fill= NA

wfsdck30

wfsdck303#

你已经用你提供的示例数据回答了你自己的问题。geom_rastergeom_tile函数有它们的使用限制。一个很大的限制是,因为它们要求数据在均匀的网格上,对于非常高的分辨率数据,机器将在绘图前尝试创建均匀的网格,因此需要39.7GB内存。通过将lon/lat值舍入到0.1x0.1网格,计算机现在可以绘制本文其他答案中的数据。尝试对数据运行以下代码块。

mapoc_temp <- mapoc_temp %>% 
  ungroup() %>% 
  mutate(Longitude = plyr::round_any(Longitude, 0.1), 
         Latitude = plyr::round_any(Latitude, 0.1))

您可以根据自己的喜好更改四舍五入的级别,但除非绝对必要,否则应避免使用太细的比例。一旦您的数据四舍五入到偶数网格中,您就可以使用以下块来绘制它们。

ggplot(mapoc_temp, aes(x = Longitude, y = Latitude)) +
  borders(fill = "grey80") +
  geom_raster(aes(fill = Temp)) +
  coord_quickmap(xlim = c(-64.5,-62.8), ylim = c(42.7,45)) +
  facet_wrap(~month)

没有必要创建任何sf形状多边形或添加任何CRS信息到绘图中。tidyverse的设计不是为了利用这些信息,在我看来,它只会使事情不必要地复杂化,以带来sf包。我希望这有帮助!

相关问题