Whilst the GMT route is the best way to go long term, in dire straits you could use the polar histogram plotting function in an old repo of mine (see instructions in the readme file there).
julia> using CircPlot
julia> θ = 2π.*rand(10_000);
julia> cplot_histogram(θ, 2π/30)
will give you this:
and you can tweak the appearance using the usual Plots keywords:
julia> cplot_histogram(θ, 2π/30, showaxis=false, grid=false, line=(:red,1), fill=(:blue, 0.5), circ=(lc=:black, fill=nothing))
giving
(Note that passing a named tuple, Dict, etc. to the circ keyword argument sets the plotting style of the outer circle.)
The docstring describes the usage:
help?> cplot_histogram
search: cplot_histogram
cplot_histogram(a, binwidth, degrees=false; azimuth=false, axial=false, weights=ones(length(a)), circ=nothing; kwargs...) -> ::Plots.Plot
Plot a polar histogram of a set of angles `a` using bins `binwidth` wide.
`binwidth` should be given in the same convention (radians or °) as the data.
If π (or 180°) is not an integer multiple of `binwidth`, then the nearest value
for which this is true is selected.
Use `azimuth=true` to plot angles clockwise from north.
Use `axial=true` if data are π-ambiguous (i.e., are orientational rather than
directional data).
`weights` must have the same length as `a`, and if passed in will scale each
point `i` by `weights[i]`.
Set options for the surrounding circle using `circ=Dict(:opt=>val)`, where
the keys and values in the Dict will be passed to Plots for plotting the circle
shape around the plot.
Additional keyword arguments `kwargs...` are passed to `Plots.plot()` when
plotting the bars.
Ideally this should be implemented as a recipe, in StatsPlots or similar, but as yet there aren’t many directional distributions or measures implemented in Distributions nor tools for these in StatsBase. There is some very basic circular statistics stuff in CircStats.
Edit: I should say by way of explanation that you will want to use the weights keyword argument and set the binwidth to match the bar locations you need, with one point per bin and the bar length specified by the respective value in weights:
cplot_histogram(theta, step(theta), weights=r, azimuth=true, circ=(lw=1, fill=nothing), showaxis=false, grid=false)
