In short, natural atomic orbitals are the orbitals you get from diagonalizing the one-particle reduced density matrix. One can then localize this set of orbitals via some maximization criteria. Natural bond orbitals are the orbitals you get when maximizing the occupancy of orbitals such that electrons occupy the space between two atoms or on a single atom. These are then interpreted as bonds and lone pairs. The definition of these criteria somewhat hard to explain as you have to define the natural hybrid orbitals first and go from there.
All localization schemes which don't begin with natural atomic orbitals define some arbitrary measure which should intuitively localize the orbitals and then maximize or minimize this measure. For instance, one might begin with the atomic orbitals and maximize electron repulsion, minimize orbital overlap on different atoms, or maximize partial charge on each atom, etc.
Orbital localization schemes and NBO orbitals both tend to result in orbitals that are localized and sort of look like the atomic orbitals chemists are familiar with. They are somewhat different from a theoretical point of view however.
All orbital localization schemes (let's just call this a definition in case it's not strictly true) are unitary transformations of the canonical Hartree-Fock (HF) orbitals or some other choice of orbitals (e.g. post-HF). This means that the energy is unchanged under this transformation.
For instance, the Edmiston-Reudenberg localization scheme iteratively maximizes the electron self-repulsion. Intuitively, this will make the orbitals compact since these integrals will be over all electrons and all orbitals. My understanding is that this method results in very good orbital localization, but it is not popular as it requires the calculation of many integrals which one would not ordinarily need to calculate.
The Pipek-Mezey localization scheme takes the same idea but maximizes the partial charges in each orbital. These charges are trivial to calculate at virtually no extra expense as they are just related to the population of each orbital. The downside is that these charges are known to be arbitrary and highly basis set dependent. Nonetheless, the results work out quite well, so this is the method I have seen used the most.
Natural Orbitals (NO) result in localized orbitals by diagonalizing the one-particle reduced density matrix (1-RDM), which can be formed from canonical HF orbitals (or other orbitals). For explanation of what that means see this Chem.se question.
NBOs are a further transformation of these NOs which variationally maximize the Lewis-like character of the orbitals. Basically, it maximally localizes the orbitals either between pairs of atoms or on a single atom. This sort of defines the dominant Lewis structure for the system.
As an aside, I think it's somewhat dangerous to take any set of orbitals too seriously. The wavefunction and orbitals are not observable, so any inference you make on the basis of the wavefunction itself is coming from something which is not observable. Maybe the orbitals lead you to make a prediction about some observable quantity, but this always feels a bit ad hoc to me.
I also have it in my mind that there are certain types of orbitals derived from natural atomic orbitals which are not unitarily related to the atomic orbital basis used to form the 1-RDM. I can't find anything that says this though. I feel like I've heard of schemes where certain off-diagonal elements of the 1-RDM are thrown out if they are quite small. This would make the transformation non-unitary, but I don't know where that idea came from, so I could be mis-remembering something.
My understanding is that NBOs have been used to provide some physical interpretation of the electronic wavefunction in simple to understand terms. I don't know how common this is today. I don't seem to see many theoretical papers which do this, but I would think it might be more common in physical organic chemistry, and I don't commonly read papers in that field.
Localization schemes are more commonly used, in my experience, to speed up certain types of post-HF calculations. This works because electron correlation between a particular pair of electrons falls off very quickly when the one-electron wavefunctions have a very small overlap. Thus, many schemes have been devised to get MP2 or CCSD energies at a cheaper cost than the full method by working in a localized basis. DLPNO-CCSD(T) is probably the most well-known method of this type. Generally, one can perform more effective integral screening when working in a localized basis as well.
