There's a lot going on here so I'll take a stab at addressing some of the issues.
Unbalanced data terminology
The "unbalanced data" terminology is misleading in my view. When we have several groups and measure the value of a response variable (e.g., pH) for different numbers of study units within each group, we say that the groups are unbalanced. In your case, the groups would be all possible combinations of SoilType and LandUse. It's not clear from your description what your study unit is (e.g., soil sample) but hopefully you have several study units per group.
Combinations of levels of SoilType and LandUse
What I think you mean by "unbalanced data" is the fact that you did not consider all possible combinations of levels of SoilType and LandUse in your study - you only considered a subset of these combinations. This would explain the error message you get from R regarding the singular gls fit. You can verify if this is the case for your data using the R command:
with(mydata, table(SoilType, LandUse))
The command will show you how often each possible combination of levels of SoilType and LandUse appears in your data. If a particular combination does not appear at all in your data, the contingency table produced by the above command will show an entry of 0 for it.
Addressing missing combinations of levels of SoilType and LandUse
If you are missing some of all possible combinations of SoilType and LandUse but you try to fit a gls model including an interaction between SoilType and LandUse, you will get the dreaded singular fit message from R.
This is because the design matrix for your model is degenerate. See this excellent post on Degenerate design matrices for more details: https://rstudio-pubs-static.s3.amazonaws.com/6311_a09169ad892f4f5499874751a5fa822d.html.
To circumvent this message, you would need to "manually" define a variable which indicates what specific combination of levels of SoilType and LandUse you are dealing with for each of the study units represented in your data. For example:
mydata$Combination <- interaction(mydata$SoilType,
mydata$LandUse)
mydata$Combination <- droplevels(mydata$Combination)
Then you can fit a gls model which does not use the weights = option like this:
gls(pH ~ Combination, data = mydata)
Here, I assume for simplicity there are no missing values in your data.
Simple example
The link on degenerate design matrices includes a simple example where Combination would be obtained by listing those combinations present in the data for two factors: f1 and f2. (In your case, f1 would be SoilType and f2 would be LandUse). In the example, f1 has 2 possible levels: a and b; f2 has 2 possible levels: A and B. All possible combinations of levels of f1 and f2 would be: aA, aB, bA, bB. However, in the data underlying that example, only the combinations aA, bA and bB are represented. If the gls() model is fitted as:
gls(y ~ 1 + Combination, data = etc.)
then R will look at the levels of Combination: aA, bA and bB. It will then set aside the first level, aA, and report 3 lines of output. The first line of output will report the estimated mean of y when Combination = aA. The second line of output will report the difference in the estimated mean value of y between Combination = bA and Combination = aA. The third line of output will report the difference in the estimated mean value of y between Combination = bB and Combination = aA. (Combination is denoted by f12 in the example).
If you compare the model gls(y ~ 1 + Combination, data = etc.) against the model gls(y ~ 1, data = etc.), you will be able to test these hypotheses:
Ho: there is no difference in the mean value of y across any of the combinations of levels of f1 and f2 represented in the study
Ha: there is a difference in the mean value of y across at least two of the combinations of levels of f1 and f2 represented in the study
The null hypothesis Ho corresponds to the model gls(y ~ 1, data = etc.), whereas the alternative hypothesis Ha corresponds to the model gls(y ~ 1 + Combination, data = etc.).
Of course, if Ho is rejected in favour of Ha, one would have to detect which specific combinations are different from each other by performing post-hoc comparisons of combinations.
Another possibility here would be that we know upfront we are only interested in comparing certain combinations - in other words, we would have some planned comparisons of combinations in mind from the study design phase: say, we would only want to compare the mean value of y for combinations bB and bA. Then we wouldn't need to test Ho against Ha - we would simply proceed to perform the comparisons of combinations of interest. This can be done by setting up and testing/estimating contrasts of combinations.
This is simple enough to do with just 2 levels per factor. But in your case, you'll have a large number of combinations, so you may either want to consolidate your LandUse definition to a smaller number of levels, or identify some specific combinations of interest, etc. Otherwise, you will have a massive number of comparisons on your hands. In any event, adjusting for multiplicity of comparisons will be a must!
Non-Constant Variability
After fitting the model gls(pH ~ Combinations, data = mydata), you can plot its residuals against SoilType to see whether there is evidence of non-constant variability across SoilType. If there is, fit a model like:
gls(pH ~ 1 + Combination,
weights = varIdent(form=~1|SoilType),
data = mydata)
Then, if you are able to fit the model without running into memory allocation issues, you would still need to consider performing multiple comparisons of combinations of levels. In particular, should you perform all possible comparisons or just a select few identified at the study planning stage?
If you can't fit the above model due to memory allocation issues, you may have no choice but to consolidate the levels of LandUse. Just use a command like:
table(mydata$LandUse)
and merge levels which are minimally represented in the data with levels that are maximally represented provided that makes substantive sense.