【问题标题】:R: how to create a heat map of averaged values from a grid and plot it with ggplot?R:如何从网格创建平均值的热图并用 ggplot 绘制它?
【发布时间】:2020-09-10 14:38:34
【问题描述】:

我有一个数据框(见下文),其中包含超过 50 000 个值,每个值都与一个位置(纬度、经度)相关联。我想计算 5° 纬度 x 5° 经度网格的每个单元格的平均值,以创建热图。最终目标是在测深地图上绘制此网格。

我查看了类似Average values of a point dataset to a grid dataset 的问题。但我无法用我自己的数据复制这些示例。可悲的是,我被困在创建网格的第一步。

我的数据如下所示:

library(sp)
library(proj4)

coordinates(data) <- c("lon", "lat")        
proj4string(data) <- CRS("+init=epsg:4326") #defined CRS to WGS 84
df<- data.frame(data)

> head(df)
         lon      lat  value
1 -48.1673562 57.71791  822.9
2 -48.7430053 57.83568 1302.3
3 -48.5662663 57.82087 1508.0
4 -48.3252052 58.29815  224.0
5 -47.1716772 58.42417   38.0
6 -46.4098311 58.67651  431.2
7 -45.8071218 58.70022  365.6
8 -45.5558936 58.46975   50.0

理想情况下,我想使用 ggplot2 在 marmap 包中的地图上绘制网格(见下文):

library(marmap)
library(ggplot2)

atlantic <- getNOAA.bathy(-80, 40, 0, 90, resolution = 25, keep = TRUE)

atl.df <- fortify(atlantic)

map <- ggplot(atl.df, aes(x=x, y=y)) +
  geom_raster(aes(fill=z), data=atl.df) +
  geom_contour(aes(z=z),
               breaks=0, #contour for continent
               colour="black", size=1) +
  scale_fill_gradientn(values = scales::rescale(c(-5000, 0, 1, 2400)),
                       colors = c("steelblue4", "#C7E0FF", "gray40", "white"))

【问题讨论】:

    标签: r ggplot2 geospatial raster spatial


    【解决方案1】:

    听起来您想将数值变量(纬度和经度)切割成偶数区间并总结每个区间内的值。以下内容对您有用吗?

    library(dplyr)
    
    df2 <- df %>%
      mutate(lon.group = cut(lon, breaks = seq(floor(min(df$lon)), ceiling(max(df$lon)), by = 5),
                             labels = seq(floor(min(df$lon)) + 2.5, ceiling(max(df$lon)), by = 5)),
             lat.group = cut(lat, breaks = seq(floor(min(df$lat)), ceiling(max(df$lat)), by = 5),
                             labels = seq(floor(min(df$lat)) + 2.5, ceiling(max(df$lat)), by = 5))) %>%
      group_by(lon.group, lat.group) %>%
      summarise(value = mean(value), .groups = "drop") %>%
      mutate(across(where(is.factor), ~as.numeric(as.character(.x))))
    

    样本数据:

    set.seed(444)
    
    n <- 10000
    df <- data.frame(lon = runif(n, min = -100, max = -50),
                     lat = runif(n, min = 30, max = 80),
                     value = runif(n, min = 0, max = 1000))
    
    > summary(df)
          lon              lat            value          
     Min.   :-99.99   Min.   :30.00   Min.   :   0.1136  
     1st Qu.:-87.55   1st Qu.:42.45   1st Qu.: 247.2377  
     Median :-75.29   Median :55.11   Median : 501.4165  
     Mean   :-75.12   Mean   :55.01   Mean   : 499.5385  
     3rd Qu.:-62.69   3rd Qu.:67.63   3rd Qu.: 748.8834  
     Max.   :-50.01   Max.   :80.00   Max.   : 999.9600 
    

    前后数据对比:

    gridExtra::grid.arrange(
      ggplot(df, 
             aes(x = lon, y = lat, colour = value)) + 
        geom_point() + 
        ggtitle("Original points"),
      ggplot(df2, 
             aes(x = lon.group, y = lat.group, fill = value)) + 
        geom_raster() + 
        ggtitle("Summarised grid"),
      nrow = 1
    )
    

    【讨论】:

    • 非常感谢@Z.Lin 的回答。这正是我一直在寻找的,因为它允许直接指定网格单元的大小。
    【解决方案2】:

    就像(几乎!)一样,有一个功能。我相信marmap::griddify() 是您正在寻找的。帮助文件指出:

    将不规则间隔的 xyz 数据转换为适合创建具有规则间隔的经度和纬度的深海对象的栅格对象。

    这是一个使用您的坐标的脚本:

    library(marmap)
    library(ggplot2)
    
    # Create fake data
    set.seed(42)
    n <- 10000
    data_irregular <- data.frame(lon = runif(n, min = -80, max = 40),
                                 lat = runif(n, min = 0, max = 90),
                                 value = runif(n, min = 0, max = 1000))
    
    # Fit data into a grid of 30 cells in longitude and 50 cells in latitude
    data_grid <- as.bathy(griddify(data_irregular, nlon = 30, nlat = 50))
    fortified_grid <- fortify(data_grid)
    
    # Get bathymetric data to plot continent contours
    atlantic <- getNOAA.bathy(-80, 40, 0, 90, resolution = 25)
    atl_df <- fortify(atlantic)
    
    # Plot with ggplot with gridded data as tiles
    map <- ggplot(atl_df, aes(x = x, y = y)) +
      geom_raster(data = fortified_grid, aes(fill = z)) +
      geom_contour(data = atl_df, aes(z = z), 
                   breaks = 0, # contour for continent
                   colour = "black", size = 1) +
      scale_fill_gradientn(values = scales::rescale(c(-5000, 0, 1, 2400)),
                           colors = c("steelblue4", "#C7E0FF", "gray40", "white")) +
      labs(x = "Longitude", y = "Latitude", fill = "Value")
    
    
    map +
      theme_bw()
    

    结果如下:

    【讨论】:

    • 谢谢@Benoit。你的回答总是很有帮助。这很好用,关于如何将结果集成到 ggplot2 中的细节非常有用。
    猜你喜欢
    • 1970-01-01
    • 2018-06-17
    • 2020-04-18
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2020-05-02
    相关资源
    最近更新 更多