While looking at the test cases for Coordinates to find out why changes to Quantity I made caused problems (#1848), I noticed this in test_distance.py:

    #always preserve the CartesianPoint's units
    c5 = 3 * u.pc + c
    assert c5.unit == c.unit

which is really odd physically: what coordinate should the 3 pc be added to?

This brought home something I worried about while looking for the problems: by using a 3,N array, how can CartesianPoints protect against unphysical mixing of x, y, and z? Or ensure that methods such as np.mean work only work along relevant axes?

Easiest may be to behave more closely to the SphericalCoordinatesBase and just have x, y, and z Quantity arrays internally, but if they are combined in a single array, my suggestion would be to use a record array as the basis instead (with dtype=[('x', 'f8'), ('y', 'f8'), ('z', 'f8')]). In either case, one needs to override methods for the relevant functions.