James Joule discovered that heat and work can be related by:
h = kw
Where k is a constant equal to 4.184 calories/joule.But really, this isn't a relationship between heat and work. This is a statement of equivalence. The constant of 4.184cal/J just tells you how the two units differ. They both describe the same property. Really, h = w. The version with k is only needed if you insist on using the traditional, different units on the left and right sides.
The same is true of e = mc^2. This does not tell you how two different fundamental properties are related, it merely tells you how two units, traditionally considered separate, actually measure the same underlying property. Really, e = m. The c^2 is only needed if you insist on using J on one side and kg on the other.
That is to say, you absolutely can have a kilogram of energy. If you have an object with a mass of one kilogram, you do have a kilogram of energy. Physicists and engineers don't typically use that unit for energy, but there's no reason they couldn't, and likewise no reason they couldn't use joules for mass.
So any time two physical quantities are related by a constant factor -- even if said constant factor has what we currently think of as "units" (and even particularly weird units like squared speed) -- they are fundamentally equivalent. I'm still mulling this over, but it's a pretty fascinating mental revelation for me.
Edit: I also want to say that I didn't find the comment I replied to rude in any way. It's direct and straightforward, but that's pretty common around here.
In natural units, c is "unitless" (AKA dimensionless), by definition.
Maybe someone has a link that goes into more depth on why all of this is okay?
When using a pseudo-Riemannian manifold one uses a metric signature where (keeping it simple, cf. "metric signature" on wikipedia) coordinates on orthogonal dimensions of can take one sign, or the opposite sign and a constant multiplier. In the Lorentzian case there will be one dimension taking one sign, and one or more taking the other; the choice of whether the solitary dimension takes a + or - sign is a matter of convention or preference. Conventionally the solitary dimension also takes the constant multiplier and is called the timelike dimension. A Lorentzian signature (minuses, 1) or (1, plusses) guarantees that one can describe paths through the manifold as null, timelike, or spacelike, and this gives one a causal structure.
Our universe can be well represented by a Lorentzian manifold (3, 1) or (1, 3), and this has been tested to exquisite precision. It does not tell us the value of the constant c, but we can determine that from tests of causal relations, the boundaries of which will be null. Alternatively, a massless pointlike object will always travel on null geodesics.
One runs into c being set to unity in systems geometrized units in relativity often; it's very handy to have mark off coordinates as e.g. "-seconds" vs "c seconds" (-,+,+,+ aka (1,3)) as long as one doesn't mess up one's dimensional analyses (which is unfortunately easy).
As a concrete example, the line element for (1,3) flat spacetime using Cartesian coordinates is dS^2 = -c^2dt^2 + dx^2 + dy^2 + dz^2. Compare the formula for Euclidean distance in the (Euclidean flat) plane between points p = (p1,p2) and q = (q1,q2) for (x,y) coordinates: d(p,q) = sqrt((q1-p1)^2 + (q2-p2)^2), which we could rewrite as dx = q1-p1, dy = q2-p2, ds = sqrt(dx^2 + dy^2) or ds^2 = dx^2 + dy^2. In Euclidean 3-space, we add another axis: ds^2 = dx^2 + dy^2 + dz^2. In Lorentzian 4-spacetime, we have to change the sign, so ds^2 = dx^2 + dy^2 + dz^2 - c^2dt^2. Note that we are not restricted to use any particular unit of distance or system of units; dx could be in metres, miles, astronomical units, light-years, gigaparsecs or practically anything else, while dt could be in seconds or fortnights or any other handy unit of time. When setting c to unity, we do need to choose appropriate units for dx (and dy and dz ...) vs dt. In SI units, that's seconds and light-seconds, but we could use another system if we wanted.
Notably we aren't restricted to Cartesian coordinates, however if we were to change to some other system of coordinates (e.g. polar ones) the line element would need to reflect that. For example, in the Lorentzian 4-spacetime case, we would write ds^2 = dr^2 + r^2 * dtheta^2 + r^2sin^2(theta)dphi^2 -c^2dt^2.
Finally, when using SI units (for example), one can see very clearly that a path taken by an object moving much slower than the speed of light is totally dominated by the amount of time between starting point and finishing point, because the value of c is large. Using the (+,-,-,-) metric signature [ds^2 = c^2dt^2 - dx^2 - dy^2 - dz^2], a large c helps make it clear that lightlike paths through spacetime are shorter, and purely timelike paths (where dx=0, dy=0, dz=0, dt != 0) have extremized length, since we don't subtract anything from cdt. This is the root of the explanation of the twin paradox in flat spacetime: the twin moving quickly compared to the speed of light takes a shorter path between together1 = (x1,y1,z1,t1) and together2 = (x2,y2,z2,t1) than the twin moving slowly compared to the speed of light. In the extreme, twin A holds x=const,y=const,z=const, whereas twin B's x coordinate is only equal to twin A's at the start and end of the journey. This holds up under any system of coordinates; we could consider r=const vs changing r in spherical coordinates, for instance.