Is there a way to see linear and surface charge density as a “special case” of volume charge density?

You can always think of line- and surface-densities as a collection of point charges. If you already have no problem in using Dirac deltas, then we can easily just extend them to more complex geometries. Weird deltas might be harder to handle, but in principle it is possible.

For example you could write, for a series of $N$ charges of value $q$ along a line at positions $x_i$ along the $x$-axis, something like

$$\rho(\vec{x})=\Sigma^N_i q \delta(x-x_i)$$

clearly this is not problematic. It's just a sum of point charges.

If you extend this reasoning to continuum space [see below for a quick derivation], you can use Dirac's delta functions that are $0$ everywhere except along a line. For example, if we choose the $x-axis$ you could write $$\rho(\vec{x})=\lambda(x) \delta(y)\delta(z)$$

where instead of the vector $\vec{x}$ I used the individual components and $\lambda(x)$ is the charge density along the $x$ axis. Those deltas are zero except where both $y=0$ and $z=0$ i.e. along the $x-axis$.

If you integrate that over the full space (integrals all go from $0$ to $\infty$), you get

$$\int\rho(\vec{x}) d\vec{x}= \int \lambda(x) \delta(y)\delta(z) dx dy dz =\int \lambda(x) dx \int \delta(y)dy\int\delta(z)dz$$ and because $\int \delta(y)dy=1$ and the same for z

$$\int\rho(\vec{x}) d\vec{x}=\int \lambda(x) dx =Q$$ which gives you the total charge as expected (which can be infinite, is $\lambda$ is taken constant, although non-physical).

In the same way, we can define $$\rho(\vec{x}) = \sigma(x, y) \delta(z)$$ which is only non-zero if $z=0$ i.e. on the full $xy$-plane.

You can then create any charge density you want, using combined delta functions. You could also define it "by hand", just saying "here is constant and here is zero" and then integrating each domain separatly if necessary, but that will give you a bit more of a problem in dealing with discontinuities, so you need to careful.

For example, of course you can also limit your charge distribution in space by defining

$$\rho(\vec{x})=\lambda(x)\delta(y)\delta(z)$$ if $|x|<L/2$ and 0 otherwise, which would give you a thin wire along the the $x$-axis of length $L$ - or you could define $\lambda(x)$ using step functions. That would give you a total charge of $Q$ along such wire which, if assumed at constant density, would give you $\lambda(x)=Q/L$.

So yes, you can use GLDF in a very general way. You just get the annoying thing of having to deal with $\delta$s.

Derivation

As an exercise, we can derive this expression. There is a bit of hand-waving but I think it sends you in the right direction.

Imagine you have charges along the $x-axis$ of value $\lambda(x)dx$ - i.e. the local density times an infinitesimal length.

We can write the density $\rho(\vec{x})$ as a sum of point charges at positions $x'$ (with $x'$ is a continous variable spanning the $x$ axis).

$$\rho(\vec{x}) = \Sigma \lambda(x')dx'\delta(\vec{x}-\vec{x'})$$

now of course, because we are dealing with infinitesimal charges, that sum actually is an integral along the line

$$\rho(\vec{x}) = \int \lambda(x')dx'\delta(\vec{x}-\vec{x'}) $$

and if we write $\delta(\vec{x}-\vec{x'})=\delta(x-x')\delta(y-y')\delta(z-z')$ we have to perform only the integral on the variable $x'$

$$\rho(\vec{x}) = \delta(z-z')\delta(y-y') \int \lambda(x')\delta(x-x')dx'$$ and using the fact that $\int f(x')\delta(x-x')dx' = f(x)$ $$\rho(\vec{x}) = \delta(z-z')\delta(y-y')\lambda(x)$$

and because $y'$ and $z'$ are $0$ as we are on the $x$-axis, then $y-y'\to y$ and $z-z'\to z$ and we get $$\rho(\vec{x}) = \delta(z)\delta(y)\lambda(x)$$