We propose a method to define quasiprobability distributions for general spin-$j$ systems of dimension $n=2j+1$, where $n$ is a prime or power of prime. The method is based on a complete set of orthonormal commuting operators related to Mutually Unbiased Bases which enable (i) a parameterisation of the density matrix and (ii) construction of measurement operators that can be physically realised. As a result we geometrically characterise the set of states for which the quasiprobability distribution is non-negative, and can be viewed as a joint distribution of classical random variables assuming values in a finite set of outcomes. The set is an $(n^2-1)$-dimensional convex polytope with $n+1$ vertices as the only pure states, $n^{n+1}$ number of higher dimensional faces, and $n^3(n+1)/2$ edges.
