采用反距离加权 (IDW) 插值的热图

Question

我有这个数据框，我需要使用反距离加权（IDW）插值。我已经查过了，但找不到办法，我想要的是一张热图，其中包含在该地方获取的水 pH 数据的 IDW 方法。

pacman::p_load(dplyr,sf,raster,sp,gstat,mapview)

dflima<-getData("GADM",country="PER", level=3) %>% st_as_sf() %>% filter(NAME_2=="Lima")

long <- c(-77.0958901385, -76.7059862868,  -77.0490377322, -77.0082420569, -77.029088411, -76.8875419285)
lat <- c(-11.702568503, -12.2926462918, -12.0510554271, -11.8780894556,  -12.1207933935,  -12.2201247629)
ph <- c(7.64, 7.66, 7.71, 8.10, 8.22, 8.29)
df<-data.frame(long,lat,ph)

我尝试制作密度图，但这不是我真正想要的。

df %>% ggplot() +
  stat_density2d(aes(x=df1$long, y=df1$lat, fill =..level..),bins=5 ,geom = "polygon")+
  scale_fill_distiller(palette = "Spectral")+
  geom_sf(data=dflima ,col = "#000000",alpha=0)

我什至进行了插值（在方形区域中），但我无法将其与地图的形状重叠在图表上

我想要类似于此图的东西：

如果您能帮助我，我将不胜感激

Answer 1

不确定您在插值方面遇到了什么问题，尽管这是一种简化的方法。这里我们将使用

sf

和

stars

栅格网格与

idw()

，结果也将是一个

stars

对象。

suppressPackageStartupMessages({
  library(geodata)
  library(ggplot2)
  library(dplyr)
  library(stars)
  library(gstat)
  library(sf)
})
geodata_path("~/geodata_dl")
dflima <- 
  gadm(country="PER", level=3) %>% 
  st_as_sf() %>% 
  filter(NAME_2=="Lima") %>% 
  st_geometry()

long <- c(-77.0958901385, -76.7059862868,  -77.0490377322, -77.0082420569, -77.029088411, -76.8875419285)
lat <- c(-11.702568503, -12.2926462918, -12.0510554271, -11.8780894556,  -12.1207933935,  -12.2201247629)
ph <- c(7.64, 7.66, 7.71, 8.10, 8.22, 8.29)
df <-data.frame(long,lat,ph)

# look for projected crs suggestions
crsuggest::suggest_crs(dflima, limit = 3)
#> # A tibble: 3 × 6
#>   crs_code crs_name                   crs_type  crs_gcs crs_units crs_proj4     
#>   <chr>    <chr>                      <chr>       <dbl> <chr>     <chr>         
#> 1 5387     Peru96 / UTM zone 18S      projected    5373 m         +proj=utm +zo…
#> 2 24892    PSAD56 / Peru central zone projected    4248 m         +proj=tmerc +…
#> 3 24878    PSAD56 / UTM zone 18S      projected    4248 m         +proj=utm +zo…

# transform points and shapes from lat/lon to EPSG:5387, this will also be the 
# projection for grid and IDW results
crs <- st_crs("EPSG:5387")
lima_3857 <- st_transform(dflima, crs)
ph_3857 <- st_as_sf(df, coords = c("long", "lat"), crs = "WGS84") %>% st_transform(crs)

# grid for IDW, cropped by lima_3857 shape (resulting raster is shaped like Lima)
# dx is grid cell size, controls raster resolution
grd <- 
  st_bbox(lima_3857) %>% 
  st_as_stars(dx = 200) %>%  
  st_crop(lima_3857) 

# interpolate
ph_idw <- idw(ph~1, ph_3857, grd)
#> [inverse distance weighted interpolation]

# plot
ggplot() +
  geom_stars(data = ph_idw, aes(fill = var1.pred)) +
  geom_sf(data=lima_3857 ,col = "#000000",alpha=0) +
  geom_sf(data = ph_3857, aes(fill = ph), shape = 21, size = 2, color = "white", stroke  = 1) +
  geom_sf_label(data = ph_3857, aes(label = ph), alpha = .8, nudge_x = 6e3, nudge_y = 5e3) +
  scale_fill_viridis_c(na.value = "transparent", name = "ph_idw")

^{创建于 2023-10-03，使用 reprex v2.0.2}