As I understand it, a natural system of units is one in which the numerical values of $c$ and $\hbar$ are unity, i.e. $c=\hbar =1$.
What I find confusing is that they are still dimensionful, i.e. $[c]=LT^{-1}$ and $[\hbar]=ML^{2}T^{-1}$. So, how can the action be dimensionless, $[S]=1$, when it has the same dimensions as $\hbar$?
It is dimensionless in the sense of mass dimension.
Setting $\hbar = c = 1$ means we only need to fix one base unit, which is usually taken to be the energy measured in $\mathrm{eV}$. Now, since $c=1$, this means that because of $E=mc^2$ becoming $E=m$ both $E$ and $m$ are measured in $\mathrm{eV}$. They might represent different dimensions (mass and energy), but they are measured in the same unit. Now, $E=h\nu$ means that inverse time is measured in $\mathrm{eV}$, so time is measured in $\mathrm{eV^{-1}}$. And so on.
Now, the "mass dimension" of a quantity is simply the power of $\mathrm{eV}$ it is measured in. Since the action is the integral of an energy against time, it has units of $\mathrm{eV}\cdot\mathrm{eV}^{-1} = \mathrm{eV}^0$, i.e. it has mass dimension zero.
You are right in that it is not "dimensionless". But having mass dimension zero means for any quantity $Q$ that there are powers of $\hbar$ and $c$ such that $\frac{Q}{\hbar^n c^m}$ is dimensionless, and since $\hbar = c = 1$, $\frac{Q}{\hbar^n c^m} = Q$, so there is no numerical difference between those quantities, and one sloppily says that $Q$ is dimensionless.
If you are somewhat worried that $\frac{Q}{\hbar^n c^m} = Q$ "looks wrong" from a dimensional analysis standpoint, then yes, that's right - the convenience in the formulae we get from $\hbar = c = 1$ comes inevitably with the loss of a large part of dimensional analysis, all that is left is the mass dimension for that.
To awnser your second question :
Is there any particular reason for wanting the action to be dimensionless? Is it simply that we want its numerical value to be independent of any system of units?
there are several reasons to it. One of them is that we don't want to bother about keeping track of the $\hbar$ and $c$'s everywhere, as they can add up pretty fast. One other, more fundamental reason is that many physical phenomenon are scale-dependent. In the natural system of units, lenght is measured in the same units as inverse energy, so looking at a phenomenon at different scales is the same as looking at it at different energies. Therefore, measuring everything in terms of energy (or mass) enables to clearly see how a theory behaves at different energy scales : this idea is one of the basic principles of the Renormalization Group method. In quantum electrodynamics for example, the mass of the electron varies according to the energy scale you're working at : this comes from the fact that the phenomenons involved in QED vary with the energy scale of the problem.
Also note that you cannot make any arbitrary choice of constants, it has to be consistent. For instance you cannot put $c=1$, $\hbar=1$ as well as the QED coupling constant $e=1$, because this would make the fine structure constant $\alpha = \frac{e^2}{\hbar c} = 1$, thus making the electromagnetic interaction much stronger than it really is ($\alpha = \frac{1}{137}$ at the atomic scale).
The convention called natural units can be separated into two parts:
Using a different set of base units. For example, instead of using the SI units $\text{kg}, \text{m}$ and $\text{s}$ as base units for a certain subspace of the unit space, the most common flavour of natural units uses $\hbar, c$ and $\text{eV}$ as base units for the same subspace. For example $1\,\text{fm}$ written in these base units equals $\frac{\hbar c}{197\,\text{MeV}}$.
Not writing down some or all of these base units. E.g., in the most common flavour of natural units, $\hbar$ and $c$ are not written down. This is usually denoted as $\hbar=1; c=1$. You could as well not write down $\text{eV}$.
Convention 1 is practical as you have frequently used constants as base units. Convention 2 spares you the trouble of writing down some of these base units and builds upon the fact that you can always reconstruct the unit from the dimension of the quantity you are considering. Note that even when using natural units, writing $\hbar = c$ is technically wrong as the unwritten units do not match.
Taking this point of view, action is not dimensionless in natural units – you just do not write down the units that indicate the dimension.
Further reading and blatant self-advertising: I wrote a didactical paper (preprint) on the topic of natural units.
