heatmap with Inverse Distance Weighting (IDW) Interpolation

Question

I have this data frame, I need to use Inverse Distance Weighting (IDW) Interpolation. I have checked but I can't find a way to do it, what I want is a heat map with the IDW method for water pH data taken in that place.

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)

I tried to make a density plot but it's not what I really want.

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)

I even did the interpolation (in a square area) but I couldn't overlay it on the graph with the shape of the map:

I want something similar to this graph:

If you can help me I would be grateful.

Answer 1

Not sure what issues you had with interpolation, though here's one simplified approach. Here we'll use sf and stars raster grid with idw(), result will also be a stars object.

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

^{Created on 2023-10-03 with reprex v2.0.2}