How can I calculate weighted standard errors and plot them in a bar plot?

There isn't (as far as I know) a built-in R function to calculate the standard error of a weighted mean, but it is fairly straightforward to calculate - with some assumptions that are probably valid in the case you describe.
See, for instance: https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Standard_error

Standard error of the weighted mean

If the elements used to calculate the weighted mean are samples from populations that all have the same variance v, then the variance of the weighted sample mean is estimated as:

var_m = v^2 * sum( wnorm^2 )   # wnorm = weights normalized to sum to 1

And the standard error of the weighted mean is equal to the square root of the variance.

sem = sqrt( var_m )

So, we need to calculate the sample variance from the weighted data.

Weighted variance

The weighted population variance (or biased sample variance) is calculated as:

pop_v = sum( w * (x-mean)^2 ) / sum( w )

However, if (as in the case you describe), we are working with samples taken from the population, rather then with the population itself, we need to make an adjustment to obtain an unbiased sample variance.
If the weights represent the frequencies of observations underlying each of the elements used to calculate the weighted mean & variance, then the adjustment is:

v = pop_v * sum( w ) / ( sum( w ) -1 )

However, this is not the case here, as the weights are the total frequenceis of 911 calls for each policeman, not the calls for each beat. So in this case the weights correspond to the reliabilities of each element, and the adjustment is:

v = pop_v * sum( w )^2 / ( sum( w )^2 - sum( w^2) ) 

weighted.var and weighted.sem functions

Putting all this together, we can define weighted.var and weighted.sem functions, similar to the base R weighted.mean function (note that several R packages, for instance "Hmisc", already include more-versatile functions to calculate the weighted variance):

weighted.var = function(x,w,type="reliability") {
    m=weighted.mean(x,w)
    if(type=="frequency"){ return( sum(w*(x-m)^2)/(sum(w)-1) ) }
    else { return( sum(w*(x-m)^2)*sum(w)/(sum(w)^2-sum(w^2)) ) }
}
weighted.sem = function(x,w,...) { return( sqrt(weighted.var(x,w,...)*sum(w^2)/sum(w)^2) ) }

applied to 911 call data in the question

In the case of the question, the elements from which we want to calculate the weighted mean and weighted sem correspond to the proportions of calls in each beat, for each policeman.
So (finally...):

props = t(apply(df,1,function(row) row[-(1:3)]/row[3]))
wmean_props = apply(props,2,function(col) weighted.mean(col,w=df[,3]))
wsem_props = apply(props,2,function(col) weighted.sem(col,w=df[,3]))

Aren't your "proportions" actually the mean of the weighted (by total) observations? Then we could simply calculate the weighted colMeans accordingly.

df2 <- df[, grep('b_', colnames(df))]

means.w <- colMeans(df2 / df$total)

For the error bars we could use the quantiles of 1 - alpha/2, i.e. for alpha==.05 we use c(.025, .975). The analytical sds would yield negative values.

q.w <- t(apply(df2 / df$total, 2, quantile, c(.025, .975)))

Now, we store the x-positions that barplot returns invisible,

# Bar plot of proportions by beat
b <- barplot(means.w, main='Distribution of Calls by Beat', 
             xlab="Nth Most Common Division", ylim=c(0,1), 
             names.arg=1:length(means.w), ylab="Percent of Total Calls")

and construct the error bars with arrows.

arrows(b, q.w[,1], b, q.w[,2], length=.02, angle=90, code=3)