The short answer is that LS means (or more modernly, estimated marginal means) are incredibly useful with experimental data. With observational data, not so much. A long-winded explanation follows.
The underlying ideas are very old (and predate SAS by at least 50 years). Look at a standard experimental design textbook -- pretty much any of them. Such a book has initial chapters on one-way CRDs, focusing especially on balanced designs. Then they move on to factorial designs, and it is a 2-way factorial I will discuss the most. For that, we randomize experimental units to combinations of the levels of two factors $A$ and $B$, and consider analyzing a model of the form
$$ Y_{ijk} = \mu_{ij} + E_{ijk}, \quad i=1,2,...,a;\quad j=1,2,...,b;\quad k=1,2,...,n_{ij} $$
where $Y_{ijk}$ is the $k$th observation at the $i$th level of $A$, $j$th level of $B$, the $\mu_{ij}$ are the cell means, and $E_{ijk}$ are random errors, typically assumed to be iid normal. Emphasis is usually placed on balanced experiments, where the $n_{ij}$ are all equal.
We proceed to talk about partitioning the total sum of squares into components attributable to marginal effects of $A$, marginal effects of $B$, perhaps interaction effects, and residual error.
If the interactions are not strong, we are attracted to an additive model whereby $\mu_{ij} = \mu + \alpha_i + \beta_j$; in that case, interpretation is simpler because the effects of $A$ and $B$ can be considered independently. Evaluating these effects can be expressed in terms of the marginal means for $A$, $\mu_{i.}=\mu+\alpha_i$, and for $B$, $\mu_{.j}=\mu+\beta_j$; and perhaps pairwise comparisons thereof. In a balanced design, these marginal means may be estimated as the marginal means of the data: $\bar Y_{i..}$ and $\bar Y_{.j.}$.
Balanced experiments lie in very foundational places in statistics, including the writings of Fisher, Youden, Snedecor, Cochran, Cox, Box, and Tukey. But what if it isn't a balanced experiment? For example, a few observations were lost (completely at random). In very old design texts, what was suggested in some cases was the "method of unweighted means", whereby we use the values of $\bar Y_{ij.}$ and construct an analysis of those cell means as if they were data from a balanced experiment with $\tilde n$ observations per cell (where, typically, $\tilde n$ is the harmonic mean of the cell sample sizes). This is far preferable to just computing marginal means of the data, because some cells receive more weight than others, which can produce Simpson's-paradox-like effects.
"Least-square means" are essentially a model-based version of unweighted means. They were developed by Walter Harvey in a technical report in 1960 and finally published as
Walter R. Harvey (1982), "Least-Squares Analysis of Discrete Data,"
Journal of Animal Science, 54(5), 1067–1071, https://doi.org/10.2527/jas1982.5451067x. They were indeed implemented in SAS in an add-on procedure called PROC HARVEY, not long after the first (1972) edition of SAS was released, and then later implemented as part of the new PROC GLM released in SAS75. You may also want to take a look at a 2014 retrospective by Weijie Cai.
Again, LS means are essentially the same idea as unweighted means, which is a very, very old idea. In LS means, we fit a model to the data and use it (in the two-way factorial case) to predict the $\mu_{ij}$; then our marginal means are estimated as equally-weighted marginal averages of these predictions, just as in unweighted-means analysis. But these LS means are linear functions of the model predictions, and hence of the regression coefficients; so they can be estimated from the model with no need for ad hoc quantities like $\tilde n$. The idea extends to other experimental designs, such as Latin squares, split-plots, etc., and even to covariance models (where typically we make predictions at the mean of each covariate). You can also generalize these ideas to experimental data analyzed with other types of models, including generalized linear models, ordinal models and multinomial models, zero-inflated models, etc. And the ideas can be applied to Bayesian models fitted via MCMC methods, simply by computing the relevant quantities from the posterior predictions.
Is this useful? Definitely yes, when the model is additive and/or the interaction effects are small. They mimic the types of analyses proposed for balanced experiments by the giants among our forebears. If, on the other hand, we do have sizeable interaction effects, then (as our forebears recommend) we simply should not even consider estimating marginal means. (In a multi-factor experiment, it may still be reasonable to construct some marginal means where we are averaging only over negligible interaction effects).
In the situation where you have data sampled from a population, then LS means or EMMs may very well be entirely inappropriate. They certainly don't estimate the marginal means of your population because the equal weighting is not appropriate. The emmeans package provides various alternative weighting schemes, and some may be suited for a particular purpose; or not.
What about type III sums of squares? What I would say is that you need to consider whether they are appropriate for an individual situation. SAS defines these in terms of estimable functions, and a given set of estimable functions corresponds to a family of contrasts that can usually be expressed as contrasts among LS means. If that family of contrasts is of interest, then the associated type III test is fine; and if not, it's not fine. And again, this has to do with the appropriateness of the LS means being considered. I personally shy away from type III tests, but that is because I don't use omnibus tests for anything but model selection.
An additional reference is Searle, Speed, and Milliken (1980), "Population marginal means in the linear model: An alternative to least squares means", The American Statistician 34(4), 216-221, doi:10.1080/00031305.1980.10483031, where they critique Harvey's method and propose the term "estimated marginal means." They also discuss instances where there are empty cells. In that part, I personally think they do it all wrong, but it depends of course on what you are trying to accomplish (if the cells are empty because those factor combinations actually cannot happen, then what they do is more defensible). In general, empty cells can be terribly problematic, and the emmeans package deals with them in terms of assessing estimability of the specified quantities.
I would also say that there is a definite anti-SAS attitude among some members of the R community (and S and S-Plus before R). This is unfortunate, because SAS is carefully developed under the direction of many top professional statisticians. It is well documented and its code is solid and reliable for what it does. There are definite reasons I prefer R to SAS myself, not the least of which is that SAS cannot quite escape some very old ruts. But I would like to encourage users to adopt a less tribal attitude about statistical software.
