code_01.R

library(sf)

## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0

library(stars)

## Loading required package: abind

set.seed(2)

# Data - points
dat = read.csv("https://geobgu.xyz/r/data/rainfall.csv")
dat = dat[sample(1:nrow(dat), 30), ]

# Data - raster
f = tempfile()
download.file("https://geobgu.xyz/r/data/rainfall.tif", f)
r = read_stars(f)
names(r) = "rainfall"

# To points
dat = st_as_sf(dat, coords = c("x_utm", "y_utm"), crs = 32636)

# Buffer
dat = st_buffer(dat, dist = 2000)

# Extract
dat$value = NA
for(i in 1:nrow(dat)) {
    r_subset = r[dat[i, ]]
    dat$value[i] = mean(r_subset[[1]], na.rm = TRUE)
}
dat$value

##  [1] 528.3299 615.0802 381.8471 304.6810 595.6562 429.1149 552.6085 556.0259
##  [9] 822.5706 407.0831 516.3507 332.4678 414.5561 220.3162      NaN 561.4758
## [17] 431.6985 622.6923 392.6890 557.3144 453.8358 429.7664      NaN 513.5844
## [25] 388.8415 511.5128 631.5563 535.3544 375.8217      NaN

# Plot
plot(st_geometry(dat), col = NA, reset = FALSE)
plot(r, add = TRUE)
plot(st_geometry(dat), add = TRUE)
text(st_coordinates(st_centroid(dat)), as.character(round(dat$value)), pos = 3, col = "red")

## Warning in st_centroid.sf(dat): st_centroid assumes attributes are constant over
## geometries of x