Extracting areal precipitation from multiple location from large netcdf files

Question

I have huge netcdf files, storing precipitation data for whole Germany with 1 km spatial and 5 min temporal resolution. I have 372 equidistant points across the country where the areal precipitation needs to be calculated within a circular area around these point with different radii (2km, 4km, 6km, 8km, ..., 18 km). The files are saved in daily resolutions and extracting what I need takes a long time. I have 22 years of data in daily files, and I need get the areal averages for 10 area sizes. So this needs to be repeated 2210365 = 80300 iterations. I ran the following code only for point.x.y[1:10,] and it tool over 15 minutes. Even with paralel processing this calculation will take months. Is there any efficient way of doing this?

I have tried brick() and extract() in R and xarray in Python.

nc.flname  = paste0(ncPath, "xxx.nc")
nc = nc_open(nc.flname)

lat_nc <- ncvar_get(nc, "lat")
lon_nc <- ncvar_get(nc, "lon")
x_nc = ncvar_get(nc, "x")
y_nc = ncvar_get(nc, "y")

Points_LatLon = read.csv(paste0(PointsPath, "StudyPoints_372_LatLon.csv"))
## This is necessary because the nc files's dimensions are not directly lat
### - lon. I find the closest cell of the grid to my desired point by 
## calculating the distance
points.row.cols = sapply(1:nrow(Points_LatLon), function(p){
  ind = (find_the_closest_pixel_ind(Points_LatLon$longitude[p], 
                                    Points_LatLon$latitude[p],
                                    lon_nc, lat_nc))
  row.col = which(lat_nc == lat_nc[ind] & lon_nc == lon_nc[ind], arr.ind = T)
  return(row.col)
})

point.x.y = cbind(x = x_nc[points.row.cols[1,]], y = y_nc[points.row.cols[2,]])

rainfall.brick = brick(nc.flname, varname ='RR') 
time_series <- extract(rainfall.brick, point.x.y, buffer = 2000, fun = mean)

I'd appreciate any solution in R, Python and/or a combination with ArcMap.