不知道插值有什么问题，但这里有一个简化的方法。这里我们将使用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