The coalescence of a neutron star with a black hole is a primary science
target of ground-based gravitational wave detectors. Constraining or measuring
the neutron star spin directly from gravitational wave observations requires
knowledge of the dependence of the emission properties of these systems on the
neutron star spin. This paper lays foundations for this task, by developing a
numerical method to construct initial data for black hole--neutron star
binaries with arbitrary spin on the neutron star. We demonstrate the robustness
of the code by constructing initial-data sets in large regions of the parameter
space. In addition to varying the neutron star spin-magnitude and
spin-direction, we also explore neutron star compactness, mass-ratio, black
hole spin, and black hole spin-direction. Specifically, we are able to
construct initial data sets with neutron stars spinning near centrifugal
break-up, and with black hole spins as large as $S_{\rm BH}/M_{\rm BH}^2=0.99$.