Assuming the observations are collected independently of each other, the easiest way I can think of is to propagate uncertainty by using simulation. The idea is to generate random vectors from the (hyper-)cube and take the average of their coordinates; do this a large number of times and collect all the values obtained (the histogram below).
Here is an R code to illustrate it.
# set up a function to handle the simulation process
gen_x <- function() {
x = c(runif(1, 16-0.5,16+0.5),
runif(1, 21-0.5,21+0.5),
runif(1, 22-0.5, 22+0.5))
return(mean(x))
}
N=1e+4
# evaluate the function N times
hh <- sapply(1:N, function(x) gen_x())
hist(hh)
abline(v = mean(c(16,21,22)))
The histogram shows the distribution of the average of 16,21,22 taking into account the uncertainty of $\pm 0.5$ and assuming independent uniform distributions centred at the given values. The vertical line shows the sample average $(16+21+22)/3$.
The above solution gives the entire distribution of the sample average. If you are only interested in the variance of the sample average, then that's just equal to $1/(12n)$. Indeed, If you denote the sample $X_1,\ldots,X_n$, setting $\bar X = n^{-1}\sum_i X_i$ the sample average and assuming $X_i$ are independent, we have
$$
\text{var}(\bar X) = n^{-2}(\text{var}(X_1)+\cdots +\text{var}(X_n)) = n^{-2}\frac{n}{12} = (12n)^{-1}.
$$
Here I have used the fact that if $X \sim \text{Unif}(a,b)$, then $\text{var}(X) = (b-a)^2/12$.
So to sum up, if your $n$ samples are independent and the measures can be thought of as $\pm 0.5$ around the value measured, then, under the uniform assumption of the measure within each interval, the variance of the average is $(12n)^{-1}$.