问题描述
使用 SF 在栅格中聚合值
我需要的是为每个栅格汇总一些指标的值。 假设我们有一些数据 - 坐标和值,我想创建一个热图。
首先,我使用简单的要素框架创建了一个网格和栅格。
现在我需要获取 df
中的每个坐标并检查它们是否在一个栅格中。然后为每个栅格计算任何其他聚合函数的平均值。
# Packages ----------------------------------------------------------------
library(raster)
library(tidyverse)
library(sf)
library(sp)
# Border Data and Grid ----------------------------------------------------
Regions <- getData("GADM",country= "CZE",level = 1)
Regions %>%
Regions <-
Regions %>% st_as_sf()
grid_spacing <- 0.25
polygony <- st_make_grid(Regions,square = T,cellsize = c(grid_spacing,grid_spacing)) %>%
st_sf()
plot(polygony,col = 'white')
plot(st_geometry(CZ),add = T)
A = st_intersection(polygony,CZ)
# Artifial Values to be use -----------------------------------------------
df <-
tibble(long = runif(500,13.27857,14),lat = runif(500,49,50),price = rnbinom(500,size = 40,0.3))
df %>%
ggplot(aes(long,lat,color = values)) +
geom_point()
A <-
A %>%
# Here Val shall be calculated as a mean of observations within each rastel grid cell
mutate(val = rnorm(n = 202))
# Since not every raster cell has observations inside some NAN will be present
A[5,'val'] = NaN
A[25,'val'] = NaN
A %>%
ggplot(aes(fill = val)) +
geom_sf()
# What I need -------------------------------------------------------------
# For each raster in A calculate the mean of price for coordinates within this raster grid.
在 sf 框架中,我基本上需要一些等效的 https://www.rdocumentation.org/packages/raster/versions/3.4-5/topics/rasterize。
但是,我想绘制看起来像在示例中我真的需要栅格具有给定网格的形状。
解决方法
要汇总网格中的点值并保留 sf,您可以执行以下操作:
library(dplyr)
library(sf)
# get data. using your example,we'll take spatial data from the raster package
Regions <- raster::getData("GADM",country= "CZE",level = 1)
# convert spdf to sf
Regions <- Regions %>% st_as_sf()
# create a grid sf object
grid_spacing <- 0.25
polygony <- st_make_grid(Regions,square = T,cellsize = c(grid_spacing,grid_spacing)) %>%
st_sf() %>%
mutate(ID = row_number()) # add a unique ID to each grid cell
# clip grid to shape of country polygons
A <- st_intersect(polygony,Regions)
# create fake data with coordinates and prices
df <- tibble(long = runif(500,13.27857,14),lat = runif(500,49,50),price = rnbinom(500,size = 40,0.3))
# convert the df to sf point layer
points <- st_as_sf(df,coords = c("long","lat"),crs = st_crs(A))
# spatially join grid to points,so that each point is assigned the grid ID into which it falls
pointsID <- st_join(points,A)
# group and summarize point values by grid ID
pointsID <- pointsID %>%
as.data.frame() %>%
group_by(ID) %>%
summarize(avg_price = mean(price))
# join aggregated values back to your grid
A<- left_join(A,pointsID,by = "ID")
plot(A["avg_price"])
,
因为您无论如何都使用光栅包
library(raster)
cze <- getData("GADM",level = 1)
r <- raster(cze,res=0.25)
df <- data.frame(long = runif(500,0.3))
r <- rasterize(df[,c("long","lat")],r,df$price,mean,background=0)
p <- as(r,"SpatialPolygonsDataFrame")
p <- crop(p,cze)
library(sf)
s <- st_as_sf(p)
plot(p)