A known issue with astropy.units is that some NumPy operations take in arguments that have units, but return arguments that do not have units even though they should. This issue requires an upstream fix in NumPy (e.g., numpy/numpy#4164). The SciPy 2017 talk on MetPy's Choice of Unit Support delves into this in great detail. An example from the talk:

>>> import numpy as np
>>> from astropy import units
>>> a = np.array([3.]) * units.m
>>> b = np.array([4.]) * units.m
>>> np.concatenate((a, b))
<Quantity [ 3., 4.] (Unit not initialised)>

There are some workarounds to some of these problems.

>>> np.isclose(500* u.km/u.s, 300 * u.km / u.s
UnitsError: Can only apply 'add' function to dimensionless quantities when other argument is not a quantity (unless the latter is all zero/infinity/nan)
>>> np.isclose(500* u.km/u.s, 300 * u.km / u.s, atol=1e-8 * u.mm / u.s)
array([False], dtype=bool)

We cannot do anything from within PlasmaPy, but will need to wait for fixes within NumPy (e.g., new support for __array_ufunc__ will improve overriding behavior for ufuncs).